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