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 #include <Teuchos_Array.hpp>
00043 
00044 #include "AnasaziConfigDefs.hpp"
00045 #include "AnasaziTypes.hpp"
00046 #include "AnasaziMultiVecTraits.hpp"
00047 #include "AnasaziOperatorTraits.hpp"
00048 
00049 namespace Anasazi {
00050 
00052   //
00053   // Implementation of the Anasazi::MultiVecTraits for Tpetra::MultiVector.
00054   //
00056 
00061   template<class Scalar, class LO, class GO, class Node>
00062   class MultiVecTraits<Scalar, Tpetra::MultiVector<Scalar,LO,GO,Node> >
00063   {
00064   public:
00065 
00066     static Teuchos::RCP<Tpetra::MultiVector<Scalar,LO,GO,Node> > Clone( const Tpetra::MultiVector<Scalar,LO,GO,Node>& mv, const int numvecs )
00067     { 
00068       return Teuchos::rcp( new Tpetra::MultiVector<Scalar,LO,GO,Node>(mv.getMap(),numvecs));
00069     }
00070 
00071     static Teuchos::RCP<Tpetra::MultiVector<Scalar,LO,GO,Node> > CloneCopy( const Tpetra::MultiVector<Scalar,LO,GO,Node>& mv )
00072     {
00073       return Teuchos::rcp( new Tpetra::MultiVector<Scalar,LO,GO,Node>( mv ) ); 
00074     }
00075 
00076     static Teuchos::RCP<Tpetra::MultiVector<Scalar,LO,GO,Node> > CloneCopy( const Tpetra::MultiVector<Scalar,LO,GO,Node>& mv, const std::vector<int>& index )
00077     { 
00078       TEST_FOR_EXCEPTION(index.size() == 0,std::runtime_error,
00079           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::Clone(mv,index): numvecs must be greater than zero.");
00080 #ifdef HAVE_TPETRA_DEBUG
00081       TEST_FOR_EXCEPTION( *std::min_element(index.begin(),index.end()) < 0, std::runtime_error,
00082           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::Clone(mv,index): indices must be >= zero.");
00083       TEST_FOR_EXCEPTION( (size_t)*std::max_element(index.begin(),index.end()) >= mv.getNumVectors(), std::runtime_error,
00084           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::Clone(mv,index): indices must be < mv.getNumVectors().");
00085 #endif
00086       for (typename std::vector<int>::size_type j=1; j<index.size(); ++j) {
00087         if (index[j] != index[j-1]+1) {
00088           // not contiguous; short circuit
00089           Teuchos::Array<size_t> stinds(index.begin(), index.end());
00090           return mv.subCopy(stinds());
00091         }
00092       }
00093       // contiguous
00094       return mv.subCopy(Teuchos::Range1D(index.front(),index.back()));
00095     }
00096 
00097     static Teuchos::RCP<Tpetra::MultiVector<Scalar,LO,GO,Node> > CloneViewNonConst( Tpetra::MultiVector<Scalar,LO,GO,Node>& mv, const std::vector<int>& index )
00098     {
00099       TEST_FOR_EXCEPTION(index.size() == 0,std::invalid_argument,
00100           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::CloneViewNonConst(mv,index): numvecs must be greater than zero.");
00101 #ifdef HAVE_TPETRA_DEBUG
00102       TEST_FOR_EXCEPTION( *std::min_element(index.begin(),index.end()) < 0, std::invalid_argument,
00103           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::CloneViewNonConst(mv,index): indices must be >= zero.");
00104       TEST_FOR_EXCEPTION( (size_t)*std::max_element(index.begin(),index.end()) >= mv.getNumVectors(), std::invalid_argument,
00105           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::CloneViewNonConst(mv,index): indices must be < mv.getNumVectors().");
00106 #endif
00107       for (typename std::vector<int>::size_type j=1; j<index.size(); ++j) {
00108         if (index[j] != index[j-1]+1) {
00109           // not contiguous; short circuit
00110           Teuchos::Array<size_t> stinds(index.begin(), index.end());
00111           return mv.subViewNonConst(stinds);
00112         }
00113       }
00114       // contiguous
00115       return mv.subViewNonConst(Teuchos::Range1D(index.front(),index.back()));
00116     }
00117 
00118     static Teuchos::RCP<const Tpetra::MultiVector<Scalar,LO,GO,Node> > CloneView(const Tpetra::MultiVector<Scalar,LO,GO,Node>& mv, const std::vector<int>& index )
00119     {
00120       TEST_FOR_EXCEPTION(index.size() == 0,std::invalid_argument,
00121           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::CloneView(mv,index): numvecs must be greater than zero.");
00122 #ifdef HAVE_TPETRA_DEBUG
00123       TEST_FOR_EXCEPTION( *std::min_element(index.begin(),index.end()) < 0, std::invalid_argument,
00124           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::CloneView(mv,index): indices must be >= zero.");
00125       TEST_FOR_EXCEPTION( (size_t)*std::max_element(index.begin(),index.end()) >= mv.getNumVectors(), std::invalid_argument,
00126           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::CloneView(mv,index): indices must be < mv.getNumVectors().");
00127 #endif
00128       for (typename std::vector<int>::size_type j=1; j<index.size(); ++j) {
00129         if (index[j] != index[j-1]+1) {
00130           // not contiguous; short circuit
00131           Teuchos::Array<size_t> stinds(index.begin(), index.end());
00132           return mv.subView(stinds);
00133         }
00134       }
00135       // contiguous
00136       return mv.subView(Teuchos::Range1D(index.front(),index.back()));
00137     }
00138 
00139     static int GetVecLength( const Tpetra::MultiVector<Scalar,LO,GO,Node>& mv )
00140     { return mv.getGlobalLength(); }
00141 
00142     static int GetNumberVecs( const Tpetra::MultiVector<Scalar,LO,GO,Node>& mv )
00143     { return mv.getNumVectors(); }
00144 
00145     static void MvTimesMatAddMv( Scalar alpha, const Tpetra::MultiVector<Scalar,LO,GO,Node>& A, 
00146                                  const Teuchos::SerialDenseMatrix<int,Scalar>& B, 
00147                                  Scalar beta, Tpetra::MultiVector<Scalar,LO,GO,Node>& mv )
00148     {
00149       // create local map
00150       Tpetra::Map<LO,GO,Node> LocalMap(B.numRows(), 0, A.getMap()->getComm(), Tpetra::LocallyReplicated);
00151       // encapsulate Teuchos::SerialDenseMatrix data in ArrayView
00152       Teuchos::ArrayView<const Scalar> Bvalues(B.values(),B.stride()*B.numCols());
00153       // create locally replicated MultiVector with a copy of this data
00154       Tpetra::MultiVector<Scalar,LO,GO,Node> B_mv(Teuchos::rcpFromRef(LocalMap),Bvalues,B.stride(),B.numCols());
00155       // multiply
00156       mv.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, alpha, A, B_mv, beta);
00157     }
00158 
00159     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 )
00160     {
00161       mv.update(alpha,A,beta,B,Teuchos::ScalarTraits<Scalar>::zero());
00162     }
00163 
00164     static void MvScale ( Tpetra::MultiVector<Scalar,LO,GO,Node>& mv, Scalar alpha )
00165     { mv.scale(alpha); }
00166 
00167     static void MvScale ( Tpetra::MultiVector<Scalar,LO,GO,Node>& mv, const std::vector<Scalar>& alphas )
00168     { mv.scale(alphas); }
00169 
00170     static void MvTransMv( Scalar alpha, const Tpetra::MultiVector<Scalar,LO,GO,Node>& A, const Tpetra::MultiVector<Scalar,LO,GO,Node>& mv, Teuchos::SerialDenseMatrix<int,Scalar>& B)
00171     { 
00172       // create local map
00173       Tpetra::Map<LO,GO,Node> LocalMap(B.numRows(), 0, A.getMap()->getComm(), Tpetra::LocallyReplicated);
00174       // create local multivector to hold the result
00175       Tpetra::MultiVector<Scalar,LO,GO,Node> B_mv(Teuchos::rcpFromRef(LocalMap),B.numCols(),true);
00176       // multiply result into local multivector
00177       B_mv.multiply(Teuchos::CONJ_TRANS,Teuchos::NO_TRANS,alpha,A,mv,Teuchos::ScalarTraits<Scalar>::zero());
00178       // create arrayview encapsulating the Teuchos::SerialDenseMatrix
00179       Teuchos::ArrayView<Scalar> av(B.values(),B.stride()*B.numCols());
00180       // extract a copy of the result into the array view (and therefore, the SerialDenseMatrix)
00181       B_mv.get1dCopy(av,B.stride());
00182     }
00183 
00184     static void MvDot( const Tpetra::MultiVector<Scalar,LO,GO,Node>& A, const Tpetra::MultiVector<Scalar,LO,GO,Node>& B, std::vector<Scalar> &dots)
00185     {
00186       TEST_FOR_EXCEPTION(A.getNumVectors() != B.getNumVectors(),std::invalid_argument,
00187           "Belos::MultiVecTraits<Scalar,Tpetra::MultiVector>::MvDot(A,B,dots): A and B must have the same number of vectors.");
00188 #ifdef HAVE_TPETRA_DEBUG
00189       TEST_FOR_EXCEPTION(dots.size() < (typename std::vector<int>::size_type)A.getNumVectors(),std::invalid_argument,
00190           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::MvDot(A,B,dots): dots must have room for all dot products.");
00191 #endif
00192       Teuchos::ArrayView<Scalar> av(dots);
00193       A.dot(B,av(0,A.getNumVectors()));
00194     }
00195 
00196     static void MvNorm(const Tpetra::MultiVector<Scalar,LO,GO,Node>& mv, std::vector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> &normvec)
00197     { 
00198 #ifdef HAVE_TPETRA_DEBUG
00199       TEST_FOR_EXCEPTION(normvec.size() < (typename std::vector<int>::size_type)mv.getNumVectors(),std::invalid_argument,
00200           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::MvNorm(mv,normvec): normvec must have room for all norms.");
00201 #endif
00202       Teuchos::ArrayView<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> av(normvec);
00203       mv.norm2(av(0,mv.getNumVectors()));
00204     }
00205 
00206     static void SetBlock( const Tpetra::MultiVector<Scalar,LO,GO,Node>& A, const std::vector<int>& index, Tpetra::MultiVector<Scalar,LO,GO,Node>& mv )
00207     {
00208 #ifdef HAVE_TPETRA_DEBUG
00209       TEST_FOR_EXCEPTION((typename std::vector<int>::size_type)A.getNumVectors() < index.size(),std::invalid_argument,
00210           "Anasazi::MultiVecTraits<Scalar,Tpetra::MultiVector>::SetBlock(A,index,mv): index must be the same size as A.");
00211 #endif
00212       Teuchos::RCP<Tpetra::MultiVector<Scalar,LO,GO,Node> > mvsub = CloneViewNonConst(mv,index);
00213       if ((typename std::vector<int>::size_type)A.getNumVectors() > index.size()) {
00214         Teuchos::RCP<const Tpetra::MultiVector<Scalar,LO,GO,Node> > Asub = A.subView(Teuchos::Range1D(0,index.size()-1));
00215         (*mvsub) = (*Asub);
00216       }
00217       else {
00218         (*mvsub) = A;
00219       }
00220       mvsub = Teuchos::null;
00221     }
00222 
00223     static void MvRandom( Tpetra::MultiVector<Scalar,LO,GO,Node>& mv )
00224     { mv.randomize(); }
00225 
00226     static void MvInit( Tpetra::MultiVector<Scalar,LO,GO,Node>& mv, Scalar alpha = Teuchos::ScalarTraits<Scalar>::zero() )
00227     { mv.putScalar(alpha); }
00228 
00229     static void MvPrint( const Tpetra::MultiVector<Scalar,LO,GO,Node>& mv, std::ostream& os )
00230     { mv.print(os); }
00231 
00232   };        
00233 
00235   //
00236   // Implementation of the Anasazi::OperatorTraits for Tpetra::Operator.
00237   //
00239 
00240   template <class Scalar, class LO, class GO, class Node> 
00241   class OperatorTraits < Scalar, Tpetra::MultiVector<Scalar,LO,GO,Node>, Tpetra::Operator<Scalar,LO,GO,Node> >
00242   {
00243   public:
00244     static void Apply ( const Tpetra::Operator<Scalar,LO,GO,Node> & Op, 
00245                         const Tpetra::MultiVector<Scalar,LO,GO,Node> & X,
00246                               Tpetra::MultiVector<Scalar,LO,GO,Node> & Y)
00247     { 
00248       Op.apply(X,Y,Teuchos::NO_TRANS);
00249     }
00250   };
00251 
00252 } // end of Anasazi namespace 
00253 
00254 #endif 
00255 // end of file ANASAZI_TPETRA_ADAPTER_HPP