Thyra Version of the Day
Thyra_EpetraThyraWrappers.cpp
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //    Thyra: Interfaces and Support for Abstract Numerical Algorithms
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 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Roscoe A. Bartlett (bartlettra@ornl.gov) 
00038 // 
00039 // ***********************************************************************
00040 // @HEADER
00041 
00042 #include "Thyra_EpetraThyraWrappers.hpp"
00043 #include "Thyra_DefaultSpmdVectorSpace.hpp"
00044 #include "Thyra_DefaultSpmdMultiVector.hpp"
00045 #include "Thyra_DefaultSpmdVector.hpp"
00046 #include "Thyra_DetachedVectorView.hpp"
00047 #include "Thyra_DetachedMultiVectorView.hpp"
00048 #include "Thyra_VectorSpaceFactoryBase.hpp"
00049 #include "Thyra_ProductVectorSpaceBase.hpp"
00050 #include "Teuchos_Assert.hpp"
00051 #include "Teuchos_dyn_cast.hpp"
00052 
00053 #include "Teuchos_DefaultSerialComm.hpp"
00054 #ifdef HAVE_MPI
00055 #include "Teuchos_DefaultMpiComm.hpp"
00056 #endif
00057 
00058 #include "Epetra_Map.h"
00059 #include "Epetra_Comm.h"
00060 #include "Epetra_SerialComm.h"
00061 #ifdef HAVE_MPI
00062 #  include "Epetra_MpiComm.h"
00063 #endif
00064 #include "Epetra_MultiVector.h"
00065 #include "Epetra_Vector.h"
00066 
00067 //
00068 // Helpers
00069 //
00070 
00071 
00072 namespace {
00073 
00074 
00075 Teuchos::RCP<const Thyra::VectorSpaceBase<double> >
00076 unwrapSingleProductVectorSpace(
00077   const Teuchos::RCP<const Thyra::VectorSpaceBase<double> > &vs_in
00078   )
00079 {
00080   using Teuchos::RCP;
00081   using Teuchos::rcp_dynamic_cast;
00082   using Thyra::ProductVectorSpaceBase;
00083   const RCP<const ProductVectorSpaceBase<double> > pvs =
00084     rcp_dynamic_cast<const ProductVectorSpaceBase<double> >(vs_in);
00085   if (nonnull(pvs)) {
00086     TEUCHOS_ASSERT_EQUALITY( pvs->numBlocks(), 1 );
00087     return pvs->getBlock(0);
00088   }
00089   return vs_in;
00090 }
00091 
00092 
00093 } // namespace
00094 
00095 
00096 //
00097 // Implementations of user function
00098 //
00099 
00100 
00101 Teuchos::RCP<const Teuchos::Comm<Thyra::Ordinal> >
00102 Thyra::create_Comm( const RCP<const Epetra_Comm> &epetraComm )
00103 {
00104   using Teuchos::rcp;
00105   using Teuchos::rcp_dynamic_cast;
00106   using Teuchos::set_extra_data;
00107 
00108   RCP<const Epetra_SerialComm>
00109     serialEpetraComm = rcp_dynamic_cast<const Epetra_SerialComm>(epetraComm);
00110   if( serialEpetraComm.get() ) {
00111     RCP<const Teuchos::SerialComm<Ordinal> >
00112       serialComm = rcp(new Teuchos::SerialComm<Ordinal>());
00113     set_extra_data( serialEpetraComm, "serialEpetraComm", Teuchos::inOutArg(serialComm) );
00114     return serialComm;
00115   }
00116 
00117 #ifdef HAVE_MPI
00118   
00119   RCP<const Epetra_MpiComm>
00120     mpiEpetraComm = rcp_dynamic_cast<const Epetra_MpiComm>(epetraComm);
00121   if( mpiEpetraComm.get() ) {
00122     RCP<const Teuchos::OpaqueWrapper<MPI_Comm> >
00123       rawMpiComm = Teuchos::opaqueWrapper(mpiEpetraComm->Comm());
00124     set_extra_data( mpiEpetraComm, "mpiEpetraComm", Teuchos::inOutArg(rawMpiComm) );
00125     RCP<const Teuchos::MpiComm<Ordinal> >
00126       mpiComm = rcp(new Teuchos::MpiComm<Ordinal>(rawMpiComm));
00127     return mpiComm;
00128   }
00129 
00130 #endif // HAVE_MPI
00131   
00132   // If you get here then the failed!
00133   return Teuchos::null;
00134 
00135 }
00136 
00137 
00138 Teuchos::RCP<const Thyra::VectorSpaceBase<double> >
00139 Thyra::create_VectorSpace(
00140   const RCP<const Epetra_Map> &epetra_map
00141   )
00142 {
00143   using Teuchos::as; using Teuchos::inoutArg;
00144 #ifdef TEUCHOS_DEBUG
00145   TEUCHOS_TEST_FOR_EXCEPTION(
00146     !epetra_map.get(), std::invalid_argument,
00147     "create_VectorSpace::initialize(...): Error!" );
00148 #endif // TEUCHOS_DEBUG
00149   RCP<const Teuchos::Comm<Ordinal> > comm =
00150     create_Comm(Teuchos::rcpFromRef(epetra_map->Comm())).assert_not_null();
00151   Teuchos::set_extra_data(epetra_map, "epetra_map", inoutArg(comm));
00152   const Ordinal localSubDim = epetra_map->NumMyElements();
00153   RCP<DefaultSpmdVectorSpace<double> > vs =
00154     defaultSpmdVectorSpace<double>(
00155       comm, localSubDim, epetra_map->NumGlobalElements(),
00156       !epetra_map->DistributedGlobal());
00157   TEUCHOS_ASSERT_EQUALITY(vs->dim(), as<Ordinal>(epetra_map->NumGlobalElements()));
00158   // NOTE: It is impossible to trigger the above exception unless
00159   // NumGlobalElemenets() overflows 'int'.  However, this is a nice sanity
00160   // check to stop the code early in case we seen an overflow in practice.
00161   // Because this assert will only likely trigger in a non-debug build, we
00162   // will unguard the assert since it is very cheap to perform.
00163   Teuchos::set_extra_data( epetra_map, "epetra_map", inoutArg(vs) );
00164   return vs;
00165 }
00166 
00167 
00168 Teuchos::RCP<const Thyra::VectorSpaceBase<double> >
00169 Thyra::create_LocallyReplicatedVectorSpace(
00170   const RCP<const VectorSpaceBase<double> > &parentSpace,
00171   const int dim
00172   )
00173 {
00174   using Teuchos::rcp_dynamic_cast;
00175 #ifdef TEUCHOS_DEBUG
00176   TEUCHOS_TEST_FOR_EXCEPT(parentSpace.get()==NULL);
00177   Teuchos::dyn_cast<const SpmdVectorSpaceBase<double> >(*parentSpace);
00178   TEUCHOS_TEST_FOR_EXCEPT(dim <= 0);
00179 #endif
00180   return parentSpace->smallVecSpcFcty()->createVecSpc(dim);
00181 }
00182 
00183 
00184 Teuchos::RCP<Thyra::VectorBase<double> >
00185 Thyra::create_Vector(
00186   const RCP<Epetra_Vector> &epetra_v,
00187   const RCP<const VectorSpaceBase<double> > &space_in
00188   )
00189 {
00190 #ifdef TEUCHOS_DEBUG
00191   TEUCHOS_TEST_FOR_EXCEPT(space_in.get()==NULL);
00192 #endif
00193   RCP<const SpmdVectorSpaceBase<double> >
00194     space = Teuchos::rcp_dynamic_cast<const SpmdVectorSpaceBase<double> >(
00195       space_in,true);
00196   if(!epetra_v.get())
00197     return Teuchos::null;
00198   // New local view of raw data
00199   double *localValues;
00200   epetra_v->ExtractView( &localValues );
00201   // Build the Vector
00202   RCP<SpmdVectorBase<double> >
00203     v = Teuchos::rcp(
00204       new DefaultSpmdVector<double>(
00205         space,
00206         Teuchos::arcp(localValues,0,epetra_v->Map().NumMyElements(),false),
00207         1
00208         )
00209       );
00210   Teuchos::set_extra_data<RCP<Epetra_Vector> >( epetra_v, "Epetra_Vector",
00211     Teuchos::inOutArg(v) );
00212   return v;
00213 }
00214 
00215 
00216 Teuchos::RCP<const Thyra::VectorBase<double> >
00217 Thyra::create_Vector(
00218   const RCP<const Epetra_Vector> &epetra_v,
00219   const RCP<const VectorSpaceBase<double> > &space_in
00220   )
00221 {
00222 #ifdef TEUCHOS_DEBUG
00223   TEUCHOS_TEST_FOR_EXCEPT(space_in.get()==NULL);
00224 #endif
00225   RCP<const SpmdVectorSpaceBase<double> >
00226     space = Teuchos::rcp_dynamic_cast<const SpmdVectorSpaceBase<double> >(
00227       space_in,true);
00228   if(!epetra_v.get())
00229     return Teuchos::null;
00230   // New local view of raw data
00231   double *localValues;
00232   epetra_v->ExtractView( &localValues );
00233   // Build the Vector
00234   RCP<const SpmdVectorBase<double> >
00235     v = Teuchos::rcp(
00236       new DefaultSpmdVector<double>(
00237         space,
00238         Teuchos::arcp(localValues,0,epetra_v->Map().NumMyElements(),false),
00239         1
00240         )
00241       );
00242   Teuchos::set_extra_data<RCP<const Epetra_Vector> >( epetra_v, "Epetra_Vector",
00243     Teuchos::inOutArg(v) );
00244   return v;
00245 }
00246 
00247 
00248 Teuchos::RCP<Thyra::MultiVectorBase<double> >
00249 Thyra::create_MultiVector(
00250   const RCP<Epetra_MultiVector> &epetra_mv,
00251   const RCP<const VectorSpaceBase<double> > &range_in,
00252   const RCP<const VectorSpaceBase<double> > &domain_in
00253   )
00254 {
00255   using Teuchos::rcp_dynamic_cast;
00256 #ifdef TEUCHOS_DEBUG
00257   TEUCHOS_TEST_FOR_EXCEPT(range_in.get()==NULL);
00258 #endif
00259   const RCP<const SpmdVectorSpaceBase<double> > range =
00260     Teuchos::rcp_dynamic_cast<const SpmdVectorSpaceBase<double> >(
00261       unwrapSingleProductVectorSpace(range_in),
00262       true
00263       );
00264   RCP<const ScalarProdVectorSpaceBase<double> > domain =
00265     Teuchos::rcp_dynamic_cast<const ScalarProdVectorSpaceBase<double> >(
00266       unwrapSingleProductVectorSpace(domain_in),
00267       true
00268       );
00269   if (!epetra_mv.get() )
00270     return Teuchos::null;
00271   if ( is_null(domain) ) {
00272     domain = rcp_dynamic_cast<const ScalarProdVectorSpaceBase<double> >(
00273       create_LocallyReplicatedVectorSpace(range,epetra_mv->NumVectors())
00274       );
00275   }
00276   // New local view of raw data
00277   double *localValues; int leadingDim;
00278   if( epetra_mv->ConstantStride() ) {
00279     epetra_mv->ExtractView( &localValues, &leadingDim );
00280   }
00281   else {
00282     TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement views of non-contiguous mult-vectors!
00283   }
00284   // Build the MultiVector
00285   RCP<SpmdMultiVectorBase<double> >
00286     mv = Teuchos::rcp(
00287       new DefaultSpmdMultiVector<double>(
00288         range,
00289         domain,
00290         Teuchos::arcp(localValues,0,leadingDim*epetra_mv->NumVectors(),false),
00291         leadingDim
00292         )
00293       );
00294   Teuchos::set_extra_data<RCP<Epetra_MultiVector> >(
00295     epetra_mv, "Epetra_MultiVector", Teuchos::inOutArg(mv) );
00296   return mv;
00297 }
00298 
00299 
00300 Teuchos::RCP<const Thyra::MultiVectorBase<double> >
00301 Thyra::create_MultiVector(
00302   const RCP<const Epetra_MultiVector> &epetra_mv,
00303   const RCP<const VectorSpaceBase<double> > &range_in,
00304   const RCP<const VectorSpaceBase<double> > &domain_in
00305   )
00306 {
00307   using Teuchos::rcp_dynamic_cast;
00308 #ifdef TEUCHOS_DEBUG
00309   TEUCHOS_TEST_FOR_EXCEPT(range_in.get()==NULL);
00310 #endif
00311   const RCP<const SpmdVectorSpaceBase<double> > range =
00312     Teuchos::rcp_dynamic_cast<const SpmdVectorSpaceBase<double> >(
00313       unwrapSingleProductVectorSpace(range_in),
00314       true
00315       );
00316   RCP<const ScalarProdVectorSpaceBase<double> > domain =
00317     Teuchos::rcp_dynamic_cast<const ScalarProdVectorSpaceBase<double> >(
00318       unwrapSingleProductVectorSpace(domain_in),
00319       true
00320       );
00321   if (!epetra_mv.get())
00322     return Teuchos::null;
00323   if ( is_null(domain) ) {
00324     domain = rcp_dynamic_cast<const ScalarProdVectorSpaceBase<double> >(
00325       create_LocallyReplicatedVectorSpace(range,epetra_mv->NumVectors())
00326       );
00327   }
00328   // New local view of raw data
00329   double *localValues; int leadingDim;
00330   if( epetra_mv->ConstantStride() ) {
00331     epetra_mv->ExtractView( &localValues, &leadingDim );
00332   }
00333   else {
00334     TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement views of non-contiguous mult-vectors!
00335   }
00336   // Build the MultiVector
00337   RCP<const SpmdMultiVectorBase<double> >
00338     mv = Teuchos::rcp(
00339       new DefaultSpmdMultiVector<double>(
00340         range,
00341         domain,
00342         Teuchos::arcp(localValues,0,leadingDim*epetra_mv->NumVectors(),false),
00343         leadingDim
00344         )
00345       );
00346   Teuchos::set_extra_data<RCP<const Epetra_MultiVector> >(
00347     epetra_mv, "Epetra_MultiVector", Teuchos::inOutArg(mv) );
00348   return mv;
00349 }
00350 
00351 
00352 Teuchos::RCP<const Epetra_Comm>
00353 Thyra::get_Epetra_Comm(const Teuchos::Comm<Ordinal>& comm_in)
00354 {
00355 
00356   using Teuchos::rcp;
00357   using Teuchos::ptrFromRef;
00358   using Teuchos::ptr_dynamic_cast;
00359   using Teuchos::SerialComm;
00360 #ifdef HAVE_MPI
00361   using Teuchos::MpiComm;
00362 #endif
00363 
00364   const Ptr<const Teuchos::Comm<Ordinal> > comm = Teuchos::ptrFromRef(comm_in);
00365 
00366   const Ptr<const SerialComm<Ordinal> > serialComm =
00367     ptr_dynamic_cast<const SerialComm<Ordinal> >(comm);
00368 
00369   RCP<const Epetra_Comm> epetraComm;
00370 
00371 #ifdef HAVE_MPI
00372 
00373   const Ptr<const MpiComm<Ordinal> > mpiComm =
00374     ptr_dynamic_cast<const MpiComm<Ordinal> >(comm);
00375 
00376   TEUCHOS_TEST_FOR_EXCEPTION(is_null(mpiComm) && is_null(serialComm),
00377     std::runtime_error, 
00378     "SPMD std::vector space has a communicator that is "
00379     "neither a serial comm nor an MPI comm");
00380 
00381   if (nonnull(mpiComm)) {
00382     epetraComm = rcp(new Epetra_MpiComm(*mpiComm->getRawMpiComm()()));
00383   }
00384   else {
00385     epetraComm = rcp(new Epetra_SerialComm());
00386   }
00387 
00388 #else
00389 
00390   TEUCHOS_TEST_FOR_EXCEPTION(is_null(serialComm), std::runtime_error, 
00391     "SPMD std::vector space has a communicator that is "
00392     "neither a serial comm nor an MPI comm");
00393 
00394   epetraComm = rcp(new Epetra_SerialComm());
00395   
00396 #endif
00397   
00398   TEUCHOS_TEST_FOR_EXCEPTION(is_null(epetraComm), std::runtime_error,
00399     "null communicator created");
00400   
00401   return epetraComm;
00402 
00403 }
00404 
00405 
00406 Teuchos::RCP<const Epetra_Map>
00407 Thyra::get_Epetra_Map(const VectorSpaceBase<double>& vs_in,
00408   const RCP<const Epetra_Comm>& comm)
00409 {
00410 
00411   using Teuchos::rcpFromRef;
00412   using Teuchos::rcpFromPtr;
00413   using Teuchos::rcp_dynamic_cast;
00414   using Teuchos::ptrFromRef;
00415   using Teuchos::ptr_dynamic_cast;
00416 
00417   const Ptr<const VectorSpaceBase<double> > vs_ptr = ptrFromRef(vs_in);
00418   
00419   const Ptr<const SpmdVectorSpaceBase<double> > spmd_vs =
00420     ptr_dynamic_cast<const SpmdVectorSpaceBase<double> >(vs_ptr);
00421 
00422   const Ptr<const ProductVectorSpaceBase<double> > &prod_vs = 
00423     ptr_dynamic_cast<const ProductVectorSpaceBase<double> >(vs_ptr);
00424 
00425   TEUCHOS_TEST_FOR_EXCEPTION( is_null(spmd_vs) && is_null(prod_vs), std::logic_error,
00426     "Error, the concrete VectorSpaceBase object of type "
00427     +Teuchos::demangleName(typeid(vs_in).name())+" does not support the"
00428     " SpmdVectorSpaceBase or the ProductVectorSpaceBase interfaces!" );
00429 
00430   const int numBlocks = (nonnull(prod_vs) ? prod_vs->numBlocks() : 1);
00431 
00432   // Get an array of SpmdVectorBase objects for the blocks
00433   
00434   Array<RCP<const SpmdVectorSpaceBase<double> > > spmd_vs_blocks;
00435   if (nonnull(prod_vs)) {
00436     for (int block_i = 0; block_i < numBlocks; ++block_i) {
00437       const RCP<const SpmdVectorSpaceBase<double> > spmd_vs_i =
00438         rcp_dynamic_cast<const SpmdVectorSpaceBase<double> >(
00439           prod_vs->getBlock(block_i), true);
00440       spmd_vs_blocks.push_back(spmd_vs_i);
00441     }
00442   }
00443   else {
00444     spmd_vs_blocks.push_back(rcpFromPtr(spmd_vs));
00445   }
00446   
00447   // Find the number of local elements, summed over all blocks
00448 
00449   int myLocalElements = 0;
00450   for (int block_i = 0; block_i < numBlocks; ++block_i) {
00451     myLocalElements += spmd_vs_blocks[block_i]->localSubDim();
00452   }
00453   
00454   // Find the GIDs owned by this processor, taken from all blocks
00455   
00456   int count=0;
00457   int blockOffset = 0;
00458   Array<int> myGIDs(myLocalElements);
00459   for (int block_i = 0; block_i < numBlocks; ++block_i) {
00460     const RCP<const SpmdVectorSpaceBase<double> > spmd_vs_i = spmd_vs_blocks[block_i];
00461     const int lowGIDInBlock = spmd_vs_i->localOffset();
00462     const int numLocalElementsInBlock = spmd_vs_i->localSubDim();
00463     for (int i=0; i < numLocalElementsInBlock; ++i, ++count) {
00464       myGIDs[count] = blockOffset + lowGIDInBlock + i;
00465     }
00466     blockOffset += spmd_vs_i->dim();
00467   }
00468   
00469   const int globalDim = vs_in.dim();
00470 
00471   return Teuchos::rcp(
00472     new Epetra_Map(globalDim, myLocalElements, &(myGIDs[0]), 0, *comm));
00473 
00474 }
00475 
00476 
00477 Teuchos::RCP<Epetra_Vector>
00478 Thyra::get_Epetra_Vector(
00479   const Epetra_Map &map,
00480   const RCP<VectorBase<double> > &v
00481   )
00482 {
00483   using Teuchos::get_optional_extra_data;
00484 #ifdef TEUCHOS_DEBUG
00485   RCP<const VectorSpaceBase<double> >
00486     epetra_vs = create_VectorSpace(Teuchos::rcp(&map,false));
00487   THYRA_ASSERT_VEC_SPACES( 
00488     "Thyra::get_Epetra_Vector(map,v)", *epetra_vs, *v->space() );
00489 #endif
00490   //
00491   // First, try to grab the Epetra_Vector straight out of the
00492   // RCP since this is the fastest way.
00493   //
00494   const Ptr<const RCP<Epetra_Vector> >
00495     epetra_v_ptr = get_optional_extra_data<RCP<Epetra_Vector> >(
00496       v,"Epetra_Vector");
00497   if(!is_null(epetra_v_ptr)) {
00498     return *epetra_v_ptr;
00499   }
00500   //
00501   // The assumption that we (rightly) make here is that if the vector spaces
00502   // are compatible, that either the multi-vectors are both in-core or the
00503   // vector spaces are both derived from SpmdVectorSpaceBase and have
00504   // compatible maps.
00505   // 
00506   const VectorSpaceBase<double>  &vs = *v->range();
00507   const SpmdVectorSpaceBase<double> *mpi_vs
00508     = dynamic_cast<const SpmdVectorSpaceBase<double>*>(&vs);
00509   const Ordinal localOffset = ( mpi_vs ? mpi_vs->localOffset() : 0 );
00510   const Ordinal localSubDim = ( mpi_vs ? mpi_vs->localSubDim() : vs.dim() );
00511   //
00512   // Here we will extract a view of the local elements in the underlying
00513   // VectorBase object.  In most cases, no data will be allocated or copied
00514   // and only some small objects will be created and a few virtual functions
00515   // will be called so the overhead should be low and fixed.
00516   //
00517   // Create a *mutable* view of the local elements, this view will be set on
00518   // the RCP that is returned.  As a result, this view will be relased
00519   // when the returned Epetra_Vector is released.
00520   //
00521   // Note that the input vector 'v' will be remembered through this detached
00522   // view!
00523   //
00524   RCP<DetachedVectorView<double> >
00525     emvv = Teuchos::rcp(
00526       new DetachedVectorView<double>(
00527         v
00528         ,Range1D(localOffset,localOffset+localSubDim-1)
00529         ,true // forceContiguous
00530         )
00531       );
00532   // Create a temporary Epetra_Vector object and give it
00533   // the above local view
00534   RCP<Epetra_Vector>
00535     epetra_v = Teuchos::rcp(
00536       new Epetra_Vector(
00537         ::View                                 // CV
00538         ,map                                   // Map
00539         ,const_cast<double*>(emvv->values())   // V
00540         )
00541       );
00542   // Give the explict view object to the above Epetra_Vector smart pointer
00543   // object.  In this way, when the client is finished with the Epetra_Vector
00544   // view the destructor from the object in emvv will automatically commit the
00545   // changes to the elements in the input v VectorBase object (reguardless of
00546   // its implementation).  This is truly an elegant result!
00547   Teuchos::set_extra_data( emvv, "emvv", Teuchos::inOutArg(epetra_v),
00548     Teuchos::PRE_DESTROY );
00549   // We are done!
00550   return epetra_v;
00551 }
00552 
00553 
00554 Teuchos::RCP<const Epetra_Vector>
00555 Thyra::get_Epetra_Vector(
00556   const Epetra_Map &map, 
00557   const RCP<const VectorBase<double> > &v
00558   )
00559 {
00560   using Teuchos::get_optional_extra_data;
00561 #ifdef TEUCHOS_DEBUG
00562   RCP<const VectorSpaceBase<double> >
00563     epetra_vs = create_VectorSpace(Teuchos::rcp(&map,false));
00564   THYRA_ASSERT_VEC_SPACES(
00565     "Thyra::get_Epetra_Vector(map,v)", *epetra_vs, *v->space() );
00566 #endif
00567   //
00568   // First, try to grab the Epetra_Vector straight out of the
00569   // RCP since this is the fastest way.
00570   //
00571   const Ptr<const RCP<const Epetra_Vector> >
00572     epetra_v_ptr = get_optional_extra_data<RCP<const Epetra_Vector> >(
00573       v,"Epetra_Vector");
00574   if(!is_null(epetra_v_ptr))
00575     return *epetra_v_ptr;
00576   const Ptr<const RCP<Epetra_Vector> >
00577     epetra_nonconst_v_ptr = get_optional_extra_data<RCP<Epetra_Vector> >(
00578       v,"Epetra_Vector");
00579   if(!is_null(epetra_nonconst_v_ptr))
00580     return *epetra_nonconst_v_ptr;
00581   //
00582   // Same as above function except as stated below
00583   //
00584   const VectorSpaceBase<double> &vs = *v->range();
00585   const SpmdVectorSpaceBase<double> *mpi_vs
00586     = dynamic_cast<const SpmdVectorSpaceBase<double>*>(&vs);
00587   const Ordinal localOffset = ( mpi_vs ? mpi_vs->localOffset() : 0 );
00588   const Ordinal localSubDim = ( mpi_vs ? mpi_vs->localSubDim() : vs.dim() );
00589   // Get an explicit *non-mutable* view of all of the elements in the multi
00590   // vector.  Note that 'v' will be remembered by this view!
00591   RCP<ConstDetachedVectorView<double> >
00592     evv = Teuchos::rcp(
00593       new ConstDetachedVectorView<double>(
00594         v
00595         ,Range1D(localOffset,localOffset+localSubDim-1)
00596         ,true // forceContiguous
00597         )
00598       );
00599   // Create a temporary Epetra_Vector object and give it
00600   // the above view
00601   RCP<Epetra_Vector>
00602     epetra_v = Teuchos::rcp(
00603       new Epetra_Vector(
00604         ::View                                  // CV
00605         ,map                                    // Map
00606         ,const_cast<double*>(evv->values())     // V
00607         )
00608       );
00609   // This next statement will cause the destructor to free the view if
00610   // needed (see above function).  Since this view is non-mutable,
00611   // only a releaseDetachedView(...) and not a commit will be called.
00612   // This is the whole reason there is a seperate implementation for
00613   // the const and non-const cases.
00614   Teuchos::set_extra_data( evv, "evv", Teuchos::inOutArg(epetra_v),
00615     Teuchos::PRE_DESTROY );
00616   // We are done!
00617   return epetra_v;
00618 }
00619 
00620 
00621 Teuchos::RCP<Epetra_MultiVector>
00622 Thyra::get_Epetra_MultiVector(
00623   const Epetra_Map &map,
00624   const RCP<MultiVectorBase<double> > &mv
00625   )
00626 {
00627   using Teuchos::get_optional_extra_data;
00628 #ifdef TEUCHOS_DEBUG
00629   RCP<const VectorSpaceBase<double> >
00630     epetra_vs = create_VectorSpace(Teuchos::rcp(&map,false));
00631   THYRA_ASSERT_VEC_SPACES(
00632     "Thyra::get_Epetra_MultiVector(map,mv)", *epetra_vs, *mv->range() );
00633 #endif
00634   //
00635   // First, try to grab the Epetra_MultiVector straight out of the
00636   // RCP since this is the fastest way.
00637   //
00638   const Ptr<const RCP<Epetra_MultiVector> >
00639     epetra_mv_ptr = get_optional_extra_data<RCP<Epetra_MultiVector> >(
00640       mv,"Epetra_MultiVector");
00641   if(!is_null(epetra_mv_ptr)) {
00642     return *epetra_mv_ptr;
00643   }
00644   //
00645   // The assumption that we (rightly) make here is that if the vector spaces
00646   // are compatible, that either the multi-vectors are both in-core or the
00647   // vector spaces are both derived from SpmdVectorSpaceBase and have
00648   // compatible maps.
00649   // 
00650   const VectorSpaceBase<double> &vs = *mv->range();
00651   const SpmdVectorSpaceBase<double> *mpi_vs
00652     = dynamic_cast<const SpmdVectorSpaceBase<double>*>(&vs);
00653   const Ordinal localOffset = ( mpi_vs ? mpi_vs->localOffset() : 0 );
00654   const Ordinal localSubDim = ( mpi_vs ? mpi_vs->localSubDim() : vs.dim() );
00655   //
00656   // Here we will extract a view of the local elements in the underlying
00657   // MultiVectorBase object.  In most cases, no data will be allocated or
00658   // copied and only some small objects will be created and a few virtual
00659   // functions will be called so the overhead should be low and fixed.
00660   //
00661   // Create a *mutable* view of the local elements, this view will be set on
00662   // the RCP that is returned.  As a result, this view will be relased
00663   // when the returned Epetra_MultiVector is released.
00664   //
00665   RCP<DetachedMultiVectorView<double> >
00666     emmvv = Teuchos::rcp(
00667       new DetachedMultiVectorView<double>(
00668         *mv
00669         ,Range1D(localOffset,localOffset+localSubDim-1)
00670         )
00671       );
00672   // Create a temporary Epetra_MultiVector object and give it
00673   // the above local view
00674   RCP<Epetra_MultiVector>
00675     epetra_mv = Teuchos::rcp(
00676       new Epetra_MultiVector(
00677         ::View                                  // CV
00678         ,map                                    // Map
00679         ,const_cast<double*>(emmvv->values())   // A
00680         ,emmvv->leadingDim()                    // MyLDA
00681         ,emmvv->numSubCols()                    // NumVectors
00682         )
00683       );
00684   // Give the explict view object to the above Epetra_MultiVector
00685   // smart pointer object.  In this way, when the client is finished
00686   // with the Epetra_MultiVector view the destructor from the object
00687   // in emmvv will automatically commit the changes to the elements in
00688   // the input mv MultiVectorBase object (reguardless of its
00689   // implementation).  This is truly an elegant result!
00690   Teuchos::set_extra_data( emmvv, "emmvv", Teuchos::inOutArg(epetra_mv),
00691     Teuchos::PRE_DESTROY );
00692   // Also set the mv itself as extra data just to be safe
00693   Teuchos::set_extra_data( mv, "mv", Teuchos::inOutArg(epetra_mv) );
00694   // We are done!
00695   return epetra_mv;
00696 }
00697 
00698 
00699 Teuchos::RCP<const Epetra_MultiVector>
00700 Thyra::get_Epetra_MultiVector(
00701   const Epetra_Map &map,
00702   const RCP<const MultiVectorBase<double> > &mv
00703   )
00704 {
00705   using Teuchos::get_optional_extra_data;
00706 #ifdef TEUCHOS_DEBUG
00707   RCP<const VectorSpaceBase<double> >
00708     epetra_vs = create_VectorSpace(Teuchos::rcp(&map,false));
00709   THYRA_ASSERT_VEC_SPACES(
00710     "Thyra::get_Epetra_MultiVector(map,mv)",
00711     *epetra_vs, *mv->range() );
00712 #endif
00713   //
00714   // First, try to grab the Epetra_MultiVector straight out of the
00715   // RCP since this is the fastest way.
00716   //
00717   const Ptr<const RCP<const Epetra_MultiVector> >
00718     epetra_mv_ptr
00719     = get_optional_extra_data<RCP<const Epetra_MultiVector> >(
00720       mv,"Epetra_MultiVector" );
00721   if(!is_null(epetra_mv_ptr)) {
00722     return *epetra_mv_ptr;
00723   }
00724   //
00725   // Same as above function except as stated below
00726   //
00727   const VectorSpaceBase<double> &vs = *mv->range();
00728   const SpmdVectorSpaceBase<double> *mpi_vs
00729     = dynamic_cast<const SpmdVectorSpaceBase<double>*>(&vs);
00730   const Ordinal localOffset = ( mpi_vs ? mpi_vs->localOffset() : 0 );
00731   const Ordinal localSubDim = ( mpi_vs ? mpi_vs->localSubDim() : vs.dim() );
00732   // Get an explicit *non-mutable* view of all of the elements in
00733   // the multi vector.
00734   RCP<ConstDetachedMultiVectorView<double> >
00735     emvv = Teuchos::rcp(
00736       new ConstDetachedMultiVectorView<double>(
00737         *mv
00738         ,Range1D(localOffset,localOffset+localSubDim-1)
00739         )
00740       );
00741   // Create a temporary Epetra_MultiVector object and give it
00742   // the above view
00743   RCP<Epetra_MultiVector>
00744     epetra_mv = Teuchos::rcp(
00745       new Epetra_MultiVector(
00746         ::View                                  // CV
00747         ,map                                    // Map
00748         ,const_cast<double*>(emvv->values())    // A
00749         ,emvv->leadingDim()                     // MyLDA
00750         ,emvv->numSubCols()                     // NumVectors
00751         )
00752       );
00753   // This next statement will cause the destructor to free the view if
00754   // needed (see above function).  Since this view is non-mutable,
00755   // only a releaseDetachedView(...) and not a commit will be called.
00756   // This is the whole reason there is a seperate implementation for
00757   // the const and non-const cases.
00758   Teuchos::set_extra_data( emvv, "emvv", Teuchos::inOutArg(epetra_mv),
00759     Teuchos::PRE_DESTROY );
00760   // Also set the mv itself as extra data just to be safe
00761   Teuchos::set_extra_data( mv, "mv", Teuchos::inOutArg(epetra_mv) );
00762   // We are done!
00763   return epetra_mv;
00764 }
00765 
00766 
00767 Teuchos::RCP<Epetra_MultiVector>
00768 Thyra::get_Epetra_MultiVector(
00769   const Epetra_Map &map,
00770   MultiVectorBase<double> &mv
00771   )
00772 {
00773   using Teuchos::rcpWithEmbeddedObj;
00774   using Teuchos::rcpFromRef;
00775   using Teuchos::ptrFromRef;
00776   using Teuchos::ptr_dynamic_cast;
00777   using Teuchos::outArg;
00778   ArrayRCP<double> mvData;
00779   Ordinal mvLeadingDim = -1;
00780   Ptr<SpmdMultiVectorBase<double> > mvSpmdMv;
00781   if (nonnull(mvSpmdMv = ptr_dynamic_cast<SpmdMultiVectorBase<double> >(ptrFromRef(mv)))) {
00782     mvSpmdMv->getNonconstLocalData(outArg(mvData), outArg(mvLeadingDim));
00783   }
00784   if (nonnull(mvData)) {
00785     return rcpWithEmbeddedObj(
00786       new Epetra_MultiVector(
00787         ::View, map, mvData.getRawPtr(), mvLeadingDim, mv.domain()->dim()
00788         ),
00789       mvData
00790       );
00791   }
00792   return ::Thyra::get_Epetra_MultiVector(map, rcpFromRef(mv));
00793 }
00794 
00795 
00796 Teuchos::RCP<const Epetra_MultiVector>
00797 Thyra::get_Epetra_MultiVector(
00798   const Epetra_Map &map,
00799   const MultiVectorBase<double> &mv
00800   )
00801 {
00802   using Teuchos::rcpWithEmbeddedObj;
00803   using Teuchos::rcpFromRef;
00804   using Teuchos::ptrFromRef;
00805   using Teuchos::ptr_dynamic_cast;
00806   using Teuchos::outArg;
00807   ArrayRCP<const double> mvData;
00808   Ordinal mvLeadingDim = -1;
00809   Ptr<const SpmdMultiVectorBase<double> > mvSpmdMv;
00810   if (nonnull(mvSpmdMv = ptr_dynamic_cast<const SpmdMultiVectorBase<double> >(ptrFromRef(mv)))) {
00811     mvSpmdMv->getLocalData(outArg(mvData), outArg(mvLeadingDim));
00812   }
00813   if (nonnull(mvData)) {
00814     return rcpWithEmbeddedObj(
00815       new Epetra_MultiVector(
00816         ::View,map, const_cast<double*>(mvData.getRawPtr()), mvLeadingDim, mv.domain()->dim()
00817         ),
00818       mvData
00819       );
00820   }
00821   return ::Thyra::get_Epetra_MultiVector(map, rcpFromRef(mv));
00822 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines