|
Tpetra Matrix/Vector Services
Version of the Day
|
|
Tpetra Matrix/Vector Services
Version of the Day
|
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 { 00173 template <class Scalar = RowMatrix<>::scalar_type, 00174 class LocalOrdinal = typename RowMatrix<Scalar>::local_ordinal_type, 00175 class GlobalOrdinal = typename RowMatrix<Scalar, LocalOrdinal>::global_ordinal_type, 00176 class Node = typename RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal>::node_type> 00177 class CrsMatrix : 00178 public RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>, 00179 public DistObject<char, LocalOrdinal,GlobalOrdinal,Node> { 00180 public: 00182 00183 00185 typedef Scalar scalar_type; 00187 typedef LocalOrdinal local_ordinal_type; 00189 typedef GlobalOrdinal global_ordinal_type; 00191 typedef Node node_type; 00192 00194 typedef Map<LocalOrdinal, GlobalOrdinal, node_type> map_type; 00195 00197 typedef Import<LocalOrdinal, GlobalOrdinal, node_type> import_type; 00198 00200 typedef Export<LocalOrdinal, GlobalOrdinal, node_type> export_type; 00201 00203 typedef CrsGraph<LocalOrdinal, GlobalOrdinal, node_type> crs_graph_type; 00204 00206 00207 00208 00226 CrsMatrix (const Teuchos::RCP<const map_type>& rowMap, 00227 size_t maxNumEntriesPerRow, 00228 ProfileType pftype = DynamicProfile, 00229 const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null); 00230 00248 CrsMatrix (const Teuchos::RCP<const map_type>& rowMap, 00249 const ArrayRCP<const size_t>& NumEntriesPerRowToAlloc, 00250 ProfileType pftype = DynamicProfile, 00251 const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null); 00252 00275 CrsMatrix (const Teuchos::RCP<const map_type>& rowMap, 00276 const Teuchos::RCP<const map_type>& colMap, 00277 size_t maxNumEntriesPerRow, 00278 ProfileType pftype = DynamicProfile, 00279 const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null); 00280 00303 CrsMatrix (const Teuchos::RCP<const map_type>& rowMap, 00304 const Teuchos::RCP<const map_type>& colMap, 00305 const ArrayRCP<const size_t>& NumEntriesPerRowToAlloc, 00306 ProfileType pftype = DynamicProfile, 00307 const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null); 00308 00333 explicit CrsMatrix (const Teuchos::RCP<const crs_graph_type>& graph, 00334 const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null); 00335 00359 CrsMatrix (const RCP<const map_type>& rowMap, 00360 const RCP<const map_type>& colMap, 00361 const ArrayRCP<size_t>& rowPointers, 00362 const ArrayRCP<LocalOrdinal>& columnIndices, 00363 const ArrayRCP<Scalar>& values, 00364 const RCP<ParameterList>& params = null); 00365 00366 // This friend declaration makes the clone() method work. 00367 template <class S2, class LO2, class GO2, class N2> 00368 friend class CrsMatrix; 00369 00394 template <class Node2> 00395 Teuchos::RCP<CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node2> > 00396 clone (const Teuchos::RCP<Node2>& node2, 00397 const Teuchos::RCP<Teuchos::ParameterList>& params = 00398 Teuchos::null) const 00399 { 00400 using Teuchos::arcp; 00401 using Teuchos::Array; 00402 using Teuchos::ArrayRCP; 00403 using Teuchos::ArrayView; 00404 using Teuchos::null; 00405 using Teuchos::ParameterList; 00406 using Teuchos::RCP; 00407 using Teuchos::rcp; 00408 using Teuchos::sublist; 00409 typedef CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node2> CrsMatrix2; 00410 typedef Map<LocalOrdinal, GlobalOrdinal, Node2> Map2; 00411 const char tfecfFuncName[] = "clone"; 00412 00413 // Get parameter values. Set them initially to their default values. 00414 bool fillCompleteClone = true; 00415 bool useLocalIndices = this->hasColMap (); 00416 ProfileType pftype = StaticProfile; 00417 if (! params.is_null ()) { 00418 fillCompleteClone = params->get ("fillComplete clone", fillCompleteClone); 00419 useLocalIndices = params->get ("Locally indexed clone", useLocalIndices); 00420 00421 bool staticProfileClone = true; 00422 staticProfileClone = params->get ("Static profile clone", staticProfileClone); 00423 pftype = staticProfileClone ? StaticProfile : DynamicProfile; 00424 } 00425 00426 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 00427 ! this->hasColMap () && useLocalIndices, std::runtime_error, 00428 ": You requested that the returned clone have local indices, but the " 00429 "the source matrix does not have a column Map yet."); 00430 00431 RCP<const Map2> clonedRowMap = 00432 this->getRowMap ()->template clone<Node2> (node2); 00433 RCP<CrsMatrix2> clonedMatrix; 00434 ArrayRCP<const size_t> numEntries; 00435 size_t numEntriesForAll = 0; 00436 if (! staticGraph_->indicesAreAllocated ()) { 00437 if (! staticGraph_->numAllocPerRow_.is_null ()) { 00438 numEntries = staticGraph_->numAllocPerRow_; 00439 } else { 00440 numEntriesForAll = staticGraph_->numAllocForAllRows_; 00441 } 00442 } else if (! staticGraph_->numRowEntries_.is_null ()) { 00443 numEntries = staticGraph_->numRowEntries_; 00444 } else if (staticGraph_->nodeNumAllocated_ == 0) { 00445 numEntriesForAll = 0; 00446 } else { 00447 // We're left with the case that we have optimized storage. 00448 // In this case, we have to construct a list of row sizes. 00449 TEUCHOS_TEST_FOR_EXCEPTION( 00450 getProfileType() != StaticProfile, std::logic_error, 00451 "Internal logic error. Please report this to Tpetra team." ) 00452 00453 const size_t numRows = this->getNodeNumRows (); 00454 numEntriesForAll = 0; 00455 ArrayRCP<size_t> numEnt; 00456 if (numRows != 0) { 00457 numEnt = arcp<size_t> (numRows); 00458 } 00459 for (size_t i = 0; i < numRows; ++i) { 00460 numEnt[i] = staticGraph_->rowPtrs_[i+1] - staticGraph_->rowPtrs_[i]; 00461 } 00462 numEntries = numEnt; 00463 } 00464 00465 RCP<ParameterList> matrixparams = 00466 params.is_null () ? null : sublist (params,"CrsMatrix"); 00467 if (useLocalIndices) { 00468 RCP<const Map2> clonedColMap = 00469 this->getColMap ()->template clone<Node2> (node2); 00470 if (numEntries.is_null ()) { 00471 clonedMatrix = rcp (new CrsMatrix2 (clonedRowMap, clonedColMap, 00472 numEntriesForAll, pftype, 00473 matrixparams)); 00474 } else { 00475 clonedMatrix = rcp (new CrsMatrix2 (clonedRowMap, clonedColMap, 00476 numEntries, pftype, 00477 matrixparams)); 00478 } 00479 } else { 00480 if (numEntries.is_null ()) { 00481 clonedMatrix = rcp (new CrsMatrix2 (clonedRowMap, numEntriesForAll, 00482 pftype, matrixparams)); 00483 } else { 00484 clonedMatrix = rcp (new CrsMatrix2 (clonedRowMap, numEntries, pftype, 00485 matrixparams)); 00486 } 00487 } 00488 // done with these 00489 numEntries = null; 00490 numEntriesForAll = 0; 00491 00492 if (useLocalIndices) { 00493 clonedMatrix->allocateValues (LocalIndices, 00494 CrsMatrix2::GraphNotYetAllocated); 00495 if (this->isLocallyIndexed ()) { 00496 ArrayView<const LocalOrdinal> linds; 00497 ArrayView<const Scalar> vals; 00498 for (LocalOrdinal lrow = clonedRowMap->getMinLocalIndex (); 00499 lrow <= clonedRowMap->getMaxLocalIndex (); 00500 ++lrow) { 00501 this->getLocalRowView (lrow, linds, vals); 00502 if (linds.size ()) { 00503 clonedMatrix->insertLocalValues (lrow, linds, vals); 00504 } 00505 } 00506 } 00507 else { // this->isGloballyIndexed() 00508 Array<LocalOrdinal> linds; 00509 Array<Scalar> vals; 00510 for (LocalOrdinal lrow = clonedRowMap->getMinLocalIndex (); 00511 lrow <= clonedRowMap->getMaxLocalIndex (); 00512 ++lrow) { 00513 size_t theNumEntries = this->getNumEntriesInLocalRow (lrow); 00514 if (theNumEntries > static_cast<size_t> (linds.size ())) { 00515 linds.resize (theNumEntries); 00516 } 00517 if (theNumEntries > static_cast<size_t> (vals.size ())) { 00518 vals.resize (theNumEntries); 00519 } 00520 this->getLocalRowCopy (clonedRowMap->getGlobalElement (lrow), 00521 linds (), vals (), theNumEntries); 00522 if (theNumEntries != 0) { 00523 clonedMatrix->insertLocalValues (lrow, linds (0, theNumEntries), 00524 vals (0, theNumEntries)); 00525 } 00526 } 00527 } 00528 } 00529 else { // useGlobalIndices 00530 clonedMatrix->allocateValues (GlobalIndices, 00531 CrsMatrix2::GraphNotYetAllocated); 00532 if (this->isGloballyIndexed ()) { 00533 ArrayView<const GlobalOrdinal> ginds; 00534 ArrayView<const Scalar> vals; 00535 for (GlobalOrdinal grow = clonedRowMap->getMinGlobalIndex (); 00536 grow <= clonedRowMap->getMaxGlobalIndex (); 00537 ++grow) { 00538 this->getGlobalRowView (grow, ginds, vals); 00539 if (ginds.size () > 0) { 00540 clonedMatrix->insertGlobalValues (grow, ginds, vals); 00541 } 00542 } 00543 } 00544 else { // this->isLocallyIndexed() 00545 Array<GlobalOrdinal> ginds; 00546 Array<Scalar> vals; 00547 for (GlobalOrdinal grow = clonedRowMap->getMinGlobalIndex (); 00548 grow <= clonedRowMap->getMaxGlobalIndex (); 00549 ++grow) { 00550 size_t theNumEntries = this->getNumEntriesInGlobalRow (grow); 00551 if (theNumEntries > static_cast<size_t> (ginds.size ())) { 00552 ginds.resize (theNumEntries); 00553 } 00554 if (theNumEntries > static_cast<size_t> (vals.size ())) { 00555 vals.resize (theNumEntries); 00556 } 00557 this->getGlobalRowCopy (grow, ginds (), vals (), theNumEntries); 00558 if (theNumEntries != 0) { 00559 clonedMatrix->insertGlobalValues (grow, ginds (0, theNumEntries), 00560 vals (0, theNumEntries)); 00561 } 00562 } 00563 } 00564 } 00565 00566 if (fillCompleteClone) { 00567 RCP<ParameterList> fillparams = 00568 params.is_null () ? Teuchos::null : sublist (params, "fillComplete"); 00569 try { 00570 RCP<const Map2> clonedRangeMap; 00571 RCP<const Map2> clonedDomainMap; 00572 if (! this->getRangeMap ().is_null () && 00573 this->getRangeMap () != clonedRowMap) { 00574 clonedRangeMap = 00575 this->getRangeMap ()->template clone<Node2> (node2); 00576 } else { 00577 clonedRangeMap = clonedRowMap; 00578 } 00579 if (! this->getDomainMap ().is_null () && 00580 this->getDomainMap () != clonedRowMap) { 00581 clonedDomainMap = 00582 this->getDomainMap ()->template clone<Node2> (node2); 00583 } else { 00584 clonedDomainMap = clonedRowMap; 00585 } 00586 clonedMatrix->fillComplete (clonedDomainMap, clonedRangeMap, 00587 fillparams); 00588 } 00589 catch (std::exception &e) { 00590 const bool caughtExceptionOnClone = true; 00591 TEUCHOS_TEST_FOR_EXCEPTION( 00592 caughtExceptionOnClone, std::runtime_error, 00593 "Tpetra::CrsMatrix::clone: Caught the following exception while " 00594 "calling fillComplete() on a clone of type " 00595 << Teuchos::typeName (*clonedMatrix) << ": " 00596 << e.what () << std::endl); 00597 } 00598 } 00599 return clonedMatrix; 00600 } 00601 00603 virtual ~CrsMatrix (); 00604 00606 00607 00608 00631 // 00677 void 00678 insertGlobalValues (const GlobalOrdinal globalRow, 00679 const ArrayView<const GlobalOrdinal>& cols, 00680 const ArrayView<const Scalar>& vals); 00681 00721 void 00722 insertLocalValues (const LocalOrdinal localRow, 00723 const ArrayView<const LocalOrdinal> &cols, 00724 const ArrayView<const Scalar> &vals); 00725 00760 LocalOrdinal 00761 replaceGlobalValues (GlobalOrdinal globalRow, 00762 const ArrayView<const GlobalOrdinal>& cols, 00763 const ArrayView<const Scalar>& vals); 00764 00795 LocalOrdinal 00796 replaceLocalValues (const LocalOrdinal localRow, 00797 const ArrayView<const LocalOrdinal>&cols, 00798 const ArrayView<const Scalar>& vals); 00799 00832 LocalOrdinal 00833 sumIntoGlobalValues (const GlobalOrdinal globalRow, 00834 const ArrayView<const GlobalOrdinal> &cols, 00835 const ArrayView<const Scalar> &vals); 00836 00851 LocalOrdinal 00852 sumIntoLocalValues (const LocalOrdinal localRow, 00853 const ArrayView<const LocalOrdinal>& cols, 00854 const ArrayView<const Scalar>& vals); 00855 00857 void setAllToScalar (const Scalar &alpha); 00858 00860 void scale (const Scalar &alpha); 00861 00863 00870 void 00871 setAllValues (const ArrayRCP<size_t>& rowPointers, 00872 const ArrayRCP<LocalOrdinal>& columnIndices, 00873 const ArrayRCP<Scalar>& values); 00874 00876 00883 void 00884 getAllValues (ArrayRCP<const size_t>& rowPointers, 00885 ArrayRCP<const LocalOrdinal>& columnIndices, 00886 ArrayRCP<const Scalar>& values) const; 00887 00889 00890 00891 00920 void globalAssemble(); 00921 00935 void resumeFill (const RCP<ParameterList>& params = null); 00936 00970 void 00971 fillComplete (const RCP<const map_type>& domainMap, 00972 const RCP<const map_type>& rangeMap, 00973 const RCP<ParameterList>& params = null); 00974 00987 void fillComplete (const RCP<ParameterList>& params = null); 00988 01000 void 01001 expertStaticFillComplete (const RCP<const map_type>& domainMap, 01002 const RCP<const map_type>& rangeMap, 01003 const RCP<const import_type>& importer = Teuchos::null, 01004 const RCP<const export_type>& exporter = Teuchos::null, 01005 const RCP<ParameterList>& params = Teuchos::null); 01006 01016 void 01017 replaceColMap (const Teuchos::RCP<const map_type>& newColMap); 01018 01100 void 01101 reindexColumns (crs_graph_type* const graph, 01102 const Teuchos::RCP<const map_type>& newColMap, 01103 const Teuchos::RCP<const import_type>& newImport = Teuchos::null, 01104 const bool sortEachRow = true); 01105 01118 void 01119 replaceDomainMapAndImporter (const Teuchos::RCP<const map_type>& newDomainMap, 01120 Teuchos::RCP<const import_type>& newImporter); 01121 01135 virtual void 01136 removeEmptyProcessesInPlace (const Teuchos::RCP<const map_type>& newMap); 01137 01139 01140 01141 01143 RCP<const Comm<int> > getComm() const; 01144 01146 RCP<node_type> getNode () const; 01147 01149 RCP<const map_type> getRowMap () const; 01150 01152 RCP<const map_type> getColMap () const; 01153 01155 RCP<const RowGraph<LocalOrdinal,GlobalOrdinal,node_type> > getGraph () const; 01156 01158 RCP<const crs_graph_type> getCrsGraph () const; 01159 01179 global_size_t getGlobalNumRows() const; 01180 01186 global_size_t getGlobalNumCols() const; 01187 01194 size_t getNodeNumRows() const; 01195 01199 size_t getNodeNumCols() const; 01200 01202 GlobalOrdinal getIndexBase() const; 01203 01205 global_size_t getGlobalNumEntries() const; 01206 01208 size_t getNodeNumEntries() const; 01209 01211 01212 size_t getNumEntriesInGlobalRow(GlobalOrdinal globalRow) const; 01213 01215 01216 size_t getNumEntriesInLocalRow(LocalOrdinal localRow) const; 01217 01219 01221 global_size_t getGlobalNumDiags() const; 01222 01224 01226 size_t getNodeNumDiags() const; 01227 01229 01231 size_t getGlobalMaxNumRowEntries() const; 01232 01234 01236 size_t getNodeMaxNumRowEntries() const; 01237 01239 bool hasColMap() const; 01240 01242 01244 bool isLowerTriangular() const; 01245 01247 01249 bool isUpperTriangular() const; 01250 01270 bool isLocallyIndexed() const; 01271 01291 bool isGloballyIndexed() const; 01292 01315 bool isFillComplete() const; 01316 01339 bool isFillActive() const; 01340 01342 01348 bool isStorageOptimized() const; 01349 01351 ProfileType getProfileType() const; 01352 01354 bool isStaticGraph() const; 01355 01367 typename ScalarTraits<Scalar>::magnitudeType getFrobeniusNorm() const; 01368 01375 virtual bool supportsRowViews() const; 01376 01378 01388 void 01389 getGlobalRowCopy (GlobalOrdinal GlobalRow, 01390 const ArrayView<GlobalOrdinal> &Indices, 01391 const ArrayView<Scalar> &Values, 01392 size_t &NumEntries) const; 01393 01395 01407 void 01408 getLocalRowCopy (LocalOrdinal LocalRow, 01409 const ArrayView<LocalOrdinal> &Indices, 01410 const ArrayView<Scalar> &Values, 01411 size_t &NumEntries) const; 01412 01426 void 01427 getGlobalRowView (GlobalOrdinal GlobalRow, 01428 ArrayView<const GlobalOrdinal> &indices, 01429 ArrayView<const Scalar> &values) const; 01430 01444 void 01445 getLocalRowView (LocalOrdinal LocalRow, 01446 ArrayView<const LocalOrdinal>& indices, 01447 ArrayView<const Scalar>& values) const; 01448 01454 void getLocalDiagCopy (Vector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>& diag) const; 01455 01493 void getLocalDiagOffsets (Teuchos::ArrayRCP<size_t>& offsets) const; 01494 01513 void 01514 getLocalDiagCopy (Vector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>& diag, 01515 const Teuchos::ArrayView<const size_t>& offsets) const; 01516 01518 void leftScale(const Vector<Scalar, LocalOrdinal, GlobalOrdinal, node_type>& x); 01519 01521 void rightScale(const Vector<Scalar, LocalOrdinal, GlobalOrdinal, node_type>& x); 01522 01524 01525 01526 01575 template <class DomainScalar, class RangeScalar> 01576 void 01577 localMultiply (const MultiVector<DomainScalar,LocalOrdinal,GlobalOrdinal,node_type>& X, 01578 MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,node_type>& Y, 01579 Teuchos::ETransp trans, 01580 RangeScalar alpha, 01581 RangeScalar beta) const; 01582 01607 template <class DomainScalar, class RangeScalar> 01608 void 01609 localGaussSeidel (const MultiVector<DomainScalar,LocalOrdinal,GlobalOrdinal,node_type> &B, 01610 MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,node_type> &X, 01611 const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type> &D, 01612 const RangeScalar& dampingFactor, 01613 const KokkosClassic::ESweepDirection direction) const; 01614 01641 template <class DomainScalar, class RangeScalar> 01642 void 01643 reorderedLocalGaussSeidel (const MultiVector<DomainScalar,LocalOrdinal,GlobalOrdinal,node_type>& B, 01644 MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,node_type>& X, 01645 const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>& D, 01646 const ArrayView<LocalOrdinal>& rowIndices, 01647 const RangeScalar& dampingFactor, 01648 const KokkosClassic::ESweepDirection direction) const; 01649 01666 template <class DomainScalar, class RangeScalar> 01667 void 01668 localSolve (const MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,node_type>& Y, 01669 MultiVector<DomainScalar,LocalOrdinal,GlobalOrdinal,node_type>& X, 01670 Teuchos::ETransp trans) const; 01671 01673 template <class T> 01674 RCP<CrsMatrix<T,LocalOrdinal,GlobalOrdinal,node_type> > convert() const; 01675 01677 01678 01679 01690 void 01691 apply (const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>& X, 01692 MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>&Y, 01693 Teuchos::ETransp mode = Teuchos::NO_TRANS, 01694 Scalar alpha = ScalarTraits<Scalar>::one(), 01695 Scalar beta = ScalarTraits<Scalar>::zero()) const; 01696 01698 bool hasTransposeApply() const; 01699 01703 RCP<const map_type> getDomainMap () const; 01704 01708 RCP<const map_type> getRangeMap () const; 01709 01711 01712 01713 01778 void 01779 gaussSeidel (const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type> &B, 01780 MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type> &X, 01781 const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type> &D, 01782 const Scalar& dampingFactor, 01783 const ESweepDirection direction, 01784 const int numSweeps) const; 01785 01852 void 01853 reorderedGaussSeidel (const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>& B, 01854 MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>& X, 01855 const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>& D, 01856 const ArrayView<LocalOrdinal>& rowIndices, 01857 const Scalar& dampingFactor, 01858 const ESweepDirection direction, 01859 const int numSweeps) const; 01860 01889 void 01890 gaussSeidelCopy (MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type> &X, 01891 const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type> &B, 01892 const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type> &D, 01893 const Scalar& dampingFactor, 01894 const ESweepDirection direction, 01895 const int numSweeps, 01896 const bool zeroInitialGuess) const; 01897 01927 void 01928 reorderedGaussSeidelCopy (MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>& X, 01929 const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>& B, 01930 const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>& D, 01931 const ArrayView<LocalOrdinal>& rowIndices, 01932 const Scalar& dampingFactor, 01933 const ESweepDirection direction, 01934 const int numSweeps, 01935 const bool zeroInitialGuess) const; 01936 01947 virtual Teuchos::RCP<RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, node_type> > 01948 add (const Scalar& alpha, 01949 const RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, node_type>& A, 01950 const Scalar& beta, 01951 const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, node_type> >& domainMap, 01952 const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, node_type> >& rangeMap, 01953 const Teuchos::RCP<Teuchos::ParameterList>& params) const; 01954 01956 01957 01958 01960 std::string description() const; 01961 01963 void 01964 describe (Teuchos::FancyOStream &out, 01965 const Teuchos::EVerbosityLevel verbLevel = 01966 Teuchos::Describable::verbLevel_default) const; 01967 01969 01970 01971 01972 virtual bool 01973 checkSizes (const SrcDistObject& source); 01974 01975 virtual void 01976 copyAndPermute (const SrcDistObject& source, 01977 size_t numSameIDs, 01978 const ArrayView<const LocalOrdinal> &permuteToLIDs, 01979 const ArrayView<const LocalOrdinal> &permuteFromLIDs); 01980 01981 virtual void 01982 packAndPrepare (const SrcDistObject& source, 01983 const Teuchos::ArrayView<const LocalOrdinal>& exportLIDs, 01984 Teuchos::Array<char>& exports, 01985 const Teuchos::ArrayView<size_t>& numPacketsPerLID, 01986 size_t& constantNumPackets, 01987 Distributor& distor); 01988 01998 void 01999 unpackAndCombine (const Teuchos::ArrayView<const LocalOrdinal> &importLIDs, 02000 const Teuchos::ArrayView<const char> &imports, 02001 const Teuchos::ArrayView<size_t> &numPacketsPerLID, 02002 size_t constantNumPackets, 02003 Distributor &distor, 02004 CombineMode combineMode); 02006 02007 02008 02017 virtual void 02018 pack (const Teuchos::ArrayView<const LocalOrdinal>& exportLIDs, 02019 Teuchos::Array<char>& exports, 02020 const Teuchos::ArrayView<size_t>& numPacketsPerLID, 02021 size_t& constantNumPackets, 02022 Distributor& distor) const; 02023 02025 private: 02026 // Friend declaration for nonmember function. 02027 template<class CrsMatrixType> 02028 friend Teuchos::RCP<CrsMatrixType> 02029 importAndFillCompleteCrsMatrix (const Teuchos::RCP<const CrsMatrixType>& sourceMatrix, 02030 const Import<typename CrsMatrixType::local_ordinal_type, 02031 typename CrsMatrixType::global_ordinal_type, 02032 typename CrsMatrixType::node_type>& importer, 02033 const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type, 02034 typename CrsMatrixType::global_ordinal_type, 02035 typename CrsMatrixType::node_type> >& domainMap, 02036 const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type, 02037 typename CrsMatrixType::global_ordinal_type, 02038 typename CrsMatrixType::node_type> >& rangeMap, 02039 const Teuchos::RCP<Teuchos::ParameterList>& params); 02040 02041 // Friend declaration for nonmember function. 02042 template<class CrsMatrixType> 02043 friend Teuchos::RCP<CrsMatrixType> 02044 exportAndFillCompleteCrsMatrix (const Teuchos::RCP<const CrsMatrixType>& sourceMatrix, 02045 const Export<typename CrsMatrixType::local_ordinal_type, 02046 typename CrsMatrixType::global_ordinal_type, 02047 typename CrsMatrixType::node_type>& exporter, 02048 const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type, 02049 typename CrsMatrixType::global_ordinal_type, 02050 typename CrsMatrixType::node_type> >& domainMap, 02051 const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type, 02052 typename CrsMatrixType::global_ordinal_type, 02053 typename CrsMatrixType::node_type> >& rangeMap, 02054 const Teuchos::RCP<Teuchos::ParameterList>& params); 02055 02056 public: 02072 void 02073 importAndFillComplete (Teuchos::RCP<CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,node_type> > & destMatrix, 02074 const import_type& importer, 02075 const Teuchos::RCP<const map_type>& domainMap = Teuchos::null, 02076 const Teuchos::RCP<const map_type>& rangeMap = Teuchos::null, 02077 const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null) const; 02078 02094 void 02095 exportAndFillComplete (Teuchos::RCP<CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,node_type> > & destMatrix, 02096 const export_type& exporter, 02097 const Teuchos::RCP<const map_type>& domainMap = Teuchos::null, 02098 const Teuchos::RCP<const map_type>& rangeMap = Teuchos::null, 02099 const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null) const; 02100 private: 02121 void 02122 transferAndFillComplete (Teuchos::RCP<CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > & destMat, 02123 const ::Tpetra::Details::Transfer<LocalOrdinal, GlobalOrdinal, Node>& rowTransfer, // either Import or Export 02124 const Teuchos::RCP<const map_type>& domainMap = Teuchos::null, 02125 const Teuchos::RCP<const map_type>& rangeMap = Teuchos::null, 02126 const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null) const; 02127 02128 // We forbid copy construction by declaring this method private 02129 // and not implementing it. 02130 CrsMatrix (const CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,node_type> &rhs); 02131 02132 // We forbid assignment (operator=) by declaring this method 02133 // private and not implementing it. 02134 CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,node_type>& 02135 operator= (const CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,node_type> &rhs); 02136 02142 void 02143 insertGlobalValuesFiltered (const GlobalOrdinal globalRow, 02144 const ArrayView<const GlobalOrdinal> &indices, 02145 const ArrayView<const Scalar> &values); 02146 02152 void 02153 insertLocalValuesFiltered (const LocalOrdinal localRow, 02154 const ArrayView<const LocalOrdinal> &indices, 02155 const ArrayView<const Scalar> &values); 02156 02164 void 02165 combineGlobalValues (const GlobalOrdinal globalRowIndex, 02166 const Teuchos::ArrayView<const GlobalOrdinal> columnIndices, 02167 const Teuchos::ArrayView<const Scalar> values, 02168 const Tpetra::CombineMode combineMode); 02169 02199 template<class BinaryFunction> 02200 LocalOrdinal 02201 transformGlobalValues (GlobalOrdinal globalRow, 02202 const Teuchos::ArrayView<const GlobalOrdinal>& indices, 02203 const Teuchos::ArrayView<const Scalar> & values, 02204 BinaryFunction f) 02205 { 02206 typedef LocalOrdinal LO; 02207 typedef GlobalOrdinal GO; 02208 using Teuchos::Array; 02209 using Teuchos::ArrayView; 02210 typedef typename ArrayView<const GO>::size_type size_type; 02211 02212 if (! isFillActive ()) { 02213 // Fill must be active in order to call this method. 02214 return Teuchos::OrdinalTraits<LO>::invalid (); 02215 } 02216 else if (values.size () != indices.size ()) { 02217 // The sizes of values and indices must match. 02218 return Teuchos::OrdinalTraits<LO>::invalid (); 02219 } 02220 02221 const LO lrow = this->getRowMap()->getLocalElement (globalRow); 02222 if (lrow == Teuchos::OrdinalTraits<LO>::invalid ()) { 02223 // We don't own the row, so we're not allowed to modify its values. 02224 return Teuchos::OrdinalTraits<LO>::invalid (); 02225 } 02226 02227 if (staticGraph_.is_null ()) { 02228 return Teuchos::OrdinalTraits<LO>::invalid (); 02229 } 02230 const crs_graph_type& graph = *staticGraph_; 02231 RowInfo rowInfo = graph.getRowInfo (lrow); 02232 if (indices.size () == 0) { 02233 return static_cast<LO> (0); 02234 } 02235 else { 02236 ArrayView<Scalar> curVals = this->getViewNonConst (rowInfo); 02237 if (isLocallyIndexed ()) { 02238 // Convert the given global indices to local indices. 02239 // 02240 // FIXME (mfh 08 Jul 2014) Why can't we ask the graph to do 02241 // that? It could do the conversions in place, so that we 02242 // wouldn't need temporary storage. 02243 const map_type& colMap = * (this->getColMap ()); 02244 const size_type numInds = indices.size (); 02245 Array<LO> lclInds (numInds); 02246 for (size_type k = 0; k < numInds; ++k) { 02247 // There is no need to filter out indices not in the 02248 // column Map. Those that aren't will be mapped to 02249 // invalid(), which the graph's transformGlobalValues() 02250 // will filter out (but not count in its return value). 02251 lclInds[k] = colMap.getLocalElement (indices[k]); 02252 } 02253 return graph.template transformLocalValues<Scalar, BinaryFunction> (rowInfo, curVals, 02254 lclInds (), values, f); 02255 } 02256 else if (isGloballyIndexed ()) { 02257 return graph.template transformGlobalValues<Scalar, BinaryFunction> (rowInfo, curVals, 02258 indices, values, f); 02259 } 02260 else { 02261 // If the graph is neither locally nor globally indexed on 02262 // the calling process, that means that the calling process 02263 // can't possibly have any entries in the owned row. Thus, 02264 // there are no entries to transform, so we return zero. 02265 return static_cast<LO> (0); 02266 } 02267 } 02268 } 02269 02270 private: 02273 void 02274 insertNonownedGlobalValues (const GlobalOrdinal globalRow, 02275 const Teuchos::ArrayView<const GlobalOrdinal>& indices, 02276 const Teuchos::ArrayView<const Scalar>& values); 02277 02278 protected: 02279 // useful typedefs 02280 typedef OrdinalTraits<LocalOrdinal> OTL; 02281 typedef ScalarTraits<Scalar> STS; 02282 typedef typename STS::magnitudeType Magnitude; 02283 typedef ScalarTraits<Magnitude> STM; 02284 typedef MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type> MV; 02285 typedef Vector<Scalar,LocalOrdinal,GlobalOrdinal,node_type> V; 02286 // This typedef stays only for the sake of backwards 02287 // compatibility. It does not follow the standard format for 02288 // typedefs in a class, nor is it a nonambiguous abbreviation (is 02289 // it RowGraph or CrsGraph?). 02290 typedef crs_graph_type Graph; 02291 02292 // mfh 30 Aug 2014: CrsMatrix used this bind_scalar business 02293 // because Chris Baker expected that some implementations of 02294 // SparseOps might not work for all Scalar types. (cuSPARSE is an 02295 // example, since NVIDIA only wrote sparse matrix-vector multiply 02296 // routines for float and double.) In cases like that, the 02297 // "other_type" typedef would point to some other, "generic" 02298 // implementation of SparseOps. 02299 // 02300 // The Kokkos refactor version of Tpetra will hide this complexity 02301 // at the KokkosLinAlg level, so that users won't have to look at 02302 // type selectors quite so much. 02303 typedef typename KokkosClassic::DefaultKernels<Scalar, LocalOrdinal, node_type>::SparseOps source_sparse_ops_type; 02304 typedef typename source_sparse_ops_type::template bind_scalar<Scalar>::other_type sparse_ops_type; 02305 typedef typename sparse_ops_type::template graph<LocalOrdinal, node_type>::graph_type local_graph_type; 02306 typedef typename sparse_ops_type::template matrix<Scalar, LocalOrdinal, node_type>::matrix_type local_matrix_type; 02307 02308 // Enums 02309 enum GraphAllocationStatus { 02310 GraphAlreadyAllocated, 02311 GraphNotYetAllocated 02312 }; 02313 02330 void allocateValues (ELocalGlobal lg, GraphAllocationStatus gas); 02331 02337 void sortEntries(); 02338 02344 void mergeRedundantEntries(); 02345 02353 void clearGlobalConstants(); 02354 02363 void computeGlobalConstants(); 02364 02377 mutable Teuchos::RCP<MV> importMV_; 02378 02391 mutable Teuchos::RCP<MV> exportMV_; 02392 02412 Teuchos::RCP<MV> 02413 getColumnMapMultiVector (const MV& X_domainMap, 02414 const bool force = false) const; 02415 02437 Teuchos::RCP<MV> 02438 getRowMapMultiVector (const MV& Y_rangeMap, 02439 const bool force = false) const; 02440 02442 void 02443 applyNonTranspose (const MV& X_in, 02444 MV& Y_in, 02445 Scalar alpha, 02446 Scalar beta) const; 02447 02449 void 02450 applyTranspose (const MV& X_in, 02451 MV& Y_in, 02452 const Teuchos::ETransp mode, 02453 Scalar alpha, 02454 Scalar beta) const; 02455 02456 // matrix data accessors 02457 ArrayView<const Scalar> getView(RowInfo rowinfo) const; 02458 ArrayView< Scalar> getViewNonConst(RowInfo rowinfo); 02459 // local Kokkos objects 02460 02466 void fillLocalMatrix (const Teuchos::RCP<Teuchos::ParameterList>& params); 02467 02473 void fillLocalGraphAndMatrix (const Teuchos::RCP<Teuchos::ParameterList>& params); 02474 02476 void checkInternalState() const; 02477 02489 02490 Teuchos::RCP<const crs_graph_type> staticGraph_; 02491 Teuchos::RCP<crs_graph_type> myGraph_; 02493 02498 Teuchos::RCP<sparse_ops_type> lclMatOps_; 02505 Teuchos::RCP<local_matrix_type> lclMatrix_; 02506 02519 02520 ArrayRCP<Scalar> values1D_; 02521 ArrayRCP<Array<Scalar> > values2D_; 02523 02525 bool fillComplete_; 02526 02554 std::map<GlobalOrdinal, std::pair<Teuchos::Array<GlobalOrdinal>, 02555 Teuchos::Array<Scalar> > > nonlocals_; 02556 02563 mutable Magnitude frobNorm_; 02564 }; // class CrsMatrix 02565 02572 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 02573 Teuchos::RCP<CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > 02574 createCrsMatrix (const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &map, 02575 size_t maxNumEntriesPerRow = 0, 02576 const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null) 02577 { 02578 using Teuchos::rcp; 02579 typedef CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> matrix_type; 02580 return rcp (new matrix_type (map, maxNumEntriesPerRow, DynamicProfile, params)); 02581 } 02582 02634 template<class CrsMatrixType> 02635 Teuchos::RCP<CrsMatrixType> 02636 importAndFillCompleteCrsMatrix (const Teuchos::RCP<const CrsMatrixType>& sourceMatrix, 02637 const Import<typename CrsMatrixType::local_ordinal_type, 02638 typename CrsMatrixType::global_ordinal_type, 02639 typename CrsMatrixType::node_type>& importer, 02640 const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type, 02641 typename CrsMatrixType::global_ordinal_type, 02642 typename CrsMatrixType::node_type> >& domainMap = Teuchos::null, 02643 const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type, 02644 typename CrsMatrixType::global_ordinal_type, 02645 typename CrsMatrixType::node_type> >& rangeMap = Teuchos::null, 02646 const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null) 02647 { 02648 Teuchos::RCP<CrsMatrixType> destMatrix; 02649 sourceMatrix->importAndFillComplete (destMatrix,importer, domainMap, rangeMap, params); 02650 return destMatrix; 02651 } 02652 02686 template<class CrsMatrixType> 02687 Teuchos::RCP<CrsMatrixType> 02688 exportAndFillCompleteCrsMatrix (const Teuchos::RCP<const CrsMatrixType>& sourceMatrix, 02689 const Export<typename CrsMatrixType::local_ordinal_type, 02690 typename CrsMatrixType::global_ordinal_type, 02691 typename CrsMatrixType::node_type>& exporter, 02692 const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type, 02693 typename CrsMatrixType::global_ordinal_type, 02694 typename CrsMatrixType::node_type> >& domainMap = Teuchos::null, 02695 const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type, 02696 typename CrsMatrixType::global_ordinal_type, 02697 typename CrsMatrixType::node_type> >& rangeMap = Teuchos::null, 02698 const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null) 02699 { 02700 Teuchos::RCP<CrsMatrixType> destMatrix; 02701 sourceMatrix->exportAndFillComplete (destMatrix,exporter, domainMap, rangeMap, params); 02702 return destMatrix; 02703 } 02704 } // namespace Tpetra 02705 02706 // Include KokkosRefactor partial specialisation if enabled 02707 #if defined(TPETRA_HAVE_KOKKOS_REFACTOR) 02708 #include "Tpetra_KokkosRefactor_CrsMatrix_decl.hpp" 02709 #endif 02710 02716 #endif // TPETRA_CRSMATRIX_DECL_HPP
1.7.6.1
