Tpetra Matrix/Vector Services Version of the Day
Tpetra_CrsMatrix_decl.hpp
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //          Tpetra: Templated Linear Algebra Services Package
00005 //                 Copyright (2008) Sandia Corporation
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ************************************************************************
00040 // @HEADER
00041 
00042 #ifndef TPETRA_CRSMATRIX_DECL_HPP
00043 #define TPETRA_CRSMATRIX_DECL_HPP
00044 
00045 #include <Kokkos_DefaultNode.hpp>
00046 #include <Kokkos_DefaultKernels.hpp>
00047 
00048 #include "Tpetra_ConfigDefs.hpp"
00049 #include "Tpetra_RowMatrix_decl.hpp"
00050 #include "Tpetra_Exceptions.hpp"
00051 #include "Tpetra_DistObject.hpp"
00052 #include "Tpetra_CrsGraph.hpp"
00053 #include "Tpetra_Vector.hpp"
00054 
00055 
00056 namespace Tpetra {
00057 
00059 
00196   template <class Scalar,
00197             class LocalOrdinal  = int,
00198             class GlobalOrdinal = LocalOrdinal,
00199             class Node          = KokkosClassic::DefaultNode::DefaultNodeType,
00200             class LocalMatOps   = typename KokkosClassic::DefaultKernels<Scalar,LocalOrdinal,Node>::SparseOps >
00201   class CrsMatrix : public RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>,
00202                     public DistObject<char, LocalOrdinal,GlobalOrdinal,Node> {
00203   public:
00205 
00206 
00208     typedef Scalar                                scalar_type;
00210     typedef LocalOrdinal                          local_ordinal_type;
00212     typedef GlobalOrdinal                         global_ordinal_type;
00214     typedef Node                                  node_type;
00218     typedef LocalMatOps   mat_vec_type;
00222     typedef LocalMatOps   mat_solve_type;
00223 
00225     typedef Map<LocalOrdinal,GlobalOrdinal,Node>  map_type;
00226 
00228 
00229 
00230 
00248     CrsMatrix (const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& rowMap,
00249                size_t maxNumEntriesPerRow,
00250                ProfileType pftype = DynamicProfile,
00251                const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null);
00252 
00270     CrsMatrix (const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& rowMap,
00271                const ArrayRCP<const size_t>& NumEntriesPerRowToAlloc,
00272                ProfileType pftype = DynamicProfile,
00273                const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null);
00274 
00297     CrsMatrix (const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& rowMap,
00298                const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& colMap,
00299                size_t maxNumEntriesPerRow,
00300                ProfileType pftype = DynamicProfile,
00301                const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null);
00302 
00325     CrsMatrix (const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& rowMap,
00326                const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& colMap,
00327                const ArrayRCP<const size_t>& NumEntriesPerRowToAlloc,
00328                ProfileType pftype = DynamicProfile,
00329                const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null);
00330 
00355     explicit CrsMatrix (const Teuchos::RCP<const CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps> >& graph,
00356                         const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null);
00357 
00380     CrsMatrix (const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& rowMap,
00381                const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& colMap,
00382                const ArrayRCP<size_t> & rowPointers,
00383                const ArrayRCP<LocalOrdinal> & columnIndices,
00384                const ArrayRCP<Scalar> & values,
00385                const RCP<ParameterList>& params = null);
00386 
00387     // This friend declaration makes the clone() method work.
00388     template <class S2, class LO2, class GO2, class N2, class LMO2>
00389     friend class CrsMatrix;
00390 
00415     template <class Node2>
00416     RCP<CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node2, typename KokkosClassic::DefaultKernels<void,LocalOrdinal,Node2>::SparseOps> >
00417     clone (const RCP<Node2>& node2, const RCP<ParameterList>& params = null) const
00418     {
00419       const char tfecfFuncName[] = "clone";
00420 
00421       // Get parameter values.  Set them initially to their default values.
00422       bool fillCompleteClone = true;
00423       bool useLocalIndices = this->hasColMap ();
00424       ProfileType pftype = StaticProfile;
00425       if (! params.is_null ()) {
00426         fillCompleteClone = params->get ("fillComplete clone", fillCompleteClone);
00427         useLocalIndices = params->get ("Locally indexed clone", useLocalIndices);
00428 
00429         bool staticProfileClone = true;
00430         staticProfileClone = params->get ("Static profile clone", staticProfileClone);
00431         pftype = staticProfileClone ? StaticProfile : DynamicProfile;
00432       }
00433 
00434       TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
00435         ! this->hasColMap () && useLocalIndices, std::runtime_error,
00436         ": You requested that the returned clone have local indices, but the "
00437         "the source matrix does not have a column Map yet.");
00438 
00439       typedef typename KokkosClassic::DefaultKernels<void,LocalOrdinal,Node2>::SparseOps LocalMatOps2;
00440       typedef CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node2, LocalMatOps2> CrsMatrix2;
00441       typedef Map<LocalOrdinal, GlobalOrdinal, Node2> Map2;
00442       RCP<const Map2> clonedRowMap = this->getRowMap ()->template clone<Node2> (node2);
00443 
00444       RCP<CrsMatrix2> clonedMatrix;
00445       ArrayRCP<const size_t> numEntries;
00446       size_t numEntriesForAll = 0;
00447       if (! staticGraph_->indicesAreAllocated ()) {
00448         if (! staticGraph_->numAllocPerRow_.is_null ()) {
00449           numEntries = staticGraph_->numAllocPerRow_;
00450         }
00451         else {
00452           numEntriesForAll = staticGraph_->numAllocForAllRows_;
00453         }
00454       }
00455       else if (! staticGraph_->numRowEntries_.is_null ()) {
00456         numEntries = staticGraph_->numRowEntries_;
00457       }
00458       else if (staticGraph_->nodeNumAllocated_ == 0) {
00459         numEntriesForAll = 0;
00460       }
00461       else {
00462         // We're left with the case that we have optimized storage.
00463         // In this case, we have to construct a list of row sizes.
00464         TEUCHOS_TEST_FOR_EXCEPTION(
00465           getProfileType() != StaticProfile, std::logic_error,
00466           "Internal logic error. Please report this to Tpetra team." )
00467 
00468         const size_t numRows = this->getNodeNumRows ();
00469         numEntriesForAll = 0;
00470         ArrayRCP<size_t> numEnt;
00471         if (numRows != 0) {
00472           numEnt = arcp<size_t> (numRows);
00473         }
00474         for (size_t i = 0; i < numRows; ++i) {
00475           numEnt[i] = staticGraph_->rowPtrs_[i+1] - staticGraph_->rowPtrs_[i];
00476         }
00477         numEntries = numEnt;
00478       }
00479 
00480       RCP<ParameterList> matrixparams =
00481         params.is_null () ? null : sublist (params,"CrsMatrix");
00482       if (useLocalIndices) {
00483         RCP<const Map2> clonedColMap =
00484           this->getColMap ()->template clone<Node2> (node2);
00485         if (numEntries.is_null ()) {
00486           clonedMatrix = rcp (new CrsMatrix2 (clonedRowMap, clonedColMap,
00487                                               numEntriesForAll, pftype,
00488                                               matrixparams));
00489         }
00490         else {
00491           clonedMatrix = rcp (new CrsMatrix2 (clonedRowMap, clonedColMap,
00492                                               numEntries, pftype,
00493                                               matrixparams));
00494         }
00495       }
00496       else {
00497         if (numEntries.is_null ()) {
00498           clonedMatrix = rcp (new CrsMatrix2 (clonedRowMap, numEntriesForAll,
00499                                               pftype, matrixparams));
00500         }
00501         else {
00502           clonedMatrix = rcp (new CrsMatrix2 (clonedRowMap, numEntries, pftype,
00503                                               matrixparams));
00504         }
00505       }
00506       // done with these
00507       numEntries = null;
00508       numEntriesForAll = 0;
00509 
00510       if (useLocalIndices) {
00511         clonedMatrix->allocateValues (LocalIndices,
00512                                       CrsMatrix2::GraphNotYetAllocated);
00513         if (this->isLocallyIndexed ()) {
00514           ArrayView<const LocalOrdinal> linds;
00515           ArrayView<const Scalar>       vals;
00516           for (LocalOrdinal lrow = clonedRowMap->getMinLocalIndex ();
00517                lrow <= clonedRowMap->getMaxLocalIndex ();
00518                ++lrow) {
00519             this->getLocalRowView (lrow, linds, vals);
00520             if (linds.size ()) {
00521               clonedMatrix->insertLocalValues (lrow, linds, vals);
00522             }
00523           }
00524         }
00525         else { // this->isGloballyIndexed()
00526           Array<LocalOrdinal> linds;
00527           Array<Scalar> vals;
00528           for (LocalOrdinal lrow = clonedRowMap->getMinLocalIndex ();
00529                lrow <= clonedRowMap->getMaxLocalIndex ();
00530                ++lrow) {
00531             size_t theNumEntries = this->getNumEntriesInLocalRow (lrow);
00532             if (theNumEntries > Teuchos::as<size_t> (linds.size ())) {
00533               linds.resize (theNumEntries);
00534             }
00535             if (theNumEntries > Teuchos::as<size_t> (vals.size ())) {
00536               vals.resize (theNumEntries);
00537             }
00538             this->getLocalRowCopy (clonedRowMap->getGlobalElement (lrow),
00539                                    linds (), vals (), theNumEntries);
00540             if (theNumEntries != 0) {
00541               clonedMatrix->insertLocalValues (lrow, linds (0, theNumEntries),
00542                                                vals (0, theNumEntries));
00543             }
00544           }
00545         }
00546       }
00547       else { // useGlobalIndices
00548         clonedMatrix->allocateValues (GlobalIndices,
00549                                       CrsMatrix2::GraphNotYetAllocated);
00550         if (this->isGloballyIndexed ()) {
00551           ArrayView<const GlobalOrdinal> ginds;
00552           ArrayView<const Scalar>         vals;
00553           for (GlobalOrdinal grow = clonedRowMap->getMinGlobalIndex ();
00554                grow <= clonedRowMap->getMaxGlobalIndex ();
00555                ++grow) {
00556             this->getGlobalRowView (grow, ginds, vals);
00557             if (ginds.size () > 0) {
00558               clonedMatrix->insertGlobalValues (grow, ginds, vals);
00559             }
00560           }
00561         }
00562         else { // this->isLocallyIndexed()
00563           Array<GlobalOrdinal> ginds;
00564           Array<Scalar> vals;
00565           for (GlobalOrdinal grow = clonedRowMap->getMinGlobalIndex ();
00566                grow <= clonedRowMap->getMaxGlobalIndex ();
00567                ++grow) {
00568             size_t theNumEntries = this->getNumEntriesInGlobalRow (grow);
00569             if (theNumEntries > Teuchos::as<size_t> (ginds.size ())) {
00570               ginds.resize (theNumEntries);
00571             }
00572             if (theNumEntries > Teuchos::as<size_t> (vals.size ())) {
00573               vals.resize (theNumEntries);
00574             }
00575             this->getGlobalRowCopy (grow, ginds (), vals (), theNumEntries);
00576             if (theNumEntries != 0) {
00577               clonedMatrix->insertGlobalValues (grow, ginds (0, theNumEntries),
00578                                                 vals (0, theNumEntries));
00579             }
00580           }
00581         }
00582       }
00583 
00584       if (fillCompleteClone) {
00585         RCP<ParameterList> fillparams =
00586           params.is_null () ? Teuchos::null : sublist (params, "fillComplete");
00587         try {
00588           RCP<const Map2> clonedRangeMap;
00589           RCP<const Map2> clonedDomainMap;
00590           if (! this->getRangeMap ().is_null () &&
00591               this->getRangeMap () != clonedRowMap) {
00592             clonedRangeMap  = this->getRangeMap ()->template clone<Node2> (node2);
00593           }
00594           else {
00595             clonedRangeMap = clonedRowMap;
00596           }
00597           if (! this->getDomainMap ().is_null () &&
00598               this->getDomainMap () != clonedRowMap) {
00599             clonedDomainMap = this->getDomainMap ()->template clone<Node2> (node2);
00600           }
00601           else {
00602             clonedDomainMap = clonedRowMap;
00603           }
00604           clonedMatrix->fillComplete (clonedDomainMap, clonedRangeMap,
00605                                       fillparams);
00606         }
00607         catch (std::exception &e) {
00608           const bool caughtExceptionOnClone = true;
00609           TEUCHOS_TEST_FOR_EXCEPTION(
00610             caughtExceptionOnClone, std::runtime_error,
00611             Teuchos::typeName (*this) << std::endl << "clone: " << std::endl <<
00612             "Caught the following exception while calling fillComplete() on a "
00613             "clone of type" << std::endl << Teuchos::typeName (*clonedMatrix)
00614             << ": " << std::endl << e.what () << std::endl);
00615         }
00616       }
00617       return clonedMatrix;
00618     }
00619 
00621     virtual ~CrsMatrix ();
00622 
00624 
00625 
00626 
00649     //
00695     void
00696     insertGlobalValues (const GlobalOrdinal globalRow,
00697                         const ArrayView<const GlobalOrdinal>& cols,
00698                         const ArrayView<const Scalar>& vals);
00699 
00739     void
00740     insertLocalValues (const LocalOrdinal localRow,
00741                        const ArrayView<const LocalOrdinal> &cols,
00742                        const ArrayView<const Scalar> &vals);
00743 
00762     void
00763     replaceGlobalValues (GlobalOrdinal globalRow,
00764                          const ArrayView<const GlobalOrdinal>& cols,
00765                          const ArrayView<const Scalar>& vals);
00766 
00780     void
00781     replaceLocalValues (LocalOrdinal localRow,
00782                         const ArrayView<const LocalOrdinal> &cols,
00783                         const ArrayView<const Scalar>       &vals);
00784 
00811     void
00812     sumIntoGlobalValues (const GlobalOrdinal globalRow,
00813                          const ArrayView<const GlobalOrdinal> &cols,
00814                          const ArrayView<const Scalar>        &vals);
00815 
00824     void
00825     sumIntoLocalValues (const LocalOrdinal localRow,
00826                         const ArrayView<const LocalOrdinal>  &cols,
00827                         const ArrayView<const Scalar>        &vals);
00828 
00830     void setAllToScalar (const Scalar &alpha);
00831 
00833     void scale (const Scalar &alpha);
00834 
00836 
00843     void
00844     setAllValues (const ArrayRCP<size_t>& rowPointers,
00845                   const ArrayRCP<LocalOrdinal>& columnIndices,
00846                   const ArrayRCP<Scalar>& values);
00847 
00849 
00850 
00851 
00880     void globalAssemble();
00881 
00895     void resumeFill (const RCP<ParameterList>& params = null);
00896 
00917     void fillComplete(const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &domainMap,
00918                       const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rangeMap,
00919                       const RCP<ParameterList> &params = null);
00920 
00933     void fillComplete(const RCP<ParameterList> &params = null);
00934 
00946     void
00947     expertStaticFillComplete (const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & domainMap,
00948                               const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & rangeMap,
00949                               const RCP<const Import<LocalOrdinal,GlobalOrdinal,Node> > &importer=Teuchos::null,
00950                               const RCP<const Export<LocalOrdinal,GlobalOrdinal,Node> > &exporter=Teuchos::null,
00951                               const RCP<ParameterList> &params=Teuchos::null);
00952 
00965     void
00966     replaceDomainMapAndImporter (const Teuchos::RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> >& newDomainMap,
00967                                  Teuchos::RCP<const Tpetra::Import<LocalOrdinal,GlobalOrdinal,Node> >& newImporter);
00968 
00982     virtual void
00983     removeEmptyProcessesInPlace (const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> >& newMap);
00984 
00986 
00987 
00988 
00990     RCP<const Comm<int> > getComm() const;
00991 
00993     RCP<Node> getNode() const;
00994 
00996     RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > getRowMap() const;
00997 
00999     RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > getColMap() const;
01000 
01002     RCP<const RowGraph<LocalOrdinal,GlobalOrdinal,Node> > getGraph() const;
01003 
01005     RCP<const CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps> > getCrsGraph() const;
01006 
01026     global_size_t getGlobalNumRows() const;
01027 
01033     global_size_t getGlobalNumCols() const;
01034 
01041     size_t getNodeNumRows() const;
01042 
01046     size_t getNodeNumCols() const;
01047 
01049     GlobalOrdinal getIndexBase() const;
01050 
01052     global_size_t getGlobalNumEntries() const;
01053 
01055     size_t getNodeNumEntries() const;
01056 
01058 
01059     size_t getNumEntriesInGlobalRow(GlobalOrdinal globalRow) const;
01060 
01062 
01063     size_t getNumEntriesInLocalRow(LocalOrdinal localRow) const;
01064 
01066 
01068     global_size_t getGlobalNumDiags() const;
01069 
01071 
01073     size_t getNodeNumDiags() const;
01074 
01076 
01078     size_t getGlobalMaxNumRowEntries() const;
01079 
01081 
01083     size_t getNodeMaxNumRowEntries() const;
01084 
01086     bool hasColMap() const;
01087 
01089 
01091     bool isLowerTriangular() const;
01092 
01094 
01096     bool isUpperTriangular() const;
01097 
01117     bool isLocallyIndexed() const;
01118 
01138     bool isGloballyIndexed() const;
01139 
01162     bool isFillComplete() const;
01163 
01186     bool isFillActive() const;
01187 
01189 
01195     bool isStorageOptimized() const;
01196 
01198     ProfileType getProfileType() const;
01199 
01201     bool isStaticGraph() const;
01202 
01204 
01210     typename ScalarTraits<Scalar>::magnitudeType getFrobeniusNorm() const;
01211 
01212 
01214     virtual bool supportsRowViews() const;
01215 
01217 
01227     void
01228     getGlobalRowCopy (GlobalOrdinal GlobalRow,
01229                       const ArrayView<GlobalOrdinal> &Indices,
01230                       const ArrayView<Scalar> &Values,
01231                       size_t &NumEntries) const;
01232 
01234 
01246     void
01247     getLocalRowCopy (LocalOrdinal LocalRow,
01248                      const ArrayView<LocalOrdinal> &Indices,
01249                      const ArrayView<Scalar> &Values,
01250                      size_t &NumEntries) const;
01251 
01253 
01262     void
01263     getGlobalRowView (GlobalOrdinal GlobalRow,
01264                       ArrayView<const GlobalOrdinal> &indices,
01265                       ArrayView<const Scalar> &values) const;
01266 
01268 
01277     void
01278     getLocalRowView (LocalOrdinal LocalRow,
01279                      ArrayView<const LocalOrdinal> &indices,
01280                      ArrayView<const Scalar> &values) const;
01281 
01283 
01285     void getLocalDiagCopy (Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &diag) const;
01286 
01324     void getLocalDiagOffsets (Teuchos::ArrayRCP<size_t>& offsets) const;
01325 
01344     void
01345     getLocalDiagCopy (Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& diag,
01346                       const Teuchos::ArrayView<const size_t>& offsets) const;
01347 
01349     void leftScale(const Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& x);
01350 
01352     void rightScale(const Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& x);
01353 
01355 
01356 
01357 
01406     template <class DomainScalar, class RangeScalar>
01407     void
01408     localMultiply (const MultiVector<DomainScalar,LocalOrdinal,GlobalOrdinal,Node>& X,
01409                    MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,Node>& Y,
01410                    Teuchos::ETransp trans,
01411                    RangeScalar alpha,
01412                    RangeScalar beta) const;
01413 
01438     template <class DomainScalar, class RangeScalar>
01439     void
01440     localGaussSeidel (const MultiVector<DomainScalar,LocalOrdinal,GlobalOrdinal,Node> &B,
01441                       MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,Node> &X,
01442                       const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &D,
01443                       const RangeScalar& dampingFactor,
01444                       const KokkosClassic::ESweepDirection direction) const;
01445 
01462     template <class DomainScalar, class RangeScalar>
01463     void
01464     localSolve (const MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,Node>& Y,
01465                 MultiVector<DomainScalar,LocalOrdinal,GlobalOrdinal,Node>& X,
01466                 Teuchos::ETransp trans) const;
01467 
01469     template <class T>
01470     RCP<CrsMatrix<T,LocalOrdinal,GlobalOrdinal,Node> > convert() const;
01471 
01473 
01474 
01475 
01486     void
01487     apply (const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& X,
01488            MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>&Y,
01489            Teuchos::ETransp mode = Teuchos::NO_TRANS,
01490            Scalar alpha = ScalarTraits<Scalar>::one(),
01491            Scalar beta = ScalarTraits<Scalar>::zero()) const;
01492 
01494     bool hasTransposeApply() const;
01495 
01499     RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > getDomainMap() const;
01500 
01504     RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > getRangeMap() const;
01505 
01507 
01508 
01509 
01574     void
01575     gaussSeidel (const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &B,
01576                  MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &X,
01577                  const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &D,
01578                  const Scalar& dampingFactor,
01579                  const ESweepDirection direction,
01580                  const int numSweeps) const;
01581 
01610     void
01611     gaussSeidelCopy (MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &X,
01612                      const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &B,
01613                      const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &D,
01614                      const Scalar& dampingFactor,
01615                      const ESweepDirection direction,
01616                      const int numSweeps,
01617                      const bool zeroInitialGuess) const;
01618 
01629     virtual Teuchos::RCP<RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
01630     add (const Scalar& alpha,
01631          const RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& A,
01632          const Scalar& beta,
01633          const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> >& domainMap,
01634          const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> >& rangeMap,
01635          const Teuchos::RCP<Teuchos::ParameterList>& params) const;
01636 
01638 
01639 
01640 
01642     std::string description() const;
01643 
01645     void
01646     describe (Teuchos::FancyOStream &out,
01647               const Teuchos::EVerbosityLevel verbLevel =
01648               Teuchos::Describable::verbLevel_default) const;
01649 
01651 
01652 
01653 
01654     virtual bool
01655     checkSizes (const SrcDistObject& source);
01656 
01657     virtual void
01658     copyAndPermute (const SrcDistObject& source,
01659                     size_t numSameIDs,
01660                     const ArrayView<const LocalOrdinal> &permuteToLIDs,
01661                     const ArrayView<const LocalOrdinal> &permuteFromLIDs);
01662 
01663     virtual void
01664     packAndPrepare (const SrcDistObject& source,
01665                     const Teuchos::ArrayView<const LocalOrdinal>& exportLIDs,
01666                     Teuchos::Array<char>& exports,
01667                     const Teuchos::ArrayView<size_t>& numPacketsPerLID,
01668                     size_t& constantNumPackets,
01669                     Distributor& distor);
01670 
01680     void
01681     unpackAndCombine (const Teuchos::ArrayView<const LocalOrdinal> &importLIDs,
01682                       const Teuchos::ArrayView<const char> &imports,
01683                       const Teuchos::ArrayView<size_t> &numPacketsPerLID,
01684                       size_t constantNumPackets,
01685                       Distributor &distor,
01686                       CombineMode combineMode);
01688 
01689 
01690 
01699     virtual void
01700     pack (const Teuchos::ArrayView<const LocalOrdinal>& exportLIDs,
01701           Teuchos::Array<char>& exports,
01702           const Teuchos::ArrayView<size_t>& numPacketsPerLID,
01703           size_t& constantNumPackets,
01704           Distributor& distor) const;
01705 
01707   private:
01708     // Friend declaration for nonmember function.
01709     template<class CrsMatrixType>
01710     friend Teuchos::RCP<CrsMatrixType>
01711     importAndFillCompleteCrsMatrix (const Teuchos::RCP<const CrsMatrixType>& sourceMatrix,
01712                                     const Import<typename CrsMatrixType::local_ordinal_type,
01713                                                  typename CrsMatrixType::global_ordinal_type,
01714                                                  typename CrsMatrixType::node_type>& importer,
01715                                     const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
01716                                                                  typename CrsMatrixType::global_ordinal_type,
01717                                                                  typename CrsMatrixType::node_type> >& domainMap,
01718                                     const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
01719                                                                  typename CrsMatrixType::global_ordinal_type,
01720                                                                  typename CrsMatrixType::node_type> >& rangeMap,
01721                                     const Teuchos::RCP<Teuchos::ParameterList>& params);
01722 
01723     // Friend declaration for nonmember function.
01724     template<class CrsMatrixType>
01725     friend Teuchos::RCP<CrsMatrixType>
01726     exportAndFillCompleteCrsMatrix (const Teuchos::RCP<const CrsMatrixType>& sourceMatrix,
01727                                     const Export<typename CrsMatrixType::local_ordinal_type,
01728                                                  typename CrsMatrixType::global_ordinal_type,
01729                                                  typename CrsMatrixType::node_type>& exporter,
01730                                     const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
01731                                                                  typename CrsMatrixType::global_ordinal_type,
01732                                                                  typename CrsMatrixType::node_type> >& domainMap,
01733                                     const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
01734                                                                  typename CrsMatrixType::global_ordinal_type,
01735                                                                  typename CrsMatrixType::node_type> >& rangeMap,
01736                                     const Teuchos::RCP<Teuchos::ParameterList>& params);
01737 
01744     Teuchos::RCP<CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node, LocalMatOps> >
01745     importAndFillComplete (const Import<LocalOrdinal, GlobalOrdinal, Node>& importer,
01746                            const Teuchos::RCP<const map_type>& domainMap,
01747                            const Teuchos::RCP<const map_type>& rangeMap,
01748                            const Teuchos::RCP<Teuchos::ParameterList>& params) const;
01749 
01756     Teuchos::RCP<CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps> >
01757     exportAndFillComplete (const Export<LocalOrdinal, GlobalOrdinal, Node>& exporter,
01758                            const Teuchos::RCP<const map_type>& domainMap = Teuchos::null,
01759                            const Teuchos::RCP<const map_type>& rangeMap = Teuchos::null,
01760                            const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null) const;
01761 
01762   private:
01763     // We forbid copy construction by declaring this method private
01764     // and not implementing it.
01765     CrsMatrix (const CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps> &rhs);
01766 
01767     // We forbid assignment (operator=) by declaring this method
01768     // private and not implementing it.
01769     CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>&
01770     operator= (const CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps> &rhs);
01771 
01777     void
01778     insertGlobalValuesFiltered (const GlobalOrdinal globalRow,
01779                                 const ArrayView<const GlobalOrdinal> &indices,
01780                                 const ArrayView<const Scalar>        &values);
01781 
01787     void
01788     insertLocalValuesFiltered (const LocalOrdinal localRow,
01789                                const ArrayView<const LocalOrdinal> &indices,
01790                                const ArrayView<const Scalar>       &values);
01791 
01799     void
01800     combineGlobalValues (const GlobalOrdinal globalRowIndex,
01801                          const Teuchos::ArrayView<const GlobalOrdinal> columnIndices,
01802                          const Teuchos::ArrayView<const Scalar> values,
01803                          const Tpetra::CombineMode combineMode);
01804 
01836     template<class BinaryFunction>
01837     void
01838     transformLocalValues (LocalOrdinal localRow,
01839                           const Teuchos::ArrayView<const LocalOrdinal>& indices,
01840                           const Teuchos::ArrayView<const Scalar>        & values,
01841                           BinaryFunction f)
01842     {
01843       typedef LocalOrdinal LO;
01844       typedef GlobalOrdinal GO;
01845       typedef Node NT;
01846       using Teuchos::Array;
01847       using Teuchos::ArrayView;
01848 
01849       TEUCHOS_TEST_FOR_EXCEPTION(
01850         ! isFillActive (),
01851         std::runtime_error,
01852         "Tpetra::CrsMatrix::transformLocalValues: Fill must be active in order "
01853         "to call this method.  That is, isFillActive() must return true.  If "
01854         "you have already called fillComplete(), you need to call resumeFill() "
01855         "before you can replace values.");
01856       TEUCHOS_TEST_FOR_EXCEPTION(
01857         values.size () != indices.size (),
01858         std::runtime_error,
01859         "Tpetra::CrsMatrix::transformLocalValues: values.size () = "
01860         << values.size () << " != indices.size () = " << indices.size ()
01861         << ".");
01862       TEUCHOS_TEST_FOR_EXCEPTION(
01863         ! this->hasColMap (),
01864         std::runtime_error,
01865         "Tpetra::CrsMatrix::transformLocalValues: We cannot transform local "
01866         "indices without a column map.");
01867       const bool isLocalRow = getRowMap ()->isNodeLocalElement (localRow);
01868       TEUCHOS_TEST_FOR_EXCEPTION(
01869         ! isLocalRow,
01870         std::runtime_error,
01871         "Tpetra::CrsMatrix::transformLocalValues: The specified local row "
01872         << localRow << " does not belong to this process "
01873         << getRowMap ()->getComm ()->getRank () << ".");
01874 
01875       RowInfo rowInfo = staticGraph_->getRowInfo (localRow);
01876       if (indices.size () > 0) {
01877         ArrayView<Scalar> curVals = this->getViewNonConst (rowInfo);
01878         if (isLocallyIndexed ()) {
01879           staticGraph_->template transformLocalValues<Scalar, BinaryFunction> (rowInfo, curVals, indices, values, f);
01880         }
01881         else if (isGloballyIndexed ()) {
01882           // Convert the given local indices to global indices.
01883           const Map<LO, GO, NT>& colMap = * (this->getColMap ());
01884           Array<GO> gindices (indices.size ());
01885           typename ArrayView<const LO>::iterator lindit = indices.begin();
01886           typename Array<GO>::iterator           gindit = gindices.begin();
01887           while (lindit != indices.end()) {
01888             // There is no need to filter out indices not in the column
01889             // Map.  Those that aren't will be mapped to invalid(),
01890             // which transformGlobalValues() will ignore.
01891             *gindit++ = colMap.getGlobalElement (*lindit++);
01892           }
01893           staticGraph_->template transformGlobalValues<Scalar, BinaryFunction> (rowInfo, curVals, gindices (), values, f);
01894         }
01895       }
01896     }
01897 
01898 
01928     template<class BinaryFunction>
01929     void
01930     transformGlobalValues (GlobalOrdinal globalRow,
01931                            const Teuchos::ArrayView<const GlobalOrdinal>& indices,
01932                            const Teuchos::ArrayView<const Scalar>        & values,
01933                            BinaryFunction f)
01934     {
01935       typedef LocalOrdinal LO;
01936       typedef GlobalOrdinal GO;
01937       typedef Node NT;
01938       using Teuchos::Array;
01939       using Teuchos::ArrayView;
01940 
01941       TEUCHOS_TEST_FOR_EXCEPTION(
01942         ! isFillActive (),
01943         std::runtime_error,
01944         "Tpetra::CrsMatrix::transformGlobalValues: Fill must be active in order "
01945         "to call this method.  That is, isFillActive() must return true.  If "
01946         "you have already called fillComplete(), you need to call resumeFill() "
01947         "before you can replace values.");
01948       TEUCHOS_TEST_FOR_EXCEPTION(
01949         values.size () != indices.size (),
01950         std::runtime_error,
01951         "Tpetra::CrsMatrix::transformGlobalValues: values.size () = "
01952         << values.size () << " != indices.size () = " << indices.size ()
01953         << ".");
01954 
01955       const LO lrow = this->getRowMap()->getLocalElement(globalRow);
01956 
01957       if (lrow == OTL::invalid()) {
01958         // FIXME (mfh 16 May 2013) We're using this exception to do
01959         // sumIntoGlobalValues for nonowned rows, so we might want to
01960         // avoid the overhead of constructing the fancy exception
01961         // message each time if we don't plan to use it.
01962 
01963         // The exception test macro doesn't let you pass an additional
01964         // argument to the exception's constructor, so we don't use it.
01965         std::ostringstream os;
01966         os << "transformGlobalValues: The given global row index "
01967            << globalRow << " is not owned by the calling process (rank "
01968            << this->getRowMap()->getComm()->getRank() << ").";
01969         throw Details::InvalidGlobalRowIndex<GO> (os.str (), globalRow);
01970       }
01971 
01972       RowInfo rowInfo = staticGraph_->getRowInfo (lrow);
01973       if (indices.size () > 0) {
01974         ArrayView<Scalar> curVals = this->getViewNonConst (rowInfo);
01975         if (isLocallyIndexed ()) {
01976           // Convert global indices to local indices.
01977           const Map<LO, GO, NT> &colMap = * (this->getColMap ());
01978           Array<LO> lindices (indices.size ());
01979           typename ArrayView<const GO>::iterator gindit = indices.begin();
01980           typename Array<LO>::iterator           lindit = lindices.begin();
01981           while (gindit != indices.end()) {
01982             // There is no need to filter out indices not in the column
01983             // Map.  Those that aren't will be mapped to invalid(),
01984             // which transformLocalValues() will ignore.
01985             *lindit++ = colMap.getLocalElement (*gindit++);
01986           }
01987           staticGraph_->template transformLocalValues<Scalar, BinaryFunction> (rowInfo, curVals, lindices (), values, f);
01988         }
01989         else if (isGloballyIndexed ()) {
01990           staticGraph_->template transformGlobalValues<Scalar, BinaryFunction> (rowInfo, curVals, indices, values, f);
01991         }
01992       }
01993     }
01994 
01995   protected:
01996     // useful typedefs
01997     typedef OrdinalTraits<LocalOrdinal>                     OTL;
01998     typedef ScalarTraits<Scalar>                            STS;
01999     typedef typename STS::magnitudeType               Magnitude;
02000     typedef ScalarTraits<Magnitude>                         STM;
02001     typedef MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>       MV;
02002     typedef Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>             V;
02003     typedef CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>  Graph;
02004     typedef typename LocalMatOps::template bind_scalar<Scalar>::other_type                    sparse_ops_type;
02005     typedef typename sparse_ops_type::template graph<LocalOrdinal,Node>::graph_type          local_graph_type;
02006     typedef typename sparse_ops_type::template matrix<Scalar,LocalOrdinal,Node>::matrix_type local_matrix_type;
02007 
02008     typedef Export<LocalOrdinal, GlobalOrdinal, Node> export_type;
02009     typedef Import<LocalOrdinal, GlobalOrdinal, Node> import_type;
02010 
02011     // Enums
02012     enum GraphAllocationStatus {
02013       GraphAlreadyAllocated,
02014       GraphNotYetAllocated
02015     };
02016 
02033     void allocateValues (ELocalGlobal lg, GraphAllocationStatus gas);
02034 
02040     void sortEntries();
02041 
02047     void mergeRedundantEntries();
02048 
02056     void clearGlobalConstants();
02057 
02066     void computeGlobalConstants();
02067 
02080     mutable Teuchos::RCP<MV> importMV_;
02081 
02094     mutable Teuchos::RCP<MV> exportMV_;
02095 
02115     Teuchos::RCP<MV>
02116     getColumnMapMultiVector (const MV& X_domainMap,
02117                              const bool force = false) const;
02118 
02140     Teuchos::RCP<MV>
02141     getRowMapMultiVector (const MV& Y_rangeMap,
02142                           const bool force = false) const;
02143 
02145     void
02146     applyNonTranspose (const MV& X_in,
02147                        MV& Y_in,
02148                        Scalar alpha,
02149                        Scalar beta) const;
02150 
02152     void
02153     applyTranspose (const MV& X_in,
02154                     MV& Y_in,
02155                     const Teuchos::ETransp mode,
02156                     Scalar alpha,
02157                     Scalar beta) const;
02158 
02159     // matrix data accessors
02160     ArrayView<const Scalar>    getView(RowInfo rowinfo) const;
02161     ArrayView<      Scalar>    getViewNonConst(RowInfo rowinfo);
02162     // local Kokkos objects
02163 
02169     void fillLocalMatrix(const RCP<ParameterList> &params);
02175     void fillLocalGraphAndMatrix(const RCP<ParameterList> &params);
02176     // debugging
02177     void checkInternalState() const;
02178 
02190 
02191     RCP<const Graph> staticGraph_;
02192     RCP<      Graph>     myGraph_;
02194 
02199     RCP<sparse_ops_type>   lclMatOps_;
02206     RCP<local_matrix_type> lclMatrix_;
02207 
02220 
02221     ArrayRCP<Scalar> values1D_;
02222     ArrayRCP<Array<Scalar> > values2D_;
02224 
02226     bool fillComplete_;
02227 
02240     std::map<GlobalOrdinal, Array<std::pair<GlobalOrdinal,Scalar> > > nonlocals_;
02241 
02248     mutable Magnitude frobNorm_;
02249   }; // class CrsMatrix
02250 
02257   template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
02258   Teuchos::RCP<CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >
02259   createCrsMatrix (const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &map,
02260                    size_t maxNumEntriesPerRow = 0,
02261                    const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null)
02262   {
02263     using Teuchos::rcp;
02264     typedef CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> matrix_type;
02265     return rcp (new matrix_type (map, maxNumEntriesPerRow, DynamicProfile, params));
02266   }
02267 
02319   template<class CrsMatrixType>
02320   Teuchos::RCP<CrsMatrixType>
02321   importAndFillCompleteCrsMatrix (const Teuchos::RCP<const CrsMatrixType>& sourceMatrix,
02322                                   const Import<typename CrsMatrixType::local_ordinal_type,
02323                                                typename CrsMatrixType::global_ordinal_type,
02324                                                typename CrsMatrixType::node_type>& importer,
02325                                   const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
02326                                                                typename CrsMatrixType::global_ordinal_type,
02327                                                                typename CrsMatrixType::node_type> >& domainMap = Teuchos::null,
02328                                   const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
02329                                                                typename CrsMatrixType::global_ordinal_type,
02330                                                                typename CrsMatrixType::node_type> >& rangeMap = Teuchos::null,
02331                                   const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null)
02332   {
02333     return sourceMatrix->importAndFillComplete (importer, domainMap, rangeMap, params);
02334   }
02335 
02369   template<class CrsMatrixType>
02370   Teuchos::RCP<CrsMatrixType>
02371   exportAndFillCompleteCrsMatrix (const Teuchos::RCP<const CrsMatrixType>& sourceMatrix,
02372                                   const Export<typename CrsMatrixType::local_ordinal_type,
02373                                                typename CrsMatrixType::global_ordinal_type,
02374                                                typename CrsMatrixType::node_type>& exporter,
02375                                   const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
02376                                                                typename CrsMatrixType::global_ordinal_type,
02377                                                                typename CrsMatrixType::node_type> >& domainMap = Teuchos::null,
02378                                   const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type,
02379                                                                typename CrsMatrixType::global_ordinal_type,
02380                                                                typename CrsMatrixType::node_type> >& rangeMap = Teuchos::null,
02381                                   const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null)
02382   {
02383     return sourceMatrix->exportAndFillComplete (exporter, domainMap, rangeMap, params);
02384   }
02385 } // namespace Tpetra
02386 
02397 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines