Anasazi Version of the Day
AnasaziTpetraAdapter.hpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //                 Anasazi: Block Eigensolvers Package
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 #ifndef ANASAZI_TPETRA_ADAPTER_HPP
00030 #define ANASAZI_TPETRA_ADAPTER_HPP
00031 
00032 #include <Kokkos_NodeTrace.hpp>
00033 
00038 // TODO: the assumption is made that the solver, multivector and operator are templated on the same scalar. this will need to be modified.
00039 
00040 #include <Tpetra_MultiVector.hpp>
00041 #include <Tpetra_Operator.hpp>
00042 #include <Teuchos_TestForException.hpp>
00043 #include <Teuchos_ScalarTraits.hpp>
00044 #include <Teuchos_Array.hpp>
00045 #include <Teuchos_DefaultSerialComm.hpp>
00046 
00047 #include "AnasaziConfigDefs.hpp"
00048 #include "AnasaziTypes.hpp"
00049 #include "AnasaziMultiVecTraits.hpp"
00050 #include "AnasaziOperatorTraits.hpp"
00051 #include <Kokkos_NodeAPIConfigDefs.hpp>
00052 
00053 #ifdef HAVE_ANASAZI_TSQR
00054 #  include "AnasaziTsqrAdaptor.hpp" 
00055 #include "TsqrAdaptor_Tpetra_MultiVector.hpp" 
00056 #include "TsqrRandomizer_Tpetra_MultiVector.hpp" 
00057 #endif // HAVE_ANASAZI_TSQR
00058 
00059 
00060 namespace Anasazi {
00061 
00063   //
00064   // Implementation of the Anasazi::MultiVecTraits for Tpetra::MultiVector.
00065   //
00067 
00072   template<class Scalar, class LO, class GO, class Node>
00073   class MultiVecTraits<Scalar, Tpetra::MultiVector<Scalar,LO,GO,Node> >
00074   {
00075   public:
00076 #ifdef HAVE_ANASAZI_TPETRA_TIMERS
00077     static Teuchos::RCP<Teuchos::Time> mvTimesMatAddMvTimer_, mvTransMvTimer_;
00078 #endif
00079 
00080     static Teuchos::RCP<Tpetra::MultiVector<Scalar,LO,GO,Node> > Clone( const Tpetra::MultiVector<Scalar,LO,GO,Node>& mv, const int numvecs )
00081     { 
00082       return Teuchos::rcp( new Tpetra::MultiVector<Scalar,LO,GO,Node>(mv.getMap(),numvecs));
00083     }
00084 
00085     static Teuchos::RCP<Tpetra::MultiVector<Scalar,LO,GO,Node> > CloneCopy( const Tpetra::MultiVector<Scalar,LO,GO,Node>& mv )
00086     {
00087       KOKKOS_NODE_TRACE("Anasazi::MVT::CloneCopy(MV)")
00088       return Teuchos::rcp( new Tpetra::MultiVector<Scalar,LO,GO,Node>( mv ) ); 
00089     }
00090 
00091     static Teuchos::RCP<Tpetra::MultiVector<Scalar,LO,GO,Node> > CloneCopy( const Tpetra::MultiVector<Scalar,LO,GO,Node>& mv, const std::vector<int>& index )
00092     { 
00093       KOKKOS_NODE_TRACE("Anasazi::MVT::CloneCopy(MV,ind)")
00094       TEST_FOR_EXCEPTION(index.size() == 0,std::runtime_error,
00095           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::Clone(mv,index): numvecs must be greater than zero.");
00096 #ifdef HAVE_TPETRA_DEBUG
00097       TEST_FOR_EXCEPTION( *std::min_element(index.begin(),index.end()) < 0, std::runtime_error,
00098           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::Clone(mv,index): indices must be >= zero.");
00099       TEST_FOR_EXCEPTION( (size_t)*std::max_element(index.begin(),index.end()) >= mv.getNumVectors(), std::runtime_error,
00100           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::Clone(mv,index): indices must be < mv.getNumVectors().");
00101 #endif
00102       for (typename std::vector<int>::size_type j=1; j<index.size(); ++j) {
00103         if (index[j] != index[j-1]+1) {
00104           // not contiguous; short circuit
00105           Teuchos::Array<size_t> stinds(index.begin(), index.end());
00106           return mv.subCopy(stinds());
00107         }
00108       }
00109       // contiguous
00110       return mv.subCopy(Teuchos::Range1D(index.front(),index.back()));
00111     }
00112 
00113     static Teuchos::RCP<Tpetra::MultiVector<Scalar,LO,GO,Node> > CloneViewNonConst( Tpetra::MultiVector<Scalar,LO,GO,Node>& mv, const std::vector<int>& index )
00114     {
00115       TEST_FOR_EXCEPTION(index.size() == 0,std::invalid_argument,
00116           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::CloneViewNonConst(mv,index): numvecs must be greater than zero.");
00117 #ifdef HAVE_TPETRA_DEBUG
00118       TEST_FOR_EXCEPTION( *std::min_element(index.begin(),index.end()) < 0, std::invalid_argument,
00119           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::CloneViewNonConst(mv,index): indices must be >= zero.");
00120       TEST_FOR_EXCEPTION( (size_t)*std::max_element(index.begin(),index.end()) >= mv.getNumVectors(), std::invalid_argument,
00121           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::CloneViewNonConst(mv,index): indices must be < mv.getNumVectors().");
00122 #endif
00123       for (typename std::vector<int>::size_type j=1; j<index.size(); ++j) {
00124         if (index[j] != index[j-1]+1) {
00125           // not contiguous; short circuit
00126           Teuchos::Array<size_t> stinds(index.begin(), index.end());
00127           return mv.subViewNonConst(stinds);
00128         }
00129       }
00130       // contiguous
00131       return mv.subViewNonConst(Teuchos::Range1D(index.front(),index.back()));
00132     }
00133 
00134     static Teuchos::RCP<const Tpetra::MultiVector<Scalar,LO,GO,Node> > CloneView(const Tpetra::MultiVector<Scalar,LO,GO,Node>& mv, const std::vector<int>& index )
00135     {
00136       TEST_FOR_EXCEPTION(index.size() == 0,std::invalid_argument,
00137           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::CloneView(mv,index): numvecs must be greater than zero.");
00138 #ifdef HAVE_TPETRA_DEBUG
00139       TEST_FOR_EXCEPTION( *std::min_element(index.begin(),index.end()) < 0, std::invalid_argument,
00140           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::CloneView(mv,index): indices must be >= zero.");
00141       TEST_FOR_EXCEPTION( (size_t)*std::max_element(index.begin(),index.end()) >= mv.getNumVectors(), std::invalid_argument,
00142           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::CloneView(mv,index): indices must be < mv.getNumVectors().");
00143 #endif
00144       for (typename std::vector<int>::size_type j=1; j<index.size(); ++j) {
00145         if (index[j] != index[j-1]+1) {
00146           // not contiguous; short circuit
00147           Teuchos::Array<size_t> stinds(index.begin(), index.end());
00148           return mv.subView(stinds);
00149         }
00150       }
00151       // contiguous
00152       return mv.subView(Teuchos::Range1D(index.front(),index.back()));
00153     }
00154 
00155     static int GetVecLength( const Tpetra::MultiVector<Scalar,LO,GO,Node>& mv )
00156     { return mv.getGlobalLength(); }
00157 
00158     static int GetNumberVecs( const Tpetra::MultiVector<Scalar,LO,GO,Node>& mv )
00159     { return mv.getNumVectors(); }
00160 
00161     static void MvTimesMatAddMv( Scalar alpha, const Tpetra::MultiVector<Scalar,LO,GO,Node>& A, 
00162                                  const Teuchos::SerialDenseMatrix<int,Scalar>& B, 
00163                                  Scalar beta, Tpetra::MultiVector<Scalar,LO,GO,Node>& mv )
00164     {
00165       KOKKOS_NODE_TRACE("Anasazi::MVT::MvTimesMatAddMv()")
00166 #ifdef HAVE_ANASAZI_TPETRA_TIMERS
00167       Teuchos::TimeMonitor lcltimer(*mvTimesMatAddMvTimer_);
00168 #endif
00169       // create local map
00170       Teuchos::SerialComm<int> scomm;
00171       Tpetra::Map<LO,GO,Node> LocalMap(B.numRows(), 0, Teuchos::rcpFromRef< const Teuchos::Comm<int> >(scomm), Tpetra::LocallyReplicated, A.getMap()->getNode());
00172       // encapsulate Teuchos::SerialDenseMatrix data in ArrayView
00173       Teuchos::ArrayView<const Scalar> Bvalues(B.values(),B.stride()*B.numCols());
00174       // create locally replicated MultiVector with a copy of this data
00175       Tpetra::MultiVector<Scalar,LO,GO,Node> B_mv(Teuchos::rcpFromRef(LocalMap),Bvalues,B.stride(),B.numCols());
00176       // multiply
00177       mv.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, alpha, A, B_mv, beta);
00178     }
00179 
00180     static void MvAddMv( Scalar alpha, const Tpetra::MultiVector<Scalar,LO,GO,Node>& A, Scalar beta, const Tpetra::MultiVector<Scalar,LO,GO,Node>& B, Tpetra::MultiVector<Scalar,LO,GO,Node>& mv )
00181     {
00182       mv.update(alpha,A,beta,B,Teuchos::ScalarTraits<Scalar>::zero());
00183     }
00184 
00185     static void MvScale ( Tpetra::MultiVector<Scalar,LO,GO,Node>& mv, Scalar alpha )
00186     { mv.scale(alpha); }
00187 
00188     static void MvScale ( Tpetra::MultiVector<Scalar,LO,GO,Node>& mv, const std::vector<Scalar>& alphas )
00189     { mv.scale(alphas); }
00190 
00191     static void MvTransMv( Scalar alpha, const Tpetra::MultiVector<Scalar,LO,GO,Node>& A, const Tpetra::MultiVector<Scalar,LO,GO,Node>& B, Teuchos::SerialDenseMatrix<int,Scalar>& C)
00192     { 
00193       KOKKOS_NODE_TRACE("Anasazi::MVT::MvTransMv()")
00194 #ifdef HAVE_ANASAZI_TPETRA_TIMERS
00195       Teuchos::TimeMonitor lcltimer(*mvTransMvTimer_);
00196 #endif
00197       // form alpha * A^H * B, then copy into SDM
00198       // we will create a multivector C_mv from a a local map
00199       // this map has a serial comm, the purpose being to short-circuit the MultiVector::reduce() call at the end of MultiVector::multiply()
00200       // otherwise, the reduced multivector data would be copied back to the GPU, only to turn around and have to get it back here.
00201       // this saves us a round trip for this data.
00202       const int numRowsC = C.numRows(),
00203                 numColsC = C.numCols(),
00204                 strideC  = C.stride();
00205       Teuchos::SerialComm<int> scomm;
00206       // create local map with serial comm
00207       Tpetra::Map<LO,GO,Node> LocalMap(numRowsC, 0, Teuchos::rcpFromRef< const Teuchos::Comm<int> >(scomm), Tpetra::LocallyReplicated, A.getMap()->getNode());
00208       // create local multivector to hold the result
00209       const bool INIT_TO_ZERO = true;
00210       Tpetra::MultiVector<Scalar,LO,GO,Node> C_mv(Teuchos::rcpFromRef(LocalMap),numColsC, INIT_TO_ZERO);
00211       // multiply result into local multivector
00212       C_mv.multiply(Teuchos::CONJ_TRANS,Teuchos::NO_TRANS,alpha,A,B,Teuchos::ScalarTraits<Scalar>::zero());
00213       // get comm
00214       Teuchos::RCP< const Teuchos::Comm<int> > pcomm = A.getMap()->getComm();
00215       // create arrayview encapsulating the Teuchos::SerialDenseMatrix
00216       Teuchos::ArrayView<Scalar> C_view(C.values(),strideC*numColsC);
00217       if (pcomm->getSize() == 1) {
00218         // no accumulation to do; simply extract the multivector data into C
00219         // extract a copy of the result into the array view (and therefore, the SerialDenseMatrix)
00220         C_mv.get1dCopy(C_view,strideC);
00221       }  
00222       else {
00223         // get a const host view of the data in C_mv
00224         Teuchos::ArrayRCP<const Scalar> C_mv_view = C_mv.get1dView();
00225         if (strideC == numRowsC) {
00226           // sumall into C
00227           Teuchos::reduceAll<int,Scalar>(*pcomm,Teuchos::REDUCE_SUM,numColsC*numRowsC,C_mv_view.getRawPtr(),C_view.getRawPtr());
00228         }
00229         else {
00230           // sumall into temp, copy into C
00231           Teuchos::Array<Scalar> destBuff(numColsC*numRowsC);
00232           Teuchos::reduceAll<int,Scalar>(*pcomm,Teuchos::REDUCE_SUM,numColsC*numRowsC,C_mv_view.getRawPtr(),destBuff.getRawPtr());
00233           for (int j=0; j < numColsC; ++j) {
00234             for (int i=0; i < numRowsC; ++i) {
00235               C_view[strideC*j+i] = destBuff[numRowsC*j+i];
00236             }
00237           }
00238         }
00239       }
00240     }
00241 
00242     static void MvDot( const Tpetra::MultiVector<Scalar,LO,GO,Node>& A, const Tpetra::MultiVector<Scalar,LO,GO,Node>& B, std::vector<Scalar> &dots)
00243     {
00244       TEST_FOR_EXCEPTION(A.getNumVectors() != B.getNumVectors(),std::invalid_argument,
00245           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::MvDot(A,B,dots): A and B must have the same number of vectors.");
00246 #ifdef HAVE_TPETRA_DEBUG
00247       TEST_FOR_EXCEPTION(dots.size() < (typename std::vector<int>::size_type)A.getNumVectors(),std::invalid_argument,
00248           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::MvDot(A,B,dots): dots must have room for all dot products.");
00249 #endif
00250       Teuchos::ArrayView<Scalar> av(dots);
00251       A.dot(B,av(0,A.getNumVectors()));
00252     }
00253 
00254     static void MvNorm(const Tpetra::MultiVector<Scalar,LO,GO,Node>& mv, std::vector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> &normvec)
00255     { 
00256 #ifdef HAVE_TPETRA_DEBUG
00257       TEST_FOR_EXCEPTION(normvec.size() < (typename std::vector<int>::size_type)mv.getNumVectors(),std::invalid_argument,
00258           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::MvNorm(mv,normvec): normvec must have room for all norms.");
00259 #endif
00260       Teuchos::ArrayView<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> av(normvec);
00261       mv.norm2(av(0,mv.getNumVectors()));
00262     }
00263 
00264     static void SetBlock( const Tpetra::MultiVector<Scalar,LO,GO,Node>& A, const std::vector<int>& index, Tpetra::MultiVector<Scalar,LO,GO,Node>& mv )
00265     {
00266       KOKKOS_NODE_TRACE("Anasazi::MVT::SetBlock()")
00267 #ifdef HAVE_TPETRA_DEBUG
00268       TEST_FOR_EXCEPTION((typename std::vector<int>::size_type)A.getNumVectors() < index.size(),std::invalid_argument,
00269           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::SetBlock(A,index,mv): index must be the same size as A.");
00270 #endif
00271       Teuchos::RCP<Tpetra::MultiVector<Scalar,LO,GO,Node> > mvsub = CloneViewNonConst(mv,index);
00272       if ((typename std::vector<int>::size_type)A.getNumVectors() > index.size()) {
00273         Teuchos::RCP<const Tpetra::MultiVector<Scalar,LO,GO,Node> > Asub = A.subView(Teuchos::Range1D(0,index.size()-1));
00274         (*mvsub) = (*Asub);
00275       }
00276       else {
00277         (*mvsub) = A;
00278       }
00279       mvsub = Teuchos::null;
00280     }
00281 
00282     static void MvRandom( Tpetra::MultiVector<Scalar,LO,GO,Node>& mv )
00283     { 
00284       KOKKOS_NODE_TRACE("Anasazi::MVT::randomize()")
00285       mv.randomize(); 
00286     }
00287 
00288     static void MvInit( Tpetra::MultiVector<Scalar,LO,GO,Node>& mv, Scalar alpha = Teuchos::ScalarTraits<Scalar>::zero() )
00289     { mv.putScalar(alpha); }
00290 
00291     static void MvPrint( const Tpetra::MultiVector<Scalar,LO,GO,Node>& mv, std::ostream& os )
00292     { mv.print(os); }
00293 
00294   };        
00295 
00296 #ifdef HAVE_ANASAZI_TSQR
00297 
00298   //
00299   // Implementation of Anasazi::TsqrAdaptor for Tpetra::MultiVector.
00300   //
00302 
00303   template< class Scalar, class LO, class GO, class Node >
00304   class TsqrAdaptor< Scalar, Tpetra::MultiVector< Scalar, LO, GO, Node > >
00305   {
00306   public:
00307     typedef TSQR::Trilinos::TsqrTpetraAdaptor< Scalar, LO, GO, Node > adaptor_type;
00308   };
00309 #endif // HAVE_ANASAZI_TSQR
00310 
00312   //
00313   // Implementation of the Anasazi::OperatorTraits for Tpetra::Operator.
00314   //
00316 
00317   template <class Scalar, class LO, class GO, class Node> 
00318   class OperatorTraits < Scalar, Tpetra::MultiVector<Scalar,LO,GO,Node>, Tpetra::Operator<Scalar,LO,GO,Node> >
00319   {
00320   public:
00321     static void Apply ( const Tpetra::Operator<Scalar,LO,GO,Node> & Op, 
00322                         const Tpetra::MultiVector<Scalar,LO,GO,Node> & X,
00323                               Tpetra::MultiVector<Scalar,LO,GO,Node> & Y)
00324     { 
00325       Op.apply(X,Y,Teuchos::NO_TRANS);
00326     }
00327   };
00328 
00329 } // end of Anasazi namespace 
00330 
00331 #endif 
00332 // end of file ANASAZI_TPETRA_ADAPTER_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends