Tpetra_CrsMatrixSolveOp_def.hpp

00001 //@HEADER
00002 // ************************************************************************
00003 // 
00004 //               Tpetra: Templated Linear Algebra Services Package 
00005 //                 Copyright (2008) 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 TPETRA_CRSMATRIXSOLVEOP_DEF_HPP
00030 #define TPETRA_CRSMATRIXSOLVEOP_DEF_HPP
00031 
00032 #include "Tpetra_CrsMatrix.hpp"
00033 
00034 namespace Tpetra {
00035 
00036   template <class OpScalar, class MatScalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatVec, class LocalMatSolve>
00037   CrsMatrixSolveOp<OpScalar,MatScalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatVec,LocalMatSolve>::CrsMatrixSolveOp(const Teuchos::RCP<const CrsMatrix<MatScalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatVec,LocalMatSolve> > &A) 
00038   : matrix_(A) {
00039   }
00040 
00041   template <class OpScalar, class MatScalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatVec, class LocalMatSolve>
00042   CrsMatrixSolveOp<OpScalar,MatScalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatVec,LocalMatSolve>::~CrsMatrixSolveOp() {
00043   }
00044 
00045 
00046   template <class OpScalar, class MatScalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatVec, class LocalMatSolve>
00047   void 
00048   CrsMatrixSolveOp<OpScalar,MatScalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatVec,LocalMatSolve>::apply(
00049               const MultiVector<OpScalar,LocalOrdinal,GlobalOrdinal,Node> & X,   
00050                     MultiVector<OpScalar,LocalOrdinal,GlobalOrdinal,Node> & Y,
00051                     Teuchos::ETransp mode, OpScalar alpha, OpScalar beta) const 
00052   {
00053     TEST_FOR_EXCEPTION(!matrix_->isFillComplete(), std::runtime_error, 
00054         Teuchos::typeName(*this) << "::apply(): underlying matrix is not fill-complete.");
00055     TEST_FOR_EXCEPTION(X.getNumVectors() != Y.getNumVectors(), std::runtime_error,
00056         Teuchos::typeName(*this) << "::apply(X,Y): X and Y must have the same number of vectors.");
00057     TEST_FOR_EXCEPTION(matrix_->isLowerTriangular() == false && matrix_->isUpperTriangular() == false, std::runtime_error,
00058         Teuchos::typeName(*this) << "::apply() requires either upper or lower triangular structure in underlying matrix.");
00059     TEST_FOR_EXCEPTION( alpha != Teuchos::ScalarTraits<OpScalar>::one() || beta != Teuchos::ScalarTraits<OpScalar>::zero(), std::runtime_error,
00060         Teuchos::typeName(*this) << "::apply(): non-trivial alpha,beta not supported at this time.");
00061     if (mode == Teuchos::NO_TRANS) {
00062       applyNonTranspose(X,Y);
00063     }
00064     else {
00065       applyTranspose(X,Y);
00066     }
00067   }
00068 
00069   template <class OpScalar, class MatScalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatVec, class LocalMatSolve>
00070   void 
00071   CrsMatrixSolveOp<OpScalar,MatScalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatVec,LocalMatSolve>::applyNonTranspose(
00072       const MultiVector<OpScalar,LocalOrdinal,GlobalOrdinal,Node> & X_in, 
00073             MultiVector<OpScalar,LocalOrdinal,GlobalOrdinal,Node> & Y_in) const 
00074   {
00075     // Solve U X = Y  or  L X = Y
00076     // X belongs to domain map, while Y belongs to range map
00077     typedef Teuchos::ScalarTraits<OpScalar> ST;
00078     using Teuchos::null;
00079 
00080     const size_t numVectors = X_in.getNumVectors();
00081     Teuchos::RCP<const Import<LocalOrdinal,GlobalOrdinal,Node> > importer = matrix_->getGraph()->getImporter();
00082     Teuchos::RCP<const Export<LocalOrdinal,GlobalOrdinal,Node> > exporter = matrix_->getGraph()->getExporter();
00083     Teuchos::RCP<const MV> X;
00084 
00085     // it is okay if X and Y reference the same data, because we can perform a triangular solve in-situ.
00086     // however, we require that column access to each is strided.
00087 
00088     // set up import/export temporary multivectors
00089     if (importer != null) {
00090       if (importMV_ != null && importMV_->getNumVectors() != numVectors) importMV_ = null;
00091       if (importMV_ == null) {
00092         importMV_ = Teuchos::rcp( new MV(matrix_->getColMap(),numVectors) );
00093       }
00094     }
00095     if (exporter != null) {
00096       if (exportMV_ != null && exportMV_->getNumVectors() != numVectors) exportMV_ = null;
00097       if (exportMV_ == null) {
00098         exportMV_ = Teuchos::rcp( new MV(matrix_->getRowMap(),numVectors) );
00099       }
00100     }
00101 
00102     // solve(NO_TRANS): RangeMap -> DomainMap
00103     // lclMatSolve_: RowMap -> ColMap
00104     // importer: DomainMap -> ColMap
00105     // exporter: RowMap -> RangeMap
00106     // 
00107     // solve = reverse(exporter)  o   lclMatSolve_  o reverse(importer)
00108     //         RangeMap   ->    RowMap     ->     ColMap         ->    DomainMap
00109     // If we have a non-trivial exporter, we must import elements that are permuted or are on other processors
00110     if (exporter != null) {
00111       exportMV_->doImport(X_in, *exporter, INSERT);
00112       X = exportMV_;
00113     }
00114     else if (X_in.isConstantStride() == false) {
00115       // cannot handle non-constant stride right now
00116       // generate a copy of X_in
00117       X = Teuchos::rcp(new MV(X_in));
00118     }
00119     else {
00120       // just temporary, so this non-owning RCP is okay
00121       X = Teuchos::rcp( &X_in, false );
00122     }
00123 
00124     // If we have a non-trivial importer, we must export elements that are permuted or belong to other processors
00125     // We will compute solution into the to-be-exported MV
00126     if (importer != null) {
00127       matrix_->template solve<OpScalar,OpScalar>(*X,*importMV_,Teuchos::NO_TRANS);
00128       // Make sure target is zero: necessary because we are adding.
00129       Y_in.putScalar(ST::zero());  
00130       Y_in.doExport(*importMV_, *importer, ADD);
00131     }
00132     // otherwise, solve into Y
00133     else {
00134       // can't solve into non-strided multivector
00135       if (Y_in.isConstantStride() == false) {
00136         // generate a strided copy of Y
00137         MV Y(Y_in);
00138         matrix_->template solve<OpScalar,OpScalar>(*X,Y,Teuchos::NO_TRANS);
00139         Y_in = Y;
00140       }
00141       else {
00142         matrix_->template solve<OpScalar,OpScalar>(*X,Y_in,Teuchos::NO_TRANS);
00143       }
00144     }
00145   }
00146 
00147 
00148   template <class OpScalar, class MatScalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatVec, class LocalMatSolve>
00149   void 
00150   CrsMatrixSolveOp<OpScalar,MatScalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatVec,LocalMatSolve>::applyTranspose(
00151         const MultiVector<OpScalar,LocalOrdinal,GlobalOrdinal,Node> & X_in, 
00152         MultiVector<OpScalar,LocalOrdinal,GlobalOrdinal,Node> &Y_in) const 
00153   {
00154     typedef Teuchos::ScalarTraits<OpScalar> ST;
00155     using Teuchos::null;
00156 
00157     const size_t numVectors = X_in.getNumVectors();
00158     Teuchos::RCP<const Import<LocalOrdinal,GlobalOrdinal,Node> > importer = matrix_->getGraph()->getImporter();
00159     Teuchos::RCP<const Export<LocalOrdinal,GlobalOrdinal,Node> > exporter = matrix_->getGraph()->getExporter();
00160     Teuchos::RCP<const MV> X;
00161 
00162     // it is okay if X and Y reference the same data, because we can perform a triangular solve in-situ.
00163     // however, we require that column access to each is strided.
00164 
00165     // set up import/export temporary multivectors
00166     if (importer != null) {
00167       if (importMV_ != null && importMV_->getNumVectors() != numVectors) importMV_ = null;
00168       if (importMV_ == null) {
00169         importMV_ = Teuchos::rcp( new MV(matrix_->getColMap(),numVectors) );
00170       }
00171     }
00172     if (exporter != null) {
00173       if (exportMV_ != null && exportMV_->getNumVectors() != numVectors) exportMV_ = null;
00174       if (exportMV_ == null) {
00175         exportMV_ = Teuchos::rcp( new MV(matrix_->getRowMap(),numVectors) );
00176       }
00177     }
00178 
00179     // solve(TRANS): DomainMap -> RangeMap
00180     // lclMatSolve_(TRANS): ColMap -> RowMap
00181     // importer: DomainMap -> ColMap
00182     // exporter: RowMap -> RangeMap
00183     // 
00184     // solve = importer o   lclMatSolve_  o  exporter
00185     //         Domainmap -> ColMap     ->      RowMap -> RangeMap
00186     // If we have a non-trivial importer, we must import elements that are permuted or are on other processors
00187     if (importer != null) {
00188       importMV_->doImport(X_in,*importer,INSERT);
00189       X = importMV_;
00190     }
00191     else if (X_in.isConstantStride() == false) {
00192       // cannot handle non-constant stride right now
00193       // generate a copy of X_in
00194       X = Teuchos::rcp(new MV(X_in));
00195     }
00196     else {
00197       // just temporary, so this non-owning RCP is okay
00198       X = Teuchos::rcp( &X_in, false );
00199     }
00200 
00201 
00202     // If we have a non-trivial exporter, we must export elements that are permuted or belong to other processors
00203     // We will compute solution into the to-be-exported MV; get a view
00204     if (exporter != null) {
00205       matrix_->template solve<OpScalar,OpScalar>(*X,*exportMV_,Teuchos::CONJ_TRANS);
00206       // Make sure target is zero: necessary because we are adding
00207       Y_in.putScalar(ST::zero());  
00208       Y_in.doExport(*importMV_, *importer, ADD);
00209     }
00210     // otherwise, solve into Y
00211     else {
00212       if (Y_in.isConstantStride() == false) {
00213         // generate a strided copy of Y
00214         MV Y(Y_in);
00215         matrix_->template solve<OpScalar,OpScalar>(*X,Y,Teuchos::CONJ_TRANS);
00216         Y_in = Y;
00217       }
00218       else {
00219         matrix_->template solve<OpScalar,OpScalar>(*X,Y_in,Teuchos::CONJ_TRANS);
00220       }
00221     }
00222   }
00223 
00224   template <class OpScalar, class MatScalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatVec, class LocalMatSolve>
00225   bool 
00226   CrsMatrixSolveOp<OpScalar,MatScalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatVec,LocalMatSolve>::hasTransposeApply() const {
00227     return true;
00228   }
00229 
00230   template <class OpScalar, class MatScalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatVec, class LocalMatSolve>
00231   const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & 
00232   CrsMatrixSolveOp<OpScalar,MatScalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatVec,LocalMatSolve>::getDomainMap() const {
00233     return matrix_->getRangeMap();
00234   }
00235 
00236   template <class OpScalar, class MatScalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatVec, class LocalMatSolve>
00237   const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & 
00238   CrsMatrixSolveOp<OpScalar,MatScalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatVec,LocalMatSolve>::getRangeMap() const {
00239     return matrix_->getDomainMap();
00240   }
00241 
00242 } // end of namespace Tpetra
00243 
00245 template <class OpScalar, class MatScalar, class LocalOrdinal, class GlobalOrdinal, class Node, class LocalMatVec, class LocalMatSolve>
00246 Teuchos::RCP< Tpetra::CrsMatrixSolveOp<OpScalar,MatScalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatVec,LocalMatSolve> >
00247 Tpetra::createCrsMatrixSolveOp(const Teuchos::RCP<const Tpetra::CrsMatrix<MatScalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatVec,LocalMatSolve> > &A) {
00248   return rcp(new Tpetra::CrsMatrixSolveOp<OpScalar,MatScalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatVec,LocalMatSolve>(A) );
00249 }
00250 
00251 //
00252 // Explicit instantiation macro
00253 //
00254 // Must be expanded from within the Tpetra namespace!
00255 //
00256 
00257 #define TPETRA_CRSMATRIX_SOLVEOP_INSTANT(OPSCALAR,MATSCALAR,LO,GO,NODE) \
00258   \
00259   template class CrsMatrixSolveOp< OPSCALAR , MATSCALAR , LO , GO , NODE >; \
00260   \
00261   template Teuchos::RCP< Tpetra::CrsMatrixSolveOp<OPSCALAR,MATSCALAR,LO,GO,NODE> > \
00262   createCrsMatrixSolveOp(const Teuchos::RCP<const Tpetra::CrsMatrix<MATSCALAR,LO,GO,NODE> > &A); \
00263   \
00264   template void CrsMatrix<MATSCALAR,LO,GO,NODE>::solve<OPSCALAR,OPSCALAR>( \
00265         const MultiVector<OPSCALAR,LO,GO,NODE> &X, \
00266               MultiVector<OPSCALAR,LO,GO,NODE> &Y, \
00267               Teuchos::ETransp mode) const; \
00268 
00269 
00270 #endif // TPETRA_CRSMATRIXSOLVEOP_DEF_HPP

Generated on Tue Jul 13 09:39:06 2010 for Tpetra Matrix/Vector Services by  doxygen 1.4.7