Thyra Package Browser (Single Doxygen Collection) Version of the Day
Thyra_EpetraThyraWrappers.cpp
Go to the documentation of this file.
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 #ifdef TEUCHOS_DEBUG
00144   TEUCHOS_TEST_FOR_EXCEPTION(
00145     !epetra_map.get(), std::invalid_argument,
00146     "create_VectorSpace::initialize(...): Error!" );
00147 #endif // TEUCHOS_DEBUG
00148   RCP<const Teuchos::Comm<Ordinal> >
00149     comm = create_Comm(Teuchos::rcp(&epetra_map->Comm(),false)).assert_not_null();
00150   Teuchos::set_extra_data( epetra_map, "epetra_map", Teuchos::inOutArg(comm) );
00151   const Ordinal localSubDim = epetra_map->NumMyElements();
00152   RCP<DefaultSpmdVectorSpace<double> > vs =
00153     defaultSpmdVectorSpace<double>(
00154       comm, localSubDim, epetra_map->NumGlobalElements());
00155 #ifndef TEUCHOS_DEBUG
00156   TEUCHOS_TEST_FOR_EXCEPTION(
00157     vs->dim() != epetra_map->NumGlobalElements(), std::logic_error
00158     ,"create_VectorSpace(...): Error, vs->dim() = "<<vs->dim()<<" != "
00159     "epetra_map->NumGlobalElements() = "<<epetra_map->NumGlobalElements()<<"!"
00160     );
00161 #endif    
00162   Teuchos::set_extra_data( epetra_map, "epetra_map", Teuchos::inOutArg(vs) );
00163   return vs;
00164 }
00165 
00166 
00167 Teuchos::RCP<const Thyra::VectorSpaceBase<double> >
00168 Thyra::create_LocallyReplicatedVectorSpace(
00169   const RCP<const VectorSpaceBase<double> > &parentSpace,
00170   const int dim
00171   )
00172 {
00173   using Teuchos::rcp_dynamic_cast;
00174 #ifdef TEUCHOS_DEBUG
00175   TEUCHOS_TEST_FOR_EXCEPT(parentSpace.get()==NULL);
00176   Teuchos::dyn_cast<const SpmdVectorSpaceBase<double> >(*parentSpace);
00177   TEUCHOS_TEST_FOR_EXCEPT(dim <= 0);
00178 #endif
00179   return parentSpace->smallVecSpcFcty()->createVecSpc(dim);
00180 }
00181 
00182 
00183 Teuchos::RCP<Thyra::VectorBase<double> >
00184 Thyra::create_Vector(
00185   const RCP<Epetra_Vector> &epetra_v,
00186   const RCP<const VectorSpaceBase<double> > &space_in
00187   )
00188 {
00189 #ifdef TEUCHOS_DEBUG
00190   TEUCHOS_TEST_FOR_EXCEPT(space_in.get()==NULL);
00191 #endif
00192   RCP<const SpmdVectorSpaceBase<double> >
00193     space = Teuchos::rcp_dynamic_cast<const SpmdVectorSpaceBase<double> >(
00194       space_in,true);
00195   if(!epetra_v.get())
00196     return Teuchos::null;
00197   // New local view of raw data
00198   double *localValues;
00199   epetra_v->ExtractView( &localValues );
00200   // Build the Vector
00201   RCP<SpmdVectorBase<double> >
00202     v = Teuchos::rcp(
00203       new DefaultSpmdVector<double>(
00204         space,
00205         Teuchos::arcp(localValues,0,epetra_v->Map().NumMyElements(),false),
00206         1
00207         )
00208       );
00209   Teuchos::set_extra_data<RCP<Epetra_Vector> >( epetra_v, "Epetra_Vector",
00210     Teuchos::inOutArg(v) );
00211   return v;
00212 }
00213 
00214 
00215 Teuchos::RCP<const Thyra::VectorBase<double> >
00216 Thyra::create_Vector(
00217   const RCP<const Epetra_Vector> &epetra_v,
00218   const RCP<const VectorSpaceBase<double> > &space_in
00219   )
00220 {
00221 #ifdef TEUCHOS_DEBUG
00222   TEUCHOS_TEST_FOR_EXCEPT(space_in.get()==NULL);
00223 #endif
00224   RCP<const SpmdVectorSpaceBase<double> >
00225     space = Teuchos::rcp_dynamic_cast<const SpmdVectorSpaceBase<double> >(
00226       space_in,true);
00227   if(!epetra_v.get())
00228     return Teuchos::null;
00229   // New local view of raw data
00230   double *localValues;
00231   epetra_v->ExtractView( &localValues );
00232   // Build the Vector
00233   RCP<const SpmdVectorBase<double> >
00234     v = Teuchos::rcp(
00235       new DefaultSpmdVector<double>(
00236         space,
00237         Teuchos::arcp(localValues,0,epetra_v->Map().NumMyElements(),false),
00238         1
00239         )
00240       );
00241   Teuchos::set_extra_data<RCP<const Epetra_Vector> >( epetra_v, "Epetra_Vector",
00242     Teuchos::inOutArg(v) );
00243   return v;
00244 }
00245 
00246 
00247 Teuchos::RCP<Thyra::MultiVectorBase<double> >
00248 Thyra::create_MultiVector(
00249   const RCP<Epetra_MultiVector> &epetra_mv,
00250   const RCP<const VectorSpaceBase<double> > &range_in,
00251   const RCP<const VectorSpaceBase<double> > &domain_in
00252   )
00253 {
00254   using Teuchos::rcp_dynamic_cast;
00255 #ifdef TEUCHOS_DEBUG
00256   TEUCHOS_TEST_FOR_EXCEPT(range_in.get()==NULL);
00257 #endif
00258   const RCP<const SpmdVectorSpaceBase<double> > range =
00259     Teuchos::rcp_dynamic_cast<const SpmdVectorSpaceBase<double> >(
00260       unwrapSingleProductVectorSpace(range_in),
00261       true
00262       );
00263   RCP<const ScalarProdVectorSpaceBase<double> > domain =
00264     Teuchos::rcp_dynamic_cast<const ScalarProdVectorSpaceBase<double> >(
00265       unwrapSingleProductVectorSpace(domain_in),
00266       true
00267       );
00268   if (!epetra_mv.get() )
00269     return Teuchos::null;
00270   if ( is_null(domain) ) {
00271     domain = rcp_dynamic_cast<const ScalarProdVectorSpaceBase<double> >(
00272       create_LocallyReplicatedVectorSpace(range,epetra_mv->NumVectors())
00273       );
00274   }
00275   // New local view of raw data
00276   double *localValues; int leadingDim;
00277   if( epetra_mv->ConstantStride() ) {
00278     epetra_mv->ExtractView( &localValues, &leadingDim );
00279   }
00280   else {
00281     TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement views of non-contiguous mult-vectors!
00282   }
00283   // Build the MultiVector
00284   RCP<SpmdMultiVectorBase<double> >
00285     mv = Teuchos::rcp(
00286       new DefaultSpmdMultiVector<double>(
00287         range,
00288         domain,
00289         Teuchos::arcp(localValues,0,leadingDim*epetra_mv->NumVectors(),false),
00290         leadingDim
00291         )
00292       );
00293   Teuchos::set_extra_data<RCP<Epetra_MultiVector> >(
00294     epetra_mv, "Epetra_MultiVector", Teuchos::inOutArg(mv) );
00295   return mv;
00296 }
00297 
00298 
00299 Teuchos::RCP<const Thyra::MultiVectorBase<double> >
00300 Thyra::create_MultiVector(
00301   const RCP<const Epetra_MultiVector> &epetra_mv,
00302   const RCP<const VectorSpaceBase<double> > &range_in,
00303   const RCP<const VectorSpaceBase<double> > &domain_in
00304   )
00305 {
00306   using Teuchos::rcp_dynamic_cast;
00307 #ifdef TEUCHOS_DEBUG
00308   TEUCHOS_TEST_FOR_EXCEPT(range_in.get()==NULL);
00309 #endif
00310   const RCP<const SpmdVectorSpaceBase<double> > range =
00311     Teuchos::rcp_dynamic_cast<const SpmdVectorSpaceBase<double> >(
00312       unwrapSingleProductVectorSpace(range_in),
00313       true
00314       );
00315   RCP<const ScalarProdVectorSpaceBase<double> > domain =
00316     Teuchos::rcp_dynamic_cast<const ScalarProdVectorSpaceBase<double> >(
00317       unwrapSingleProductVectorSpace(domain_in),
00318       true
00319       );
00320   if (!epetra_mv.get())
00321     return Teuchos::null;
00322   if ( is_null(domain) ) {
00323     domain = rcp_dynamic_cast<const ScalarProdVectorSpaceBase<double> >(
00324       create_LocallyReplicatedVectorSpace(range,epetra_mv->NumVectors())
00325       );
00326   }
00327   // New local view of raw data
00328   double *localValues; int leadingDim;
00329   if( epetra_mv->ConstantStride() ) {
00330     epetra_mv->ExtractView( &localValues, &leadingDim );
00331   }
00332   else {
00333     TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement views of non-contiguous mult-vectors!
00334   }
00335   // Build the MultiVector
00336   RCP<const SpmdMultiVectorBase<double> >
00337     mv = Teuchos::rcp(
00338       new DefaultSpmdMultiVector<double>(
00339         range,
00340         domain,
00341         Teuchos::arcp(localValues,0,leadingDim*epetra_mv->NumVectors(),false),
00342         leadingDim
00343         )
00344       );
00345   Teuchos::set_extra_data<RCP<const Epetra_MultiVector> >(
00346     epetra_mv, "Epetra_MultiVector", Teuchos::inOutArg(mv) );
00347   return mv;
00348 }
00349 
00350 
00351 Teuchos::RCP<const Epetra_Comm>
00352 Thyra::get_Epetra_Comm(const Teuchos::Comm<Ordinal>& comm_in)
00353 {
00354 
00355   using Teuchos::rcp;
00356   using Teuchos::ptrFromRef;
00357   using Teuchos::ptr_dynamic_cast;
00358   using Teuchos::SerialComm;
00359 #ifdef HAVE_MPI
00360   using Teuchos::MpiComm;
00361 #endif
00362 
00363   const Ptr<const Teuchos::Comm<Ordinal> > comm = Teuchos::ptrFromRef(comm_in);
00364 
00365   const Ptr<const SerialComm<Ordinal> > serialComm =
00366     ptr_dynamic_cast<const SerialComm<Ordinal> >(comm);
00367 
00368   RCP<const Epetra_Comm> epetraComm;
00369 
00370 #ifdef HAVE_MPI
00371 
00372   const Ptr<const MpiComm<Ordinal> > mpiComm =
00373     ptr_dynamic_cast<const MpiComm<Ordinal> >(comm);
00374 
00375   TEUCHOS_TEST_FOR_EXCEPTION(is_null(mpiComm) && is_null(serialComm),
00376     std::runtime_error, 
00377     "SPMD std::vector space has a communicator that is "
00378     "neither a serial comm nor an MPI comm");
00379 
00380   if (nonnull(mpiComm)) {
00381     epetraComm = rcp(new Epetra_MpiComm(*mpiComm->getRawMpiComm()()));
00382   }
00383   else {
00384     epetraComm = rcp(new Epetra_SerialComm());
00385   }
00386 
00387 #else
00388 
00389   TEUCHOS_TEST_FOR_EXCEPTION(is_null(serialComm), std::runtime_error, 
00390     "SPMD std::vector space has a communicator that is "
00391     "neither a serial comm nor an MPI comm");
00392 
00393   epetraComm = rcp(new Epetra_SerialComm());
00394   
00395 #endif
00396   
00397   TEUCHOS_TEST_FOR_EXCEPTION(is_null(epetraComm), std::runtime_error,
00398     "null communicator created");
00399   
00400   return epetraComm;
00401 
00402 }
00403 
00404 
00405 Teuchos::RCP<const Epetra_Map>
00406 Thyra::get_Epetra_Map(const VectorSpaceBase<double>& vs_in,
00407   const RCP<const Epetra_Comm>& comm)
00408 {
00409 
00410   using Teuchos::rcpFromRef;
00411   using Teuchos::rcpFromPtr;
00412   using Teuchos::rcp_dynamic_cast;
00413   using Teuchos::ptrFromRef;
00414   using Teuchos::ptr_dynamic_cast;
00415 
00416   const Ptr<const VectorSpaceBase<double> > vs_ptr = ptrFromRef(vs_in);
00417   
00418   const Ptr<const SpmdVectorSpaceBase<double> > spmd_vs =
00419     ptr_dynamic_cast<const SpmdVectorSpaceBase<double> >(vs_ptr);
00420 
00421   const Ptr<const ProductVectorSpaceBase<double> > &prod_vs = 
00422     ptr_dynamic_cast<const ProductVectorSpaceBase<double> >(vs_ptr);
00423 
00424   TEUCHOS_TEST_FOR_EXCEPTION( is_null(spmd_vs) && is_null(prod_vs), std::logic_error,
00425     "Error, the concrete VectorSpaceBase object of type "
00426     +Teuchos::demangleName(typeid(vs_in).name())+" does not support the"
00427     " SpmdVectorSpaceBase or the ProductVectorSpaceBase interfaces!" );
00428 
00429   const int numBlocks = (nonnull(prod_vs) ? prod_vs->numBlocks() : 1);
00430 
00431   // Get an array of SpmdVectorBase objects for the blocks
00432   
00433   Array<RCP<const SpmdVectorSpaceBase<double> > > spmd_vs_blocks;
00434   if (nonnull(prod_vs)) {
00435     for (int block_i = 0; block_i < numBlocks; ++block_i) {
00436       const RCP<const SpmdVectorSpaceBase<double> > spmd_vs_i =
00437         rcp_dynamic_cast<const SpmdVectorSpaceBase<double> >(
00438           prod_vs->getBlock(block_i), true);
00439       spmd_vs_blocks.push_back(spmd_vs_i);
00440     }
00441   }
00442   else {
00443     spmd_vs_blocks.push_back(rcpFromPtr(spmd_vs));
00444   }
00445   
00446   // Find the number of local elements, summed over all blocks
00447 
00448   int myLocalElements = 0;
00449   for (int block_i = 0; block_i < numBlocks; ++block_i) {
00450     myLocalElements += spmd_vs_blocks[block_i]->localSubDim();
00451   }
00452   
00453   // Find the GIDs owned by this processor, taken from all blocks
00454   
00455   int count=0;
00456   int blockOffset = 0;
00457   Array<int> myGIDs(myLocalElements);
00458   for (int block_i = 0; block_i < numBlocks; ++block_i) {
00459     const RCP<const SpmdVectorSpaceBase<double> > spmd_vs_i = spmd_vs_blocks[block_i];
00460     const int lowGIDInBlock = spmd_vs_i->localOffset();
00461     const int numLocalElementsInBlock = spmd_vs_i->localSubDim();
00462     for (int i=0; i < numLocalElementsInBlock; ++i, ++count) {
00463       myGIDs[count] = blockOffset + lowGIDInBlock + i;
00464     }
00465     blockOffset += spmd_vs_i->dim();
00466   }
00467   
00468   const int globalDim = vs_in.dim();
00469 
00470   return Teuchos::rcp(
00471     new Epetra_Map(globalDim, myLocalElements, &(myGIDs[0]), 0, *comm));
00472 
00473 }
00474 
00475 
00476 Teuchos::RCP<Epetra_Vector>
00477 Thyra::get_Epetra_Vector(
00478   const Epetra_Map &map,
00479   const RCP<VectorBase<double> > &v
00480   )
00481 {
00482   using Teuchos::get_optional_extra_data;
00483 #ifdef TEUCHOS_DEBUG
00484   RCP<const VectorSpaceBase<double> >
00485     epetra_vs = create_VectorSpace(Teuchos::rcp(&map,false));
00486   THYRA_ASSERT_VEC_SPACES( 
00487     "Thyra::get_Epetra_Vector(map,v)", *epetra_vs, *v->space() );
00488 #endif
00489   //
00490   // First, try to grab the Epetra_Vector straight out of the
00491   // RCP since this is the fastest way.
00492   //
00493   const Ptr<const RCP<Epetra_Vector> >
00494     epetra_v_ptr = get_optional_extra_data<RCP<Epetra_Vector> >(
00495       v,"Epetra_Vector");
00496   if(!is_null(epetra_v_ptr)) {
00497     return *epetra_v_ptr;
00498   }
00499   //
00500   // The assumption that we (rightly) make here is that if the vector spaces
00501   // are compatible, that either the multi-vectors are both in-core or the
00502   // vector spaces are both derived from SpmdVectorSpaceBase and have
00503   // compatible maps.
00504   // 
00505   const VectorSpaceBase<double>  &vs = *v->range();
00506   const SpmdVectorSpaceBase<double> *mpi_vs
00507     = dynamic_cast<const SpmdVectorSpaceBase<double>*>(&vs);
00508   const Ordinal localOffset = ( mpi_vs ? mpi_vs->localOffset() : 0 );
00509   const Ordinal localSubDim = ( mpi_vs ? mpi_vs->localSubDim() : vs.dim() );
00510   //
00511   // Here we will extract a view of the local elements in the underlying
00512   // VectorBase object.  In most cases, no data will be allocated or copied
00513   // and only some small objects will be created and a few virtual functions
00514   // will be called so the overhead should be low and fixed.
00515   //
00516   // Create a *mutable* view of the local elements, this view will be set on
00517   // the RCP that is returned.  As a result, this view will be relased
00518   // when the returned Epetra_Vector is released.
00519   //
00520   // Note that the input vector 'v' will be remembered through this detached
00521   // view!
00522   //
00523   RCP<DetachedVectorView<double> >
00524     emvv = Teuchos::rcp(
00525       new DetachedVectorView<double>(
00526         v
00527         ,Range1D(localOffset,localOffset+localSubDim-1)
00528         ,true // forceContiguous
00529         )
00530       );
00531   // Create a temporary Epetra_Vector object and give it
00532   // the above local view
00533   RCP<Epetra_Vector>
00534     epetra_v = Teuchos::rcp(
00535       new Epetra_Vector(
00536         ::View                                 // CV
00537         ,map                                   // Map
00538         ,const_cast<double*>(emvv->values())   // V
00539         )
00540       );
00541   // Give the explict view object to the above Epetra_Vector smart pointer
00542   // object.  In this way, when the client is finished with the Epetra_Vector
00543   // view the destructor from the object in emvv will automatically commit the
00544   // changes to the elements in the input v VectorBase object (reguardless of
00545   // its implementation).  This is truly an elegant result!
00546   Teuchos::set_extra_data( emvv, "emvv", Teuchos::inOutArg(epetra_v),
00547     Teuchos::PRE_DESTROY );
00548   // We are done!
00549   return epetra_v;
00550 }
00551 
00552 
00553 Teuchos::RCP<const Epetra_Vector>
00554 Thyra::get_Epetra_Vector(
00555   const Epetra_Map &map, 
00556   const RCP<const VectorBase<double> > &v
00557   )
00558 {
00559   using Teuchos::get_optional_extra_data;
00560 #ifdef TEUCHOS_DEBUG
00561   RCP<const VectorSpaceBase<double> >
00562     epetra_vs = create_VectorSpace(Teuchos::rcp(&map,false));
00563   THYRA_ASSERT_VEC_SPACES(
00564     "Thyra::get_Epetra_Vector(map,v)", *epetra_vs, *v->space() );
00565 #endif
00566   //
00567   // First, try to grab the Epetra_Vector straight out of the
00568   // RCP since this is the fastest way.
00569   //
00570   const Ptr<const RCP<const Epetra_Vector> >
00571     epetra_v_ptr = get_optional_extra_data<RCP<const Epetra_Vector> >(
00572       v,"Epetra_Vector");
00573   if(!is_null(epetra_v_ptr))
00574     return *epetra_v_ptr;
00575   const Ptr<const RCP<Epetra_Vector> >
00576     epetra_nonconst_v_ptr = get_optional_extra_data<RCP<Epetra_Vector> >(
00577       v,"Epetra_Vector");
00578   if(!is_null(epetra_nonconst_v_ptr))
00579     return *epetra_nonconst_v_ptr;
00580   //
00581   // Same as above function except as stated below
00582   //
00583   const VectorSpaceBase<double> &vs = *v->range();
00584   const SpmdVectorSpaceBase<double> *mpi_vs
00585     = dynamic_cast<const SpmdVectorSpaceBase<double>*>(&vs);
00586   const Ordinal localOffset = ( mpi_vs ? mpi_vs->localOffset() : 0 );
00587   const Ordinal localSubDim = ( mpi_vs ? mpi_vs->localSubDim() : vs.dim() );
00588   // Get an explicit *non-mutable* view of all of the elements in the multi
00589   // vector.  Note that 'v' will be remembered by this view!
00590   RCP<ConstDetachedVectorView<double> >
00591     evv = Teuchos::rcp(
00592       new ConstDetachedVectorView<double>(
00593         v
00594         ,Range1D(localOffset,localOffset+localSubDim-1)
00595         ,true // forceContiguous
00596         )
00597       );
00598   // Create a temporary Epetra_Vector object and give it
00599   // the above view
00600   RCP<Epetra_Vector>
00601     epetra_v = Teuchos::rcp(
00602       new Epetra_Vector(
00603         ::View                                  // CV
00604         ,map                                    // Map
00605         ,const_cast<double*>(evv->values())     // V
00606         )
00607       );
00608   // This next statement will cause the destructor to free the view if
00609   // needed (see above function).  Since this view is non-mutable,
00610   // only a releaseDetachedView(...) and not a commit will be called.
00611   // This is the whole reason there is a seperate implementation for
00612   // the const and non-const cases.
00613   Teuchos::set_extra_data( evv, "evv", Teuchos::inOutArg(epetra_v),
00614     Teuchos::PRE_DESTROY );
00615   // We are done!
00616   return epetra_v;
00617 }
00618 
00619 
00620 Teuchos::RCP<Epetra_MultiVector>
00621 Thyra::get_Epetra_MultiVector(
00622   const Epetra_Map &map,
00623   const RCP<MultiVectorBase<double> > &mv
00624   )
00625 {
00626   using Teuchos::get_optional_extra_data;
00627 #ifdef TEUCHOS_DEBUG
00628   RCP<const VectorSpaceBase<double> >
00629     epetra_vs = create_VectorSpace(Teuchos::rcp(&map,false));
00630   THYRA_ASSERT_VEC_SPACES(
00631     "Thyra::get_Epetra_MultiVector(map,mv)", *epetra_vs, *mv->range() );
00632 #endif
00633   //
00634   // First, try to grab the Epetra_MultiVector straight out of the
00635   // RCP since this is the fastest way.
00636   //
00637   const Ptr<const RCP<Epetra_MultiVector> >
00638     epetra_mv_ptr = get_optional_extra_data<RCP<Epetra_MultiVector> >(
00639       mv,"Epetra_MultiVector");
00640   if(!is_null(epetra_mv_ptr)) {
00641     return *epetra_mv_ptr;
00642   }
00643   //
00644   // The assumption that we (rightly) make here is that if the vector spaces
00645   // are compatible, that either the multi-vectors are both in-core or the
00646   // vector spaces are both derived from SpmdVectorSpaceBase and have
00647   // compatible maps.
00648   // 
00649   const VectorSpaceBase<double> &vs = *mv->range();
00650   const SpmdVectorSpaceBase<double> *mpi_vs
00651     = dynamic_cast<const SpmdVectorSpaceBase<double>*>(&vs);
00652   const Ordinal localOffset = ( mpi_vs ? mpi_vs->localOffset() : 0 );
00653   const Ordinal localSubDim = ( mpi_vs ? mpi_vs->localSubDim() : vs.dim() );
00654   //
00655   // Here we will extract a view of the local elements in the underlying
00656   // MultiVectorBase object.  In most cases, no data will be allocated or
00657   // copied and only some small objects will be created and a few virtual
00658   // functions will be called so the overhead should be low and fixed.
00659   //
00660   // Create a *mutable* view of the local elements, this view will be set on
00661   // the RCP that is returned.  As a result, this view will be relased
00662   // when the returned Epetra_MultiVector is released.
00663   //
00664   RCP<DetachedMultiVectorView<double> >
00665     emmvv = Teuchos::rcp(
00666       new DetachedMultiVectorView<double>(
00667         *mv
00668         ,Range1D(localOffset,localOffset+localSubDim-1)
00669         )
00670       );
00671   // Create a temporary Epetra_MultiVector object and give it
00672   // the above local view
00673   RCP<Epetra_MultiVector>
00674     epetra_mv = Teuchos::rcp(
00675       new Epetra_MultiVector(
00676         ::View                                  // CV
00677         ,map                                    // Map
00678         ,const_cast<double*>(emmvv->values())   // A
00679         ,emmvv->leadingDim()                    // MyLDA
00680         ,emmvv->numSubCols()                    // NumVectors
00681         )
00682       );
00683   // Give the explict view object to the above Epetra_MultiVector
00684   // smart pointer object.  In this way, when the client is finished
00685   // with the Epetra_MultiVector view the destructor from the object
00686   // in emmvv will automatically commit the changes to the elements in
00687   // the input mv MultiVectorBase object (reguardless of its
00688   // implementation).  This is truly an elegant result!
00689   Teuchos::set_extra_data( emmvv, "emmvv", Teuchos::inOutArg(epetra_mv),
00690     Teuchos::PRE_DESTROY );
00691   // Also set the mv itself as extra data just to be safe
00692   Teuchos::set_extra_data( mv, "mv", Teuchos::inOutArg(epetra_mv) );
00693   // We are done!
00694   return epetra_mv;
00695 }
00696 
00697 
00698 Teuchos::RCP<const Epetra_MultiVector>
00699 Thyra::get_Epetra_MultiVector(
00700   const Epetra_Map &map,
00701   const RCP<const MultiVectorBase<double> > &mv
00702   )
00703 {
00704   using Teuchos::get_optional_extra_data;
00705 #ifdef TEUCHOS_DEBUG
00706   RCP<const VectorSpaceBase<double> >
00707     epetra_vs = create_VectorSpace(Teuchos::rcp(&map,false));
00708   THYRA_ASSERT_VEC_SPACES(
00709     "Thyra::get_Epetra_MultiVector(map,mv)",
00710     *epetra_vs, *mv->range() );
00711 #endif
00712   //
00713   // First, try to grab the Epetra_MultiVector straight out of the
00714   // RCP since this is the fastest way.
00715   //
00716   const Ptr<const RCP<const Epetra_MultiVector> >
00717     epetra_mv_ptr
00718     = get_optional_extra_data<RCP<const Epetra_MultiVector> >(
00719       mv,"Epetra_MultiVector" );
00720   if(!is_null(epetra_mv_ptr)) {
00721     return *epetra_mv_ptr;
00722   }
00723   //
00724   // Same as above function except as stated below
00725   //
00726   const VectorSpaceBase<double> &vs = *mv->range();
00727   const SpmdVectorSpaceBase<double> *mpi_vs
00728     = dynamic_cast<const SpmdVectorSpaceBase<double>*>(&vs);
00729   const Ordinal localOffset = ( mpi_vs ? mpi_vs->localOffset() : 0 );
00730   const Ordinal localSubDim = ( mpi_vs ? mpi_vs->localSubDim() : vs.dim() );
00731   // Get an explicit *non-mutable* view of all of the elements in
00732   // the multi vector.
00733   RCP<ConstDetachedMultiVectorView<double> >
00734     emvv = Teuchos::rcp(
00735       new ConstDetachedMultiVectorView<double>(
00736         *mv
00737         ,Range1D(localOffset,localOffset+localSubDim-1)
00738         )
00739       );
00740   // Create a temporary Epetra_MultiVector object and give it
00741   // the above view
00742   RCP<Epetra_MultiVector>
00743     epetra_mv = Teuchos::rcp(
00744       new Epetra_MultiVector(
00745         ::View                                  // CV
00746         ,map                                    // Map
00747         ,const_cast<double*>(emvv->values())    // A
00748         ,emvv->leadingDim()                     // MyLDA
00749         ,emvv->numSubCols()                     // NumVectors
00750         )
00751       );
00752   // This next statement will cause the destructor to free the view if
00753   // needed (see above function).  Since this view is non-mutable,
00754   // only a releaseDetachedView(...) and not a commit will be called.
00755   // This is the whole reason there is a seperate implementation for
00756   // the const and non-const cases.
00757   Teuchos::set_extra_data( emvv, "emvv", Teuchos::inOutArg(epetra_mv),
00758     Teuchos::PRE_DESTROY );
00759   // Also set the mv itself as extra data just to be safe
00760   Teuchos::set_extra_data( mv, "mv", Teuchos::inOutArg(epetra_mv) );
00761   // We are done!
00762   return epetra_mv;
00763 }
00764 
00765 
00766 Teuchos::RCP<Epetra_MultiVector>
00767 Thyra::get_Epetra_MultiVector(
00768   const Epetra_Map &map,
00769   MultiVectorBase<double> &mv
00770   )
00771 {
00772   using Teuchos::rcpWithEmbeddedObj;
00773   using Teuchos::rcpFromRef;
00774   using Teuchos::outArg;
00775   ArrayRCP<double> mvData;
00776   Ordinal mvLeadingDim = -1;
00777   SpmdMultiVectorBase<double> *mvSpmdMv = 0;
00778   SpmdVectorBase<double> *mvSpmdV = 0;
00779   if ((mvSpmdMv = dynamic_cast<SpmdMultiVectorBase<double>*>(&mv))) {
00780     mvSpmdMv->getNonconstLocalData(outArg(mvData), outArg(mvLeadingDim));
00781   }
00782   else if ((mvSpmdV = dynamic_cast<SpmdVectorBase<double>*>(&mv))) {
00783     mvSpmdV->getNonconstLocalData(outArg(mvData));
00784     mvLeadingDim = mvSpmdV->spmdSpace()->localSubDim();
00785   }
00786   if (nonnull(mvData)) {
00787     return rcpWithEmbeddedObj(
00788       new Epetra_MultiVector(
00789         ::View, map, mvData.getRawPtr(), mvLeadingDim, mv.domain()->dim()
00790         ),
00791       mvData
00792       );
00793   }
00794   return ::Thyra::get_Epetra_MultiVector(map, rcpFromRef(mv));
00795 }
00796 
00797 
00798 Teuchos::RCP<const Epetra_MultiVector>
00799 Thyra::get_Epetra_MultiVector(
00800   const Epetra_Map &map,
00801   const MultiVectorBase<double> &mv
00802   )
00803 {
00804   using Teuchos::rcpWithEmbeddedObj;
00805   using Teuchos::rcpFromRef;
00806   using Teuchos::outArg;
00807   ArrayRCP<const double> mvData;
00808   Ordinal mvLeadingDim = -1;
00809   const SpmdMultiVectorBase<double> *mvSpmdMv = 0;
00810   const SpmdVectorBase<double> *mvSpmdV = 0;
00811   if ((mvSpmdMv = dynamic_cast<const SpmdMultiVectorBase<double>*>(&mv))) {
00812     mvSpmdMv->getLocalData(outArg(mvData), outArg(mvLeadingDim));
00813   }
00814   else if ((mvSpmdV = dynamic_cast<const SpmdVectorBase<double>*>(&mv))) {
00815     mvSpmdV->getLocalData(outArg(mvData));
00816     mvLeadingDim = mvSpmdV->spmdSpace()->localSubDim();
00817   }
00818   if (nonnull(mvData)) {
00819     return rcpWithEmbeddedObj(
00820       new Epetra_MultiVector(
00821         ::View,map, const_cast<double*>(mvData.getRawPtr()), mvLeadingDim, mv.domain()->dim()
00822         ),
00823       mvData
00824       );
00825   }
00826   return ::Thyra::get_Epetra_MultiVector(map, rcpFromRef(mv));
00827 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines