Amesos2 - Direct Sparse Solver Interfaces Version of the Day
Amesos2_EpetraMultiVecAdapter_def.hpp
Go to the documentation of this file.
00001 // @HEADER
00002 //
00003 // ***********************************************************************
00004 //
00005 //           Amesos2: Templated Direct Sparse Solver Package 
00006 //                  Copyright 2011 Sandia Corporation
00007 //
00008 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00009 // the U.S. Government retains certain rights in this software.
00010 //
00011 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00039 //
00040 // ***********************************************************************
00041 //
00042 // @HEADER
00043 
00053 #ifndef AMESOS2_EPETRA_MULTIVEC_ADAPTER_DEF_HPP
00054 #define AMESOS2_EPETRA_MULTIVEC_ADAPTER_DEF_HPP
00055 
00056 #include <Teuchos_as.hpp>
00057 
00058 #include <Epetra_SerialComm.h>
00059 #ifdef HAVE_MPI
00060 #include <Epetra_MpiComm.h>
00061 #endif
00062 #include <Epetra_LocalMap.h>
00063 #include <Epetra_Import.h>
00064 #include <Epetra_Export.h>
00065 
00066 #include "Amesos2_EpetraMultiVecAdapter_decl.hpp"
00067 #include "Amesos2_Util.hpp"
00068 
00069 namespace Amesos2 {
00070 
00071 MultiVecAdapter<Epetra_MultiVector>::MultiVecAdapter(const MultiVecAdapter<multivec_t>& adapter)
00072   : mv_(adapter.mv_)
00073   , mv_map_(adapter.mv_map_)
00074 { }
00075 
00076 MultiVecAdapter<Epetra_MultiVector>::MultiVecAdapter(const Teuchos::RCP<multivec_t>& m)
00077   : mv_(m)
00078 {
00079   mv_map_ = Teuchos::rcpFromRef(mv_->Map());
00080 }
00081   
00082 
00083 bool MultiVecAdapter<Epetra_MultiVector>::isLocallyIndexed() const
00084 {
00085   return !mv_->DistributedGlobal();
00086 }
00087 
00088 bool MultiVecAdapter<Epetra_MultiVector>::isGloballyIndexed() const
00089 {
00090   return mv_->DistributedGlobal();
00091 }
00092 
00093 
00094 const Teuchos::RCP<const Teuchos::Comm<int> >
00095 MultiVecAdapter<Epetra_MultiVector>::getComm() const
00096 {
00097   return Util::to_teuchos_comm(Teuchos::rcpFromRef(mv_->Comm()));
00098 }
00099 
00100 
00101 size_t MultiVecAdapter<Epetra_MultiVector>::getLocalLength() const
00102 {
00103   return Teuchos::as<size_t>(mv_->MyLength());
00104 }
00105 
00106 
00107 size_t MultiVecAdapter<Epetra_MultiVector>::getLocalNumVectors() const
00108 {
00109   return Teuchos::as<size_t>(mv_->NumVectors());
00110 }
00111 
00112 
00113 MultiVecAdapter<Epetra_MultiVector>::global_size_t
00114 MultiVecAdapter<Epetra_MultiVector>::getGlobalLength() const
00115 {
00116   return Teuchos::as<global_size_t>(mv_->GlobalLength());
00117 }
00118 
00119 
00120 size_t MultiVecAdapter<Epetra_MultiVector>::getGlobalNumVectors() const
00121 {
00122   return Teuchos::as<size_t>(mv_->NumVectors());
00123 }
00124 
00125 
00126 size_t MultiVecAdapter<Epetra_MultiVector>::getStride() const
00127 {
00128   return Teuchos::as<size_t>(mv_->Stride());
00129 }
00130 
00131 
00132 bool MultiVecAdapter<Epetra_MultiVector>::isConstantStride() const
00133 {
00134   return mv_->ConstantStride();
00135 }
00136 
00137 
00138 Teuchos::RCP<const Tpetra::Vector<MultiVecAdapter<Epetra_MultiVector>::scalar_t,
00139                                   MultiVecAdapter<Epetra_MultiVector>::local_ordinal_t,
00140                                   MultiVecAdapter<Epetra_MultiVector>::global_ordinal_t,
00141                                   MultiVecAdapter<Epetra_MultiVector>::node_t> >
00142 MultiVecAdapter<Epetra_MultiVector>::getVector( size_t j ) const
00143 {
00144   using Teuchos::RCP;
00145   using Teuchos::ArrayRCP;
00146   using Tpetra::MultiVector;
00147 
00148   typedef scalar_t st;
00149   typedef local_ordinal_t lot;
00150   typedef global_ordinal_t got;
00151   typedef node_t nt;
00152   
00153   RCP<MultiVector<st,lot,got,nt> > vector = rcp(new MultiVector<st,lot,got,nt>(this->getMap(),1));
00154 
00155   // Copy vector contents into Tpetra multi-vector
00156   ArrayRCP<st> it = vector->getDataNonConst(0);
00157   double* vector_data = mv_->operator[](Teuchos::as<int>(j)); // values from j^th vector
00158   Tpetra::global_size_t size = vector->getGlobalLength();
00159 
00160   for( Tpetra::global_size_t i = 0; i < size; ++i ){
00161     *it = vector_data[i];
00162   }
00163 
00164   return vector->getVector(j);
00165 }
00166 
00167 
00168 // Implementation is essentially the same as getVector()
00169 Teuchos::RCP<Tpetra::Vector<MultiVecAdapter<Epetra_MultiVector>::scalar_t,
00170                             MultiVecAdapter<Epetra_MultiVector>::local_ordinal_t,
00171                             MultiVecAdapter<Epetra_MultiVector>::global_ordinal_t,
00172                             MultiVecAdapter<Epetra_MultiVector>::node_t> >
00173 MultiVecAdapter<Epetra_MultiVector>::getVectorNonConst( size_t j )
00174 {
00175   using Teuchos::RCP;
00176   using Teuchos::ArrayRCP;
00177   using Tpetra::MultiVector;
00178 
00179   typedef scalar_t st;
00180   typedef local_ordinal_t lot;
00181   typedef global_ordinal_t got;
00182   typedef node_t nt;
00183   
00184   RCP<MultiVector<st,lot,got,nt> > vector = rcp(new MultiVector<st,lot,got,nt>(this->getMap(),1));
00185 
00186   // Copy vector contents into Tpetra multi-vector
00187   ArrayRCP<st> it = vector->getDataNonConst(0);
00188   double* vector_data = mv_->operator[](Teuchos::as<int>(j)); // values from j^th vector
00189   Tpetra::global_size_t size = vector->getGlobalLength();
00190 
00191   for( Tpetra::global_size_t i = 0; i < size; ++i ){
00192     *it = vector_data[i];
00193   }
00194 
00195   return vector->getVectorNonConst(j);
00196 }
00197 
00198 
00199 void MultiVecAdapter<Epetra_MultiVector>::get1dCopy(
00200   const Teuchos::ArrayView<MultiVecAdapter<Epetra_MultiVector>::scalar_t>& av,
00201   size_t lda,
00202   Teuchos::Ptr<
00203     const Tpetra::Map<MultiVecAdapter<Epetra_MultiVector>::local_ordinal_t,
00204                       MultiVecAdapter<Epetra_MultiVector>::global_ordinal_t,
00205                       MultiVecAdapter<Epetra_MultiVector>::node_t> > distribution_map ) const
00206 {
00207   using Teuchos::rcpFromPtr;
00208   using Teuchos::as;
00209 
00210   const size_t num_vecs = getGlobalNumVectors();
00211   
00212 #ifdef HAVE_AMESOS2_DEBUG
00213   const size_t requested_vector_length = distribution_map->getNodeNumElements();
00214   TEUCHOS_TEST_FOR_EXCEPTION( lda < requested_vector_length,
00215           std::invalid_argument,
00216           "Given stride is not large enough for local vector length" );
00217   TEUCHOS_TEST_FOR_EXCEPTION( as<size_t>(av.size()) < (num_vecs-1) * lda + requested_vector_length,
00218           std::invalid_argument,
00219           "MultiVector storage not large enough given leading dimension "
00220           "and number of vectors" );
00221 #endif
00222 
00223   Epetra_Map e_dist_map
00224     = *Util::tpetra_map_to_epetra_map<local_ordinal_t,
00225                                       global_ordinal_t,
00226                                       global_size_t,
00227                                       node_t>(*distribution_map);
00228 
00229   multivec_t redist_mv(e_dist_map, as<int>(num_vecs));
00230   const Epetra_Import importer(e_dist_map, *mv_map_); // Note, target/source order is reversed in Tpetra
00231   redist_mv.Import(*mv_, importer, Insert);
00232 
00233   // Finally, do copy
00234   redist_mv.ExtractCopy(av.getRawPtr(), lda);
00235 }
00236 
00237 
00238 Teuchos::ArrayRCP<MultiVecAdapter<Epetra_MultiVector>::scalar_t>
00239 MultiVecAdapter<Epetra_MultiVector>::get1dViewNonConst(bool local)
00240 {
00241   // TEUCHOS_TEST_FOR_EXCEPTION( !this->isConstantStride(),
00242   //   std::logic_error,
00243   //   "get1dViewNonConst() : can only get 1d view if stride is constant");
00244 
00245   // if( local ){
00246   //   TEUCHOS_TEST_FOR_EXCEPTION(
00247   //     true,
00248   //     std::logic_error,
00249   //     "Amesos2::MultiVecAdapter<Epetra_MultiVector> : 1D views not yet supported for local-local Epetra multi-vectors");
00250     
00251   //   // localize();
00252   //   // /* Use the global element list returned by
00253   //   //  * mv_->getMap()->getNodeElementList() to get a subCopy of mv_ which we
00254   //   //  * assign to l_l_mv_, then return get1dViewNonConst() of l_l_mv_
00255   //   //  */
00256   //   // l_l_mv_ = Teuchos::null;
00257     
00258   //   // Teuchos::Array<GlobalOrdinal> nodeElements_go(mv_->Map().NumMyElements());
00259   //   // mv_->Map().MyGlobalElements(nodeElements_go.getRawPtr());
00260   //   // Teuchos::Array<size_t> nodeElements_st(nodeElements_go.size());
00261 
00262   //   // // Convert the node element to a list of size_t type objects
00263   //   // typename Teuchos::ArrayView<const GlobalOrdinal>::iterator it_go;
00264   //   // Teuchos::Array<size_t>::iterator it_st = nodeElements_st.begin();
00265   //   // for( it_go = nodeElements_go.begin(); it_go != nodeElements_go.end(); ++it_go ){
00266   //   //   *(it_st++) = Teuchos::as<size_t>(*it_go);
00267   //   // }
00268 
00269   //   // l_l_mv_ = l_mv_->subViewNonConst(nodeElements_st);
00270 
00271   //   // return(l_l_mv_->get1dViewNonConst());
00272   // } else {
00273   //   scalar_t* values;
00274   //   int lda;
00275 
00276   //   if( !isLocal() ){
00277   //     this->localize();
00278   //     l_mv_->ExtractView(&values, &lda);
00279   //   } else {
00280   //     mv_->ExtractView(&values, &lda);
00281   //   }
00282 
00283   //   TEUCHOS_TEST_FOR_EXCEPTION( lda != Teuchos::as<int>(this->getStride()),
00284   //     std::logic_error,
00285   //     "Stride reported during extraction not consistent with what multivector reports");
00286 
00287   //   return Teuchos::arcp(values,0,lda*this->getGlobalNumVectors(),false);
00288   // }
00289   return Teuchos::null;
00290 }
00291 
00292 
00293 void
00294 MultiVecAdapter<Epetra_MultiVector>::put1dData(
00295   const Teuchos::ArrayView<const MultiVecAdapter<Epetra_MultiVector>::scalar_t>& new_data,
00296   size_t lda,
00297   Teuchos::Ptr<
00298     const Tpetra::Map<MultiVecAdapter<Epetra_MultiVector>::local_ordinal_t,
00299                       MultiVecAdapter<Epetra_MultiVector>::global_ordinal_t,
00300                       MultiVecAdapter<Epetra_MultiVector>::node_t> > source_map)
00301 {
00302   using Teuchos::rcpFromPtr;
00303   using Teuchos::as;
00304 
00305   const size_t num_vecs  = getGlobalNumVectors();
00306   // TODO: check that the following const_cast is safe
00307   double* data_ptr = const_cast<double*>(new_data.getRawPtr());
00308   const Epetra_BlockMap e_source_map
00309     = *Util::tpetra_map_to_epetra_map<local_ordinal_t,global_ordinal_t,global_size_t,node_t>(*source_map);
00310   const multivec_t source_mv(Copy, e_source_map, data_ptr, as<int>(lda), as<int>(num_vecs));
00311   const Epetra_Import importer(*mv_map_, e_source_map);
00312   
00313   mv_->Import(source_mv, importer, Insert);
00314 }
00315 
00316 
00317 std::string MultiVecAdapter<Epetra_MultiVector>::description() const
00318 {
00319   std::ostringstream oss;
00320   oss << "Amesos2 adapter wrapping: Epetra_MultiVector";
00321   return oss.str();
00322 }
00323 
00324 
00325 void MultiVecAdapter<Epetra_MultiVector>::describe(
00326   Teuchos::FancyOStream& os,
00327   const Teuchos::EVerbosityLevel verbLevel=
00328   Teuchos::Describable::verbLevel_default) const
00329 {
00330   // TODO: implement!
00331 }
00332 
00333 
00334 Teuchos::RCP<const Tpetra::Map<MultiVecAdapter<Epetra_MultiVector>::local_ordinal_t,
00335                                MultiVecAdapter<Epetra_MultiVector>::global_ordinal_t,
00336                                MultiVecAdapter<Epetra_MultiVector>::node_t> >
00337 MultiVecAdapter<Epetra_MultiVector>::getMap() const
00338 {
00339   return Util::epetra_map_to_tpetra_map<local_ordinal_t,global_ordinal_t,global_size_t,node_t>(*mv_map_);
00340 }
00341 
00342 
00343 const char* MultiVecAdapter<Epetra_MultiVector>::name
00344 = "Amesos2 adapter for Epetra_MultiVector";
00345 
00346 
00347 } // end namespace Amesos2
00348 
00349 #endif // AMESOS2_EPETRA_MULTIVEC_ADAPTER_DEF_HPP