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 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_CRSMATRIX_DECL_HPP
00030 #define TPETRA_CRSMATRIX_DECL_HPP
00031 
00032 // TODO: row-wise insertion of entries in globalAssemble() may be more efficient
00033 
00034 // TODO: add typeglobs: CrsMatrix<Scalar,typeglob>
00035 // TODO: add template (template) parameter for nonlocal container (this will be part of typeglob)
00036 
00037 #include <Kokkos_DefaultNode.hpp>
00038 #include <Kokkos_DefaultKernels.hpp>
00039 #include <Kokkos_CrsMatrix.hpp>
00040 
00041 #include "Tpetra_ConfigDefs.hpp"
00042 #include "Tpetra_RowMatrix.hpp"
00043 #include "Tpetra_DistObject.hpp"
00044 #include "Tpetra_CrsGraph.hpp"
00045 #include "Tpetra_Vector.hpp"
00046 #include "Tpetra_CrsMatrixMultiplyOp_decl.hpp"
00047 
00048 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00049 namespace Tpetra {
00050   // struct for i,j,v triplets
00051   template <class Ordinal, class Scalar>
00052   struct CrsIJV {
00053     CrsIJV();
00054     CrsIJV(Ordinal row, Ordinal col, const Scalar &val);
00055     Ordinal i,j;
00056     Scalar  v;
00057   };
00058 }
00059 
00060 namespace Teuchos {
00061   // SerializationTraits specialization for CrsIJV, using DirectSerialization
00062   template <typename Ordinal, typename Scalar>
00063   class SerializationTraits<int,Tpetra::CrsIJV<Ordinal,Scalar> >
00064   : public DirectSerializationTraits<int,Tpetra::CrsIJV<Ordinal,Scalar> >
00065   {};
00066 }
00067 
00068 namespace std {
00069   template <class Ordinal, class Scalar>
00070   bool operator<(const Tpetra::CrsIJV<Ordinal,Scalar> &ijv1, const Tpetra::CrsIJV<Ordinal,Scalar> &ijv2);
00071 }
00072 #endif
00073 
00074 namespace Tpetra {
00075 
00077 
00103   template <class Scalar, 
00104             class LocalOrdinal  = int, 
00105             class GlobalOrdinal = LocalOrdinal, 
00106             class Node          = Kokkos::DefaultNode::DefaultNodeType, 
00107             class LocalMatOps   = typename Kokkos::DefaultKernels<Scalar,LocalOrdinal,Node>::SparseOps >
00108   class CrsMatrix : public RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>,
00109                     public DistObject<char, LocalOrdinal,GlobalOrdinal,Node> {
00110     public:
00111       typedef Scalar        scalar_type;
00112       typedef LocalOrdinal  local_ordinal_type;
00113       typedef GlobalOrdinal global_ordinal_type;
00114       typedef Node          node_type;
00115       // backwards compatibility defines both of these
00116       typedef LocalMatOps   mat_vec_type;
00117       typedef LocalMatOps   mat_solve_type;
00118 
00120 
00121 
00123       CrsMatrix(const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rowMap, size_t maxNumEntriesPerRow, ProfileType pftype = DynamicProfile);
00124 
00126       CrsMatrix(const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rowMap, const ArrayRCP<const size_t> &NumEntriesPerRowToAlloc, ProfileType pftype = DynamicProfile);
00127 
00129 
00131       CrsMatrix(const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rowMap, const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &colMap, size_t maxNumEntriesPerRow, ProfileType pftype = DynamicProfile);
00132 
00134 
00136       CrsMatrix(const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rowMap, const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &colMap, const ArrayRCP<const size_t> &NumEntriesPerRowToAlloc, ProfileType pftype = DynamicProfile);
00137 
00139       explicit CrsMatrix(const RCP<const CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps> > &graph);
00140 
00141       // !Destructor.
00142       virtual ~CrsMatrix();
00143 
00145 
00147 
00148 
00150 
00159       void insertGlobalValues(GlobalOrdinal globalRow, const ArrayView<const GlobalOrdinal> &cols, const ArrayView<const Scalar> &vals);
00160 
00162 
00173       void insertLocalValues(LocalOrdinal localRow, const ArrayView<const LocalOrdinal> &cols, const ArrayView<const Scalar> &vals);
00174 
00176 
00181       void replaceGlobalValues(GlobalOrdinal globalRow, 
00182                                const ArrayView<const GlobalOrdinal> &cols,
00183                                const ArrayView<const Scalar>        &vals);
00184 
00186 
00188       void replaceLocalValues(LocalOrdinal localRow, 
00189                               const ArrayView<const LocalOrdinal> &cols,
00190                               const ArrayView<const Scalar>       &vals);
00191 
00193 
00198       void sumIntoGlobalValues(GlobalOrdinal globalRow, 
00199                                const ArrayView<const GlobalOrdinal> &cols,
00200                                const ArrayView<const Scalar>        &vals);
00201 
00202 
00204 
00209       void sumIntoLocalValues(LocalOrdinal globalRow, 
00210                               const ArrayView<const LocalOrdinal>  &cols,
00211                               const ArrayView<const Scalar>        &vals); 
00212 
00214       void setAllToScalar(const Scalar &alpha);
00215 
00217       void scale(const Scalar &alpha);
00218 
00220     
00221 
00223 
00224 
00226       void globalAssemble();
00227 
00236       void resumeFill();
00237 
00249       void fillComplete(const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &domainMap, const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rangeMap, OptimizeOption os = DoOptimizeStorage);
00250 
00264       void fillComplete(OptimizeOption os = DoOptimizeStorage);
00265 
00267 
00269 
00270 
00272       const RCP<const Comm<int> > & getComm() const;
00273 
00275       RCP<Node> getNode() const;
00276 
00278       const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & getRowMap() const;
00279 
00281       const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & getColMap() const;
00282 
00284       RCP<const RowGraph<LocalOrdinal,GlobalOrdinal,Node> > getGraph() const;
00285 
00287       RCP<const CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps> > getCrsGraph() const;
00288 
00309       global_size_t getGlobalNumRows() const;
00310 
00312 
00315       global_size_t getGlobalNumCols() const;
00316 
00318       size_t getNodeNumRows() const;
00319 
00321 
00323       size_t getNodeNumCols() const;
00324 
00326       GlobalOrdinal getIndexBase() const;
00327 
00329       global_size_t getGlobalNumEntries() const;
00330 
00332       size_t getNodeNumEntries() const;
00333 
00335 
00336       size_t getNumEntriesInGlobalRow(GlobalOrdinal globalRow) const;
00337 
00339 
00340       size_t getNumEntriesInLocalRow(LocalOrdinal localRow) const;
00341 
00343 
00345       global_size_t getGlobalNumDiags() const;
00346 
00348 
00350       size_t getNodeNumDiags() const;
00351 
00353 
00355       size_t getGlobalMaxNumRowEntries() const;
00356 
00358 
00360       size_t getNodeMaxNumRowEntries() const;
00361 
00363       bool hasColMap() const; 
00364 
00366 
00368       bool isLowerTriangular() const;
00369 
00371 
00373       bool isUpperTriangular() const;
00374 
00376       bool isLocallyIndexed() const;
00377 
00379       bool isGloballyIndexed() const;
00380 
00382       bool isFillComplete() const;
00383 
00385       bool isFillActive() const;
00386 
00388 
00394       bool isStorageOptimized() const;
00395 
00397       ProfileType getProfileType() const;
00398 
00400       bool isStaticGraph() const;
00401 
00403 
00409       typename ScalarTraits<Scalar>::magnitudeType getFrobeniusNorm() const;
00410 
00412 
00422       void getGlobalRowCopy(GlobalOrdinal GlobalRow,
00423                             const ArrayView<GlobalOrdinal> &Indices,
00424                             const ArrayView<Scalar> &Values,
00425                             size_t &NumEntries
00426                             ) const;
00427 
00429 
00441       void getLocalRowCopy(LocalOrdinal LocalRow, 
00442                            const ArrayView<LocalOrdinal> &Indices, 
00443                            const ArrayView<Scalar> &Values,
00444                            size_t &NumEntries
00445                            ) const;
00446 
00448 
00457       void getGlobalRowView(GlobalOrdinal GlobalRow, ArrayView<const GlobalOrdinal> &indices, ArrayView<const Scalar> &values) const;
00458 
00460 
00469       void getLocalRowView(LocalOrdinal LocalRow, ArrayView<const LocalOrdinal> &indices, ArrayView<const Scalar> &values) const;
00470 
00472 
00474       void getLocalDiagCopy(Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &diag) const;
00475 
00477       void leftScale(const Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& x);
00478 
00480       void rightScale(const Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& x);
00481 
00483 
00485 
00486 
00488 
00498       template <class DomainScalar, class RangeScalar>
00499       void localMultiply(const MultiVector<DomainScalar,LocalOrdinal,GlobalOrdinal,Node> & X, MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,Node> &Y, Teuchos::ETransp trans, RangeScalar alpha, RangeScalar beta) const;
00500 
00502 
00508       template <class DomainScalar, class RangeScalar>
00509       void localSolve(const MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,Node> & Y, MultiVector<DomainScalar,LocalOrdinal,GlobalOrdinal,Node> &X, Teuchos::ETransp trans) const;
00510           
00512 
00514       template <class T>
00515       RCP<CrsMatrix<T,LocalOrdinal,GlobalOrdinal,Node> > convert() const;
00516           
00518 
00520 
00521 
00523 
00526       void apply(const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> & X, MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &Y, 
00527                  Teuchos::ETransp mode = Teuchos::NO_TRANS,
00528                  Scalar alpha = ScalarTraits<Scalar>::one(),
00529                  Scalar beta = ScalarTraits<Scalar>::zero()) const;
00530 
00532       bool hasTransposeApply() const;
00533 
00536       const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & getDomainMap() const;
00537 
00540       const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & getRangeMap() const;
00541 
00543 
00545 
00546 
00548       std::string description() const;
00549 
00551       void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const;
00552 
00554 
00556 
00557 
00558       bool checkSizes(const DistObject<char, LocalOrdinal,GlobalOrdinal,Node>& source);
00559 
00560       void copyAndPermute(const DistObject<char, LocalOrdinal,GlobalOrdinal,Node>& source,
00561                           size_t numSameIDs,
00562                           const ArrayView<const LocalOrdinal> &permuteToLIDs,
00563                           const ArrayView<const LocalOrdinal> &permuteFromLIDs);
00564 
00565       void packAndPrepare(const DistObject<char, LocalOrdinal,GlobalOrdinal,Node>& source,
00566                           const ArrayView<const LocalOrdinal> &exportLIDs,
00567                           Array<char> &exports,
00568                           const ArrayView<size_t> & numPacketsPerLID,
00569                           size_t& constantNumPackets,
00570                           Distributor &distor);
00571 
00572       void unpackAndCombine(const ArrayView<const LocalOrdinal> &importLIDs,
00573                             const ArrayView<const char> &imports,
00574                             const ArrayView<size_t> &numPacketsPerLID,
00575                             size_t constantNumPackets,
00576                             Distributor &distor,
00577                             CombineMode CM);
00578 
00580 
00582 
00583 
00592       TPETRA_DEPRECATED void optimizeStorage();
00593 
00595       TPETRA_DEPRECATED void getGlobalRowView(GlobalOrdinal GlobalRow, ArrayRCP<const GlobalOrdinal> &indices, ArrayRCP<const Scalar> &values) const;
00596 
00598       TPETRA_DEPRECATED void getLocalRowView(LocalOrdinal LocalRow, ArrayRCP<const LocalOrdinal> &indices, ArrayRCP<const Scalar> &values) const;
00599 
00601       template <class DomainScalar, class RangeScalar>
00602       TPETRA_DEPRECATED 
00603       void multiply(const MultiVector<DomainScalar,LocalOrdinal,GlobalOrdinal,Node> & X, MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,Node> &Y, Teuchos::ETransp trans, RangeScalar alpha, RangeScalar beta) const;
00604 
00606       template <class DomainScalar, class RangeScalar>
00607       TPETRA_DEPRECATED 
00608       void solve(const MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,Node> & Y, MultiVector<DomainScalar,LocalOrdinal,GlobalOrdinal,Node> &X, Teuchos::ETransp trans) const;
00609 
00611 
00612     private:
00613       // copy constructor disabled
00614       CrsMatrix(const CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps> &Source);
00615       // operator= disabled
00616       CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps> & operator=(const CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps> &rhs);
00617     protected:
00618       // useful typedefs
00619       typedef OrdinalTraits<LocalOrdinal>                     LOT;
00620       typedef OrdinalTraits<GlobalOrdinal>                    GOT;
00621       typedef ScalarTraits<Scalar>                             ST;
00622       typedef typename ST::magnitudeType                Magnitude;
00623       typedef MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>       MV;
00624       typedef Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>             V;
00625       typedef CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>  Graph;
00626       // Enums
00627       enum GraphAllocationStatus {
00628         GraphAlreadyAllocated,
00629         GraphNotYetAllocated
00630       };
00631       // Allocation
00632       void allocateValues(ELocalGlobal lg, GraphAllocationStatus gas);
00633       // Sorting and merging
00634       void sortEntries();
00635       void mergeRedundantEntries();
00636       // global consts
00637       void clearGlobalConstants();
00638       void computeGlobalConstants();
00639       // matrix data accessors
00640       ArrayView<const Scalar>    getView(RowInfo rowinfo) const;
00641       ArrayView<      Scalar>    getViewNonConst(RowInfo rowinfo);
00642       // local Kokkos objects
00643       void pushToLocalMatrix();
00644       void pullFromLocalMatrix();
00645       void fillLocalMatrix(OptimizeOption os);
00646       void fillLocalSparseOps();
00647       // debugging
00648       void checkInternalState() const;
00649 
00650       // Two graph pointers needed in order to maintain const-correctness:
00651       // staticGraph_ is a graph passed to the constructor. We are not allowed to modify it. it is always a valid pointer.
00652       // myGraph_     is a graph created here. We are allowed to modify it. if myGraph_ != null, then staticGraph_ = myGraph_
00653       RCP<const Graph> staticGraph_;
00654       RCP<      Graph>     myGraph_;
00655 
00656       Kokkos::CrsMatrix<Scalar,LocalOrdinal,Node,LocalMatOps> lclMatrix_;
00657       typename LocalMatOps::template rebind<Scalar>::other    lclMatOps_;
00658 
00659       // matrix values. before allocation, both are null.
00660       // after allocation, one is null.
00661       // 1D == StaticAllocation, 2D == DynamicAllocation
00662       // The allocation always matches that of graph_, as the graph does the allocation for the matrix.
00663       ArrayRCP<Scalar>                       values1D_;
00664       ArrayRCP<ArrayRCP<Scalar> >            values2D_;
00665       // TODO: these could be allocated at resumeFill() and de-allocated at fillComplete() to make for very fast getView()/getViewNonConst()
00666       // ArrayRCP< typedef ArrayRCP<const Scalar>::iterator > rowPtrs_;
00667       // ArrayRCP< typedef ArrayRCP<      Scalar>::iterator > rowPtrsNC_;
00668 
00669       bool fillComplete_;
00670 
00671       // non-local data
00672       std::map<GlobalOrdinal, Array<std::pair<GlobalOrdinal,Scalar> > > nonlocals_;
00673 
00674       // a wrapper around multiply, for use in apply; it contains a non-owning RCP to *this, therefore, it is not allowed 
00675       // to persist past the destruction of *this. therefore, WE MAY NOT SHARE THIS POINTER.
00676       RCP< const CrsMatrixMultiplyOp<Scalar,Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps> > sameScalarMultiplyOp_;
00677 
00678       // cached frobenius norm: -ST::one() means invalid
00679       mutable Magnitude frobNorm_;
00680 
00681   }; // class CrsMatrix
00682 
00689   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00690   RCP<CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
00691   createCrsMatrix(const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &map, size_t maxNumEntriesPerRow = 0)
00692   {
00693     RCP<CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > ret;
00694     ret = rcp(new CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>(map, maxNumEntriesPerRow, DynamicProfile) );
00695     return ret;
00696   }
00697 
00698 } // namespace Tpetra
00699 
00710 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines