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 #include <Kokkos_CrsMatrix.hpp>
00053 
00054 #include "Tpetra_ConfigDefs.hpp"
00055 #include "Tpetra_RowMatrix.hpp"
00056 #include "Tpetra_DistObject.hpp"
00057 #include "Tpetra_CrsGraph.hpp"
00058 #include "Tpetra_Vector.hpp"
00059 #include "Tpetra_CrsMatrixMultiplyOp_decl.hpp"
00060 
00061 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00062 namespace Tpetra {
00063   // struct for i,j,v triplets
00064   template <class Ordinal, class Scalar>
00065   struct CrsIJV {
00066     CrsIJV();
00067     CrsIJV(Ordinal row, Ordinal col, const Scalar &val);
00068     Ordinal i,j;
00069     Scalar  v;
00070   };
00071 }
00072 
00073 namespace Teuchos {
00074   // SerializationTraits specialization for CrsIJV, using DirectSerialization
00075   template <typename Ordinal, typename Scalar>
00076   class SerializationTraits<int,Tpetra::CrsIJV<Ordinal,Scalar> >
00077   : public DirectSerializationTraits<int,Tpetra::CrsIJV<Ordinal,Scalar> >
00078   {};
00079 }
00080 
00081 namespace std {
00082   template <class Ordinal, class Scalar>
00083   bool operator<(const Tpetra::CrsIJV<Ordinal,Scalar> &ijv1, const Tpetra::CrsIJV<Ordinal,Scalar> &ijv2);
00084 }
00085 #endif
00086 
00087 namespace Tpetra {
00088 
00090 
00116   template <class Scalar, 
00117             class LocalOrdinal  = int, 
00118             class GlobalOrdinal = LocalOrdinal, 
00119             class Node          = Kokkos::DefaultNode::DefaultNodeType, 
00120             class LocalMatOps   = typename Kokkos::DefaultKernels<Scalar,LocalOrdinal,Node>::SparseOps >
00121   class CrsMatrix : public RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>,
00122                     public DistObject<char, LocalOrdinal,GlobalOrdinal,Node> {
00123     public:
00124       typedef Scalar        scalar_type;
00125       typedef LocalOrdinal  local_ordinal_type;
00126       typedef GlobalOrdinal global_ordinal_type;
00127       typedef Node          node_type;
00128       // backwards compatibility defines both of these
00129       typedef LocalMatOps   mat_vec_type;
00130       typedef LocalMatOps   mat_solve_type;
00131 
00133 
00134 
00136       CrsMatrix(const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rowMap, size_t maxNumEntriesPerRow, ProfileType pftype = DynamicProfile);
00137 
00139       CrsMatrix(const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rowMap, const ArrayRCP<const size_t> &NumEntriesPerRowToAlloc, ProfileType pftype = DynamicProfile);
00140 
00142 
00144       CrsMatrix(const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rowMap, const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &colMap, size_t maxNumEntriesPerRow, ProfileType pftype = DynamicProfile);
00145 
00147 
00149       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);
00150 
00152       explicit CrsMatrix(const RCP<const CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps> > &graph);
00153 
00154       // !Destructor.
00155       virtual ~CrsMatrix();
00156 
00158 
00160 
00161 
00163 
00172       void insertGlobalValues(GlobalOrdinal globalRow, const ArrayView<const GlobalOrdinal> &cols, const ArrayView<const Scalar> &vals);
00173 
00175 
00186       void insertLocalValues(LocalOrdinal localRow, const ArrayView<const LocalOrdinal> &cols, const ArrayView<const Scalar> &vals);
00187 
00189 
00194       void replaceGlobalValues(GlobalOrdinal globalRow, 
00195                                const ArrayView<const GlobalOrdinal> &cols,
00196                                const ArrayView<const Scalar>        &vals);
00197 
00199 
00201       void replaceLocalValues(LocalOrdinal localRow, 
00202                               const ArrayView<const LocalOrdinal> &cols,
00203                               const ArrayView<const Scalar>       &vals);
00204 
00206 
00211       void sumIntoGlobalValues(GlobalOrdinal globalRow, 
00212                                const ArrayView<const GlobalOrdinal> &cols,
00213                                const ArrayView<const Scalar>        &vals);
00214 
00215 
00217 
00222       void sumIntoLocalValues(LocalOrdinal globalRow, 
00223                               const ArrayView<const LocalOrdinal>  &cols,
00224                               const ArrayView<const Scalar>        &vals); 
00225 
00227       void setAllToScalar(const Scalar &alpha);
00228 
00230       void scale(const Scalar &alpha);
00231 
00233     
00234 
00236 
00237 
00239       void globalAssemble();
00240 
00249       void resumeFill();
00250 
00262       void fillComplete(const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &domainMap, const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rangeMap, OptimizeOption os = DoOptimizeStorage);
00263 
00277       void fillComplete(OptimizeOption os = DoOptimizeStorage);
00278 
00280 
00282 
00283 
00285       const RCP<const Comm<int> > & getComm() const;
00286 
00288       RCP<Node> getNode() const;
00289 
00291       const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & getRowMap() const;
00292 
00294       const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & getColMap() const;
00295 
00297       RCP<const RowGraph<LocalOrdinal,GlobalOrdinal,Node> > getGraph() const;
00298 
00300       RCP<const CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps> > getCrsGraph() const;
00301 
00322       global_size_t getGlobalNumRows() const;
00323 
00325 
00328       global_size_t getGlobalNumCols() const;
00329 
00331       size_t getNodeNumRows() const;
00332 
00334 
00336       size_t getNodeNumCols() const;
00337 
00339       GlobalOrdinal getIndexBase() const;
00340 
00342       global_size_t getGlobalNumEntries() const;
00343 
00345       size_t getNodeNumEntries() const;
00346 
00348 
00349       size_t getNumEntriesInGlobalRow(GlobalOrdinal globalRow) const;
00350 
00352 
00353       size_t getNumEntriesInLocalRow(LocalOrdinal localRow) const;
00354 
00356 
00358       global_size_t getGlobalNumDiags() const;
00359 
00361 
00363       size_t getNodeNumDiags() const;
00364 
00366 
00368       size_t getGlobalMaxNumRowEntries() const;
00369 
00371 
00373       size_t getNodeMaxNumRowEntries() const;
00374 
00376       bool hasColMap() const; 
00377 
00379 
00381       bool isLowerTriangular() const;
00382 
00384 
00386       bool isUpperTriangular() const;
00387 
00389       bool isLocallyIndexed() const;
00390 
00392       bool isGloballyIndexed() const;
00393 
00395       bool isFillComplete() const;
00396 
00398       bool isFillActive() const;
00399 
00401 
00407       bool isStorageOptimized() const;
00408 
00410       ProfileType getProfileType() const;
00411 
00413       bool isStaticGraph() const;
00414 
00416 
00422       typename ScalarTraits<Scalar>::magnitudeType getFrobeniusNorm() const;
00423 
00425 
00435       void getGlobalRowCopy(GlobalOrdinal GlobalRow,
00436                             const ArrayView<GlobalOrdinal> &Indices,
00437                             const ArrayView<Scalar> &Values,
00438                             size_t &NumEntries
00439                             ) const;
00440 
00442 
00454       void getLocalRowCopy(LocalOrdinal LocalRow, 
00455                            const ArrayView<LocalOrdinal> &Indices, 
00456                            const ArrayView<Scalar> &Values,
00457                            size_t &NumEntries
00458                            ) const;
00459 
00461 
00470       void getGlobalRowView(GlobalOrdinal GlobalRow, ArrayView<const GlobalOrdinal> &indices, ArrayView<const Scalar> &values) const;
00471 
00473 
00482       void getLocalRowView(LocalOrdinal LocalRow, ArrayView<const LocalOrdinal> &indices, ArrayView<const Scalar> &values) const;
00483 
00485 
00487       void getLocalDiagCopy(Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &diag) const;
00488 
00490       void leftScale(const Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& x);
00491 
00493       void rightScale(const Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& x);
00494 
00496 
00498 
00499 
00501 
00511       template <class DomainScalar, class RangeScalar>
00512       void localMultiply(const MultiVector<DomainScalar,LocalOrdinal,GlobalOrdinal,Node> & X, MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,Node> &Y, Teuchos::ETransp trans, RangeScalar alpha, RangeScalar beta) const;
00513 
00515 
00521       template <class DomainScalar, class RangeScalar>
00522       void localSolve(const MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,Node> & Y, MultiVector<DomainScalar,LocalOrdinal,GlobalOrdinal,Node> &X, Teuchos::ETransp trans) const;
00523           
00525 
00527       template <class T>
00528       RCP<CrsMatrix<T,LocalOrdinal,GlobalOrdinal,Node> > convert() const;
00529           
00531 
00533 
00534 
00536 
00539       void apply(const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> & X, MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &Y, 
00540                  Teuchos::ETransp mode = Teuchos::NO_TRANS,
00541                  Scalar alpha = ScalarTraits<Scalar>::one(),
00542                  Scalar beta = ScalarTraits<Scalar>::zero()) const;
00543 
00545       bool hasTransposeApply() const;
00546 
00549       const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & getDomainMap() const;
00550 
00553       const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & getRangeMap() const;
00554 
00556 
00558 
00559 
00561       std::string description() const;
00562 
00564       void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const;
00565 
00567 
00569 
00570 
00571       bool checkSizes(const DistObject<char, LocalOrdinal,GlobalOrdinal,Node>& source);
00572 
00573       void copyAndPermute(const DistObject<char, LocalOrdinal,GlobalOrdinal,Node>& source,
00574                           size_t numSameIDs,
00575                           const ArrayView<const LocalOrdinal> &permuteToLIDs,
00576                           const ArrayView<const LocalOrdinal> &permuteFromLIDs);
00577 
00578       void packAndPrepare(const DistObject<char, LocalOrdinal,GlobalOrdinal,Node>& source,
00579                           const ArrayView<const LocalOrdinal> &exportLIDs,
00580                           Array<char> &exports,
00581                           const ArrayView<size_t> & numPacketsPerLID,
00582                           size_t& constantNumPackets,
00583                           Distributor &distor);
00584 
00585       void unpackAndCombine(const ArrayView<const LocalOrdinal> &importLIDs,
00586                             const ArrayView<const char> &imports,
00587                             const ArrayView<size_t> &numPacketsPerLID,
00588                             size_t constantNumPackets,
00589                             Distributor &distor,
00590                             CombineMode CM);
00591 
00593 
00595 
00596 
00605       TPETRA_DEPRECATED void optimizeStorage();
00606 
00608       TPETRA_DEPRECATED void getGlobalRowView(GlobalOrdinal GlobalRow, ArrayRCP<const GlobalOrdinal> &indices, ArrayRCP<const Scalar> &values) const;
00609 
00611       TPETRA_DEPRECATED void getLocalRowView(LocalOrdinal LocalRow, ArrayRCP<const LocalOrdinal> &indices, ArrayRCP<const Scalar> &values) const;
00612 
00614       template <class DomainScalar, class RangeScalar>
00615       TPETRA_DEPRECATED 
00616       void multiply(const MultiVector<DomainScalar,LocalOrdinal,GlobalOrdinal,Node> & X, MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,Node> &Y, Teuchos::ETransp trans, RangeScalar alpha, RangeScalar beta) const;
00617 
00619       template <class DomainScalar, class RangeScalar>
00620       TPETRA_DEPRECATED 
00621       void solve(const MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,Node> & Y, MultiVector<DomainScalar,LocalOrdinal,GlobalOrdinal,Node> &X, Teuchos::ETransp trans) const;
00622 
00624 
00625     private:
00626       // copy constructor disabled
00627       CrsMatrix(const CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps> &Source);
00628       // operator= disabled
00629       CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps> & operator=(const CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps> &rhs);
00630     protected:
00631       // useful typedefs
00632       typedef OrdinalTraits<LocalOrdinal>                     LOT;
00633       typedef OrdinalTraits<GlobalOrdinal>                    GOT;
00634       typedef ScalarTraits<Scalar>                             ST;
00635       typedef typename ST::magnitudeType                Magnitude;
00636       typedef MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>       MV;
00637       typedef Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>             V;
00638       typedef CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>  Graph;
00639       // Enums
00640       enum GraphAllocationStatus {
00641         GraphAlreadyAllocated,
00642         GraphNotYetAllocated
00643       };
00644       // Allocation
00645       void allocateValues(ELocalGlobal lg, GraphAllocationStatus gas);
00646       // Sorting and merging
00647       void sortEntries();
00648       void mergeRedundantEntries();
00649       // global consts
00650       void clearGlobalConstants();
00651       void computeGlobalConstants();
00652       // matrix data accessors
00653       ArrayView<const Scalar>    getView(RowInfo rowinfo) const;
00654       ArrayView<      Scalar>    getViewNonConst(RowInfo rowinfo);
00655       // local Kokkos objects
00656       void pushToLocalMatrix();
00657       void pullFromLocalMatrix();
00658       void fillLocalMatrix(OptimizeOption os);
00659       void fillLocalSparseOps();
00660       // debugging
00661       void checkInternalState() const;
00662 
00663       // Two graph pointers needed in order to maintain const-correctness:
00664       // staticGraph_ is a graph passed to the constructor. We are not allowed to modify it. it is always a valid pointer.
00665       // myGraph_     is a graph created here. We are allowed to modify it. if myGraph_ != null, then staticGraph_ = myGraph_
00666       RCP<const Graph> staticGraph_;
00667       RCP<      Graph>     myGraph_;
00668 
00669       Kokkos::CrsMatrix<Scalar,LocalOrdinal,Node,LocalMatOps> lclMatrix_;
00670       typename LocalMatOps::template rebind<Scalar>::other    lclMatOps_;
00671 
00672       // matrix values. before allocation, both are null.
00673       // after allocation, one is null.
00674       // 1D == StaticAllocation, 2D == DynamicAllocation
00675       // The allocation always matches that of graph_, as the graph does the allocation for the matrix.
00676       ArrayRCP<Scalar>                       values1D_;
00677       ArrayRCP<ArrayRCP<Scalar> >            values2D_;
00678       // TODO: these could be allocated at resumeFill() and de-allocated at fillComplete() to make for very fast getView()/getViewNonConst()
00679       // ArrayRCP< typedef ArrayRCP<const Scalar>::iterator > rowPtrs_;
00680       // ArrayRCP< typedef ArrayRCP<      Scalar>::iterator > rowPtrsNC_;
00681 
00682       bool fillComplete_;
00683 
00684       // non-local data
00685       std::map<GlobalOrdinal, Array<std::pair<GlobalOrdinal,Scalar> > > nonlocals_;
00686 
00687       // a wrapper around multiply, for use in apply; it contains a non-owning RCP to *this, therefore, it is not allowed 
00688       // to persist past the destruction of *this. therefore, WE MAY NOT SHARE THIS POINTER.
00689       RCP< const CrsMatrixMultiplyOp<Scalar,Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps> > sameScalarMultiplyOp_;
00690 
00691       // cached frobenius norm: -ST::one() means invalid
00692       mutable Magnitude frobNorm_;
00693 
00694   }; // class CrsMatrix
00695 
00702   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
00703   RCP<CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
00704   createCrsMatrix(const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &map, size_t maxNumEntriesPerRow = 0)
00705   {
00706     RCP<CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > ret;
00707     ret = rcp(new CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>(map, maxNumEntriesPerRow, DynamicProfile) );
00708     return ret;
00709   }
00710 
00711 } // namespace Tpetra
00712 
00723 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines