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 
00176   template <class Scalar        = double,
00177             class LocalOrdinal  = int,
00178             class GlobalOrdinal = LocalOrdinal,
00179             class Node          = KokkosClassic::DefaultNode::DefaultNodeType,
00180             class LocalMatOps   = typename KokkosClassic::DefaultKernels<Scalar,LocalOrdinal,Node>::SparseOps >
00181   class CrsMatrix : public RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>,
00182                     public DistObject<char, LocalOrdinal,GlobalOrdinal,Node> {
00183   public:
00185 
00186 
00188     typedef Scalar scalar_type;
00190     typedef LocalOrdinal local_ordinal_type;
00192     typedef GlobalOrdinal global_ordinal_type;
00194     typedef Node node_type;
00195 
00199     typedef LocalMatOps mat_vec_type;
00203     typedef LocalMatOps mat_solve_type;
00204 
00206     typedef Map<LocalOrdinal, GlobalOrdinal, node_type> map_type;
00207 
00209     typedef Import<LocalOrdinal, GlobalOrdinal, node_type> import_type;
00210 
00212     typedef Export<LocalOrdinal, GlobalOrdinal, node_type> export_type;
00213 
00215     typedef CrsGraph<LocalOrdinal, GlobalOrdinal, node_type, LocalMatOps> crs_graph_type;
00216 
00218 
00219 
00220 
00238     CrsMatrix (const Teuchos::RCP<const map_type>& rowMap,
00239                size_t maxNumEntriesPerRow,
00240                ProfileType pftype = DynamicProfile,
00241                const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null);
00242 
00260     CrsMatrix (const Teuchos::RCP<const map_type>& rowMap,
00261                const ArrayRCP<const size_t>& NumEntriesPerRowToAlloc,
00262                ProfileType pftype = DynamicProfile,
00263                const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null);
00264 
00287     CrsMatrix (const Teuchos::RCP<const map_type>& rowMap,
00288                const Teuchos::RCP<const map_type>& colMap,
00289                size_t maxNumEntriesPerRow,
00290                ProfileType pftype = DynamicProfile,
00291                const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null);
00292 
00315     CrsMatrix (const Teuchos::RCP<const map_type>& rowMap,
00316                const Teuchos::RCP<const map_type>& colMap,
00317                const ArrayRCP<const size_t>& NumEntriesPerRowToAlloc,
00318                ProfileType pftype = DynamicProfile,
00319                const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null);
00320 
00345     explicit CrsMatrix (const Teuchos::RCP<const crs_graph_type>& graph,
00346                         const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null);
00347 
00370     CrsMatrix (const RCP<const map_type>& rowMap,
00371                const RCP<const map_type>& colMap,
00372                const ArrayRCP<size_t>& rowPointers,
00373                const ArrayRCP<LocalOrdinal>& columnIndices,
00374                const ArrayRCP<Scalar>& values,
00375                const RCP<ParameterList>& params = null);
00376 
00377     // This friend declaration makes the clone() method work.
00378     template <class S2, class LO2, class GO2, class N2, class LMO2>
00379     friend class CrsMatrix;
00380 
00405     template <class Node2>
00406     RCP<CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node2, typename KokkosClassic::DefaultKernels<void,LocalOrdinal,Node2>::SparseOps> >
00407     clone (const RCP<Node2>& node2, const RCP<ParameterList>& params = null) const
00408     {
00409       const char tfecfFuncName[] = "clone";
00410 
00411       // Get parameter values.  Set them initially to their default values.
00412       bool fillCompleteClone = true;
00413       bool useLocalIndices = this->hasColMap ();
00414       ProfileType pftype = StaticProfile;
00415       if (! params.is_null ()) {
00416         fillCompleteClone = params->get ("fillComplete clone", fillCompleteClone);
00417         useLocalIndices = params->get ("Locally indexed clone", useLocalIndices);
00418 
00419         bool staticProfileClone = true;
00420         staticProfileClone = params->get ("Static profile clone", staticProfileClone);
00421         pftype = staticProfileClone ? StaticProfile : DynamicProfile;
00422       }
00423 
00424       TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00425         ! this->hasColMap () && useLocalIndices, std::runtime_error,
00426         ": You requested that the returned clone have local indices, but the "
00427         "the source matrix does not have a column Map yet.");
00428 
00429       typedef typename KokkosClassic::DefaultKernels<void,LocalOrdinal,Node2>::SparseOps LocalMatOps2;
00430       typedef CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node2, LocalMatOps2> CrsMatrix2;
00431       typedef Map<LocalOrdinal, GlobalOrdinal, Node2> Map2;
00432       RCP<const Map2> clonedRowMap = this->getRowMap ()->template clone<Node2> (node2);
00433 
00434       RCP<CrsMatrix2> clonedMatrix;
00435       ArrayRCP<const size_t> numEntries;
00436       size_t numEntriesForAll = 0;
00437       if (! staticGraph_->indicesAreAllocated ()) {
00438         if (! staticGraph_->numAllocPerRow_.is_null ()) {
00439           numEntries = staticGraph_->numAllocPerRow_;
00440         }
00441         else {
00442           numEntriesForAll = staticGraph_->numAllocForAllRows_;
00443         }
00444       }
00445       else if (! staticGraph_->numRowEntries_.is_null ()) {
00446         numEntries = staticGraph_->numRowEntries_;
00447       }
00448       else if (staticGraph_->nodeNumAllocated_ == 0) {
00449         numEntriesForAll = 0;
00450       }
00451       else {
00452         // We're left with the case that we have optimized storage.
00453         // In this case, we have to construct a list of row sizes.
00454         TEUCHOS_TEST_FOR_EXCEPTION(
00455           getProfileType() != StaticProfile, std::logic_error,
00456           "Internal logic error. Please report this to Tpetra team." )
00457 
00458         const size_t numRows = this->getNodeNumRows ();
00459         numEntriesForAll = 0;
00460         ArrayRCP<size_t> numEnt;
00461         if (numRows != 0) {
00462           numEnt = arcp<size_t> (numRows);
00463         }
00464         for (size_t i = 0; i < numRows; ++i) {
00465           numEnt[i] = staticGraph_->rowPtrs_[i+1] - staticGraph_->rowPtrs_[i];
00466         }
00467         numEntries = numEnt;
00468       }
00469 
00470       RCP<ParameterList> matrixparams =
00471         params.is_null () ? null : sublist (params,"CrsMatrix");
00472       if (useLocalIndices) {
00473         RCP<const Map2> clonedColMap =
00474           this->getColMap ()->template clone<Node2> (node2);
00475         if (numEntries.is_null ()) {
00476           clonedMatrix = rcp (new CrsMatrix2 (clonedRowMap, clonedColMap,
00477                                               numEntriesForAll, pftype,
00478                                               matrixparams));
00479         }
00480         else {
00481           clonedMatrix = rcp (new CrsMatrix2 (clonedRowMap, clonedColMap,
00482                                               numEntries, pftype,
00483                                               matrixparams));
00484         }
00485       }
00486       else {
00487         if (numEntries.is_null ()) {
00488           clonedMatrix = rcp (new CrsMatrix2 (clonedRowMap, numEntriesForAll,
00489                                               pftype, matrixparams));
00490         }
00491         else {
00492           clonedMatrix = rcp (new CrsMatrix2 (clonedRowMap, numEntries, pftype,
00493                                               matrixparams));
00494         }
00495       }
00496       // done with these
00497       numEntries = null;
00498       numEntriesForAll = 0;
00499 
00500       if (useLocalIndices) {
00501         clonedMatrix->allocateValues (LocalIndices,
00502                                       CrsMatrix2::GraphNotYetAllocated);
00503         if (this->isLocallyIndexed ()) {
00504           ArrayView<const LocalOrdinal> linds;
00505           ArrayView<const Scalar>       vals;
00506           for (LocalOrdinal lrow = clonedRowMap->getMinLocalIndex ();
00507                lrow <= clonedRowMap->getMaxLocalIndex ();
00508                ++lrow) {
00509             this->getLocalRowView (lrow, linds, vals);
00510             if (linds.size ()) {
00511               clonedMatrix->insertLocalValues (lrow, linds, vals);
00512             }
00513           }
00514         }
00515         else { // this->isGloballyIndexed()
00516           Array<LocalOrdinal> linds;
00517           Array<Scalar> vals;
00518           for (LocalOrdinal lrow = clonedRowMap->getMinLocalIndex ();
00519                lrow <= clonedRowMap->getMaxLocalIndex ();
00520                ++lrow) {
00521             size_t theNumEntries = this->getNumEntriesInLocalRow (lrow);
00522             if (theNumEntries > Teuchos::as<size_t> (linds.size ())) {
00523               linds.resize (theNumEntries);
00524             }
00525             if (theNumEntries > Teuchos::as<size_t> (vals.size ())) {
00526               vals.resize (theNumEntries);
00527             }
00528             this->getLocalRowCopy (clonedRowMap->getGlobalElement (lrow),
00529                                    linds (), vals (), theNumEntries);
00530             if (theNumEntries != 0) {
00531               clonedMatrix->insertLocalValues (lrow, linds (0, theNumEntries),
00532                                                vals (0, theNumEntries));
00533             }
00534           }
00535         }
00536       }
00537       else { // useGlobalIndices
00538         clonedMatrix->allocateValues (GlobalIndices,
00539                                       CrsMatrix2::GraphNotYetAllocated);
00540         if (this->isGloballyIndexed ()) {
00541           ArrayView<const GlobalOrdinal> ginds;
00542           ArrayView<const Scalar>         vals;
00543           for (GlobalOrdinal grow = clonedRowMap->getMinGlobalIndex ();
00544                grow <= clonedRowMap->getMaxGlobalIndex ();
00545                ++grow) {
00546             this->getGlobalRowView (grow, ginds, vals);
00547             if (ginds.size () > 0) {
00548               clonedMatrix->insertGlobalValues (grow, ginds, vals);
00549             }
00550           }
00551         }
00552         else { // this->isLocallyIndexed()
00553           Array<GlobalOrdinal> ginds;
00554           Array<Scalar> vals;
00555           for (GlobalOrdinal grow = clonedRowMap->getMinGlobalIndex ();
00556                grow <= clonedRowMap->getMaxGlobalIndex ();
00557                ++grow) {
00558             size_t theNumEntries = this->getNumEntriesInGlobalRow (grow);
00559             if (theNumEntries > Teuchos::as<size_t> (ginds.size ())) {
00560               ginds.resize (theNumEntries);
00561             }
00562             if (theNumEntries > Teuchos::as<size_t> (vals.size ())) {
00563               vals.resize (theNumEntries);
00564             }
00565             this->getGlobalRowCopy (grow, ginds (), vals (), theNumEntries);
00566             if (theNumEntries != 0) {
00567               clonedMatrix->insertGlobalValues (grow, ginds (0, theNumEntries),
00568                                                 vals (0, theNumEntries));
00569             }
00570           }
00571         }
00572       }
00573 
00574       if (fillCompleteClone) {
00575         RCP<ParameterList> fillparams =
00576           params.is_null () ? Teuchos::null : sublist (params, "fillComplete");
00577         try {
00578           RCP<const Map2> clonedRangeMap;
00579           RCP<const Map2> clonedDomainMap;
00580           if (! this->getRangeMap ().is_null () &&
00581               this->getRangeMap () != clonedRowMap) {
00582             clonedRangeMap  = this->getRangeMap ()->template clone<Node2> (node2);
00583           }
00584           else {
00585             clonedRangeMap = clonedRowMap;
00586           }
00587           if (! this->getDomainMap ().is_null () &&
00588               this->getDomainMap () != clonedRowMap) {
00589             clonedDomainMap = this->getDomainMap ()->template clone<Node2> (node2);
00590           }
00591           else {
00592             clonedDomainMap = clonedRowMap;
00593           }
00594           clonedMatrix->fillComplete (clonedDomainMap, clonedRangeMap,
00595                                       fillparams);
00596         }
00597         catch (std::exception &e) {
00598           const bool caughtExceptionOnClone = true;
00599           TEUCHOS_TEST_FOR_EXCEPTION(
00600             caughtExceptionOnClone, std::runtime_error,
00601             Teuchos::typeName (*this) << std::endl << "clone: " << std::endl <<
00602             "Caught the following exception while calling fillComplete() on a "
00603             "clone of type" << std::endl << Teuchos::typeName (*clonedMatrix)
00604             << ": " << std::endl << e.what () << std::endl);
00605         }
00606       }
00607       return clonedMatrix;
00608     }
00609 
00611     virtual ~CrsMatrix ();
00612 
00614 
00615 
00616 
00639     //
00685     void
00686     insertGlobalValues (const GlobalOrdinal globalRow,
00687                         const ArrayView<const GlobalOrdinal>& cols,
00688                         const ArrayView<const Scalar>& vals);
00689 
00729     void
00730     insertLocalValues (const LocalOrdinal localRow,
00731                        const ArrayView<const LocalOrdinal> &cols,
00732                        const ArrayView<const Scalar> &vals);
00733 
00768     LocalOrdinal
00769     replaceGlobalValues (GlobalOrdinal globalRow,
00770                          const ArrayView<const GlobalOrdinal>& cols,
00771                          const ArrayView<const Scalar>& vals);
00772 
00803     LocalOrdinal
00804     replaceLocalValues (const LocalOrdinal localRow,
00805                         const ArrayView<const LocalOrdinal>&cols,
00806                         const ArrayView<const Scalar>& vals);
00807 
00840     LocalOrdinal
00841     sumIntoGlobalValues (const GlobalOrdinal globalRow,
00842                          const ArrayView<const GlobalOrdinal> &cols,
00843                          const ArrayView<const Scalar>        &vals);
00844 
00859     LocalOrdinal
00860     sumIntoLocalValues (const LocalOrdinal localRow,
00861                         const ArrayView<const LocalOrdinal>& cols,
00862                         const ArrayView<const Scalar>& vals);
00863 
00865     void setAllToScalar (const Scalar &alpha);
00866 
00868     void scale (const Scalar &alpha);
00869 
00871 
00878     void
00879     setAllValues (const ArrayRCP<size_t>& rowPointers,
00880                   const ArrayRCP<LocalOrdinal>& columnIndices,
00881                   const ArrayRCP<Scalar>& values);
00882 
00884 
00891     void
00892     getAllValues (ArrayRCP<const size_t>& rowPointers,
00893                   ArrayRCP<const LocalOrdinal>& columnIndices,
00894                   ArrayRCP<const Scalar>& values) const;
00895 
00897 
00898 
00899 
00928     void globalAssemble();
00929 
00943     void resumeFill (const RCP<ParameterList>& params = null);
00944 
00978     void
00979     fillComplete (const RCP<const map_type>& domainMap,
00980                   const RCP<const map_type>& rangeMap,
00981                   const RCP<ParameterList>& params = null);
00982 
00995     void fillComplete (const RCP<ParameterList>& params = null);
00996 
01008     void
01009     expertStaticFillComplete (const RCP<const map_type>& domainMap,
01010                               const RCP<const map_type>& rangeMap,
01011                               const RCP<const import_type>& importer = Teuchos::null,
01012                               const RCP<const export_type>& exporter = Teuchos::null,
01013                               const RCP<ParameterList>& params = Teuchos::null);
01014 
01024     void
01025     replaceColMap (const Teuchos::RCP<const map_type>& newColMap);
01026 
01039     void
01040     replaceDomainMapAndImporter (const Teuchos::RCP<const map_type>& newDomainMap,
01041                                  Teuchos::RCP<const import_type>& newImporter);
01042 
01056     virtual void
01057     removeEmptyProcessesInPlace (const Teuchos::RCP<const map_type>& newMap);
01058 
01060 
01061 
01062 
01064     RCP<const Comm<int> > getComm() const;
01065 
01067     RCP<node_type> getNode () const;
01068 
01070     RCP<const map_type> getRowMap () const;
01071 
01073     RCP<const map_type> getColMap () const;
01074 
01076     RCP<const RowGraph<LocalOrdinal,GlobalOrdinal,node_type> > getGraph () const;
01077 
01079     RCP<const crs_graph_type> getCrsGraph () const;
01080 
01100     global_size_t getGlobalNumRows() const;
01101 
01107     global_size_t getGlobalNumCols() const;
01108 
01115     size_t getNodeNumRows() const;
01116 
01120     size_t getNodeNumCols() const;
01121 
01123     GlobalOrdinal getIndexBase() const;
01124 
01126     global_size_t getGlobalNumEntries() const;
01127 
01129     size_t getNodeNumEntries() const;
01130 
01132 
01133     size_t getNumEntriesInGlobalRow(GlobalOrdinal globalRow) const;
01134 
01136 
01137     size_t getNumEntriesInLocalRow(LocalOrdinal localRow) const;
01138 
01140 
01142     global_size_t getGlobalNumDiags() const;
01143 
01145 
01147     size_t getNodeNumDiags() const;
01148 
01150 
01152     size_t getGlobalMaxNumRowEntries() const;
01153 
01155 
01157     size_t getNodeMaxNumRowEntries() const;
01158 
01160     bool hasColMap() const;
01161 
01163 
01165     bool isLowerTriangular() const;
01166 
01168 
01170     bool isUpperTriangular() const;
01171 
01191     bool isLocallyIndexed() const;
01192 
01212     bool isGloballyIndexed() const;
01213 
01236     bool isFillComplete() const;
01237 
01260     bool isFillActive() const;
01261 
01263 
01269     bool isStorageOptimized() const;
01270 
01272     ProfileType getProfileType() const;
01273 
01275     bool isStaticGraph() const;
01276 
01288     typename ScalarTraits<Scalar>::magnitudeType getFrobeniusNorm() const;
01289 
01296     virtual bool supportsRowViews() const;
01297 
01299 
01309     void
01310     getGlobalRowCopy (GlobalOrdinal GlobalRow,
01311                       const ArrayView<GlobalOrdinal> &Indices,
01312                       const ArrayView<Scalar> &Values,
01313                       size_t &NumEntries) const;
01314 
01316 
01328     void
01329     getLocalRowCopy (LocalOrdinal LocalRow,
01330                      const ArrayView<LocalOrdinal> &Indices,
01331                      const ArrayView<Scalar> &Values,
01332                      size_t &NumEntries) const;
01333 
01347     void
01348     getGlobalRowView (GlobalOrdinal GlobalRow,
01349                       ArrayView<const GlobalOrdinal> &indices,
01350                       ArrayView<const Scalar> &values) const;
01351 
01365     void
01366     getLocalRowView (LocalOrdinal LocalRow,
01367                      ArrayView<const LocalOrdinal>& indices,
01368                      ArrayView<const Scalar>& values) const;
01369 
01375     void getLocalDiagCopy (Vector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>& diag) const;
01376 
01414     void getLocalDiagOffsets (Teuchos::ArrayRCP<size_t>& offsets) const;
01415 
01434     void
01435     getLocalDiagCopy (Vector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>& diag,
01436                       const Teuchos::ArrayView<const size_t>& offsets) const;
01437 
01439     void leftScale(const Vector<Scalar, LocalOrdinal, GlobalOrdinal, node_type>& x);
01440 
01442     void rightScale(const Vector<Scalar, LocalOrdinal, GlobalOrdinal, node_type>& x);
01443 
01445 
01446 
01447 
01496     template <class DomainScalar, class RangeScalar>
01497     void
01498     localMultiply (const MultiVector<DomainScalar,LocalOrdinal,GlobalOrdinal,node_type>& X,
01499                    MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,node_type>& Y,
01500                    Teuchos::ETransp trans,
01501                    RangeScalar alpha,
01502                    RangeScalar beta) const;
01503 
01528     template <class DomainScalar, class RangeScalar>
01529     void
01530     localGaussSeidel (const MultiVector<DomainScalar,LocalOrdinal,GlobalOrdinal,node_type> &B,
01531                       MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,node_type> &X,
01532                       const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type> &D,
01533                       const RangeScalar& dampingFactor,
01534                       const KokkosClassic::ESweepDirection direction) const;
01535 
01562     template <class DomainScalar, class RangeScalar>
01563     void
01564     reorderedLocalGaussSeidel (const MultiVector<DomainScalar,LocalOrdinal,GlobalOrdinal,node_type>& B,
01565                                MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,node_type>& X,
01566                                const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>& D,
01567                                const ArrayView<LocalOrdinal>& rowIndices,
01568                                const RangeScalar& dampingFactor,
01569                                const KokkosClassic::ESweepDirection direction) const;
01570 
01587     template <class DomainScalar, class RangeScalar>
01588     void
01589     localSolve (const MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,node_type>& Y,
01590                 MultiVector<DomainScalar,LocalOrdinal,GlobalOrdinal,node_type>& X,
01591                 Teuchos::ETransp trans) const;
01592 
01594     template <class T>
01595     RCP<CrsMatrix<T,LocalOrdinal,GlobalOrdinal,node_type> > convert() const;
01596 
01598 
01599 
01600 
01611     void
01612     apply (const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>& X,
01613            MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>&Y,
01614            Teuchos::ETransp mode = Teuchos::NO_TRANS,
01615            Scalar alpha = ScalarTraits<Scalar>::one(),
01616            Scalar beta = ScalarTraits<Scalar>::zero()) const;
01617 
01619     bool hasTransposeApply() const;
01620 
01624     RCP<const map_type> getDomainMap () const;
01625 
01629     RCP<const map_type> getRangeMap () const;
01630 
01632 
01633 
01634 
01699     void
01700     gaussSeidel (const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type> &B,
01701                  MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type> &X,
01702                  const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type> &D,
01703                  const Scalar& dampingFactor,
01704                  const ESweepDirection direction,
01705                  const int numSweeps) const;
01706 
01773     void
01774     reorderedGaussSeidel (const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>& B,
01775                           MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>& X,
01776                           const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>& D,
01777                           const ArrayView<LocalOrdinal>& rowIndices,
01778                           const Scalar& dampingFactor,
01779                           const ESweepDirection direction,
01780                           const int numSweeps) const;
01781 
01810     void
01811     gaussSeidelCopy (MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type> &X,
01812                      const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type> &B,
01813                      const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type> &D,
01814                      const Scalar& dampingFactor,
01815                      const ESweepDirection direction,
01816                      const int numSweeps,
01817                      const bool zeroInitialGuess) const;
01818 
01848     void
01849     reorderedGaussSeidelCopy (MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>& X,
01850                               const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>& B,
01851                               const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>& D,
01852                               const ArrayView<LocalOrdinal>& rowIndices,
01853                               const Scalar& dampingFactor,
01854                               const ESweepDirection direction,
01855                               const int numSweeps,
01856                               const bool zeroInitialGuess) const;
01857 
01868     virtual Teuchos::RCP<RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, node_type> >
01869     add (const Scalar& alpha,
01870          const RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, node_type>& A,
01871          const Scalar& beta,
01872          const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, node_type> >& domainMap,
01873          const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, node_type> >& rangeMap,
01874          const Teuchos::RCP<Teuchos::ParameterList>& params) const;
01875 
01877 
01878 
01879 
01881     std::string description() const;
01882 
01884     void
01885     describe (Teuchos::FancyOStream &out,
01886               const Teuchos::EVerbosityLevel verbLevel =
01887               Teuchos::Describable::verbLevel_default) const;
01888 
01890 
01891 
01892 
01893     virtual bool
01894     checkSizes (const SrcDistObject& source);
01895 
01896     virtual void
01897     copyAndPermute (const SrcDistObject& source,
01898                     size_t numSameIDs,
01899                     const ArrayView<const LocalOrdinal> &permuteToLIDs,
01900                     const ArrayView<const LocalOrdinal> &permuteFromLIDs);
01901 
01902     virtual void
01903     packAndPrepare (const SrcDistObject& source,
01904                     const Teuchos::ArrayView<const LocalOrdinal>& exportLIDs,
01905                     Teuchos::Array<char>& exports,
01906                     const Teuchos::ArrayView<size_t>& numPacketsPerLID,
01907                     size_t& constantNumPackets,
01908                     Distributor& distor);
01909 
01919     void
01920     unpackAndCombine (const Teuchos::ArrayView<const LocalOrdinal> &importLIDs,
01921                       const Teuchos::ArrayView<const char> &imports,
01922                       const Teuchos::ArrayView<size_t> &numPacketsPerLID,
01923                       size_t constantNumPackets,
01924                       Distributor &distor,
01925                       CombineMode combineMode);
01927 
01928 
01929 
01938     virtual void
01939     pack (const Teuchos::ArrayView<const LocalOrdinal>& exportLIDs,
01940           Teuchos::Array<char>& exports,
01941           const Teuchos::ArrayView<size_t>& numPacketsPerLID,
01942           size_t& constantNumPackets,
01943           Distributor& distor) const;
01944 
01946   private:
01947     // Friend declaration for nonmember function.
01948     template<class CrsMatrixType>
01949     friend Teuchos::RCP<CrsMatrixType>
01950     importAndFillCompleteCrsMatrix (const Teuchos::RCP<const CrsMatrixType>& sourceMatrix,
01951                                     const Import<typename CrsMatrixType::local_ordinal_type,
01952                                                  typename CrsMatrixType::global_ordinal_type,
01953                                                  typename CrsMatrixType::node_type>& importer,
01954                                     const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
01955                                                                  typename CrsMatrixType::global_ordinal_type,
01956                                                                  typename CrsMatrixType::node_type> >& domainMap,
01957                                     const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
01958                                                                  typename CrsMatrixType::global_ordinal_type,
01959                                                                  typename CrsMatrixType::node_type> >& rangeMap,
01960                                     const Teuchos::RCP<Teuchos::ParameterList>& params);
01961 
01962     // Friend declaration for nonmember function.
01963     template<class CrsMatrixType>
01964     friend Teuchos::RCP<CrsMatrixType>
01965     exportAndFillCompleteCrsMatrix (const Teuchos::RCP<const CrsMatrixType>& sourceMatrix,
01966                                     const Export<typename CrsMatrixType::local_ordinal_type,
01967                                                  typename CrsMatrixType::global_ordinal_type,
01968                                                  typename CrsMatrixType::node_type>& exporter,
01969                                     const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
01970                                                                  typename CrsMatrixType::global_ordinal_type,
01971                                                                  typename CrsMatrixType::node_type> >& domainMap,
01972                                     const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
01973                                                                  typename CrsMatrixType::global_ordinal_type,
01974                                                                  typename CrsMatrixType::node_type> >& rangeMap,
01975                                     const Teuchos::RCP<Teuchos::ParameterList>& params);
01976 
01977   public:
01993     void
01994     importAndFillComplete (Teuchos::RCP<CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,node_type,LocalMatOps> > & destMatrix,
01995                            const import_type& importer,
01996                            const Teuchos::RCP<const map_type>& domainMap = Teuchos::null,
01997                            const Teuchos::RCP<const map_type>& rangeMap = Teuchos::null,
01998                            const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null) const;
01999 
02015     void
02016     exportAndFillComplete (Teuchos::RCP<CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,node_type,LocalMatOps> > & destMatrix,
02017                            const export_type& exporter,
02018                            const Teuchos::RCP<const map_type>& domainMap = Teuchos::null,
02019                            const Teuchos::RCP<const map_type>& rangeMap = Teuchos::null,
02020                            const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null) const;
02021   private:
02042     void
02043     transferAndFillComplete (Teuchos::RCP<CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps> > & destMat,
02044                              const ::Tpetra::Details::Transfer<LocalOrdinal, GlobalOrdinal, Node>& rowTransfer, // either Import or Export
02045                              const Teuchos::RCP<const map_type>& domainMap = Teuchos::null,
02046                              const Teuchos::RCP<const map_type>& rangeMap = Teuchos::null,
02047                              const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null) const;
02048 
02049     // We forbid copy construction by declaring this method private
02050     // and not implementing it.
02051     CrsMatrix (const CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,node_type,LocalMatOps> &rhs);
02052 
02053     // We forbid assignment (operator=) by declaring this method
02054     // private and not implementing it.
02055     CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,node_type,LocalMatOps>&
02056     operator= (const CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,node_type,LocalMatOps> &rhs);
02057 
02063     void
02064     insertGlobalValuesFiltered (const GlobalOrdinal globalRow,
02065                                 const ArrayView<const GlobalOrdinal> &indices,
02066                                 const ArrayView<const Scalar>        &values);
02067 
02073     void
02074     insertLocalValuesFiltered (const LocalOrdinal localRow,
02075                                const ArrayView<const LocalOrdinal> &indices,
02076                                const ArrayView<const Scalar>       &values);
02077 
02085     void
02086     combineGlobalValues (const GlobalOrdinal globalRowIndex,
02087                          const Teuchos::ArrayView<const GlobalOrdinal> columnIndices,
02088                          const Teuchos::ArrayView<const Scalar> values,
02089                          const Tpetra::CombineMode combineMode);
02090 
02120     template<class BinaryFunction>
02121     LocalOrdinal
02122     transformGlobalValues (GlobalOrdinal globalRow,
02123                            const Teuchos::ArrayView<const GlobalOrdinal>& indices,
02124                            const Teuchos::ArrayView<const Scalar>        & values,
02125                            BinaryFunction f)
02126     {
02127       typedef LocalOrdinal LO;
02128       typedef GlobalOrdinal GO;
02129       using Teuchos::Array;
02130       using Teuchos::ArrayView;
02131       typedef typename ArrayView<const GO>::size_type size_type;
02132 
02133       if (! isFillActive ()) {
02134         // Fill must be active in order to call this method.
02135         return Teuchos::OrdinalTraits<LO>::invalid ();
02136       }
02137       else if (values.size () != indices.size ()) {
02138         // The sizes of values and indices must match.
02139         return Teuchos::OrdinalTraits<LO>::invalid ();
02140       }
02141 
02142       const LO lrow = this->getRowMap()->getLocalElement (globalRow);
02143       if (lrow == Teuchos::OrdinalTraits<LO>::invalid ()) {
02144         // We don't own the row, so we're not allowed to modify its values.
02145         return Teuchos::OrdinalTraits<LO>::invalid ();
02146       }
02147 
02148       if (staticGraph_.is_null ()) {
02149         return Teuchos::OrdinalTraits<LO>::invalid ();
02150       }
02151       const crs_graph_type& graph = *staticGraph_;
02152       RowInfo rowInfo = graph.getRowInfo (lrow);
02153       if (indices.size () == 0) {
02154         return static_cast<LO> (0);
02155       }
02156       else {
02157         ArrayView<Scalar> curVals = this->getViewNonConst (rowInfo);
02158         if (isLocallyIndexed ()) {
02159           // Convert the given global indices to local indices.
02160           //
02161           // FIXME (mfh 08 Jul 2014) Why can't we ask the graph to do
02162           // that?  It could do the conversions in place, so that we
02163           // wouldn't need temporary storage.
02164           const map_type& colMap = * (this->getColMap ());
02165           const size_type numInds = indices.size ();
02166           Array<LO> lclInds (numInds);
02167           for (size_type k = 0; k < numInds; ++k) {
02168             // There is no need to filter out indices not in the
02169             // column Map.  Those that aren't will be mapped to
02170             // invalid(), which the graph's transformGlobalValues()
02171             // will filter out (but not count in its return value).
02172             lclInds[k] = colMap.getLocalElement (indices[k]);
02173           }
02174           return graph.template transformLocalValues<Scalar, BinaryFunction> (rowInfo, curVals,
02175                                                                               lclInds (), values, f);
02176         }
02177         else if (isGloballyIndexed ()) {
02178           return graph.template transformGlobalValues<Scalar, BinaryFunction> (rowInfo, curVals,
02179                                                                                indices, values, f);
02180         }
02181         else {
02182           // If the graph is neither locally nor globally indexed on
02183           // the calling process, that means that the calling process
02184           // can't possibly have any entries in the owned row.  Thus,
02185           // there are no entries to transform, so we return zero.
02186           return static_cast<LO> (0);
02187         }
02188       }
02189     }
02190 
02191   private:
02194     void
02195     insertNonownedGlobalValues (const GlobalOrdinal globalRow,
02196                                 const Teuchos::ArrayView<const GlobalOrdinal>& indices,
02197                                 const Teuchos::ArrayView<const Scalar>& values);
02198 
02199   protected:
02200     // useful typedefs
02201     typedef OrdinalTraits<LocalOrdinal>                     OTL;
02202     typedef ScalarTraits<Scalar>                            STS;
02203     typedef typename STS::magnitudeType               Magnitude;
02204     typedef ScalarTraits<Magnitude>                         STM;
02205     typedef MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type> MV;
02206     typedef Vector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>      V;
02207     typedef crs_graph_type Graph;
02208     typedef typename LocalMatOps::template bind_scalar<Scalar>::other_type                    sparse_ops_type;
02209     typedef typename sparse_ops_type::template graph<LocalOrdinal,node_type>::graph_type          local_graph_type;
02210     typedef typename sparse_ops_type::template matrix<Scalar,LocalOrdinal,node_type>::matrix_type local_matrix_type;
02211 
02212     // Enums
02213     enum GraphAllocationStatus {
02214       GraphAlreadyAllocated,
02215       GraphNotYetAllocated
02216     };
02217 
02234     void allocateValues (ELocalGlobal lg, GraphAllocationStatus gas);
02235 
02241     void sortEntries();
02242 
02248     void mergeRedundantEntries();
02249 
02257     void clearGlobalConstants();
02258 
02267     void computeGlobalConstants();
02268 
02281     mutable Teuchos::RCP<MV> importMV_;
02282 
02295     mutable Teuchos::RCP<MV> exportMV_;
02296 
02316     Teuchos::RCP<MV>
02317     getColumnMapMultiVector (const MV& X_domainMap,
02318                              const bool force = false) const;
02319 
02341     Teuchos::RCP<MV>
02342     getRowMapMultiVector (const MV& Y_rangeMap,
02343                           const bool force = false) const;
02344 
02346     void
02347     applyNonTranspose (const MV& X_in,
02348                        MV& Y_in,
02349                        Scalar alpha,
02350                        Scalar beta) const;
02351 
02353     void
02354     applyTranspose (const MV& X_in,
02355                     MV& Y_in,
02356                     const Teuchos::ETransp mode,
02357                     Scalar alpha,
02358                     Scalar beta) const;
02359 
02360     // matrix data accessors
02361     ArrayView<const Scalar>    getView(RowInfo rowinfo) const;
02362     ArrayView<      Scalar>    getViewNonConst(RowInfo rowinfo);
02363     // local Kokkos objects
02364 
02370     void fillLocalMatrix(const RCP<ParameterList> &params);
02376     void fillLocalGraphAndMatrix(const RCP<ParameterList> &params);
02377     // debugging
02378     void checkInternalState() const;
02379 
02391 
02392     RCP<const Graph> staticGraph_;
02393     RCP<      Graph>     myGraph_;
02395 
02400     RCP<sparse_ops_type>   lclMatOps_;
02407     RCP<local_matrix_type> lclMatrix_;
02408 
02421 
02422     ArrayRCP<Scalar> values1D_;
02423     ArrayRCP<Array<Scalar> > values2D_;
02425 
02427     bool fillComplete_;
02428 
02456     std::map<GlobalOrdinal, std::pair<Teuchos::Array<GlobalOrdinal>,
02457                                       Teuchos::Array<Scalar> > > nonlocals_;
02458 
02465     mutable Magnitude frobNorm_;
02466   }; // class CrsMatrix
02467 
02474   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
02475   Teuchos::RCP<CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
02476   createCrsMatrix (const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &map,
02477                    size_t maxNumEntriesPerRow = 0,
02478                    const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null)
02479   {
02480     using Teuchos::rcp;
02481     typedef CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> matrix_type;
02482     return rcp (new matrix_type (map, maxNumEntriesPerRow, DynamicProfile, params));
02483   }
02484 
02536   template<class CrsMatrixType>
02537   Teuchos::RCP<CrsMatrixType>
02538   importAndFillCompleteCrsMatrix (const Teuchos::RCP<const CrsMatrixType>& sourceMatrix,
02539                                   const Import<typename CrsMatrixType::local_ordinal_type,
02540                                                typename CrsMatrixType::global_ordinal_type,
02541                                                typename CrsMatrixType::node_type>& importer,
02542                                   const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
02543                                                                typename CrsMatrixType::global_ordinal_type,
02544                                                                typename CrsMatrixType::node_type> >& domainMap = Teuchos::null,
02545                                   const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
02546                                                                typename CrsMatrixType::global_ordinal_type,
02547                                                                typename CrsMatrixType::node_type> >& rangeMap = Teuchos::null,
02548                                   const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null)
02549   {
02550     Teuchos::RCP<CrsMatrixType> destMatrix;
02551     sourceMatrix->importAndFillComplete (destMatrix,importer, domainMap, rangeMap, params);
02552     return destMatrix;
02553   }
02554 
02588   template<class CrsMatrixType>
02589   Teuchos::RCP<CrsMatrixType>
02590   exportAndFillCompleteCrsMatrix (const Teuchos::RCP<const CrsMatrixType>& sourceMatrix,
02591                                   const Export<typename CrsMatrixType::local_ordinal_type,
02592                                                typename CrsMatrixType::global_ordinal_type,
02593                                                typename CrsMatrixType::node_type>& exporter,
02594                                   const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
02595                                                                typename CrsMatrixType::global_ordinal_type,
02596                                                                typename CrsMatrixType::node_type> >& domainMap = Teuchos::null,
02597                                   const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
02598                                                                typename CrsMatrixType::global_ordinal_type,
02599                                                                typename CrsMatrixType::node_type> >& rangeMap = Teuchos::null,
02600                                   const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null)
02601   {
02602     Teuchos::RCP<CrsMatrixType> destMatrix;
02603     sourceMatrix->exportAndFillComplete (destMatrix,exporter, domainMap, rangeMap, params);
02604     return destMatrix;
02605   }
02606 } // namespace Tpetra
02607 
02608 // Include KokkosRefactor partial specialisation if enabled
02609 #if defined(TPETRA_HAVE_KOKKOS_REFACTOR)
02610 #include "Tpetra_KokkosRefactor_CrsMatrix_decl.hpp"
02611 #endif
02612 
02623 #endif // TPETRA_CRSMATRIX_DECL_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines