Anasazi: AnasaziTpetraAdapter.hpp Source 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 
00036 // TODO: the assumption is made that the solver, multivector and operator are templated on the same scalar. this will need to be modified.
00037 
00038 #include <Tpetra_MultiVector.hpp>
00039 #include <Tpetra_Operator.hpp>
00040 #include <Teuchos_TestForException.hpp>
00041 #include <Teuchos_ScalarTraits.hpp>
00042 
00043 #include "AnasaziConfigDefs.hpp"
00044 #include "AnasaziTypes.hpp"
00045 #include "AnasaziMultiVecTraits.hpp"
00046 #include "AnasaziOperatorTraits.hpp"
00047 
00048 namespace Anasazi {
00049  
00051   //
00052   // Implementation of the Anasazi::MultiVecTraits for Tpetra::MultiVector
00053   //
00055 
00061   template <class Scalar, class LO, class GO, class Node>
00062   class MultiVecTraits<Scalar, Tpetra::MultiVector<Scalar,LO,GO,Node> >
00063   {
00064   public:
00065 
00067 
00068 
00073     static Teuchos::RCP<Tpetra::MultiVector<Scalar,LO,GO,Node> > Clone( const Tpetra::MultiVector<Scalar,LO,GO,Node>& mv, const int numvecs )
00074     {
00075 #ifdef TPETRA_DEBUG
00076       TEST_FOR_EXCEPTION(numvecs <= 0,std::invalid_argument,
00077           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::Clone(mv,numvecs): numvecs must be greater than zero.");
00078 #endif
00079       return Teuchos::rcp( new Tpetra::MultiVector<Scalar,LO,GO,Node>(mv.getMap(), numvecs) ); 
00080     }
00081 
00086     static Teuchos::RCP<Tpetra::MultiVector<Scalar,LO,GO,Node> > CloneCopy( const Tpetra::MultiVector<Scalar,LO,GO,Node>& mv )
00087     {
00088       return Teuchos::rcp( new Tpetra::MultiVector<Scalar,LO,GO,Node>( mv ) ); 
00089     }
00090 
00096     static Teuchos::RCP<Tpetra::MultiVector<Scalar,LO,GO,Node> > CloneCopy( const Tpetra::MultiVector<Scalar,LO,GO,Node>& mv, const std::vector<int>& index )
00097     {
00098       TEST_FOR_EXCEPTION(index.size() == 0,std::runtime_error,
00099           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::Clone(mv,index): numvecs must be greater than zero.");
00100 #ifdef TPETRA_DEBUG
00101       TEST_FOR_EXCEPTION( *std::min_element(index.begin(),index.end()) < 0, std::runtime_error,
00102           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::Clone(mv,index): indices must be >= zero.");
00103       TEST_FOR_EXCEPTION( *std::max_element(index.begin(),index.end()) >= mv.getNumVectors(), std::runtime_error,
00104           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::Clone(mv,index): indices must be < mv.getNumVectors().");
00105 #endif
00106       for (typename std::vector<int>::size_type j=1; j<index.size(); ++j) {
00107         if (index[j] != index[j-1]+1) {
00108           // not contiguous; short circuit
00109           Teuchos::Array<size_t> inds(index.begin(), index.end());
00110           return mv.subCopy(inds());
00111         }
00112       }
00113       // contiguous
00114       return mv.subCopy(Teuchos::Range1D(index.front(),index.back()));
00115     }
00116 
00122     static Teuchos::RCP<Tpetra::MultiVector<Scalar,LO,GO,Node> > CloneView( Tpetra::MultiVector<Scalar,LO,GO,Node>& mv, const std::vector<int>& index )
00123     {
00124       TEST_FOR_EXCEPTION(index.size() == 0,std::invalid_argument,
00125           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::CloneView(mv,index): numvecs must be greater than zero.");
00126 #ifdef TPETRA_DEBUG
00127       TEST_FOR_EXCEPTION( *std::min_element(index.begin(),index.end()) < 0, std::invalid_argument,
00128           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::CloneView(mv,index): indices must be >= zero.");
00129       TEST_FOR_EXCEPTION( *std::max_element(index.begin(),index.end()) >= mv.getNumVectors(), std::invalid_argument,
00130           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::CloneView(mv,index): indices must be < mv.getNumVectors().");
00131 #endif
00132       for (typename std::vector<int>::size_type j=1; j<index.size(); ++j) {
00133         if (index[j] != index[j-1]+1) {
00134           // not contiguous; short circuit
00135           Teuchos::Array<size_t> inds(index.begin(), index.end());
00136           return mv.subViewNonConst(inds());
00137         }
00138       }
00139       // contiguous
00140       return mv.subViewNonConst(Teuchos::Range1D(index.front(),index.back()));
00141     }
00142 
00148     static Teuchos::RCP<const Tpetra::MultiVector<Scalar,LO,GO,Node> > CloneView( const Tpetra::MultiVector<Scalar,LO,GO,Node>& mv, const std::vector<int>& index )
00149     {
00150       TEST_FOR_EXCEPTION(index.size() == 0,std::invalid_argument,
00151           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::CloneView(mv,index): numvecs must be greater than zero.");
00152 #ifdef TPETRA_DEBUG
00153       TEST_FOR_EXCEPTION( *std::min_element(index.begin(),index.end()) < 0, std::invalid_argument,
00154           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::CloneView(mv,index): indices must be >= zero.");
00155       TEST_FOR_EXCEPTION( *std::max_element(index.begin(),index.end()) >= mv.getNumVectors(), std::invalid_argument,
00156           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::CloneView(mv,index): indices must be < mv.getNumVectors().");
00157 #endif
00158       for (typename std::vector<int>::size_type j=1; j<index.size(); ++j) {
00159         if (index[j] != index[j-1]+1) {
00160           // not contiguous; short circuit
00161           Teuchos::Array<size_t> inds(index.begin(), index.end());
00162           return mv.subView(inds());
00163         }
00164       }
00165       // contiguous
00166       return mv.subView(Teuchos::Range1D(index.front(),index.back()));
00167     }
00168 
00170 
00172 
00173 
00175     static int GetVecLength( const Tpetra::MultiVector<Scalar,LO,GO,Node>& mv )
00176     { return mv.getGlobalLength(); }
00177 
00179     static int GetNumberVecs( const Tpetra::MultiVector<Scalar,LO,GO,Node>& mv )
00180     { return mv.getNumVectors(); }
00182 
00184 
00185 
00188     static void MvTimesMatAddMv( const Scalar alpha, const Tpetra::MultiVector<Scalar,LO,GO,Node>& A, 
00189                                  const Teuchos::SerialDenseMatrix<int,Scalar>& B, 
00190                                  const Scalar beta, Tpetra::MultiVector<Scalar,LO,GO,Node>& mv )
00191     {
00192       Tpetra::Map<LO,GO> LocalMap(B.numRows(), static_cast<GO>(0), A.getMap()->getComm(), Tpetra::LocallyReplicated, A.getMap()->getNode());
00193       Teuchos::ArrayView<const Scalar> Bvalues(B.values(),B.stride()*B.numCols());
00194       Tpetra::MultiVector<Scalar,LO,GO,Node> B_mv(rcpFromRef(LocalMap),Bvalues,B.stride(),B.numCols());
00195       mv.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, alpha, A, B_mv, beta);
00196     }
00197 
00200     static void MvAddMv( const Scalar alpha, const Tpetra::MultiVector<Scalar,LO,GO,Node>& A, const Scalar beta, const Tpetra::MultiVector<Scalar,LO,GO,Node>& B, Tpetra::MultiVector<Scalar,LO,GO,Node>& mv )
00201     {
00202       mv.update(alpha,A,beta,B,Teuchos::ScalarTraits<Scalar>::zero());
00203     }
00204 
00207     static void MvScale ( Tpetra::MultiVector<Scalar,LO,GO,Node>& mv, Scalar alpha )
00208     { mv.scale(alpha); }
00209 
00212     static void MvScale ( Tpetra::MultiVector<Scalar,LO,GO,Node>& mv, const std::vector<Scalar>& alphas )
00213     { mv.scale(alphas); }
00214 
00217     static void MvTransMv(const Scalar alpha, const Tpetra::MultiVector<Scalar,LO,GO,Node>& A, const Tpetra::MultiVector<Scalar,LO,GO,Node>& mv, Teuchos::SerialDenseMatrix<int,Scalar>& B)
00218     { 
00219       Tpetra::Map<LO,GO> LocalMap(B.numRows(), static_cast<GO>(0), A.getMap()->getComm(), Tpetra::LocallyReplicated, A.getMap()->getNode());
00220       Tpetra::MultiVector<Scalar,LO,GO,Node> B_mv(Teuchos::rcpFromRef(LocalMap),B.numCols(),true);
00221       B_mv.multiply(Teuchos::CONJ_TRANS,Teuchos::NO_TRANS,alpha,A,mv,Teuchos::ScalarTraits<Scalar>::zero());
00222       Teuchos::ArrayView<Scalar> av(B.values(),B.stride()*B.numCols());
00223       B_mv.get1dCopy(av,B.stride());
00224     }
00225 
00228     static void MvDot( const Tpetra::MultiVector<Scalar,LO,GO,Node>& A, const Tpetra::MultiVector<Scalar,LO,GO,Node>& B, std::vector<Scalar> &dots)
00229     {
00230       TEST_FOR_EXCEPTION(A.getNumVectors() != B.getNumVectors(),std::invalid_argument,
00231           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::MvDot(A,B,dots): A and B must have the same number of vectors.");
00232 #ifdef TPETRA_DEBUG
00233       TEST_FOR_EXCEPTION(dots.size() < (typename std::vector<int>::size_type)A.getNumVectors(),std::invalid_argument,
00234           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::MvDot(A,B,dots): dots must have room for all dot products.");
00235 #endif
00236       Teuchos::ArrayView<Scalar> av(dots);
00237       A.dot(B,av(0,A.getNumVectors()));
00238     }
00239 
00241 
00243 
00244 
00248     static void MvNorm(const Tpetra::MultiVector<Scalar,LO,GO,Node>& mv, std::vector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> &normvec)
00249     { 
00250 #ifdef TPETRA_DEBUG
00251       TEST_FOR_EXCEPTION(normvec.size() < (typename std::vector<int>::size_type)mv.getNumVectors(),std::invalid_argument,
00252           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::MvNorm(mv,normvec): normvec must have room for all norms.");
00253 #endif
00254       Teuchos::ArrayView<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> av(normvec);
00255       mv.norm2(av(0,mv.getNumVectors()));
00256     }
00257 
00259     
00261 
00262 
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 #ifdef TPETRA_DEBUG
00267       TEST_FOR_EXCEPTION((typename std::vector<int>::size_type)A.getNumVectors() < index.size(),std::invalid_argument,
00268           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::SetBlock(A,index,mv): index must be the same size as A.");
00269 #endif
00270       Teuchos::Array<size_t> inds(index.begin(), index.end());
00271       Teuchos::RCP<Tpetra::MultiVector<Scalar,LO,GO,Node> > mvsub = mv.subViewNonConst(inds());
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 
00284     static void MvRandom( Tpetra::MultiVector<Scalar,LO,GO,Node>& mv )
00285     { mv.randomize(); }
00286 
00289     static void MvInit( Tpetra::MultiVector<Scalar,LO,GO,Node>& mv, Scalar alpha = Teuchos::ScalarTraits<Scalar>::zero() )
00290     { mv.putScalar(alpha); }
00291 
00293 
00295 
00296 
00299     static void MvPrint( const Tpetra::MultiVector<Scalar,LO,GO,Node>& mv, std::ostream& os )
00300     { mv.print(os); }
00301 
00303   };        
00304 
00306   //
00307   // Implementation of the Anasazi::OperatorTraits for Tpetra::Operator.
00308   //
00310 
00311   template <class Scalar, class LO, class GO, class Node>
00312   class OperatorTraits < Scalar, Tpetra::MultiVector<Scalar,LO,GO,Node>, Tpetra::Operator<Scalar,LO,GO,Node> >
00313   {
00314   public:
00315     static void Apply ( const Tpetra::Operator<Scalar,LO,GO,Node> & Op, 
00316                         const Tpetra::MultiVector<Scalar,LO,GO,Node> & X,
00317                               Tpetra::MultiVector<Scalar,LO,GO,Node> & Y)
00318     { 
00319       Op.apply(X,Y,Teuchos::NO_TRANS);
00320     }
00321   };
00322 
00323 } // end of Anasazi namespace 
00324 
00325 #endif 
00326 // end of file ANASAZI_TPETRA_ADAPTER_HPP