Tpetra Matrix/Vector Services Version of the Day
Tpetra_CrsMatrix_decl.hpp
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //          Tpetra: Templated Linear Algebra Services Package
00005 //                 Copyright (2008) Sandia Corporation
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ************************************************************************
00040 // @HEADER
00041 
00042 #ifndef TPETRA_CRSMATRIX_DECL_HPP
00043 #define TPETRA_CRSMATRIX_DECL_HPP
00044 
00045 // TODO: row-wise insertion of entries in globalAssemble() may be more efficient
00046 
00047 // TODO: add typeglobs: CrsMatrix<Scalar,typeglob>
00048 // TODO: add template (template) parameter for nonlocal container (this will be part of typeglob)
00049 
00050 #include <Kokkos_DefaultNode.hpp>
00051 #include <Kokkos_DefaultKernels.hpp>
00052 
00053 #include "Tpetra_ConfigDefs.hpp"
00054 #include "Tpetra_RowMatrix.hpp"
00055 #include "Tpetra_DistObject.hpp"
00056 #include "Tpetra_CrsGraph.hpp"
00057 #include "Tpetra_Vector.hpp"
00058 #include "Tpetra_CrsMatrixMultiplyOp_decl.hpp"
00059 
00060 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00061 namespace Tpetra {
00062   // struct for i,j,v triplets
00063   template <class Ordinal, class Scalar>
00064   struct CrsIJV {
00065     CrsIJV();
00066     CrsIJV(Ordinal row, Ordinal col, const Scalar &val);
00067     Ordinal i,j;
00068     Scalar  v;
00069   };
00070 }
00071 
00072 namespace Teuchos {
00073   // SerializationTraits specialization for CrsIJV, using DirectSerialization
00074   template <typename Ordinal, typename Scalar>
00075   class SerializationTraits<int,Tpetra::CrsIJV<Ordinal,Scalar> >
00076   : public DirectSerializationTraits<int,Tpetra::CrsIJV<Ordinal,Scalar> >
00077   {};
00078 }
00079 
00080 namespace std {
00081   template <class Ordinal, class Scalar>
00082   bool operator<(const Tpetra::CrsIJV<Ordinal,Scalar> &ijv1, const Tpetra::CrsIJV<Ordinal,Scalar> &ijv2);
00083 }
00084 #endif
00085 
00086 namespace Tpetra {
00087 
00089 
00159   template <class Scalar,
00160             class LocalOrdinal  = int,
00161             class GlobalOrdinal = LocalOrdinal,
00162             class Node          = Kokkos::DefaultNode::DefaultNodeType,
00163             class LocalMatOps   = typename Kokkos::DefaultKernels<Scalar,LocalOrdinal,Node>::SparseOps >
00164   class CrsMatrix : public RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>,
00165                     public DistObject<char, LocalOrdinal,GlobalOrdinal,Node> {
00166   public:
00167     typedef Scalar                                scalar_type;
00168     typedef LocalOrdinal                          local_ordinal_type;
00169     typedef GlobalOrdinal                         global_ordinal_type;
00170     typedef Node                                  node_type;
00171     typedef Map<LocalOrdinal,GlobalOrdinal,Node>  map_type;
00172     // backwards compatibility defines both of these
00173     typedef LocalMatOps   mat_vec_type;
00174     typedef LocalMatOps   mat_solve_type;
00175 
00177 
00178 
00196     CrsMatrix (const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& rowMap,
00197                size_t maxNumEntriesPerRow,
00198                ProfileType pftype = DynamicProfile,
00199                const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null);
00200 
00218     CrsMatrix (const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& rowMap,
00219                const ArrayRCP<const size_t>& NumEntriesPerRowToAlloc,
00220                ProfileType pftype = DynamicProfile,
00221                const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null);
00222 
00245     CrsMatrix (const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& rowMap,
00246                const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& colMap,
00247                size_t maxNumEntriesPerRow,
00248                ProfileType pftype = DynamicProfile,
00249                const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null);
00250 
00273     CrsMatrix (const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& rowMap,
00274                const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& colMap,
00275                const ArrayRCP<const size_t>& NumEntriesPerRowToAlloc,
00276                ProfileType pftype = DynamicProfile,
00277                const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null);
00278 
00299     explicit CrsMatrix (const Teuchos::RCP<const CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps> >& graph,
00300                         const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null);
00301 
00303     virtual ~CrsMatrix();
00304 
00306 
00307 
00308 
00310 
00319     void insertGlobalValues(GlobalOrdinal globalRow, const ArrayView<const GlobalOrdinal> &cols, const ArrayView<const Scalar> &vals);
00320 
00322 
00333     void insertLocalValues(LocalOrdinal localRow, const ArrayView<const LocalOrdinal> &cols, const ArrayView<const Scalar> &vals);
00334 
00336 
00341     void replaceGlobalValues(GlobalOrdinal globalRow,
00342                              const ArrayView<const GlobalOrdinal> &cols,
00343                              const ArrayView<const Scalar>        &vals);
00344 
00346 
00348     void replaceLocalValues(LocalOrdinal localRow,
00349                             const ArrayView<const LocalOrdinal> &cols,
00350                             const ArrayView<const Scalar>       &vals);
00351 
00353 
00358     void sumIntoGlobalValues(GlobalOrdinal globalRow,
00359                              const ArrayView<const GlobalOrdinal> &cols,
00360                              const ArrayView<const Scalar>        &vals);
00361 
00362 
00364 
00369     void sumIntoLocalValues(LocalOrdinal globalRow,
00370                             const ArrayView<const LocalOrdinal>  &cols,
00371                             const ArrayView<const Scalar>        &vals);
00372 
00374     void setAllToScalar(const Scalar &alpha);
00375 
00377     void scale(const Scalar &alpha);
00378 
00380 
00381 
00382 
00389     void globalAssemble();
00390 
00399     void resumeFill(const RCP<ParameterList> &params = null);
00400 
00411     void fillComplete(const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &domainMap,
00412                       const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rangeMap,
00413                       const RCP<ParameterList> &params = null);
00414 
00427     void fillComplete(const RCP<ParameterList> &params = null);
00428 
00430 
00431 
00432 
00434     const RCP<const Comm<int> > & getComm() const;
00435 
00437     RCP<Node> getNode() const;
00438 
00440     const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & getRowMap() const;
00441 
00443     const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & getColMap() const;
00444 
00446     RCP<const RowGraph<LocalOrdinal,GlobalOrdinal,Node> > getGraph() const;
00447 
00449     RCP<const CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps> > getCrsGraph() const;
00450 
00470     global_size_t getGlobalNumRows() const;
00471 
00476     global_size_t getGlobalNumCols() const;
00477 
00479     size_t getNodeNumRows() const;
00480 
00482 
00484     size_t getNodeNumCols() const;
00485 
00487     GlobalOrdinal getIndexBase() const;
00488 
00490     global_size_t getGlobalNumEntries() const;
00491 
00493     size_t getNodeNumEntries() const;
00494 
00496 
00497     size_t getNumEntriesInGlobalRow(GlobalOrdinal globalRow) const;
00498 
00500 
00501     size_t getNumEntriesInLocalRow(LocalOrdinal localRow) const;
00502 
00504 
00506     global_size_t getGlobalNumDiags() const;
00507 
00509 
00511     size_t getNodeNumDiags() const;
00512 
00514 
00516     size_t getGlobalMaxNumRowEntries() const;
00517 
00519 
00521     size_t getNodeMaxNumRowEntries() const;
00522 
00524     bool hasColMap() const;
00525 
00527 
00529     bool isLowerTriangular() const;
00530 
00532 
00534     bool isUpperTriangular() const;
00535 
00537     bool isLocallyIndexed() const;
00538 
00540     bool isGloballyIndexed() const;
00541 
00543     bool isFillComplete() const;
00544 
00546     bool isFillActive() const;
00547 
00549 
00555     bool isStorageOptimized() const;
00556 
00558     ProfileType getProfileType() const;
00559 
00561     bool isStaticGraph() const;
00562 
00564 
00570     typename ScalarTraits<Scalar>::magnitudeType getFrobeniusNorm() const;
00571 
00573 
00583     void getGlobalRowCopy(GlobalOrdinal GlobalRow,
00584                           const ArrayView<GlobalOrdinal> &Indices,
00585                           const ArrayView<Scalar> &Values,
00586                           size_t &NumEntries
00587                           ) const;
00588 
00590 
00602     void getLocalRowCopy(LocalOrdinal LocalRow,
00603                          const ArrayView<LocalOrdinal> &Indices,
00604                          const ArrayView<Scalar> &Values,
00605                          size_t &NumEntries
00606                          ) const;
00607 
00609 
00618     void getGlobalRowView(GlobalOrdinal GlobalRow, ArrayView<const GlobalOrdinal> &indices, ArrayView<const Scalar> &values) const;
00619 
00621 
00630     void getLocalRowView(LocalOrdinal LocalRow, ArrayView<const LocalOrdinal> &indices, ArrayView<const Scalar> &values) const;
00631 
00633 
00635     void getLocalDiagCopy(Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &diag) const;
00636 
00638     void leftScale(const Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& x);
00639 
00641     void rightScale(const Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& x);
00642 
00644 
00645 
00646 
00648 
00658     template <class DomainScalar, class RangeScalar>
00659     void
00660     localMultiply (const MultiVector<DomainScalar,LocalOrdinal,GlobalOrdinal,Node>& X,
00661                    MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,Node>& Y,
00662                    Teuchos::ETransp trans,
00663                    RangeScalar alpha,
00664                    RangeScalar beta) const;
00665 
00682     template <class DomainScalar, class RangeScalar>
00683     void
00684     localSolve (const MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,Node>& Y,
00685                 MultiVector<DomainScalar,LocalOrdinal,GlobalOrdinal,Node>& X,
00686                 Teuchos::ETransp trans) const;
00687 
00689     template <class T>
00690     RCP<CrsMatrix<T,LocalOrdinal,GlobalOrdinal,Node> > convert() const;
00691 
00693 
00694 
00695 
00697 
00700     void apply(const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> & X, MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &Y,
00701                Teuchos::ETransp mode = Teuchos::NO_TRANS,
00702                Scalar alpha = ScalarTraits<Scalar>::one(),
00703                Scalar beta = ScalarTraits<Scalar>::zero()) const;
00704 
00706     bool hasTransposeApply() const;
00707 
00710     const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & getDomainMap() const;
00711 
00714     const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & getRangeMap() const;
00715 
00717 
00718 
00719 
00721     std::string description() const;
00722 
00724     void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const;
00725 
00727 
00728 
00729 
00730     bool checkSizes(const DistObject<char, LocalOrdinal,GlobalOrdinal,Node>& source);
00731 
00732     void copyAndPermute(const DistObject<char, LocalOrdinal,GlobalOrdinal,Node>& source,
00733                         size_t numSameIDs,
00734                         const ArrayView<const LocalOrdinal> &permuteToLIDs,
00735                         const ArrayView<const LocalOrdinal> &permuteFromLIDs);
00736 
00737     void packAndPrepare(const DistObject<char, LocalOrdinal,GlobalOrdinal,Node>& source,
00738                         const ArrayView<const LocalOrdinal> &exportLIDs,
00739                         Array<char> &exports,
00740                         const ArrayView<size_t> & numPacketsPerLID,
00741                         size_t& constantNumPackets,
00742                         Distributor &distor);
00743 
00753     void
00754     unpackAndCombine (const Teuchos::ArrayView<const LocalOrdinal> &importLIDs,
00755                       const Teuchos::ArrayView<const char> &imports,
00756                       const Teuchos::ArrayView<size_t> &numPacketsPerLID,
00757                       size_t constantNumPackets,
00758                       Distributor &distor,
00759                       CombineMode combineMode);
00761 
00762 
00763 
00772     TPETRA_DEPRECATED void optimizeStorage();
00773 
00775     TPETRA_DEPRECATED void getGlobalRowView(GlobalOrdinal GlobalRow, ArrayRCP<const GlobalOrdinal> &indices, ArrayRCP<const Scalar> &values) const;
00776 
00778     TPETRA_DEPRECATED void getLocalRowView(LocalOrdinal LocalRow, ArrayRCP<const LocalOrdinal> &indices, ArrayRCP<const Scalar> &values) const;
00779 
00781     template <class DomainScalar, class RangeScalar>
00782     TPETRA_DEPRECATED
00783     void multiply(const MultiVector<DomainScalar,LocalOrdinal,GlobalOrdinal,Node> & X, MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,Node> &Y, Teuchos::ETransp trans, RangeScalar alpha, RangeScalar beta) const;
00784 
00786     template <class DomainScalar, class RangeScalar>
00787     TPETRA_DEPRECATED
00788     void solve(const MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,Node> & Y, MultiVector<DomainScalar,LocalOrdinal,GlobalOrdinal,Node> &X, Teuchos::ETransp trans) const;
00789 
00791     TPETRA_DEPRECATED void fillComplete(const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &domainMap, const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rangeMap, OptimizeOption os);
00792 
00794     TPETRA_DEPRECATED void fillComplete(OptimizeOption os);
00795 
00797 
00798   private:
00799     // We forbid copy construction by declaring this method private
00800     // and not implementing it.
00801     CrsMatrix (const CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps> &rhs);
00802 
00803     // We forbid assignment (operator=) by declaring this method
00804     // private and not implementing it.
00805     CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>&
00806     operator= (const CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps> &rhs);
00807 
00815     void
00816     combineGlobalValues (const GlobalOrdinal globalRowIndex,
00817                          const Teuchos::ArrayView<const GlobalOrdinal> columnIndices,
00818                          const Teuchos::ArrayView<const Scalar> values,
00819                          const Tpetra::CombineMode combineMode);
00820 
00835     template<class BinaryFunction>
00836     void
00837     transformGlobalValues (GlobalOrdinal globalRow,
00838                            const Teuchos::ArrayView<const GlobalOrdinal>& indices,
00839                            const Teuchos::ArrayView<const Scalar>        & values,
00840                            BinaryFunction f)
00841     {
00842       typedef LocalOrdinal LO;
00843       typedef GlobalOrdinal GO;
00844       typedef Node NT;
00845       using Teuchos::Array;
00846       using Teuchos::ArrayView;
00847 
00848       TEUCHOS_TEST_FOR_EXCEPTION(values.size() != indices.size(),
00849         std::invalid_argument, "transformGlobalValues: values.size() = "
00850         << values.size() << " != indices.size() = " << indices.size() << ".");
00851 
00852       const LO lrow = this->getRowMap()->getLocalElement(globalRow);
00853 
00854       TEUCHOS_TEST_FOR_EXCEPTION(lrow == LOT::invalid(), std::invalid_argument,
00855         "transformGlobalValues: The given global row index " << globalRow
00856         << " is not owned by the calling process (rank "
00857         << this->getRowMap()->getComm()->getRank() << ").");
00858 
00859       RowInfo rowInfo = staticGraph_->getRowInfo(lrow);
00860       if (indices.size() > 0) {
00861         if (isLocallyIndexed()) {
00862           // Convert global indices to local indices.
00863           const Map<LO, GO, NT> &colMap = *(this->getColMap());
00864           Array<LO> lindices (indices.size());
00865           typename ArrayView<const GO>::iterator gindit = indices.begin();
00866           typename Array<LO>::iterator           lindit = lindices.begin();
00867           while (gindit != indices.end()) {
00868             // No need to filter before asking the column Map to
00869             // convert GID->LID.  If the GID doesn't exist in the
00870             // column Map, the GID will be mapped to invalid(), which
00871             // will not be found in the graph.
00872             *lindit++ = colMap.getLocalElement(*gindit++);
00873           }
00874           typename Graph::SLocalGlobalViews inds_view;
00875           inds_view.linds = lindices();
00876           staticGraph_->template transformValues<LocalIndices>(rowInfo, inds_view, this->getViewNonConst(rowInfo).begin(), values.begin(), f);
00877         }
00878         else if (isGloballyIndexed()) {
00879           typename Graph::SLocalGlobalViews inds_view;
00880           inds_view.ginds = indices;
00881           staticGraph_->template transformValues<GlobalIndices>(rowInfo, inds_view, this->getViewNonConst(rowInfo).begin(), values.begin(), f);
00882         }
00883       }
00884     }
00885 
00886   protected:
00887     // useful typedefs
00888     typedef OrdinalTraits<LocalOrdinal>                     LOT;
00889     typedef OrdinalTraits<GlobalOrdinal>                    GOT;
00890     typedef ScalarTraits<Scalar>                             ST;
00891     typedef typename ST::magnitudeType                Magnitude;
00892     typedef MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>       MV;
00893     typedef Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>             V;
00894     typedef CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>  Graph;
00895     typedef typename LocalMatOps::template bind_scalar<Scalar>::other_type                    sparse_ops_type;
00896     typedef typename sparse_ops_type::template graph<LocalOrdinal,Node>::graph_type          local_graph_type;
00897     typedef typename sparse_ops_type::template matrix<Scalar,LocalOrdinal,Node>::matrix_type local_matrix_type;
00898     // Enums
00899     enum GraphAllocationStatus {
00900       GraphAlreadyAllocated,
00901       GraphNotYetAllocated
00902     };
00903 
00920     void allocateValues (ELocalGlobal lg, GraphAllocationStatus gas);
00921 
00927     void sortEntries();
00928 
00934     void mergeRedundantEntries();
00935 
00943     void clearGlobalConstants();
00944 
00953     void computeGlobalConstants();
00954 
00955     // matrix data accessors
00956     ArrayView<const Scalar>    getView(RowInfo rowinfo) const;
00957     ArrayView<      Scalar>    getViewNonConst(RowInfo rowinfo);
00958     // local Kokkos objects
00959     void fillLocalMatrix(const RCP<ParameterList> &params);
00960     void fillLocalGraphAndMatrix(const RCP<ParameterList> &params);
00961     // debugging
00962     void checkInternalState() const;
00963 
00964     // Two graph pointers needed in order to maintain const-correctness:
00965     // staticGraph_ is a graph passed to the constructor. We are not allowed to modify it. it is always a valid pointer.
00966     // myGraph_     is a graph created here. We are allowed to modify it. if myGraph_ != null, then staticGraph_ = myGraph_
00967     RCP<const Graph> staticGraph_;
00968     RCP<      Graph>     myGraph_;
00969 
00970     RCP<sparse_ops_type>   lclMatOps_;
00971     RCP<local_matrix_type> lclMatrix_;
00972 
00973     // matrix values. before allocation, both are null.
00974     // after allocation, one is null.
00975     // 1D == StaticAllocation, 2D == DynamicAllocation
00976     // The allocation always matches that of graph_, as the graph does the allocation for the matrix.
00977     ArrayRCP<Scalar>                       values1D_;
00978     ArrayRCP<ArrayRCP<Scalar> >            values2D_;
00979     // TODO: these could be allocated at resumeFill() and de-allocated at fillComplete() to make for very fast getView()/getViewNonConst()
00980     // ArrayRCP< typedef ArrayRCP<const Scalar>::iterator > rowPtrs_;
00981     // ArrayRCP< typedef ArrayRCP<      Scalar>::iterator > rowPtrsNC_;
00982 
00984     bool fillComplete_;
00985 
00998     std::map<GlobalOrdinal, Array<std::pair<GlobalOrdinal,Scalar> > > nonlocals_;
00999 
01000     // a wrapper around multiply, for use in apply; it contains a non-owning RCP to *this, therefore, it is not allowed
01001     // to persist past the destruction of *this. therefore, WE MAY NOT SHARE THIS POINTER.
01002     RCP< const CrsMatrixMultiplyOp<Scalar,Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps> > sameScalarMultiplyOp_;
01003 
01004     // cached frobenius norm: -ST::one() means invalid
01005     mutable Magnitude frobNorm_;
01006 
01007   }; // class CrsMatrix
01008 
01015   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
01016   Teuchos::RCP<CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
01017   createCrsMatrix (const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &map,
01018                    size_t maxNumEntriesPerRow = 0,
01019                    const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null)
01020   {
01021     using Teuchos::rcp;
01022     typedef CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> matrix_type;
01023 
01024     return rcp (new matrix_type (map, maxNumEntriesPerRow, DynamicProfile, params));
01025   }
01026 
01076   template<class CrsMatrixType>
01077   Teuchos::RCP<CrsMatrixType>
01078   importAndFillCompleteCrsMatrix (const Teuchos::RCP<const CrsMatrixType>& sourceMatrix,
01079                                   const Import<typename CrsMatrixType::local_ordinal_type,
01080                                                typename CrsMatrixType::global_ordinal_type,
01081                                                typename CrsMatrixType::node_type>& importer,
01082                                   const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
01083                                                                typename CrsMatrixType::global_ordinal_type,
01084                                                                typename CrsMatrixType::node_type> >& domainMap = Teuchos::null,
01085                                   const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
01086                                                                typename CrsMatrixType::global_ordinal_type,
01087                                                                typename CrsMatrixType::node_type> >& rangeMap = Teuchos::null,
01088                                   const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null)
01089   {
01090     using Teuchos::as;
01091     using Teuchos::RCP;
01092     using Teuchos::rcp;
01093     typedef Map<typename CrsMatrixType::local_ordinal_type,
01094       typename CrsMatrixType::global_ordinal_type,
01095       typename CrsMatrixType::node_type> map_type;
01096 
01097     // FIXME (mfh 11 Apr 2012) The current implementation of this
01098     // method doesn't actually fuse the Import with fillComplete().
01099     // This will change in the future.
01100     RCP<CrsMatrixType> destMat =
01101       rcp (new CrsMatrixType (importer.getTargetMap (),
01102                               as<size_t> (0),
01103                               DynamicProfile,
01104                               params));
01105     destMat->doImport (*sourceMatrix, importer, INSERT);
01106 
01107     // Use the source matrix's domain Map as the default.
01108     RCP<const map_type> theDomainMap =
01109       domainMap.is_null () ? sourceMatrix->getDomainMap () : domainMap;
01110     // Use the source matrix's range Map as the default.
01111     RCP<const map_type> theRangeMap =
01112       rangeMap.is_null () ? sourceMatrix->getRangeMap () : rangeMap;
01113 
01114     destMat->fillComplete (theDomainMap, theRangeMap);
01115     return destMat;
01116   }
01117 
01151   template<class CrsMatrixType>
01152   Teuchos::RCP<CrsMatrixType>
01153   exportAndFillCompleteCrsMatrix (const Teuchos::RCP<const CrsMatrixType>& sourceMatrix,
01154                                   const Export<typename CrsMatrixType::local_ordinal_type,
01155                                                typename CrsMatrixType::global_ordinal_type,
01156                                                typename CrsMatrixType::node_type>& exporter,
01157                                   const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
01158                                                                typename CrsMatrixType::global_ordinal_type,
01159                                                                typename CrsMatrixType::node_type> >& domainMap = Teuchos::null,
01160                                   const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
01161                                                                typename CrsMatrixType::global_ordinal_type,
01162                                                                typename CrsMatrixType::node_type> >& rangeMap = Teuchos::null,
01163                                   const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null)
01164   {
01165     using Teuchos::as;
01166     using Teuchos::RCP;
01167     using Teuchos::rcp;
01168     typedef Map<typename CrsMatrixType::local_ordinal_type,
01169       typename CrsMatrixType::global_ordinal_type,
01170       typename CrsMatrixType::node_type> map_type;
01171 
01172     // FIXME (mfh 11 Apr 2012) The current implementation of this
01173     // method doesn't actually fuse the Export with fillComplete().
01174     // This will change in the future.
01175     RCP<CrsMatrixType> destMat =
01176       rcp (new CrsMatrixType (exporter.getTargetMap (),
01177                               as<size_t> (0),
01178                               DynamicProfile,
01179                               params));
01180     destMat->doExport (*sourceMatrix, exporter, INSERT);
01181 
01182     // Use the source matrix's domain Map as the default.
01183     RCP<const map_type> theDomainMap =
01184       domainMap.is_null () ? sourceMatrix->getDomainMap () : domainMap;
01185     // Use the source matrix's range Map as the default.
01186     RCP<const map_type> theRangeMap =
01187       rangeMap.is_null () ? sourceMatrix->getRangeMap () : rangeMap;
01188 
01189     destMat->fillComplete (theDomainMap, theRangeMap);
01190     return destMat;
01191   }
01192 } // namespace Tpetra
01193 
01204 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines