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 Import<LocalOrdinal, GlobalOrdinal, node_type> import_type;
00233 
00235     typedef Export<LocalOrdinal, GlobalOrdinal, node_type> export_type;
00236 
00238     typedef CrsGraph<LocalOrdinal, GlobalOrdinal, node_type, LocalMatOps> crs_graph_type;
00239 
00241 
00242 
00243 
00261     CrsMatrix (const Teuchos::RCP<const map_type>& rowMap,
00262                size_t maxNumEntriesPerRow,
00263                ProfileType pftype = DynamicProfile,
00264                const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null);
00265 
00283     CrsMatrix (const Teuchos::RCP<const map_type>& rowMap,
00284                const ArrayRCP<const size_t>& NumEntriesPerRowToAlloc,
00285                ProfileType pftype = DynamicProfile,
00286                const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null);
00287 
00310     CrsMatrix (const Teuchos::RCP<const map_type>& rowMap,
00311                const Teuchos::RCP<const map_type>& colMap,
00312                size_t maxNumEntriesPerRow,
00313                ProfileType pftype = DynamicProfile,
00314                const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null);
00315 
00338     CrsMatrix (const Teuchos::RCP<const map_type>& rowMap,
00339                const Teuchos::RCP<const map_type>& colMap,
00340                const ArrayRCP<const size_t>& NumEntriesPerRowToAlloc,
00341                ProfileType pftype = DynamicProfile,
00342                const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null);
00343 
00368     explicit CrsMatrix (const Teuchos::RCP<const crs_graph_type>& graph,
00369                         const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null);
00370 
00393     CrsMatrix (const RCP<const map_type>& rowMap,
00394                const RCP<const map_type>& colMap,
00395                const ArrayRCP<size_t>& rowPointers,
00396                const ArrayRCP<LocalOrdinal>& columnIndices,
00397                const ArrayRCP<Scalar>& values,
00398                const RCP<ParameterList>& params = null);
00399 
00400     // This friend declaration makes the clone() method work.
00401     template <class S2, class LO2, class GO2, class N2, class LMO2>
00402     friend class CrsMatrix;
00403 
00428     template <class Node2>
00429     RCP<CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node2, typename KokkosClassic::DefaultKernels<void,LocalOrdinal,Node2>::SparseOps> >
00430     clone (const RCP<Node2>& node2, const RCP<ParameterList>& params = null) const
00431     {
00432       const char tfecfFuncName[] = "clone";
00433 
00434       // Get parameter values.  Set them initially to their default values.
00435       bool fillCompleteClone = true;
00436       bool useLocalIndices = this->hasColMap ();
00437       ProfileType pftype = StaticProfile;
00438       if (! params.is_null ()) {
00439         fillCompleteClone = params->get ("fillComplete clone", fillCompleteClone);
00440         useLocalIndices = params->get ("Locally indexed clone", useLocalIndices);
00441 
00442         bool staticProfileClone = true;
00443         staticProfileClone = params->get ("Static profile clone", staticProfileClone);
00444         pftype = staticProfileClone ? StaticProfile : DynamicProfile;
00445       }
00446 
00447       TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00448         ! this->hasColMap () && useLocalIndices, std::runtime_error,
00449         ": You requested that the returned clone have local indices, but the "
00450         "the source matrix does not have a column Map yet.");
00451 
00452       typedef typename KokkosClassic::DefaultKernels<void,LocalOrdinal,Node2>::SparseOps LocalMatOps2;
00453       typedef CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node2, LocalMatOps2> CrsMatrix2;
00454       typedef Map<LocalOrdinal, GlobalOrdinal, Node2> Map2;
00455       RCP<const Map2> clonedRowMap = this->getRowMap ()->template clone<Node2> (node2);
00456 
00457       RCP<CrsMatrix2> clonedMatrix;
00458       ArrayRCP<const size_t> numEntries;
00459       size_t numEntriesForAll = 0;
00460       if (! staticGraph_->indicesAreAllocated ()) {
00461         if (! staticGraph_->numAllocPerRow_.is_null ()) {
00462           numEntries = staticGraph_->numAllocPerRow_;
00463         }
00464         else {
00465           numEntriesForAll = staticGraph_->numAllocForAllRows_;
00466         }
00467       }
00468       else if (! staticGraph_->numRowEntries_.is_null ()) {
00469         numEntries = staticGraph_->numRowEntries_;
00470       }
00471       else if (staticGraph_->nodeNumAllocated_ == 0) {
00472         numEntriesForAll = 0;
00473       }
00474       else {
00475         // We're left with the case that we have optimized storage.
00476         // In this case, we have to construct a list of row sizes.
00477         TEUCHOS_TEST_FOR_EXCEPTION(
00478           getProfileType() != StaticProfile, std::logic_error,
00479           "Internal logic error. Please report this to Tpetra team." )
00480 
00481         const size_t numRows = this->getNodeNumRows ();
00482         numEntriesForAll = 0;
00483         ArrayRCP<size_t> numEnt;
00484         if (numRows != 0) {
00485           numEnt = arcp<size_t> (numRows);
00486         }
00487         for (size_t i = 0; i < numRows; ++i) {
00488           numEnt[i] = staticGraph_->rowPtrs_[i+1] - staticGraph_->rowPtrs_[i];
00489         }
00490         numEntries = numEnt;
00491       }
00492 
00493       RCP<ParameterList> matrixparams =
00494         params.is_null () ? null : sublist (params,"CrsMatrix");
00495       if (useLocalIndices) {
00496         RCP<const Map2> clonedColMap =
00497           this->getColMap ()->template clone<Node2> (node2);
00498         if (numEntries.is_null ()) {
00499           clonedMatrix = rcp (new CrsMatrix2 (clonedRowMap, clonedColMap,
00500                                               numEntriesForAll, pftype,
00501                                               matrixparams));
00502         }
00503         else {
00504           clonedMatrix = rcp (new CrsMatrix2 (clonedRowMap, clonedColMap,
00505                                               numEntries, pftype,
00506                                               matrixparams));
00507         }
00508       }
00509       else {
00510         if (numEntries.is_null ()) {
00511           clonedMatrix = rcp (new CrsMatrix2 (clonedRowMap, numEntriesForAll,
00512                                               pftype, matrixparams));
00513         }
00514         else {
00515           clonedMatrix = rcp (new CrsMatrix2 (clonedRowMap, numEntries, pftype,
00516                                               matrixparams));
00517         }
00518       }
00519       // done with these
00520       numEntries = null;
00521       numEntriesForAll = 0;
00522 
00523       if (useLocalIndices) {
00524         clonedMatrix->allocateValues (LocalIndices,
00525                                       CrsMatrix2::GraphNotYetAllocated);
00526         if (this->isLocallyIndexed ()) {
00527           ArrayView<const LocalOrdinal> linds;
00528           ArrayView<const Scalar>       vals;
00529           for (LocalOrdinal lrow = clonedRowMap->getMinLocalIndex ();
00530                lrow <= clonedRowMap->getMaxLocalIndex ();
00531                ++lrow) {
00532             this->getLocalRowView (lrow, linds, vals);
00533             if (linds.size ()) {
00534               clonedMatrix->insertLocalValues (lrow, linds, vals);
00535             }
00536           }
00537         }
00538         else { // this->isGloballyIndexed()
00539           Array<LocalOrdinal> linds;
00540           Array<Scalar> vals;
00541           for (LocalOrdinal lrow = clonedRowMap->getMinLocalIndex ();
00542                lrow <= clonedRowMap->getMaxLocalIndex ();
00543                ++lrow) {
00544             size_t theNumEntries = this->getNumEntriesInLocalRow (lrow);
00545             if (theNumEntries > Teuchos::as<size_t> (linds.size ())) {
00546               linds.resize (theNumEntries);
00547             }
00548             if (theNumEntries > Teuchos::as<size_t> (vals.size ())) {
00549               vals.resize (theNumEntries);
00550             }
00551             this->getLocalRowCopy (clonedRowMap->getGlobalElement (lrow),
00552                                    linds (), vals (), theNumEntries);
00553             if (theNumEntries != 0) {
00554               clonedMatrix->insertLocalValues (lrow, linds (0, theNumEntries),
00555                                                vals (0, theNumEntries));
00556             }
00557           }
00558         }
00559       }
00560       else { // useGlobalIndices
00561         clonedMatrix->allocateValues (GlobalIndices,
00562                                       CrsMatrix2::GraphNotYetAllocated);
00563         if (this->isGloballyIndexed ()) {
00564           ArrayView<const GlobalOrdinal> ginds;
00565           ArrayView<const Scalar>         vals;
00566           for (GlobalOrdinal grow = clonedRowMap->getMinGlobalIndex ();
00567                grow <= clonedRowMap->getMaxGlobalIndex ();
00568                ++grow) {
00569             this->getGlobalRowView (grow, ginds, vals);
00570             if (ginds.size () > 0) {
00571               clonedMatrix->insertGlobalValues (grow, ginds, vals);
00572             }
00573           }
00574         }
00575         else { // this->isLocallyIndexed()
00576           Array<GlobalOrdinal> ginds;
00577           Array<Scalar> vals;
00578           for (GlobalOrdinal grow = clonedRowMap->getMinGlobalIndex ();
00579                grow <= clonedRowMap->getMaxGlobalIndex ();
00580                ++grow) {
00581             size_t theNumEntries = this->getNumEntriesInGlobalRow (grow);
00582             if (theNumEntries > Teuchos::as<size_t> (ginds.size ())) {
00583               ginds.resize (theNumEntries);
00584             }
00585             if (theNumEntries > Teuchos::as<size_t> (vals.size ())) {
00586               vals.resize (theNumEntries);
00587             }
00588             this->getGlobalRowCopy (grow, ginds (), vals (), theNumEntries);
00589             if (theNumEntries != 0) {
00590               clonedMatrix->insertGlobalValues (grow, ginds (0, theNumEntries),
00591                                                 vals (0, theNumEntries));
00592             }
00593           }
00594         }
00595       }
00596 
00597       if (fillCompleteClone) {
00598         RCP<ParameterList> fillparams =
00599           params.is_null () ? Teuchos::null : sublist (params, "fillComplete");
00600         try {
00601           RCP<const Map2> clonedRangeMap;
00602           RCP<const Map2> clonedDomainMap;
00603           if (! this->getRangeMap ().is_null () &&
00604               this->getRangeMap () != clonedRowMap) {
00605             clonedRangeMap  = this->getRangeMap ()->template clone<Node2> (node2);
00606           }
00607           else {
00608             clonedRangeMap = clonedRowMap;
00609           }
00610           if (! this->getDomainMap ().is_null () &&
00611               this->getDomainMap () != clonedRowMap) {
00612             clonedDomainMap = this->getDomainMap ()->template clone<Node2> (node2);
00613           }
00614           else {
00615             clonedDomainMap = clonedRowMap;
00616           }
00617           clonedMatrix->fillComplete (clonedDomainMap, clonedRangeMap,
00618                                       fillparams);
00619         }
00620         catch (std::exception &e) {
00621           const bool caughtExceptionOnClone = true;
00622           TEUCHOS_TEST_FOR_EXCEPTION(
00623             caughtExceptionOnClone, std::runtime_error,
00624             Teuchos::typeName (*this) << std::endl << "clone: " << std::endl <<
00625             "Caught the following exception while calling fillComplete() on a "
00626             "clone of type" << std::endl << Teuchos::typeName (*clonedMatrix)
00627             << ": " << std::endl << e.what () << std::endl);
00628         }
00629       }
00630       return clonedMatrix;
00631     }
00632 
00634     virtual ~CrsMatrix ();
00635 
00637 
00638 
00639 
00662     //
00708     void
00709     insertGlobalValues (const GlobalOrdinal globalRow,
00710                         const ArrayView<const GlobalOrdinal>& cols,
00711                         const ArrayView<const Scalar>& vals);
00712 
00752     void
00753     insertLocalValues (const LocalOrdinal localRow,
00754                        const ArrayView<const LocalOrdinal> &cols,
00755                        const ArrayView<const Scalar> &vals);
00756 
00791     LocalOrdinal
00792     replaceGlobalValues (GlobalOrdinal globalRow,
00793                          const ArrayView<const GlobalOrdinal>& cols,
00794                          const ArrayView<const Scalar>& vals);
00795 
00826     LocalOrdinal
00827     replaceLocalValues (const LocalOrdinal localRow,
00828                         const ArrayView<const LocalOrdinal>&cols,
00829                         const ArrayView<const Scalar>& vals);
00830 
00863     LocalOrdinal
00864     sumIntoGlobalValues (const GlobalOrdinal globalRow,
00865                          const ArrayView<const GlobalOrdinal> &cols,
00866                          const ArrayView<const Scalar>        &vals);
00867 
00882     LocalOrdinal
00883     sumIntoLocalValues (const LocalOrdinal localRow,
00884                         const ArrayView<const LocalOrdinal>& cols,
00885                         const ArrayView<const Scalar>& vals);
00886 
00888     void setAllToScalar (const Scalar &alpha);
00889 
00891     void scale (const Scalar &alpha);
00892 
00894 
00901     void
00902     setAllValues (const ArrayRCP<size_t>& rowPointers,
00903                   const ArrayRCP<LocalOrdinal>& columnIndices,
00904                   const ArrayRCP<Scalar>& values);
00905 
00907 
00914     void
00915     getAllValues (ArrayRCP<const size_t>& rowPointers,
00916                   ArrayRCP<const LocalOrdinal>& columnIndices,
00917                   ArrayRCP<const Scalar>& values) const;
00918 
00920 
00921 
00922 
00951     void globalAssemble();
00952 
00966     void resumeFill (const RCP<ParameterList>& params = null);
00967 
01001     void
01002     fillComplete (const RCP<const map_type>& domainMap,
01003                   const RCP<const map_type>& rangeMap,
01004                   const RCP<ParameterList>& params = null);
01005 
01018     void fillComplete (const RCP<ParameterList>& params = null);
01019 
01031     void
01032     expertStaticFillComplete (const RCP<const map_type>& domainMap,
01033                               const RCP<const map_type>& rangeMap,
01034                               const RCP<const import_type>& importer = Teuchos::null,
01035                               const RCP<const export_type>& exporter = Teuchos::null,
01036                               const RCP<ParameterList>& params = Teuchos::null);
01037 
01047     void
01048     replaceColMap (const Teuchos::RCP<const map_type>& newColMap);
01049 
01062     void
01063     replaceDomainMapAndImporter (const Teuchos::RCP<const map_type>& newDomainMap,
01064                                  Teuchos::RCP<const import_type>& newImporter);
01065 
01079     virtual void
01080     removeEmptyProcessesInPlace (const Teuchos::RCP<const map_type>& newMap);
01081 
01083 
01084 
01085 
01087     RCP<const Comm<int> > getComm() const;
01088 
01090     RCP<node_type> getNode () const;
01091 
01093     RCP<const map_type> getRowMap () const;
01094 
01096     RCP<const map_type> getColMap () const;
01097 
01099     RCP<const RowGraph<LocalOrdinal,GlobalOrdinal,node_type> > getGraph () const;
01100 
01102     RCP<const crs_graph_type> getCrsGraph () const;
01103 
01123     global_size_t getGlobalNumRows() const;
01124 
01130     global_size_t getGlobalNumCols() const;
01131 
01138     size_t getNodeNumRows() const;
01139 
01143     size_t getNodeNumCols() const;
01144 
01146     GlobalOrdinal getIndexBase() const;
01147 
01149     global_size_t getGlobalNumEntries() const;
01150 
01152     size_t getNodeNumEntries() const;
01153 
01155 
01156     size_t getNumEntriesInGlobalRow(GlobalOrdinal globalRow) const;
01157 
01159 
01160     size_t getNumEntriesInLocalRow(LocalOrdinal localRow) const;
01161 
01163 
01165     global_size_t getGlobalNumDiags() const;
01166 
01168 
01170     size_t getNodeNumDiags() const;
01171 
01173 
01175     size_t getGlobalMaxNumRowEntries() const;
01176 
01178 
01180     size_t getNodeMaxNumRowEntries() const;
01181 
01183     bool hasColMap() const;
01184 
01186 
01188     bool isLowerTriangular() const;
01189 
01191 
01193     bool isUpperTriangular() const;
01194 
01214     bool isLocallyIndexed() const;
01215 
01235     bool isGloballyIndexed() const;
01236 
01259     bool isFillComplete() const;
01260 
01283     bool isFillActive() const;
01284 
01286 
01292     bool isStorageOptimized() const;
01293 
01295     ProfileType getProfileType() const;
01296 
01298     bool isStaticGraph() const;
01299 
01311     typename ScalarTraits<Scalar>::magnitudeType getFrobeniusNorm() const;
01312 
01319     virtual bool supportsRowViews() const;
01320 
01322 
01332     void
01333     getGlobalRowCopy (GlobalOrdinal GlobalRow,
01334                       const ArrayView<GlobalOrdinal> &Indices,
01335                       const ArrayView<Scalar> &Values,
01336                       size_t &NumEntries) const;
01337 
01339 
01351     void
01352     getLocalRowCopy (LocalOrdinal LocalRow,
01353                      const ArrayView<LocalOrdinal> &Indices,
01354                      const ArrayView<Scalar> &Values,
01355                      size_t &NumEntries) const;
01356 
01370     void
01371     getGlobalRowView (GlobalOrdinal GlobalRow,
01372                       ArrayView<const GlobalOrdinal> &indices,
01373                       ArrayView<const Scalar> &values) const;
01374 
01388     void
01389     getLocalRowView (LocalOrdinal LocalRow,
01390                      ArrayView<const LocalOrdinal>& indices,
01391                      ArrayView<const Scalar>& values) const;
01392 
01398     void getLocalDiagCopy (Vector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>& diag) const;
01399 
01437     void getLocalDiagOffsets (Teuchos::ArrayRCP<size_t>& offsets) const;
01438 
01457     void
01458     getLocalDiagCopy (Vector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>& diag,
01459                       const Teuchos::ArrayView<const size_t>& offsets) const;
01460 
01462     void leftScale(const Vector<Scalar, LocalOrdinal, GlobalOrdinal, node_type>& x);
01463 
01465     void rightScale(const Vector<Scalar, LocalOrdinal, GlobalOrdinal, node_type>& x);
01466 
01468 
01469 
01470 
01519     template <class DomainScalar, class RangeScalar>
01520     void
01521     localMultiply (const MultiVector<DomainScalar,LocalOrdinal,GlobalOrdinal,node_type>& X,
01522                    MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,node_type>& Y,
01523                    Teuchos::ETransp trans,
01524                    RangeScalar alpha,
01525                    RangeScalar beta) const;
01526 
01551     template <class DomainScalar, class RangeScalar>
01552     void
01553     localGaussSeidel (const MultiVector<DomainScalar,LocalOrdinal,GlobalOrdinal,node_type> &B,
01554                       MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,node_type> &X,
01555                       const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type> &D,
01556                       const RangeScalar& dampingFactor,
01557                       const KokkosClassic::ESweepDirection direction) const;
01558 
01585     template <class DomainScalar, class RangeScalar>
01586     void
01587     reorderedLocalGaussSeidel (const MultiVector<DomainScalar,LocalOrdinal,GlobalOrdinal,node_type>& B,
01588                                MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,node_type>& X,
01589                                const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>& D,
01590                                const ArrayView<LocalOrdinal>& rowIndices,
01591                                const RangeScalar& dampingFactor,
01592                                const KokkosClassic::ESweepDirection direction) const;
01593 
01610     template <class DomainScalar, class RangeScalar>
01611     void
01612     localSolve (const MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,node_type>& Y,
01613                 MultiVector<DomainScalar,LocalOrdinal,GlobalOrdinal,node_type>& X,
01614                 Teuchos::ETransp trans) const;
01615 
01617     template <class T>
01618     RCP<CrsMatrix<T,LocalOrdinal,GlobalOrdinal,node_type> > convert() const;
01619 
01621 
01622 
01623 
01634     void
01635     apply (const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>& X,
01636            MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>&Y,
01637            Teuchos::ETransp mode = Teuchos::NO_TRANS,
01638            Scalar alpha = ScalarTraits<Scalar>::one(),
01639            Scalar beta = ScalarTraits<Scalar>::zero()) const;
01640 
01642     bool hasTransposeApply() const;
01643 
01647     RCP<const map_type> getDomainMap () const;
01648 
01652     RCP<const map_type> getRangeMap () const;
01653 
01655 
01656 
01657 
01722     void
01723     gaussSeidel (const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type> &B,
01724                  MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type> &X,
01725                  const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type> &D,
01726                  const Scalar& dampingFactor,
01727                  const ESweepDirection direction,
01728                  const int numSweeps) const;
01729 
01796     void
01797     reorderedGaussSeidel (const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>& B,
01798                           MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>& X,
01799                           const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>& D,
01800                           const ArrayView<LocalOrdinal>& rowIndices,
01801                           const Scalar& dampingFactor,
01802                           const ESweepDirection direction,
01803                           const int numSweeps) const;
01804 
01833     void
01834     gaussSeidelCopy (MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type> &X,
01835                      const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type> &B,
01836                      const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type> &D,
01837                      const Scalar& dampingFactor,
01838                      const ESweepDirection direction,
01839                      const int numSweeps,
01840                      const bool zeroInitialGuess) const;
01841 
01871     void
01872     reorderedGaussSeidelCopy (MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>& X,
01873                               const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>& B,
01874                               const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>& D,
01875                               const ArrayView<LocalOrdinal>& rowIndices,
01876                               const Scalar& dampingFactor,
01877                               const ESweepDirection direction,
01878                               const int numSweeps,
01879                               const bool zeroInitialGuess) const;
01880 
01891     virtual Teuchos::RCP<RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, node_type> >
01892     add (const Scalar& alpha,
01893          const RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, node_type>& A,
01894          const Scalar& beta,
01895          const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, node_type> >& domainMap,
01896          const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, node_type> >& rangeMap,
01897          const Teuchos::RCP<Teuchos::ParameterList>& params) const;
01898 
01900 
01901 
01902 
01904     std::string description() const;
01905 
01907     void
01908     describe (Teuchos::FancyOStream &out,
01909               const Teuchos::EVerbosityLevel verbLevel =
01910               Teuchos::Describable::verbLevel_default) const;
01911 
01913 
01914 
01915 
01916     virtual bool
01917     checkSizes (const SrcDistObject& source);
01918 
01919     virtual void
01920     copyAndPermute (const SrcDistObject& source,
01921                     size_t numSameIDs,
01922                     const ArrayView<const LocalOrdinal> &permuteToLIDs,
01923                     const ArrayView<const LocalOrdinal> &permuteFromLIDs);
01924 
01925     virtual void
01926     packAndPrepare (const SrcDistObject& source,
01927                     const Teuchos::ArrayView<const LocalOrdinal>& exportLIDs,
01928                     Teuchos::Array<char>& exports,
01929                     const Teuchos::ArrayView<size_t>& numPacketsPerLID,
01930                     size_t& constantNumPackets,
01931                     Distributor& distor);
01932 
01942     void
01943     unpackAndCombine (const Teuchos::ArrayView<const LocalOrdinal> &importLIDs,
01944                       const Teuchos::ArrayView<const char> &imports,
01945                       const Teuchos::ArrayView<size_t> &numPacketsPerLID,
01946                       size_t constantNumPackets,
01947                       Distributor &distor,
01948                       CombineMode combineMode);
01950 
01951 
01952 
01961     virtual void
01962     pack (const Teuchos::ArrayView<const LocalOrdinal>& exportLIDs,
01963           Teuchos::Array<char>& exports,
01964           const Teuchos::ArrayView<size_t>& numPacketsPerLID,
01965           size_t& constantNumPackets,
01966           Distributor& distor) const;
01967 
01969   private:
01970     // Friend declaration for nonmember function.
01971     template<class CrsMatrixType>
01972     friend Teuchos::RCP<CrsMatrixType>
01973     importAndFillCompleteCrsMatrix (const Teuchos::RCP<const CrsMatrixType>& sourceMatrix,
01974                                     const Import<typename CrsMatrixType::local_ordinal_type,
01975                                                  typename CrsMatrixType::global_ordinal_type,
01976                                                  typename CrsMatrixType::node_type>& importer,
01977                                     const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
01978                                                                  typename CrsMatrixType::global_ordinal_type,
01979                                                                  typename CrsMatrixType::node_type> >& domainMap,
01980                                     const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
01981                                                                  typename CrsMatrixType::global_ordinal_type,
01982                                                                  typename CrsMatrixType::node_type> >& rangeMap,
01983                                     const Teuchos::RCP<Teuchos::ParameterList>& params);
01984 
01985     // Friend declaration for nonmember function.
01986     template<class CrsMatrixType>
01987     friend Teuchos::RCP<CrsMatrixType>
01988     exportAndFillCompleteCrsMatrix (const Teuchos::RCP<const CrsMatrixType>& sourceMatrix,
01989                                     const Export<typename CrsMatrixType::local_ordinal_type,
01990                                                  typename CrsMatrixType::global_ordinal_type,
01991                                                  typename CrsMatrixType::node_type>& exporter,
01992                                     const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
01993                                                                  typename CrsMatrixType::global_ordinal_type,
01994                                                                  typename CrsMatrixType::node_type> >& domainMap,
01995                                     const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
01996                                                                  typename CrsMatrixType::global_ordinal_type,
01997                                                                  typename CrsMatrixType::node_type> >& rangeMap,
01998                                     const Teuchos::RCP<Teuchos::ParameterList>& params);
01999 
02000   public:
02016     void
02017     importAndFillComplete (Teuchos::RCP<CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,node_type,LocalMatOps> > & destMatrix,
02018                            const import_type& importer,
02019                            const Teuchos::RCP<const map_type>& domainMap = Teuchos::null,
02020                            const Teuchos::RCP<const map_type>& rangeMap = Teuchos::null,
02021                            const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null) const;
02022 
02038     void
02039     exportAndFillComplete (Teuchos::RCP<CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,node_type,LocalMatOps> > & destMatrix,
02040                            const export_type& exporter,
02041                            const Teuchos::RCP<const map_type>& domainMap = Teuchos::null,
02042                            const Teuchos::RCP<const map_type>& rangeMap = Teuchos::null,
02043                            const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null) const;
02044   private:
02065     void
02066     transferAndFillComplete (Teuchos::RCP<CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps> > & destMat,
02067                              const ::Tpetra::Details::Transfer<LocalOrdinal, GlobalOrdinal, Node>& rowTransfer, // either Import or Export
02068                              const Teuchos::RCP<const map_type>& domainMap = Teuchos::null,
02069                              const Teuchos::RCP<const map_type>& rangeMap = Teuchos::null,
02070                              const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null) const;
02071 
02072     // We forbid copy construction by declaring this method private
02073     // and not implementing it.
02074     CrsMatrix (const CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,node_type,LocalMatOps> &rhs);
02075 
02076     // We forbid assignment (operator=) by declaring this method
02077     // private and not implementing it.
02078     CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,node_type,LocalMatOps>&
02079     operator= (const CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,node_type,LocalMatOps> &rhs);
02080 
02086     void
02087     insertGlobalValuesFiltered (const GlobalOrdinal globalRow,
02088                                 const ArrayView<const GlobalOrdinal> &indices,
02089                                 const ArrayView<const Scalar>        &values);
02090 
02096     void
02097     insertLocalValuesFiltered (const LocalOrdinal localRow,
02098                                const ArrayView<const LocalOrdinal> &indices,
02099                                const ArrayView<const Scalar>       &values);
02100 
02108     void
02109     combineGlobalValues (const GlobalOrdinal globalRowIndex,
02110                          const Teuchos::ArrayView<const GlobalOrdinal> columnIndices,
02111                          const Teuchos::ArrayView<const Scalar> values,
02112                          const Tpetra::CombineMode combineMode);
02113 
02143     template<class BinaryFunction>
02144     LocalOrdinal
02145     transformGlobalValues (GlobalOrdinal globalRow,
02146                            const Teuchos::ArrayView<const GlobalOrdinal>& indices,
02147                            const Teuchos::ArrayView<const Scalar>        & values,
02148                            BinaryFunction f)
02149     {
02150       typedef LocalOrdinal LO;
02151       typedef GlobalOrdinal GO;
02152       using Teuchos::Array;
02153       using Teuchos::ArrayView;
02154       typedef typename ArrayView<const GO>::size_type size_type;
02155 
02156       if (! isFillActive ()) {
02157         // Fill must be active in order to call this method.
02158         return Teuchos::OrdinalTraits<LO>::invalid ();
02159       }
02160       else if (values.size () != indices.size ()) {
02161         // The sizes of values and indices must match.
02162         return Teuchos::OrdinalTraits<LO>::invalid ();
02163       }
02164 
02165       const LO lrow = this->getRowMap()->getLocalElement (globalRow);
02166       if (lrow == Teuchos::OrdinalTraits<LO>::invalid ()) {
02167         // We don't own the row, so we're not allowed to modify its values.
02168         return Teuchos::OrdinalTraits<LO>::invalid ();
02169       }
02170 
02171       if (staticGraph_.is_null ()) {
02172         return Teuchos::OrdinalTraits<LO>::invalid ();
02173       }
02174       const crs_graph_type& graph = *staticGraph_;
02175       RowInfo rowInfo = graph.getRowInfo (lrow);
02176       if (indices.size () == 0) {
02177         return static_cast<LO> (0);
02178       }
02179       else {
02180         ArrayView<Scalar> curVals = this->getViewNonConst (rowInfo);
02181         if (isLocallyIndexed ()) {
02182           // Convert the given global indices to local indices.
02183           //
02184           // FIXME (mfh 08 Jul 2014) Why can't we ask the graph to do
02185           // that?  It could do the conversions in place, so that we
02186           // wouldn't need temporary storage.
02187           const map_type& colMap = * (this->getColMap ());
02188           const size_type numInds = indices.size ();
02189           Array<LO> lclInds (numInds);
02190           for (size_type k = 0; k < numInds; ++k) {
02191             // There is no need to filter out indices not in the
02192             // column Map.  Those that aren't will be mapped to
02193             // invalid(), which the graph's transformGlobalValues()
02194             // will filter out (but not count in its return value).
02195             lclInds[k] = colMap.getLocalElement (indices[k]);
02196           }
02197           return graph.template transformLocalValues<Scalar, BinaryFunction> (rowInfo, curVals,
02198                                                                               lclInds (), values, f);
02199         }
02200         else if (isGloballyIndexed ()) {
02201           return graph.template transformGlobalValues<Scalar, BinaryFunction> (rowInfo, curVals,
02202                                                                                indices, values, f);
02203         }
02204         else {
02205           // If the graph is neither locally nor globally indexed on
02206           // the calling process, that means that the calling process
02207           // can't possibly have any entries in the owned row.  Thus,
02208           // there are no entries to transform, so we return zero.
02209           return static_cast<LO> (0);
02210         }
02211       }
02212     }
02213 
02214   private:
02217     void
02218     insertNonownedGlobalValues (const GlobalOrdinal globalRow,
02219                                 const Teuchos::ArrayView<const GlobalOrdinal>& indices,
02220                                 const Teuchos::ArrayView<const Scalar>& values);
02221 
02222   protected:
02223     // useful typedefs
02224     typedef OrdinalTraits<LocalOrdinal>                     OTL;
02225     typedef ScalarTraits<Scalar>                            STS;
02226     typedef typename STS::magnitudeType               Magnitude;
02227     typedef ScalarTraits<Magnitude>                         STM;
02228     typedef MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type> MV;
02229     typedef Vector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>      V;
02230     typedef crs_graph_type Graph;
02231     typedef typename LocalMatOps::template bind_scalar<Scalar>::other_type                    sparse_ops_type;
02232     typedef typename sparse_ops_type::template graph<LocalOrdinal,node_type>::graph_type          local_graph_type;
02233     typedef typename sparse_ops_type::template matrix<Scalar,LocalOrdinal,node_type>::matrix_type local_matrix_type;
02234 
02235     // Enums
02236     enum GraphAllocationStatus {
02237       GraphAlreadyAllocated,
02238       GraphNotYetAllocated
02239     };
02240 
02257     void allocateValues (ELocalGlobal lg, GraphAllocationStatus gas);
02258 
02264     void sortEntries();
02265 
02271     void mergeRedundantEntries();
02272 
02280     void clearGlobalConstants();
02281 
02290     void computeGlobalConstants();
02291 
02304     mutable Teuchos::RCP<MV> importMV_;
02305 
02318     mutable Teuchos::RCP<MV> exportMV_;
02319 
02339     Teuchos::RCP<MV>
02340     getColumnMapMultiVector (const MV& X_domainMap,
02341                              const bool force = false) const;
02342 
02364     Teuchos::RCP<MV>
02365     getRowMapMultiVector (const MV& Y_rangeMap,
02366                           const bool force = false) const;
02367 
02369     void
02370     applyNonTranspose (const MV& X_in,
02371                        MV& Y_in,
02372                        Scalar alpha,
02373                        Scalar beta) const;
02374 
02376     void
02377     applyTranspose (const MV& X_in,
02378                     MV& Y_in,
02379                     const Teuchos::ETransp mode,
02380                     Scalar alpha,
02381                     Scalar beta) const;
02382 
02383     // matrix data accessors
02384     ArrayView<const Scalar>    getView(RowInfo rowinfo) const;
02385     ArrayView<      Scalar>    getViewNonConst(RowInfo rowinfo);
02386     // local Kokkos objects
02387 
02393     void fillLocalMatrix(const RCP<ParameterList> &params);
02399     void fillLocalGraphAndMatrix(const RCP<ParameterList> &params);
02400     // debugging
02401     void checkInternalState() const;
02402 
02414 
02415     RCP<const Graph> staticGraph_;
02416     RCP<      Graph>     myGraph_;
02418 
02423     RCP<sparse_ops_type>   lclMatOps_;
02430     RCP<local_matrix_type> lclMatrix_;
02431 
02444 
02445     ArrayRCP<Scalar> values1D_;
02446     ArrayRCP<Array<Scalar> > values2D_;
02448 
02450     bool fillComplete_;
02451 
02479     std::map<GlobalOrdinal, std::pair<Teuchos::Array<GlobalOrdinal>,
02480                                       Teuchos::Array<Scalar> > > nonlocals_;
02481 
02488     mutable Magnitude frobNorm_;
02489   }; // class CrsMatrix
02490 
02497   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
02498   Teuchos::RCP<CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
02499   createCrsMatrix (const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &map,
02500                    size_t maxNumEntriesPerRow = 0,
02501                    const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null)
02502   {
02503     using Teuchos::rcp;
02504     typedef CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> matrix_type;
02505     return rcp (new matrix_type (map, maxNumEntriesPerRow, DynamicProfile, params));
02506   }
02507 
02559   template<class CrsMatrixType>
02560   Teuchos::RCP<CrsMatrixType>
02561   importAndFillCompleteCrsMatrix (const Teuchos::RCP<const CrsMatrixType>& sourceMatrix,
02562                                   const Import<typename CrsMatrixType::local_ordinal_type,
02563                                                typename CrsMatrixType::global_ordinal_type,
02564                                                typename CrsMatrixType::node_type>& importer,
02565                                   const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
02566                                                                typename CrsMatrixType::global_ordinal_type,
02567                                                                typename CrsMatrixType::node_type> >& domainMap = Teuchos::null,
02568                                   const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
02569                                                                typename CrsMatrixType::global_ordinal_type,
02570                                                                typename CrsMatrixType::node_type> >& rangeMap = Teuchos::null,
02571                                   const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null)
02572   {
02573     Teuchos::RCP<CrsMatrixType> destMatrix;
02574     sourceMatrix->importAndFillComplete (destMatrix,importer, domainMap, rangeMap, params);
02575     return destMatrix;
02576   }
02577 
02611   template<class CrsMatrixType>
02612   Teuchos::RCP<CrsMatrixType>
02613   exportAndFillCompleteCrsMatrix (const Teuchos::RCP<const CrsMatrixType>& sourceMatrix,
02614                                   const Export<typename CrsMatrixType::local_ordinal_type,
02615                                                typename CrsMatrixType::global_ordinal_type,
02616                                                typename CrsMatrixType::node_type>& exporter,
02617                                   const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
02618                                                                typename CrsMatrixType::global_ordinal_type,
02619                                                                typename CrsMatrixType::node_type> >& domainMap = Teuchos::null,
02620                                   const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
02621                                                                typename CrsMatrixType::global_ordinal_type,
02622                                                                typename CrsMatrixType::node_type> >& rangeMap = Teuchos::null,
02623                                   const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null)
02624   {
02625     Teuchos::RCP<CrsMatrixType> destMatrix;
02626     sourceMatrix->exportAndFillComplete (destMatrix,exporter, domainMap, rangeMap, params);
02627     return destMatrix;
02628   }
02629 } // namespace Tpetra
02630 
02631 // Include KokkosRefactor partial specialisation if enabled
02632 #if defined(TPETRA_HAVE_KOKKOS_REFACTOR)
02633 #include "Tpetra_KokkosRefactor_CrsMatrix_decl.hpp"
02634 #endif
02635 
02646 #endif // TPETRA_CRSMATRIX_DECL_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines