Thyra_EpetraThyraWrappers.cpp

Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //               Thyra: Trilinos Solver Framework Core
00005 //                 Copyright (2004) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #include "Thyra_EpetraThyraWrappers.hpp"
00030 #include "Thyra_DefaultSpmdVectorSpace.hpp"
00031 #include "Thyra_DefaultSpmdMultiVector.hpp"
00032 #include "Thyra_DefaultSpmdVector.hpp"
00033 #include "Thyra_DetachedVectorView.hpp"
00034 #include "Thyra_DetachedMultiVectorView.hpp"
00035 #include "Teuchos_dyn_cast.hpp"
00036 
00037 #include "Teuchos_DefaultSerialComm.hpp"
00038 #ifdef HAVE_MPI
00039 #include "Teuchos_DefaultMpiComm.hpp"
00040 #endif
00041 
00042 #include "Epetra_Map.h"
00043 #include "Epetra_Comm.h"
00044 #include "Epetra_SerialComm.h"
00045 #ifdef HAVE_MPI
00046 #  include "Epetra_MpiComm.h"
00047 #endif
00048 #include "Epetra_MultiVector.h"
00049 #include "Epetra_Vector.h"
00050 
00051 Teuchos::RefCountPtr<const Teuchos::Comm<Thyra::Index> >
00052 Thyra::create_Comm( const Teuchos::RefCountPtr<const Epetra_Comm> &epetraComm )
00053 {
00054   using Teuchos::RefCountPtr;
00055   using Teuchos::rcp;
00056   using Teuchos::rcp_dynamic_cast;
00057   using Teuchos::set_extra_data;
00058 
00059   RefCountPtr<const Epetra_SerialComm>
00060     serialEpetraComm = rcp_dynamic_cast<const Epetra_SerialComm>(epetraComm);
00061   if( serialEpetraComm.get() ) {
00062     RefCountPtr<const Teuchos::SerialComm<Index> >
00063       serialComm = rcp(new Teuchos::SerialComm<Index>());
00064     set_extra_data( serialEpetraComm, "serialEpetraComm", &serialComm );
00065     return serialComm;
00066   }
00067 
00068 #ifdef HAVE_MPI
00069   
00070   RefCountPtr<const Epetra_MpiComm>
00071     mpiEpetraComm = rcp_dynamic_cast<const Epetra_MpiComm>(epetraComm);
00072   if( mpiEpetraComm.get() ) {
00073     RefCountPtr<const Teuchos::OpaqueWrapper<MPI_Comm> >
00074       rawMpiComm = Teuchos::opaqueWrapper(mpiEpetraComm->Comm());
00075     set_extra_data( mpiEpetraComm, "mpiEpetraComm", &rawMpiComm );
00076     RefCountPtr<const Teuchos::MpiComm<Index> >
00077       mpiComm = rcp(new Teuchos::MpiComm<Index>(rawMpiComm));
00078     return mpiComm;
00079   }
00080 
00081 #endif // HAVE_MPI
00082   
00083   // If you get here then the failed!
00084   return Teuchos::null;
00085 
00086 }
00087 
00088 Teuchos::RefCountPtr<const Thyra::SpmdVectorSpaceDefaultBase<double> >
00089 Thyra::create_VectorSpace(
00090   const Teuchos::RefCountPtr<const Epetra_Map> &epetra_map
00091   )
00092 {
00093 #ifdef TEUCHOS_DEBUG
00094   TEST_FOR_EXCEPTION( !epetra_map.get(), std::invalid_argument, "create_VectorSpace::initialize(...): Error!" );
00095 #endif // TEUCHOS_DEBUG
00096   Teuchos::RefCountPtr<const Teuchos::Comm<Index> >
00097     comm = create_Comm(Teuchos::rcp(&epetra_map->Comm(),false)).assert_not_null();
00098   Teuchos::set_extra_data( epetra_map, "epetra_map", &comm );
00099   const Index localSubDim = epetra_map->NumMyElements();
00100   Teuchos::RefCountPtr<DefaultSpmdVectorSpace<double> >
00101     vs = Teuchos::rcp(
00102       new DefaultSpmdVectorSpace<double>(
00103         comm
00104         ,localSubDim
00105         ,epetra_map->NumGlobalElements()
00106         )
00107       );
00108 #ifndef TEUCHOS_DEBUG
00109       TEST_FOR_EXCEPTION(
00110         vs->dim() != epetra_map->NumGlobalElements(), std::logic_error
00111         ,"create_VectorSpace(...): Error, vs->dim() = "<<vs->dim()<<" != "
00112         "epetra_map->NumGlobalElements() = "<<epetra_map->NumGlobalElements()<<"!"
00113         );
00114 #endif    
00115   Teuchos::set_extra_data( epetra_map, "epetra_map", &vs );
00116   return vs;
00117 }
00118 
00119 Teuchos::RefCountPtr<Thyra::SpmdVectorBase<double> >
00120 Thyra::create_Vector(
00121   const Teuchos::RefCountPtr<Epetra_Vector>                                &epetra_v
00122   ,const Teuchos::RefCountPtr<const SpmdVectorSpaceBase<double> >          &space
00123   )
00124 {
00125 #ifdef TEUCHOS_DEBUG
00126   TEST_FOR_EXCEPT(space.get()==NULL);
00127 #endif
00128   if(!epetra_v.get()) return Teuchos::null;
00129   // New local view of raw data
00130   double *localValues;
00131   epetra_v->ExtractView( &localValues );
00132   // Build the Vector
00133   Teuchos::RefCountPtr<SpmdVectorBase<double> >
00134     v = Teuchos::rcp(new DefaultSpmdVector<double>(space,Teuchos::rcp(localValues,false),1));
00135   Teuchos::set_extra_data( epetra_v, "Epetra_Vector", &v );
00136   return v;
00137 }
00138 
00139 Teuchos::RefCountPtr<const Thyra::SpmdVectorBase<double> >
00140 Thyra::create_Vector(
00141   const Teuchos::RefCountPtr<const Epetra_Vector>                          &epetra_v
00142   ,const Teuchos::RefCountPtr<const SpmdVectorSpaceBase<double> >           &space
00143   )
00144 {
00145 #ifdef TEUCHOS_DEBUG
00146   TEST_FOR_EXCEPT(space.get()==NULL);
00147 #endif
00148   if(!epetra_v.get()) return Teuchos::null;
00149   // New local view of raw data
00150   double *localValues;
00151   epetra_v->ExtractView( &localValues );
00152   // Build the Vector
00153   Teuchos::RefCountPtr<const SpmdVectorBase<double> >
00154     v = Teuchos::rcp(new DefaultSpmdVector<double>(space,Teuchos::rcp(localValues,false),1));
00155   Teuchos::set_extra_data( epetra_v, "Epetra_Vector", &v );
00156   return v;
00157 }
00158 
00159 Teuchos::RefCountPtr<Thyra::SpmdMultiVectorBase<double> >
00160 Thyra::create_MultiVector(
00161   const Teuchos::RefCountPtr<Epetra_MultiVector>                           &epetra_mv
00162   ,const Teuchos::RefCountPtr<const SpmdVectorSpaceBase<double> >           &range
00163   ,const Teuchos::RefCountPtr<const ScalarProdVectorSpaceBase<double> >    &domain
00164   )
00165 {
00166 #ifdef TEUCHOS_DEBUG
00167   TEST_FOR_EXCEPT(range.get()==NULL);
00168   TEST_FOR_EXCEPT(domain.get()==NULL);
00169 #endif
00170   if(!epetra_mv.get()) return Teuchos::null;
00171   // New local view of raw data
00172   double *localValues; int leadingDim;
00173   if( epetra_mv->ConstantStride() ) {
00174     epetra_mv->ExtractView( &localValues, &leadingDim );
00175   }
00176   else {
00177     TEST_FOR_EXCEPT(true); // ToDo: Implement!
00178   }
00179   // Build the MultiVector
00180   Teuchos::RefCountPtr<SpmdMultiVectorBase<double> >
00181     mv = Teuchos::rcp(new DefaultSpmdMultiVector<double>(range,domain,Teuchos::rcp(localValues,false),leadingDim));
00182   Teuchos::set_extra_data( epetra_mv, "Epetra_MultiVector", &mv );
00183   return mv;
00184 }
00185 
00186 Teuchos::RefCountPtr<const Thyra::SpmdMultiVectorBase<double> >
00187 Thyra::create_MultiVector(
00188   const Teuchos::RefCountPtr<const Epetra_MultiVector>                     &epetra_mv
00189   ,const Teuchos::RefCountPtr<const SpmdVectorSpaceBase<double> >           &range
00190   ,const Teuchos::RefCountPtr<const ScalarProdVectorSpaceBase<double> >    &domain
00191   )
00192 {
00193 #ifdef TEUCHOS_DEBUG
00194   TEST_FOR_EXCEPT(range.get()==NULL);
00195   TEST_FOR_EXCEPT(domain.get()==NULL);
00196 #endif
00197   if(!epetra_mv.get()) return Teuchos::null;
00198   // New local view of raw data
00199   double *localValues; int leadingDim;
00200   if( epetra_mv->ConstantStride() ) {
00201     epetra_mv->ExtractView( &localValues, &leadingDim );
00202   }
00203   else {
00204     TEST_FOR_EXCEPT(true); // ToDo: Implement!
00205   }
00206   // Build the MultiVector
00207   Teuchos::RefCountPtr<const SpmdMultiVectorBase<double> >
00208     mv = Teuchos::rcp(new DefaultSpmdMultiVector<double>(range,domain,Teuchos::rcp(localValues,false),leadingDim));
00209   Teuchos::set_extra_data( epetra_mv, "Epetra_MultiVector", &mv );
00210   return mv;
00211 }
00212 
00213 Teuchos::RefCountPtr<Epetra_Vector>
00214 Thyra::get_Epetra_Vector(
00215   const Epetra_Map                                    &map
00216   ,const Teuchos::RefCountPtr<VectorBase<double> >    &v
00217   )
00218 {
00219 #ifdef TEUCHOS_DEBUG
00220   Teuchos::RefCountPtr<const VectorSpaceBase<double> >
00221     epetra_vs = create_VectorSpace(Teuchos::rcp(&map,false));
00222   THYRA_ASSERT_VEC_SPACES( "Thyra::get_Epetra_Vector(map,v)", *epetra_vs, *v->space() );
00223 #endif
00224   //
00225   // First, try to grab the Epetra_Vector straight out of the
00226   // RCP since this is the fastest way.
00227   //
00228   const Teuchos::RefCountPtr<Epetra_Vector>
00229     *epetra_v_ptr = Teuchos::get_optional_extra_data<Teuchos::RefCountPtr<Epetra_Vector> >(v,"Epetra_Vector");
00230   if(epetra_v_ptr) {
00231     return *epetra_v_ptr;
00232   }
00233   //
00234   // The assumption that we (rightly) make here is that if the vector spaces
00235   // are compatible, that either the multi-vectors are both in-core or the
00236   // vector spaces are both derived from SpmdVectorSpaceBase and have
00237   // compatible maps.
00238   // 
00239   const VectorSpaceBase<double>  &vs = *v->range();
00240   const SpmdVectorSpaceBase<double> *mpi_vs = dynamic_cast<const SpmdVectorSpaceBase<double>*>(&vs);
00241   const Index localOffset = ( mpi_vs ? mpi_vs->localOffset() : 0 );
00242   const Index localSubDim = ( mpi_vs ? mpi_vs->localSubDim() : vs.dim() );
00243   //
00244   // Here we will extract a view of the local elements in the underlying
00245   // VectorBase object.  In most cases, no data will be allocated or copied
00246   // and only some small objects will be created and a few virtual functions
00247   // will be called so the overhead should be low and fixed.
00248   //
00249   // Create a *mutable* view of the local elements, this view will be set on
00250   // the RefCountPtr that is returned.  As a result, this view will be relased
00251   // when the returned Epetra_Vector is released.
00252   //
00253   Teuchos::RefCountPtr<DetachedVectorView<double> >
00254     emvv = Teuchos::rcp(
00255       new DetachedVectorView<double>(
00256         v
00257         ,Range1D(localOffset,localOffset+localSubDim-1)
00258         ,true // forceContiguous
00259         )
00260       );
00261   // Create a temporary Epetra_Vector object and give it
00262   // the above local view
00263   Teuchos::RefCountPtr<Epetra_Vector>
00264     epetra_v = Teuchos::rcp(
00265       new Epetra_Vector(
00266         ::View                                 // CV
00267         ,map                                   // Map
00268         ,const_cast<double*>(emvv->values())   // V
00269         )
00270       );
00271   // Give the explict view object to the above Epetra_Vector
00272   // smart pointer object.  In this way, when the client is finished
00273   // with the Epetra_Vector view the destructor from the object
00274   // in emvv will automatically commit the changes to the elements in
00275   // the input v MultiVectorBase object (reguardless of its
00276   // implementation).  This is truly an elegant result!
00277   Teuchos::set_extra_data( emvv, "emvv", &epetra_v, Teuchos::PRE_DESTROY );
00278   // We are done!
00279   return epetra_v;
00280 }
00281 
00282 Teuchos::RefCountPtr<const Epetra_Vector>
00283 Thyra::get_Epetra_Vector(
00284   const Epetra_Map                                         &map 
00285   ,const Teuchos::RefCountPtr<const VectorBase<double> >   &v
00286   )
00287 {
00288 #ifdef TEUCHOS_DEBUG
00289   Teuchos::RefCountPtr<const VectorSpaceBase<double> >
00290     epetra_vs = create_VectorSpace(Teuchos::rcp(&map,false));
00291   THYRA_ASSERT_VEC_SPACES( "Thyra::get_Epetra_Vector(map,v)", *epetra_vs, *v->space() );
00292 #endif
00293   //
00294   // First, try to grab the Epetra_Vector straight out of the
00295   // RCP since this is the fastest way.
00296   //
00297   const Teuchos::RefCountPtr<const Epetra_Vector>
00298     *epetra_v_ptr = Teuchos::get_optional_extra_data<Teuchos::RefCountPtr<const Epetra_Vector> >(v,"Epetra_Vector");
00299   if(epetra_v_ptr)
00300     return *epetra_v_ptr;
00301   const Teuchos::RefCountPtr<Epetra_Vector>
00302     *epetra_nonconst_v_ptr = Teuchos::get_optional_extra_data<Teuchos::RefCountPtr<Epetra_Vector> >(v,"Epetra_Vector");
00303   if(epetra_nonconst_v_ptr)
00304     return *epetra_nonconst_v_ptr;
00305   //
00306   // Same as above function except as stated below
00307   //
00308   const VectorSpaceBase<double> &vs = *v->range();
00309   const SpmdVectorSpaceBase<double> *mpi_vs = dynamic_cast<const SpmdVectorSpaceBase<double>*>(&vs);
00310   const Index localOffset = ( mpi_vs ? mpi_vs->localOffset() : 0 );
00311   const Index localSubDim = ( mpi_vs ? mpi_vs->localSubDim() : vs.dim() );
00312   // Get an explicit *non-mutable* view of all of the elements in
00313   // the multi vector.
00314   Teuchos::RefCountPtr<ConstDetachedVectorView<double> >
00315     evv = Teuchos::rcp(
00316       new ConstDetachedVectorView<double>(
00317         *v
00318         ,Range1D(localOffset,localOffset+localSubDim-1)
00319         ,true // forceContiguous
00320         )
00321       );
00322   // Create a temporary Epetra_Vector object and give it
00323   // the above view
00324   Teuchos::RefCountPtr<Epetra_Vector>
00325     epetra_v = Teuchos::rcp(
00326       new Epetra_Vector(
00327         ::View                                  // CV
00328         ,map                                    // Map
00329         ,const_cast<double*>(evv->values())     // V
00330         )
00331       );
00332   // This next statement will cause the destructor to free the view if
00333   // needed (see above function).  Since this view is non-mutable,
00334   // only a releaseDetachedView(...) and not a commit will be called.
00335   // This is the whole reason there is a seperate implementation for
00336   // the const and non-const cases.
00337   Teuchos::set_extra_data( evv, "evv", &epetra_v, Teuchos::PRE_DESTROY );
00338   // Also set the v itself as extra data just to be safe
00339   Teuchos::set_extra_data( v, "v", &epetra_v );
00340   // We are done!
00341   return epetra_v;
00342 }
00343 
00344 Teuchos::RefCountPtr<Epetra_MultiVector>
00345 Thyra::get_Epetra_MultiVector(
00346   const Epetra_Map                                         &map
00347   ,const Teuchos::RefCountPtr<MultiVectorBase<double> >    &mv
00348   )
00349 {
00350 #ifdef TEUCHOS_DEBUG
00351   Teuchos::RefCountPtr<const VectorSpaceBase<double> >
00352     epetra_vs = create_VectorSpace(Teuchos::rcp(&map,false));
00353   THYRA_ASSERT_VEC_SPACES( "Thyra::get_Epetra_MultiVector(map,mv)", *epetra_vs, *mv->range() );
00354 #endif
00355   //
00356   // First, try to grab the Epetra_MultiVector straight out of the
00357   // RCP since this is the fastest way.
00358   //
00359   const Teuchos::RefCountPtr<Epetra_MultiVector>
00360     *epetra_mv_ptr = Teuchos::get_optional_extra_data<Teuchos::RefCountPtr<Epetra_MultiVector> >(mv,"Epetra_MultiVector");
00361   if(epetra_mv_ptr) {
00362     return *epetra_mv_ptr;
00363   }
00364   //
00365   // The assumption that we (rightly) make here is that if the vector spaces
00366   // are compatible, that either the multi-vectors are both in-core or the
00367   // vector spaces are both derived from SpmdVectorSpaceBase and have
00368   // compatible maps.
00369   // 
00370   const VectorSpaceBase<double> &vs = *mv->range();
00371   const SpmdVectorSpaceBase<double> *mpi_vs = dynamic_cast<const SpmdVectorSpaceBase<double>*>(&vs);
00372   const Index localOffset = ( mpi_vs ? mpi_vs->localOffset() : 0 );
00373   const Index localSubDim = ( mpi_vs ? mpi_vs->localSubDim() : vs.dim() );
00374   //
00375   // Here we will extract a view of the local elements in the underlying
00376   // MultiVectorBase object.  In most cases, no data will be allocated or
00377   // copied and only some small objects will be created and a few virtual
00378   // functions will be called so the overhead should be low and fixed.
00379   //
00380   // Create a *mutable* view of the local elements, this view will be set on
00381   // the RefCountPtr that is returned.  As a result, this view will be relased
00382   // when the returned Epetra_MultiVector is released.
00383   //
00384   Teuchos::RefCountPtr<DetachedMultiVectorView<double> >
00385     emmvv = Teuchos::rcp(
00386       new DetachedMultiVectorView<double>(
00387         *mv
00388         ,Range1D(localOffset,localOffset+localSubDim-1)
00389         )
00390       );
00391   // Create a temporary Epetra_MultiVector object and give it
00392   // the above local view
00393   Teuchos::RefCountPtr<Epetra_MultiVector>
00394     epetra_mv = Teuchos::rcp(
00395       new Epetra_MultiVector(
00396         ::View                                  // CV
00397         ,map                                    // Map
00398         ,const_cast<double*>(emmvv->values())   // A
00399         ,emmvv->leadingDim()                    // MyLDA
00400         ,emmvv->numSubCols()                    // NumVectors
00401         )
00402       );
00403   // Give the explict view object to the above Epetra_MultiVector
00404   // smart pointer object.  In this way, when the client is finished
00405   // with the Epetra_MultiVector view the destructor from the object
00406   // in emmvv will automatically commit the changes to the elements in
00407   // the input mv MultiVectorBase object (reguardless of its
00408   // implementation).  This is truly an elegant result!
00409   Teuchos::set_extra_data( emmvv, "emmvv", &epetra_mv, Teuchos::PRE_DESTROY );
00410   // Also set the mv itself as extra data just to be safe
00411   Teuchos::set_extra_data( mv, "mv", &epetra_mv );
00412   // We are done!
00413   return epetra_mv;
00414 }
00415 
00416 Teuchos::RefCountPtr<const Epetra_MultiVector>
00417 Thyra::get_Epetra_MultiVector(
00418   const Epetra_Map                                              &map
00419   ,const Teuchos::RefCountPtr<const MultiVectorBase<double> >   &mv
00420   )
00421 {
00422 #ifdef TEUCHOS_DEBUG
00423   Teuchos::RefCountPtr<const VectorSpaceBase<double> >
00424     epetra_vs = create_VectorSpace(Teuchos::rcp(&map,false));
00425   THYRA_ASSERT_VEC_SPACES( "Thyra::get_Epetra_MultiVector(map,mv)", *epetra_vs, *mv->range() );
00426 #endif
00427   //
00428   // First, try to grab the Epetra_MultiVector straight out of the
00429   // RCP since this is the fastest way.
00430   //
00431   const Teuchos::RefCountPtr<const Epetra_MultiVector>
00432     *epetra_mv_ptr = Teuchos::get_optional_extra_data<Teuchos::RefCountPtr<const Epetra_MultiVector> >(mv,"Epetra_MultiVector");
00433   if(epetra_mv_ptr) {
00434     return *epetra_mv_ptr;
00435   }
00436   //
00437   // Same as above function except as stated below
00438   //
00439   const VectorSpaceBase<double> &vs = *mv->range();
00440   const SpmdVectorSpaceBase<double> *mpi_vs = dynamic_cast<const SpmdVectorSpaceBase<double>*>(&vs);
00441   const Index localOffset = ( mpi_vs ? mpi_vs->localOffset() : 0 );
00442   const Index localSubDim = ( mpi_vs ? mpi_vs->localSubDim() : vs.dim() );
00443   // Get an explicit *non-mutable* view of all of the elements in
00444   // the multi vector.
00445   Teuchos::RefCountPtr<ConstDetachedMultiVectorView<double> >
00446     emvv = Teuchos::rcp(
00447       new ConstDetachedMultiVectorView<double>(
00448         *mv
00449         ,Range1D(localOffset,localOffset+localSubDim-1)
00450         )
00451       );
00452   // Create a temporary Epetra_MultiVector object and give it
00453   // the above view
00454   Teuchos::RefCountPtr<Epetra_MultiVector>
00455     epetra_mv = Teuchos::rcp(
00456       new Epetra_MultiVector(
00457         ::View                                  // CV
00458         ,map                                    // Map
00459         ,const_cast<double*>(emvv->values())    // A
00460         ,emvv->leadingDim()                     // MyLDA
00461         ,emvv->numSubCols()                     // NumVectors
00462         )
00463       );
00464   // This next statement will cause the destructor to free the view if
00465   // needed (see above function).  Since this view is non-mutable,
00466   // only a releaseDetachedView(...) and not a commit will be called.
00467   // This is the whole reason there is a seperate implementation for
00468   // the const and non-const cases.
00469   Teuchos::set_extra_data( emvv, "emvv", &epetra_mv, Teuchos::PRE_DESTROY );
00470   // Also set the mv itself as extra data just to be safe
00471   Teuchos::set_extra_data( mv, "mv", &epetra_mv );
00472   // We are done!
00473   return epetra_mv;
00474 }

Generated on Thu Sep 18 12:37:58 2008 for Epetra Package Browser (Single Doxygen Collection) by doxygen 1.3.9.1