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