Tpetra Matrix/Vector Services Version of the Day
Tpetra_CrsMatrix_decl.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 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 
00047 
00048 #include <Kokkos_DefaultNode.hpp>
00049 #include <Kokkos_DefaultKernels.hpp>
00050 
00051 #include "Tpetra_ConfigDefs.hpp"
00052 #include "Tpetra_RowMatrix_decl.hpp"
00053 #include "Tpetra_Exceptions.hpp"
00054 #include "Tpetra_DistObject.hpp"
00055 #include "Tpetra_CrsGraph.hpp"
00056 #include "Tpetra_Vector.hpp"
00057 
00058 
00059 namespace Tpetra {
00060 
00062 
00199   template <class Scalar,
00200             class LocalOrdinal  = int,
00201             class GlobalOrdinal = LocalOrdinal,
00202             class Node          = KokkosClassic::DefaultNode::DefaultNodeType,
00203             class LocalMatOps   = typename KokkosClassic::DefaultKernels<Scalar,LocalOrdinal,Node>::SparseOps >
00204   class CrsMatrix : public RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>,
00205                     public DistObject<char, LocalOrdinal,GlobalOrdinal,Node> {
00206   public:
00208 
00209 
00211     typedef Scalar scalar_type;
00213     typedef LocalOrdinal local_ordinal_type;
00215     typedef GlobalOrdinal global_ordinal_type;
00217     typedef Node node_type;
00218 
00222     typedef LocalMatOps mat_vec_type;
00226     typedef LocalMatOps mat_solve_type;
00227 
00229     typedef Map<LocalOrdinal, GlobalOrdinal, node_type> map_type;
00230 
00232     typedef CrsGraph<LocalOrdinal, GlobalOrdinal, node_type, LocalMatOps> crs_graph_type;
00233 
00235 
00236 
00237 
00255     CrsMatrix (const Teuchos::RCP<const map_type>& rowMap,
00256                size_t maxNumEntriesPerRow,
00257                ProfileType pftype = DynamicProfile,
00258                const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null);
00259 
00277     CrsMatrix (const Teuchos::RCP<const map_type>& rowMap,
00278                const ArrayRCP<const size_t>& NumEntriesPerRowToAlloc,
00279                ProfileType pftype = DynamicProfile,
00280                const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null);
00281 
00304     CrsMatrix (const Teuchos::RCP<const map_type>& rowMap,
00305                const Teuchos::RCP<const map_type>& colMap,
00306                size_t maxNumEntriesPerRow,
00307                ProfileType pftype = DynamicProfile,
00308                const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null);
00309 
00332     CrsMatrix (const Teuchos::RCP<const map_type>& rowMap,
00333                const Teuchos::RCP<const map_type>& colMap,
00334                const ArrayRCP<const size_t>& NumEntriesPerRowToAlloc,
00335                ProfileType pftype = DynamicProfile,
00336                const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null);
00337 
00362     explicit CrsMatrix (const Teuchos::RCP<const crs_graph_type>& graph,
00363                         const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null);
00364 
00387     CrsMatrix (const RCP<const map_type>& rowMap,
00388                const RCP<const map_type>& colMap,
00389                const ArrayRCP<size_t>& rowPointers,
00390                const ArrayRCP<LocalOrdinal>& columnIndices,
00391                const ArrayRCP<Scalar>& values,
00392                const RCP<ParameterList>& params = null);
00393 
00394     // This friend declaration makes the clone() method work.
00395     template <class S2, class LO2, class GO2, class N2, class LMO2>
00396     friend class CrsMatrix;
00397 
00422     template <class Node2>
00423     RCP<CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node2, typename KokkosClassic::DefaultKernels<void,LocalOrdinal,Node2>::SparseOps> >
00424     clone (const RCP<Node2>& node2, const RCP<ParameterList>& params = null) const
00425     {
00426       const char tfecfFuncName[] = "clone";
00427 
00428       // Get parameter values.  Set them initially to their default values.
00429       bool fillCompleteClone = true;
00430       bool useLocalIndices = this->hasColMap ();
00431       ProfileType pftype = StaticProfile;
00432       if (! params.is_null ()) {
00433         fillCompleteClone = params->get ("fillComplete clone", fillCompleteClone);
00434         useLocalIndices = params->get ("Locally indexed clone", useLocalIndices);
00435 
00436         bool staticProfileClone = true;
00437         staticProfileClone = params->get ("Static profile clone", staticProfileClone);
00438         pftype = staticProfileClone ? StaticProfile : DynamicProfile;
00439       }
00440 
00441       TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00442         ! this->hasColMap () && useLocalIndices, std::runtime_error,
00443         ": You requested that the returned clone have local indices, but the "
00444         "the source matrix does not have a column Map yet.");
00445 
00446       typedef typename KokkosClassic::DefaultKernels<void,LocalOrdinal,Node2>::SparseOps LocalMatOps2;
00447       typedef CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node2, LocalMatOps2> CrsMatrix2;
00448       typedef Map<LocalOrdinal, GlobalOrdinal, Node2> Map2;
00449       RCP<const Map2> clonedRowMap = this->getRowMap ()->template clone<Node2> (node2);
00450 
00451       RCP<CrsMatrix2> clonedMatrix;
00452       ArrayRCP<const size_t> numEntries;
00453       size_t numEntriesForAll = 0;
00454       if (! staticGraph_->indicesAreAllocated ()) {
00455         if (! staticGraph_->numAllocPerRow_.is_null ()) {
00456           numEntries = staticGraph_->numAllocPerRow_;
00457         }
00458         else {
00459           numEntriesForAll = staticGraph_->numAllocForAllRows_;
00460         }
00461       }
00462       else if (! staticGraph_->numRowEntries_.is_null ()) {
00463         numEntries = staticGraph_->numRowEntries_;
00464       }
00465       else if (staticGraph_->nodeNumAllocated_ == 0) {
00466         numEntriesForAll = 0;
00467       }
00468       else {
00469         // We're left with the case that we have optimized storage.
00470         // In this case, we have to construct a list of row sizes.
00471         TEUCHOS_TEST_FOR_EXCEPTION(
00472           getProfileType() != StaticProfile, std::logic_error,
00473           "Internal logic error. Please report this to Tpetra team." )
00474 
00475         const size_t numRows = this->getNodeNumRows ();
00476         numEntriesForAll = 0;
00477         ArrayRCP<size_t> numEnt;
00478         if (numRows != 0) {
00479           numEnt = arcp<size_t> (numRows);
00480         }
00481         for (size_t i = 0; i < numRows; ++i) {
00482           numEnt[i] = staticGraph_->rowPtrs_[i+1] - staticGraph_->rowPtrs_[i];
00483         }
00484         numEntries = numEnt;
00485       }
00486 
00487       RCP<ParameterList> matrixparams =
00488         params.is_null () ? null : sublist (params,"CrsMatrix");
00489       if (useLocalIndices) {
00490         RCP<const Map2> clonedColMap =
00491           this->getColMap ()->template clone<Node2> (node2);
00492         if (numEntries.is_null ()) {
00493           clonedMatrix = rcp (new CrsMatrix2 (clonedRowMap, clonedColMap,
00494                                               numEntriesForAll, pftype,
00495                                               matrixparams));
00496         }
00497         else {
00498           clonedMatrix = rcp (new CrsMatrix2 (clonedRowMap, clonedColMap,
00499                                               numEntries, pftype,
00500                                               matrixparams));
00501         }
00502       }
00503       else {
00504         if (numEntries.is_null ()) {
00505           clonedMatrix = rcp (new CrsMatrix2 (clonedRowMap, numEntriesForAll,
00506                                               pftype, matrixparams));
00507         }
00508         else {
00509           clonedMatrix = rcp (new CrsMatrix2 (clonedRowMap, numEntries, pftype,
00510                                               matrixparams));
00511         }
00512       }
00513       // done with these
00514       numEntries = null;
00515       numEntriesForAll = 0;
00516 
00517       if (useLocalIndices) {
00518         clonedMatrix->allocateValues (LocalIndices,
00519                                       CrsMatrix2::GraphNotYetAllocated);
00520         if (this->isLocallyIndexed ()) {
00521           ArrayView<const LocalOrdinal> linds;
00522           ArrayView<const Scalar>       vals;
00523           for (LocalOrdinal lrow = clonedRowMap->getMinLocalIndex ();
00524                lrow <= clonedRowMap->getMaxLocalIndex ();
00525                ++lrow) {
00526             this->getLocalRowView (lrow, linds, vals);
00527             if (linds.size ()) {
00528               clonedMatrix->insertLocalValues (lrow, linds, vals);
00529             }
00530           }
00531         }
00532         else { // this->isGloballyIndexed()
00533           Array<LocalOrdinal> linds;
00534           Array<Scalar> vals;
00535           for (LocalOrdinal lrow = clonedRowMap->getMinLocalIndex ();
00536                lrow <= clonedRowMap->getMaxLocalIndex ();
00537                ++lrow) {
00538             size_t theNumEntries = this->getNumEntriesInLocalRow (lrow);
00539             if (theNumEntries > Teuchos::as<size_t> (linds.size ())) {
00540               linds.resize (theNumEntries);
00541             }
00542             if (theNumEntries > Teuchos::as<size_t> (vals.size ())) {
00543               vals.resize (theNumEntries);
00544             }
00545             this->getLocalRowCopy (clonedRowMap->getGlobalElement (lrow),
00546                                    linds (), vals (), theNumEntries);
00547             if (theNumEntries != 0) {
00548               clonedMatrix->insertLocalValues (lrow, linds (0, theNumEntries),
00549                                                vals (0, theNumEntries));
00550             }
00551           }
00552         }
00553       }
00554       else { // useGlobalIndices
00555         clonedMatrix->allocateValues (GlobalIndices,
00556                                       CrsMatrix2::GraphNotYetAllocated);
00557         if (this->isGloballyIndexed ()) {
00558           ArrayView<const GlobalOrdinal> ginds;
00559           ArrayView<const Scalar>         vals;
00560           for (GlobalOrdinal grow = clonedRowMap->getMinGlobalIndex ();
00561                grow <= clonedRowMap->getMaxGlobalIndex ();
00562                ++grow) {
00563             this->getGlobalRowView (grow, ginds, vals);
00564             if (ginds.size () > 0) {
00565               clonedMatrix->insertGlobalValues (grow, ginds, vals);
00566             }
00567           }
00568         }
00569         else { // this->isLocallyIndexed()
00570           Array<GlobalOrdinal> ginds;
00571           Array<Scalar> vals;
00572           for (GlobalOrdinal grow = clonedRowMap->getMinGlobalIndex ();
00573                grow <= clonedRowMap->getMaxGlobalIndex ();
00574                ++grow) {
00575             size_t theNumEntries = this->getNumEntriesInGlobalRow (grow);
00576             if (theNumEntries > Teuchos::as<size_t> (ginds.size ())) {
00577               ginds.resize (theNumEntries);
00578             }
00579             if (theNumEntries > Teuchos::as<size_t> (vals.size ())) {
00580               vals.resize (theNumEntries);
00581             }
00582             this->getGlobalRowCopy (grow, ginds (), vals (), theNumEntries);
00583             if (theNumEntries != 0) {
00584               clonedMatrix->insertGlobalValues (grow, ginds (0, theNumEntries),
00585                                                 vals (0, theNumEntries));
00586             }
00587           }
00588         }
00589       }
00590 
00591       if (fillCompleteClone) {
00592         RCP<ParameterList> fillparams =
00593           params.is_null () ? Teuchos::null : sublist (params, "fillComplete");
00594         try {
00595           RCP<const Map2> clonedRangeMap;
00596           RCP<const Map2> clonedDomainMap;
00597           if (! this->getRangeMap ().is_null () &&
00598               this->getRangeMap () != clonedRowMap) {
00599             clonedRangeMap  = this->getRangeMap ()->template clone<Node2> (node2);
00600           }
00601           else {
00602             clonedRangeMap = clonedRowMap;
00603           }
00604           if (! this->getDomainMap ().is_null () &&
00605               this->getDomainMap () != clonedRowMap) {
00606             clonedDomainMap = this->getDomainMap ()->template clone<Node2> (node2);
00607           }
00608           else {
00609             clonedDomainMap = clonedRowMap;
00610           }
00611           clonedMatrix->fillComplete (clonedDomainMap, clonedRangeMap,
00612                                       fillparams);
00613         }
00614         catch (std::exception &e) {
00615           const bool caughtExceptionOnClone = true;
00616           TEUCHOS_TEST_FOR_EXCEPTION(
00617             caughtExceptionOnClone, std::runtime_error,
00618             Teuchos::typeName (*this) << std::endl << "clone: " << std::endl <<
00619             "Caught the following exception while calling fillComplete() on a "
00620             "clone of type" << std::endl << Teuchos::typeName (*clonedMatrix)
00621             << ": " << std::endl << e.what () << std::endl);
00622         }
00623       }
00624       return clonedMatrix;
00625     }
00626 
00628     virtual ~CrsMatrix ();
00629 
00631 
00632 
00633 
00656     //
00702     void
00703     insertGlobalValues (const GlobalOrdinal globalRow,
00704                         const ArrayView<const GlobalOrdinal>& cols,
00705                         const ArrayView<const Scalar>& vals);
00706 
00746     void
00747     insertLocalValues (const LocalOrdinal localRow,
00748                        const ArrayView<const LocalOrdinal> &cols,
00749                        const ArrayView<const Scalar> &vals);
00750 
00769     void
00770     replaceGlobalValues (GlobalOrdinal globalRow,
00771                          const ArrayView<const GlobalOrdinal>& cols,
00772                          const ArrayView<const Scalar>& vals);
00773 
00787     void
00788     replaceLocalValues (LocalOrdinal localRow,
00789                         const ArrayView<const LocalOrdinal> &cols,
00790                         const ArrayView<const Scalar>       &vals);
00791 
00818     void
00819     sumIntoGlobalValues (const GlobalOrdinal globalRow,
00820                          const ArrayView<const GlobalOrdinal> &cols,
00821                          const ArrayView<const Scalar>        &vals);
00822 
00831     void
00832     sumIntoLocalValues (const LocalOrdinal localRow,
00833                         const ArrayView<const LocalOrdinal>  &cols,
00834                         const ArrayView<const Scalar>        &vals);
00835 
00837     void setAllToScalar (const Scalar &alpha);
00838 
00840     void scale (const Scalar &alpha);
00841 
00843 
00850     void
00851     setAllValues (const ArrayRCP<size_t>& rowPointers,
00852                   const ArrayRCP<LocalOrdinal>& columnIndices,
00853                   const ArrayRCP<Scalar>& values);
00854 
00856 
00863     void
00864     getAllValues (ArrayRCP<const size_t>& rowPointers,
00865                   ArrayRCP<const LocalOrdinal>& columnIndices,
00866                   ArrayRCP<const Scalar>& values) const;
00867 
00869 
00870 
00871 
00900     void globalAssemble();
00901 
00915     void resumeFill (const RCP<ParameterList>& params = null);
00916 
00937     void
00938     fillComplete (const RCP<const map_type>& domainMap,
00939                   const RCP<const map_type>& rangeMap,
00940                   const RCP<ParameterList>& params = null);
00941 
00954     void fillComplete (const RCP<ParameterList>& params = null);
00955 
00967     void
00968     expertStaticFillComplete (const RCP<const map_type>& domainMap,
00969                               const RCP<const map_type>& rangeMap,
00970                               const RCP<const Import<LocalOrdinal,GlobalOrdinal,node_type> >& importer = Teuchos::null,
00971                               const RCP<const Export<LocalOrdinal,GlobalOrdinal,node_type> >& exporter = Teuchos::null,
00972                               const RCP<ParameterList>& params = Teuchos::null);
00973 
00983     void
00984     replaceColMap (const Teuchos::RCP<const map_type>& newColMap);
00985 
00998     void
00999     replaceDomainMapAndImporter (const Teuchos::RCP<const map_type>& newDomainMap,
01000                                  Teuchos::RCP<const Tpetra::Import<LocalOrdinal,GlobalOrdinal,node_type> >& newImporter);
01001 
01015     virtual void
01016     removeEmptyProcessesInPlace (const Teuchos::RCP<const map_type>& newMap);
01017 
01019 
01020 
01021 
01023     RCP<const Comm<int> > getComm() const;
01024 
01026     RCP<node_type> getNode () const;
01027 
01029     RCP<const map_type> getRowMap () const;
01030 
01032     RCP<const map_type> getColMap () const;
01033 
01035     RCP<const RowGraph<LocalOrdinal,GlobalOrdinal,node_type> > getGraph () const;
01036 
01038     RCP<const crs_graph_type> getCrsGraph () const;
01039 
01059     global_size_t getGlobalNumRows() const;
01060 
01066     global_size_t getGlobalNumCols() const;
01067 
01074     size_t getNodeNumRows() const;
01075 
01079     size_t getNodeNumCols() const;
01080 
01082     GlobalOrdinal getIndexBase() const;
01083 
01085     global_size_t getGlobalNumEntries() const;
01086 
01088     size_t getNodeNumEntries() const;
01089 
01091 
01092     size_t getNumEntriesInGlobalRow(GlobalOrdinal globalRow) const;
01093 
01095 
01096     size_t getNumEntriesInLocalRow(LocalOrdinal localRow) const;
01097 
01099 
01101     global_size_t getGlobalNumDiags() const;
01102 
01104 
01106     size_t getNodeNumDiags() const;
01107 
01109 
01111     size_t getGlobalMaxNumRowEntries() const;
01112 
01114 
01116     size_t getNodeMaxNumRowEntries() const;
01117 
01119     bool hasColMap() const;
01120 
01122 
01124     bool isLowerTriangular() const;
01125 
01127 
01129     bool isUpperTriangular() const;
01130 
01150     bool isLocallyIndexed() const;
01151 
01171     bool isGloballyIndexed() const;
01172 
01195     bool isFillComplete() const;
01196 
01219     bool isFillActive() const;
01220 
01222 
01228     bool isStorageOptimized() const;
01229 
01231     ProfileType getProfileType() const;
01232 
01234     bool isStaticGraph() const;
01235 
01237 
01243     typename ScalarTraits<Scalar>::magnitudeType getFrobeniusNorm() const;
01244 
01245 
01247     virtual bool supportsRowViews() const;
01248 
01250 
01260     void
01261     getGlobalRowCopy (GlobalOrdinal GlobalRow,
01262                       const ArrayView<GlobalOrdinal> &Indices,
01263                       const ArrayView<Scalar> &Values,
01264                       size_t &NumEntries) const;
01265 
01267 
01279     void
01280     getLocalRowCopy (LocalOrdinal LocalRow,
01281                      const ArrayView<LocalOrdinal> &Indices,
01282                      const ArrayView<Scalar> &Values,
01283                      size_t &NumEntries) const;
01284 
01286 
01295     void
01296     getGlobalRowView (GlobalOrdinal GlobalRow,
01297                       ArrayView<const GlobalOrdinal> &indices,
01298                       ArrayView<const Scalar> &values) const;
01299 
01301 
01310     void
01311     getLocalRowView (LocalOrdinal LocalRow,
01312                      ArrayView<const LocalOrdinal> &indices,
01313                      ArrayView<const Scalar> &values) const;
01314 
01316 
01318     void getLocalDiagCopy (Vector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>& diag) const;
01319 
01357     void getLocalDiagOffsets (Teuchos::ArrayRCP<size_t>& offsets) const;
01358 
01377     void
01378     getLocalDiagCopy (Vector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>& diag,
01379                       const Teuchos::ArrayView<const size_t>& offsets) const;
01380 
01382     void leftScale(const Vector<Scalar, LocalOrdinal, GlobalOrdinal, node_type>& x);
01383 
01385     void rightScale(const Vector<Scalar, LocalOrdinal, GlobalOrdinal, node_type>& x);
01386 
01388 
01389 
01390 
01439     template <class DomainScalar, class RangeScalar>
01440     void
01441     localMultiply (const MultiVector<DomainScalar,LocalOrdinal,GlobalOrdinal,node_type>& X,
01442                    MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,node_type>& Y,
01443                    Teuchos::ETransp trans,
01444                    RangeScalar alpha,
01445                    RangeScalar beta) const;
01446 
01471     template <class DomainScalar, class RangeScalar>
01472     void
01473     localGaussSeidel (const MultiVector<DomainScalar,LocalOrdinal,GlobalOrdinal,node_type> &B,
01474                       MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,node_type> &X,
01475                       const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type> &D,
01476                       const RangeScalar& dampingFactor,
01477                       const KokkosClassic::ESweepDirection direction) const;
01478 
01505     template <class DomainScalar, class RangeScalar>
01506     void
01507     reorderedLocalGaussSeidel (const MultiVector<DomainScalar,LocalOrdinal,GlobalOrdinal,node_type>& B,
01508                                MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,node_type>& X,
01509                                const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>& D,
01510                                const ArrayView<LocalOrdinal>& rowIndices,
01511                                const RangeScalar& dampingFactor,
01512                                const KokkosClassic::ESweepDirection direction) const;
01513 
01530     template <class DomainScalar, class RangeScalar>
01531     void
01532     localSolve (const MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,node_type>& Y,
01533                 MultiVector<DomainScalar,LocalOrdinal,GlobalOrdinal,node_type>& X,
01534                 Teuchos::ETransp trans) const;
01535 
01537     template <class T>
01538     RCP<CrsMatrix<T,LocalOrdinal,GlobalOrdinal,node_type> > convert() const;
01539 
01541 
01542 
01543 
01554     void
01555     apply (const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>& X,
01556            MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>&Y,
01557            Teuchos::ETransp mode = Teuchos::NO_TRANS,
01558            Scalar alpha = ScalarTraits<Scalar>::one(),
01559            Scalar beta = ScalarTraits<Scalar>::zero()) const;
01560 
01562     bool hasTransposeApply() const;
01563 
01567     RCP<const map_type> getDomainMap () const;
01568 
01572     RCP<const map_type> getRangeMap () const;
01573 
01575 
01576 
01577 
01642     void
01643     gaussSeidel (const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type> &B,
01644                  MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type> &X,
01645                  const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type> &D,
01646                  const Scalar& dampingFactor,
01647                  const ESweepDirection direction,
01648                  const int numSweeps) const;
01649 
01716     void
01717     reorderedGaussSeidel (const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>& B,
01718                           MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>& X,
01719                           const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>& D,
01720                           const ArrayView<LocalOrdinal>& rowIndices,
01721                           const Scalar& dampingFactor,
01722                           const ESweepDirection direction,
01723                           const int numSweeps) const;
01724 
01753     void
01754     gaussSeidelCopy (MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type> &X,
01755                      const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type> &B,
01756                      const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type> &D,
01757                      const Scalar& dampingFactor,
01758                      const ESweepDirection direction,
01759                      const int numSweeps,
01760                      const bool zeroInitialGuess) const;
01761 
01791     void
01792     reorderedGaussSeidelCopy (MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>& X,
01793                               const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>& B,
01794                               const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>& D,
01795                               const ArrayView<LocalOrdinal>& rowIndices,
01796                               const Scalar& dampingFactor,
01797                               const ESweepDirection direction,
01798                               const int numSweeps,
01799                               const bool zeroInitialGuess) const;
01800 
01811     virtual Teuchos::RCP<RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, node_type> >
01812     add (const Scalar& alpha,
01813          const RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, node_type>& A,
01814          const Scalar& beta,
01815          const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, node_type> >& domainMap,
01816          const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, node_type> >& rangeMap,
01817          const Teuchos::RCP<Teuchos::ParameterList>& params) const;
01818 
01820 
01821 
01822 
01824     std::string description() const;
01825 
01827     void
01828     describe (Teuchos::FancyOStream &out,
01829               const Teuchos::EVerbosityLevel verbLevel =
01830               Teuchos::Describable::verbLevel_default) const;
01831 
01833 
01834 
01835 
01836     virtual bool
01837     checkSizes (const SrcDistObject& source);
01838 
01839     virtual void
01840     copyAndPermute (const SrcDistObject& source,
01841                     size_t numSameIDs,
01842                     const ArrayView<const LocalOrdinal> &permuteToLIDs,
01843                     const ArrayView<const LocalOrdinal> &permuteFromLIDs);
01844 
01845     virtual void
01846     packAndPrepare (const SrcDistObject& source,
01847                     const Teuchos::ArrayView<const LocalOrdinal>& exportLIDs,
01848                     Teuchos::Array<char>& exports,
01849                     const Teuchos::ArrayView<size_t>& numPacketsPerLID,
01850                     size_t& constantNumPackets,
01851                     Distributor& distor);
01852 
01862     void
01863     unpackAndCombine (const Teuchos::ArrayView<const LocalOrdinal> &importLIDs,
01864                       const Teuchos::ArrayView<const char> &imports,
01865                       const Teuchos::ArrayView<size_t> &numPacketsPerLID,
01866                       size_t constantNumPackets,
01867                       Distributor &distor,
01868                       CombineMode combineMode);
01870 
01871 
01872 
01881     virtual void
01882     pack (const Teuchos::ArrayView<const LocalOrdinal>& exportLIDs,
01883           Teuchos::Array<char>& exports,
01884           const Teuchos::ArrayView<size_t>& numPacketsPerLID,
01885           size_t& constantNumPackets,
01886           Distributor& distor) const;
01887 
01889   private:
01890     // Friend declaration for nonmember function.
01891     template<class CrsMatrixType>
01892     friend Teuchos::RCP<CrsMatrixType>
01893     importAndFillCompleteCrsMatrix (const Teuchos::RCP<const CrsMatrixType>& sourceMatrix,
01894                                     const Import<typename CrsMatrixType::local_ordinal_type,
01895                                                  typename CrsMatrixType::global_ordinal_type,
01896                                                  typename CrsMatrixType::node_type>& importer,
01897                                     const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
01898                                                                  typename CrsMatrixType::global_ordinal_type,
01899                                                                  typename CrsMatrixType::node_type> >& domainMap,
01900                                     const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
01901                                                                  typename CrsMatrixType::global_ordinal_type,
01902                                                                  typename CrsMatrixType::node_type> >& rangeMap,
01903                                     const Teuchos::RCP<Teuchos::ParameterList>& params);
01904 
01905     // Friend declaration for nonmember function.
01906     template<class CrsMatrixType>
01907     friend Teuchos::RCP<CrsMatrixType>
01908     exportAndFillCompleteCrsMatrix (const Teuchos::RCP<const CrsMatrixType>& sourceMatrix,
01909                                     const Export<typename CrsMatrixType::local_ordinal_type,
01910                                                  typename CrsMatrixType::global_ordinal_type,
01911                                                  typename CrsMatrixType::node_type>& exporter,
01912                                     const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
01913                                                                  typename CrsMatrixType::global_ordinal_type,
01914                                                                  typename CrsMatrixType::node_type> >& domainMap,
01915                                     const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
01916                                                                  typename CrsMatrixType::global_ordinal_type,
01917                                                                  typename CrsMatrixType::node_type> >& rangeMap,
01918                                     const Teuchos::RCP<Teuchos::ParameterList>& params);
01919 
01926     Teuchos::RCP<CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, node_type, LocalMatOps> >
01927     importAndFillComplete (const Import<LocalOrdinal, GlobalOrdinal, node_type>& importer,
01928                            const Teuchos::RCP<const map_type>& domainMap,
01929                            const Teuchos::RCP<const map_type>& rangeMap,
01930                            const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null) const;
01931 
01938     Teuchos::RCP<CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,node_type,LocalMatOps> >
01939     exportAndFillComplete (const Export<LocalOrdinal, GlobalOrdinal, node_type>& exporter,
01940                            const Teuchos::RCP<const map_type>& domainMap = Teuchos::null,
01941                            const Teuchos::RCP<const map_type>& rangeMap = Teuchos::null,
01942                            const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null) const;
01943 
01949     template<class TransferType>
01950     Teuchos::RCP<CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,node_type,LocalMatOps> >
01951     transferAndFillComplete (const TransferType& rowTransfer,
01952                              const Teuchos::RCP<const map_type>& domainMap = Teuchos::null,
01953                              const Teuchos::RCP<const map_type>& rangeMap = Teuchos::null,
01954                              const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null) const;
01955 
01956 
01957     // We forbid copy construction by declaring this method private
01958     // and not implementing it.
01959     CrsMatrix (const CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,node_type,LocalMatOps> &rhs);
01960 
01961     // We forbid assignment (operator=) by declaring this method
01962     // private and not implementing it.
01963     CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,node_type,LocalMatOps>&
01964     operator= (const CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,node_type,LocalMatOps> &rhs);
01965 
01971     void
01972     insertGlobalValuesFiltered (const GlobalOrdinal globalRow,
01973                                 const ArrayView<const GlobalOrdinal> &indices,
01974                                 const ArrayView<const Scalar>        &values);
01975 
01981     void
01982     insertLocalValuesFiltered (const LocalOrdinal localRow,
01983                                const ArrayView<const LocalOrdinal> &indices,
01984                                const ArrayView<const Scalar>       &values);
01985 
01993     void
01994     combineGlobalValues (const GlobalOrdinal globalRowIndex,
01995                          const Teuchos::ArrayView<const GlobalOrdinal> columnIndices,
01996                          const Teuchos::ArrayView<const Scalar> values,
01997                          const Tpetra::CombineMode combineMode);
01998 
02030     template<class BinaryFunction>
02031     void
02032     transformLocalValues (LocalOrdinal localRow,
02033                           const Teuchos::ArrayView<const LocalOrdinal>& indices,
02034                           const Teuchos::ArrayView<const Scalar>        & values,
02035                           BinaryFunction f)
02036     {
02037       typedef LocalOrdinal LO;
02038       typedef GlobalOrdinal GO;
02039       typedef node_type NT;
02040       using Teuchos::Array;
02041       using Teuchos::ArrayView;
02042 
02043       TEUCHOS_TEST_FOR_EXCEPTION(
02044         ! isFillActive (),
02045         std::runtime_error,
02046         "Tpetra::CrsMatrix::transformLocalValues: Fill must be active in order "
02047         "to call this method.  That is, isFillActive() must return true.  If "
02048         "you have already called fillComplete(), you need to call resumeFill() "
02049         "before you can replace values.");
02050       TEUCHOS_TEST_FOR_EXCEPTION(
02051         values.size () != indices.size (),
02052         std::runtime_error,
02053         "Tpetra::CrsMatrix::transformLocalValues: values.size () = "
02054         << values.size () << " != indices.size () = " << indices.size ()
02055         << ".");
02056       TEUCHOS_TEST_FOR_EXCEPTION(
02057         ! this->hasColMap (),
02058         std::runtime_error,
02059         "Tpetra::CrsMatrix::transformLocalValues: We cannot transform local "
02060         "indices without a column map.");
02061       const bool isLocalRow = getRowMap ()->isNodeLocalElement (localRow);
02062       TEUCHOS_TEST_FOR_EXCEPTION(
02063         ! isLocalRow,
02064         std::runtime_error,
02065         "Tpetra::CrsMatrix::transformLocalValues: The specified local row "
02066         << localRow << " does not belong to this process "
02067         << getRowMap ()->getComm ()->getRank () << ".");
02068 
02069       RowInfo rowInfo = staticGraph_->getRowInfo (localRow);
02070       if (indices.size () > 0) {
02071         ArrayView<Scalar> curVals = this->getViewNonConst (rowInfo);
02072         if (isLocallyIndexed ()) {
02073           staticGraph_->template transformLocalValues<Scalar, BinaryFunction> (rowInfo, curVals, indices, values, f);
02074         }
02075         else if (isGloballyIndexed ()) {
02076           // Convert the given local indices to global indices.
02077           const Map<LO, GO, NT>& colMap = * (this->getColMap ());
02078           Array<GO> gindices (indices.size ());
02079           typename ArrayView<const LO>::iterator lindit = indices.begin();
02080           typename Array<GO>::iterator           gindit = gindices.begin();
02081           while (lindit != indices.end()) {
02082             // There is no need to filter out indices not in the column
02083             // Map.  Those that aren't will be mapped to invalid(),
02084             // which transformGlobalValues() will ignore.
02085             *gindit++ = colMap.getGlobalElement (*lindit++);
02086           }
02087           staticGraph_->template transformGlobalValues<Scalar, BinaryFunction> (rowInfo, curVals, gindices (), values, f);
02088         }
02089       }
02090     }
02091 
02092 
02122     template<class BinaryFunction>
02123     void
02124     transformGlobalValues (GlobalOrdinal globalRow,
02125                            const Teuchos::ArrayView<const GlobalOrdinal>& indices,
02126                            const Teuchos::ArrayView<const Scalar>        & values,
02127                            BinaryFunction f)
02128     {
02129       typedef LocalOrdinal LO;
02130       typedef GlobalOrdinal GO;
02131       typedef node_type NT;
02132       using Teuchos::Array;
02133       using Teuchos::ArrayView;
02134 
02135       TEUCHOS_TEST_FOR_EXCEPTION(
02136         ! isFillActive (),
02137         std::runtime_error,
02138         "Tpetra::CrsMatrix::transformGlobalValues: Fill must be active in order "
02139         "to call this method.  That is, isFillActive() must return true.  If "
02140         "you have already called fillComplete(), you need to call resumeFill() "
02141         "before you can replace values.");
02142       TEUCHOS_TEST_FOR_EXCEPTION(
02143         values.size () != indices.size (),
02144         std::runtime_error,
02145         "Tpetra::CrsMatrix::transformGlobalValues: values.size () = "
02146         << values.size () << " != indices.size () = " << indices.size ()
02147         << ".");
02148 
02149       const LO lrow = this->getRowMap()->getLocalElement(globalRow);
02150 
02151       if (lrow == OTL::invalid()) {
02152         // FIXME (mfh 16 May 2013) We're using this exception to do
02153         // sumIntoGlobalValues for nonowned rows, so we might want to
02154         // avoid the overhead of constructing the fancy exception
02155         // message each time if we don't plan to use it.
02156 
02157         // The exception test macro doesn't let you pass an additional
02158         // argument to the exception's constructor, so we don't use it.
02159         std::ostringstream os;
02160         os << "transformGlobalValues: The given global row index "
02161            << globalRow << " is not owned by the calling process (rank "
02162            << this->getRowMap()->getComm()->getRank() << ").";
02163         throw Details::InvalidGlobalRowIndex<GO> (os.str (), globalRow);
02164       }
02165 
02166       RowInfo rowInfo = staticGraph_->getRowInfo (lrow);
02167       if (indices.size () > 0) {
02168         ArrayView<Scalar> curVals = this->getViewNonConst (rowInfo);
02169         if (isLocallyIndexed ()) {
02170           // Convert global indices to local indices.
02171           const Map<LO, GO, NT> &colMap = * (this->getColMap ());
02172           Array<LO> lindices (indices.size ());
02173           typename ArrayView<const GO>::iterator gindit = indices.begin();
02174           typename Array<LO>::iterator           lindit = lindices.begin();
02175           while (gindit != indices.end()) {
02176             // There is no need to filter out indices not in the column
02177             // Map.  Those that aren't will be mapped to invalid(),
02178             // which transformLocalValues() will ignore.
02179             *lindit++ = colMap.getLocalElement (*gindit++);
02180           }
02181           staticGraph_->template transformLocalValues<Scalar, BinaryFunction> (rowInfo, curVals, lindices (), values, f);
02182         }
02183         else if (isGloballyIndexed ()) {
02184           staticGraph_->template transformGlobalValues<Scalar, BinaryFunction> (rowInfo, curVals, indices, values, f);
02185         }
02186       }
02187     }
02188 
02189   private:
02192     void
02193     insertNonownedGlobalValues (const GlobalOrdinal globalRow,
02194                                 const Teuchos::ArrayView<const GlobalOrdinal>& indices,
02195                                 const Teuchos::ArrayView<const Scalar>& values);
02196 
02197   protected:
02198     // useful typedefs
02199     typedef OrdinalTraits<LocalOrdinal>                     OTL;
02200     typedef ScalarTraits<Scalar>                            STS;
02201     typedef typename STS::magnitudeType               Magnitude;
02202     typedef ScalarTraits<Magnitude>                         STM;
02203     typedef MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type> MV;
02204     typedef Vector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>      V;
02205     typedef crs_graph_type Graph;
02206     typedef typename LocalMatOps::template bind_scalar<Scalar>::other_type                    sparse_ops_type;
02207     typedef typename sparse_ops_type::template graph<LocalOrdinal,node_type>::graph_type          local_graph_type;
02208     typedef typename sparse_ops_type::template matrix<Scalar,LocalOrdinal,node_type>::matrix_type local_matrix_type;
02209 
02210     typedef Export<LocalOrdinal, GlobalOrdinal, node_type> export_type;
02211     typedef Import<LocalOrdinal, GlobalOrdinal, node_type> import_type;
02212 
02213     // Enums
02214     enum GraphAllocationStatus {
02215       GraphAlreadyAllocated,
02216       GraphNotYetAllocated
02217     };
02218 
02235     void allocateValues (ELocalGlobal lg, GraphAllocationStatus gas);
02236 
02242     void sortEntries();
02243 
02249     void mergeRedundantEntries();
02250 
02258     void clearGlobalConstants();
02259 
02268     void computeGlobalConstants();
02269 
02282     mutable Teuchos::RCP<MV> importMV_;
02283 
02296     mutable Teuchos::RCP<MV> exportMV_;
02297 
02317     Teuchos::RCP<MV>
02318     getColumnMapMultiVector (const MV& X_domainMap,
02319                              const bool force = false) const;
02320 
02342     Teuchos::RCP<MV>
02343     getRowMapMultiVector (const MV& Y_rangeMap,
02344                           const bool force = false) const;
02345 
02347     void
02348     applyNonTranspose (const MV& X_in,
02349                        MV& Y_in,
02350                        Scalar alpha,
02351                        Scalar beta) const;
02352 
02354     void
02355     applyTranspose (const MV& X_in,
02356                     MV& Y_in,
02357                     const Teuchos::ETransp mode,
02358                     Scalar alpha,
02359                     Scalar beta) const;
02360 
02361     // matrix data accessors
02362     ArrayView<const Scalar>    getView(RowInfo rowinfo) const;
02363     ArrayView<      Scalar>    getViewNonConst(RowInfo rowinfo);
02364     // local Kokkos objects
02365 
02371     void fillLocalMatrix(const RCP<ParameterList> &params);
02377     void fillLocalGraphAndMatrix(const RCP<ParameterList> &params);
02378     // debugging
02379     void checkInternalState() const;
02380 
02392 
02393     RCP<const Graph> staticGraph_;
02394     RCP<      Graph>     myGraph_;
02396 
02401     RCP<sparse_ops_type>   lclMatOps_;
02408     RCP<local_matrix_type> lclMatrix_;
02409 
02422 
02423     ArrayRCP<Scalar> values1D_;
02424     ArrayRCP<Array<Scalar> > values2D_;
02426 
02428     bool fillComplete_;
02429 
02457     std::map<GlobalOrdinal, std::pair<Teuchos::Array<GlobalOrdinal>,
02458                                       Teuchos::Array<Scalar> > > nonlocals_;
02459 
02466     mutable Magnitude frobNorm_;
02467   }; // class CrsMatrix
02468 
02475   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
02476   Teuchos::RCP<CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
02477   createCrsMatrix (const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &map,
02478                    size_t maxNumEntriesPerRow = 0,
02479                    const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null)
02480   {
02481     using Teuchos::rcp;
02482     typedef CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> matrix_type;
02483     return rcp (new matrix_type (map, maxNumEntriesPerRow, DynamicProfile, params));
02484   }
02485 
02537   template<class CrsMatrixType>
02538   Teuchos::RCP<CrsMatrixType>
02539   importAndFillCompleteCrsMatrix (const Teuchos::RCP<const CrsMatrixType>& sourceMatrix,
02540                                   const Import<typename CrsMatrixType::local_ordinal_type,
02541                                                typename CrsMatrixType::global_ordinal_type,
02542                                                typename CrsMatrixType::node_type>& importer,
02543                                   const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
02544                                                                typename CrsMatrixType::global_ordinal_type,
02545                                                                typename CrsMatrixType::node_type> >& domainMap = Teuchos::null,
02546                                   const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
02547                                                                typename CrsMatrixType::global_ordinal_type,
02548                                                                typename CrsMatrixType::node_type> >& rangeMap = Teuchos::null,
02549                                   const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null)
02550   {
02551     return sourceMatrix->importAndFillComplete (importer, domainMap, rangeMap, params);
02552   }
02553 
02587   template<class CrsMatrixType>
02588   Teuchos::RCP<CrsMatrixType>
02589   exportAndFillCompleteCrsMatrix (const Teuchos::RCP<const CrsMatrixType>& sourceMatrix,
02590                                   const Export<typename CrsMatrixType::local_ordinal_type,
02591                                                typename CrsMatrixType::global_ordinal_type,
02592                                                typename CrsMatrixType::node_type>& exporter,
02593                                   const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
02594                                                                typename CrsMatrixType::global_ordinal_type,
02595                                                                typename CrsMatrixType::node_type> >& domainMap = Teuchos::null,
02596                                   const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
02597                                                                typename CrsMatrixType::global_ordinal_type,
02598                                                                typename CrsMatrixType::node_type> >& rangeMap = Teuchos::null,
02599                                   const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null)
02600   {
02601     return sourceMatrix->exportAndFillComplete (exporter, domainMap, rangeMap, params);
02602   }
02603 } // namespace Tpetra
02604 
02605 // Include KokkosRefactor partial specialisation if enabled
02606 #if defined(TPETRA_HAVE_KOKKOS_REFACTOR)
02607 #include "Tpetra_KokkosRefactor_CrsMatrix_decl.hpp"
02608 #endif
02609 
02620 #endif // TPETRA_CRSMATRIX_DECL_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines