|
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_KOKKOSREFACTOR_CRSMATRIX_DECL_HPP 00043 #define TPETRA_KOKKOSREFACTOR_CRSMATRIX_DECL_HPP 00044 00045 // This file gets included by Tpetra_CrsMatrix_decl.hpp, 00046 // so it inherits all of that file's includes. 00047 00048 #include <KokkosCompat_ClassicNodeAPI_Wrapper.hpp> 00049 #include <Kokkos_CrsMatrix.hpp> 00050 00051 00052 namespace Tpetra { 00053 00059 template <class Scalar, 00060 class LocalOrdinal, 00061 class GlobalOrdinal, 00062 class DeviceType> 00063 class CrsMatrix<Scalar, 00064 LocalOrdinal, 00065 GlobalOrdinal, 00066 Kokkos::Compat::KokkosDeviceWrapperNode<DeviceType> > : 00067 public RowMatrix<Scalar, 00068 LocalOrdinal, 00069 GlobalOrdinal, 00070 Kokkos::Compat::KokkosDeviceWrapperNode<DeviceType> >, 00071 public DistObject<char, 00072 LocalOrdinal, 00073 GlobalOrdinal, 00074 Kokkos::Compat::KokkosDeviceWrapperNode<DeviceType> > 00075 { 00076 public: 00078 00079 00081 typedef Scalar scalar_type; 00083 typedef LocalOrdinal local_ordinal_type; 00085 typedef GlobalOrdinal global_ordinal_type; 00087 typedef DeviceType device_type; 00093 typedef Kokkos::Compat::KokkosDeviceWrapperNode<DeviceType> node_type; 00094 00096 typedef Map<LocalOrdinal, GlobalOrdinal, node_type> map_type; 00097 00099 typedef Import<LocalOrdinal, GlobalOrdinal, node_type> import_type; 00100 00102 typedef Export<LocalOrdinal, GlobalOrdinal, node_type> export_type; 00103 00105 typedef CrsGraph<LocalOrdinal, GlobalOrdinal, node_type> crs_graph_type; 00106 00107 typedef typename crs_graph_type::t_RowPtrs t_RowPtrs; 00108 typedef typename crs_graph_type::t_LocalOrdinal_1D t_LocalOrdinal_1D; 00109 typedef Kokkos::View<Scalar*, device_type> t_ValuesType; 00110 typedef Kokkos::CrsMatrix<Scalar, LocalOrdinal, device_type, void, size_t> k_local_matrix_type; 00111 00113 00114 00115 00133 CrsMatrix (const Teuchos::RCP<const map_type>& rowMap, 00134 size_t maxNumEntriesPerRow, 00135 ProfileType pftype = DynamicProfile, 00136 const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null); 00137 00155 CrsMatrix (const Teuchos::RCP<const map_type>& rowMap, 00156 const ArrayRCP<const size_t>& NumEntriesPerRowToAlloc, 00157 ProfileType pftype = DynamicProfile, 00158 const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null); 00159 00182 CrsMatrix (const Teuchos::RCP<const map_type>& rowMap, 00183 const Teuchos::RCP<const map_type>& colMap, 00184 size_t maxNumEntriesPerRow, 00185 ProfileType pftype = DynamicProfile, 00186 const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null); 00187 00210 CrsMatrix (const Teuchos::RCP<const map_type>& rowMap, 00211 const Teuchos::RCP<const map_type>& colMap, 00212 const ArrayRCP<const size_t>& NumEntriesPerRowToAlloc, 00213 ProfileType pftype = DynamicProfile, 00214 const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null); 00215 00240 explicit CrsMatrix (const Teuchos::RCP<const crs_graph_type>& graph, 00241 const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null); 00242 00266 CrsMatrix (const RCP<const map_type>& rowMap, 00267 const RCP<const map_type>& colMap, 00268 const t_RowPtrs & rowPointers, 00269 const t_LocalOrdinal_1D & columnIndices, 00270 const t_ValuesType & values, 00271 const RCP<ParameterList>& params = null); 00272 00295 CrsMatrix (const RCP<const map_type>& rowMap, 00296 const RCP<const map_type>& colMap, 00297 const ArrayRCP<size_t>& rowPointers, 00298 const ArrayRCP<LocalOrdinal>& columnIndices, 00299 const ArrayRCP<Scalar>& values, 00300 const RCP<ParameterList>& params = null); 00301 00322 CrsMatrix (const RCP<const map_type>& rowMap, 00323 const RCP<const map_type>& colMap, 00324 const k_local_matrix_type& lclMatrix, 00325 const RCP<Teuchos::ParameterList>& params = null); 00326 00327 00328 // This friend declaration makes the clone() method work. 00329 template <class S2, class LO2, class GO2, class N2> 00330 friend class CrsMatrix; 00331 00356 template <class Node2> 00357 Teuchos::RCP<CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node2> > 00358 clone (const Teuchos::RCP<Node2>& node2, 00359 const Teuchos::RCP<Teuchos::ParameterList>& params = null) const 00360 { 00361 using Teuchos::ArrayRCP; 00362 using Teuchos::null; 00363 using Teuchos::ParameterList; 00364 using Teuchos::RCP; 00365 using Teuchos::rcp; 00366 using Teuchos::sublist; 00367 typedef CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node2> CrsMatrix2; 00368 typedef Map<LocalOrdinal, GlobalOrdinal, Node2> Map2; 00369 const char tfecfFuncName[] = "clone"; 00370 00371 // Get parameter values. Set them initially to their default values. 00372 bool fillCompleteClone = true; 00373 bool useLocalIndices = this->hasColMap (); 00374 ProfileType pftype = StaticProfile; 00375 if (! params.is_null ()) { 00376 fillCompleteClone = params->get ("fillComplete clone", fillCompleteClone); 00377 useLocalIndices = params->get ("Locally indexed clone", useLocalIndices); 00378 00379 bool staticProfileClone = true; 00380 staticProfileClone = params->get ("Static profile clone", staticProfileClone); 00381 pftype = staticProfileClone ? StaticProfile : DynamicProfile; 00382 } 00383 00384 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 00385 ! this->hasColMap () && useLocalIndices, std::runtime_error, 00386 ": You requested that the returned clone have local indices, but the " 00387 "the source matrix does not have a column Map yet."); 00388 00389 RCP<const Map2> clonedRowMap = this->getRowMap ()->template clone<Node2> (node2); 00390 00391 // Get an upper bound on the number of entries per row. 00392 RCP<CrsMatrix2> clonedMatrix; 00393 ArrayRCP<const size_t> numEntriesPerRow; 00394 size_t numEntriesForAll = 0; 00395 bool boundSameForAllLocalRows = false; 00396 staticGraph_->getNumEntriesPerLocalRowUpperBound (numEntriesPerRow, 00397 numEntriesForAll, 00398 boundSameForAllLocalRows); 00399 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 00400 numEntriesForAll != 0 && 00401 static_cast<size_t> (numEntriesPerRow.size ()) != 0, 00402 std::logic_error, ": getNumEntriesPerLocalRowUpperBound returned a " 00403 "nonzero numEntriesForAll = " << numEntriesForAll << " , as well as a " 00404 "numEntriesPerRow array of nonzero length " << numEntriesPerRow.size () 00405 << ". This should never happen. Please report this bug to the Tpetra " 00406 "developers."); 00407 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 00408 numEntriesForAll != 0 && ! boundSameForAllLocalRows, 00409 std::logic_error, ": getNumEntriesPerLocalRowUpperBound returned a " 00410 "nonzero numEntriesForAll = " << numEntriesForAll << " , but claims " 00411 "(via its third output value) that the upper bound is not the same for " 00412 "all rows. This should never happen. Please report this bug to the " 00413 "Tpetra developers."); 00414 TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC( 00415 numEntriesPerRow.size () != 0 && boundSameForAllLocalRows, 00416 std::logic_error, ": getNumEntriesPerLocalRowUpperBound returned a " 00417 "numEntriesPerRow array of nonzero length " << numEntriesPerRow.size () 00418 << ", but claims (via its third output value) that the upper bound is " 00419 "not the same for all rows. This should never happen. Please report " 00420 "this bug to the Tpetra developers."); 00421 00422 RCP<ParameterList> matParams = 00423 params.is_null () ? null : sublist (params,"CrsMatrix"); 00424 if (useLocalIndices) { 00425 RCP<const Map2> clonedColMap = 00426 this->getColMap ()->template clone<Node2> (node2); 00427 if (numEntriesPerRow.is_null ()) { 00428 clonedMatrix = rcp (new CrsMatrix2 (clonedRowMap, clonedColMap, 00429 numEntriesForAll, pftype, 00430 matParams)); 00431 } 00432 else { 00433 clonedMatrix = rcp (new CrsMatrix2 (clonedRowMap, clonedColMap, 00434 numEntriesPerRow, pftype, 00435 matParams)); 00436 } 00437 } 00438 else { 00439 if (numEntriesPerRow.is_null ()) { 00440 clonedMatrix = rcp (new CrsMatrix2 (clonedRowMap, numEntriesForAll, 00441 pftype, matParams)); 00442 } 00443 else { 00444 clonedMatrix = rcp (new CrsMatrix2 (clonedRowMap, numEntriesPerRow, 00445 pftype, matParams)); 00446 } 00447 } 00448 // done with these 00449 numEntriesPerRow = Teuchos::null; 00450 numEntriesForAll = 0; 00451 00452 if (useLocalIndices) { 00453 clonedMatrix->allocateValues (LocalIndices, 00454 CrsMatrix2::GraphNotYetAllocated); 00455 if (this->isLocallyIndexed ()) { 00456 ArrayView<const LocalOrdinal> linds; 00457 ArrayView<const Scalar> vals; 00458 for (LocalOrdinal lrow = clonedRowMap->getMinLocalIndex (); 00459 lrow <= clonedRowMap->getMaxLocalIndex (); 00460 ++lrow) { 00461 this->getLocalRowView (lrow, linds, vals); 00462 if (linds.size ()) { 00463 clonedMatrix->insertLocalValues (lrow, linds, vals); 00464 } 00465 } 00466 } 00467 else { // this->isGloballyIndexed() 00468 Array<LocalOrdinal> linds; 00469 Array<Scalar> vals; 00470 for (LocalOrdinal lrow = clonedRowMap->getMinLocalIndex (); 00471 lrow <= clonedRowMap->getMaxLocalIndex (); 00472 ++lrow) { 00473 size_t theNumEntries = this->getNumEntriesInLocalRow (lrow); 00474 if (theNumEntries > Teuchos::as<size_t> (linds.size ())) { 00475 linds.resize (theNumEntries); 00476 } 00477 if (theNumEntries > Teuchos::as<size_t> (vals.size ())) { 00478 vals.resize (theNumEntries); 00479 } 00480 this->getLocalRowCopy (clonedRowMap->getGlobalElement (lrow), 00481 linds (), vals (), theNumEntries); 00482 if (theNumEntries != 0) { 00483 clonedMatrix->insertLocalValues (lrow, linds (0, theNumEntries), 00484 vals (0, theNumEntries)); 00485 } 00486 } 00487 } 00488 } 00489 else { // useGlobalIndices 00490 clonedMatrix->allocateValues (GlobalIndices, 00491 CrsMatrix2::GraphNotYetAllocated); 00492 if (this->isGloballyIndexed ()) { 00493 ArrayView<const GlobalOrdinal> ginds; 00494 ArrayView<const Scalar> vals; 00495 for (GlobalOrdinal grow = clonedRowMap->getMinGlobalIndex (); 00496 grow <= clonedRowMap->getMaxGlobalIndex (); 00497 ++grow) { 00498 this->getGlobalRowView (grow, ginds, vals); 00499 if (ginds.size () > 0) { 00500 clonedMatrix->insertGlobalValues (grow, ginds, vals); 00501 } 00502 } 00503 } 00504 else { // this->isLocallyIndexed() 00505 Array<GlobalOrdinal> ginds; 00506 Array<Scalar> vals; 00507 for (GlobalOrdinal grow = clonedRowMap->getMinGlobalIndex (); 00508 grow <= clonedRowMap->getMaxGlobalIndex (); 00509 ++grow) { 00510 size_t theNumEntries = this->getNumEntriesInGlobalRow (grow); 00511 if (theNumEntries > Teuchos::as<size_t> (ginds.size ())) { 00512 ginds.resize (theNumEntries); 00513 } 00514 if (theNumEntries > Teuchos::as<size_t> (vals.size ())) { 00515 vals.resize (theNumEntries); 00516 } 00517 this->getGlobalRowCopy (grow, ginds (), vals (), theNumEntries); 00518 if (theNumEntries != 0) { 00519 clonedMatrix->insertGlobalValues (grow, ginds (0, theNumEntries), 00520 vals (0, theNumEntries)); 00521 } 00522 } 00523 } 00524 } 00525 00526 if (fillCompleteClone) { 00527 RCP<ParameterList> fillparams = 00528 params.is_null () ? Teuchos::null : sublist (params, "fillComplete"); 00529 try { 00530 RCP<const Map2> clonedRangeMap; 00531 RCP<const Map2> clonedDomainMap; 00532 if (! this->getRangeMap ().is_null () && 00533 this->getRangeMap () != clonedRowMap) { 00534 clonedRangeMap = this->getRangeMap ()->template clone<Node2> (node2); 00535 } 00536 else { 00537 clonedRangeMap = clonedRowMap; 00538 } 00539 if (! this->getDomainMap ().is_null () && 00540 this->getDomainMap () != clonedRowMap) { 00541 clonedDomainMap = this->getDomainMap ()->template clone<Node2> (node2); 00542 } 00543 else { 00544 clonedDomainMap = clonedRowMap; 00545 } 00546 clonedMatrix->fillComplete (clonedDomainMap, clonedRangeMap, 00547 fillparams); 00548 } 00549 catch (std::exception &e) { 00550 const bool caughtExceptionOnClone = true; 00551 TEUCHOS_TEST_FOR_EXCEPTION( 00552 caughtExceptionOnClone, std::runtime_error, 00553 Teuchos::typeName (*this) << std::endl << "clone: " << std::endl << 00554 "Caught the following exception while calling fillComplete() on a " 00555 "clone of type" << std::endl << Teuchos::typeName (*clonedMatrix) 00556 << ": " << std::endl << e.what () << std::endl); 00557 } 00558 } 00559 return clonedMatrix; 00560 } 00561 00563 virtual ~CrsMatrix (); 00564 00566 00567 00568 00591 // 00637 void 00638 insertGlobalValues (const GlobalOrdinal globalRow, 00639 const ArrayView<const GlobalOrdinal>& cols, 00640 const ArrayView<const Scalar>& vals); 00641 00681 void 00682 insertLocalValues (const LocalOrdinal localRow, 00683 const ArrayView<const LocalOrdinal> &cols, 00684 const ArrayView<const Scalar> &vals); 00685 00720 LocalOrdinal 00721 replaceGlobalValues (GlobalOrdinal globalRow, 00722 const ArrayView<const GlobalOrdinal>& cols, 00723 const ArrayView<const Scalar>& vals); 00724 00755 LocalOrdinal 00756 replaceLocalValues (const LocalOrdinal localRow, 00757 const ArrayView<const LocalOrdinal>& cols, 00758 const ArrayView<const Scalar>& vals); 00759 00792 LocalOrdinal 00793 sumIntoGlobalValues (const GlobalOrdinal globalRow, 00794 const ArrayView<const GlobalOrdinal> &cols, 00795 const ArrayView<const Scalar> &vals); 00796 00811 LocalOrdinal 00812 sumIntoLocalValues (const LocalOrdinal localRow, 00813 const ArrayView<const LocalOrdinal>& cols, 00814 const ArrayView<const Scalar>& vals); 00815 00817 void setAllToScalar (const Scalar &alpha); 00818 00820 void scale (const Scalar &alpha); 00821 00823 00830 void 00831 setAllValues (const t_RowPtrs& rowPointers, 00832 const t_LocalOrdinal_1D& columnIndices, 00833 const t_ValuesType& values); 00834 00836 00847 void 00848 setAllValues (const ArrayRCP<size_t>& rowPointers, 00849 const ArrayRCP<LocalOrdinal>& columnIndices, 00850 const ArrayRCP<Scalar>& values); 00851 00852 void 00853 getAllValues (ArrayRCP<const size_t>& rowPointers, 00854 ArrayRCP<const LocalOrdinal>& columnIndices, 00855 ArrayRCP<const Scalar>& values) const; 00856 00858 00859 00860 00889 void globalAssemble(); 00890 00904 void resumeFill (const RCP<ParameterList>& params = null); 00905 00939 void 00940 fillComplete (const RCP<const map_type>& domainMap, 00941 const RCP<const map_type>& rangeMap, 00942 const RCP<ParameterList>& params = null); 00943 00956 void fillComplete (const RCP<ParameterList>& params = null); 00957 00969 void 00970 expertStaticFillComplete (const RCP<const map_type>& domainMap, 00971 const RCP<const map_type>& rangeMap, 00972 const RCP<const import_type>& importer = Teuchos::null, 00973 const RCP<const export_type>& exporter = Teuchos::null, 00974 const RCP<ParameterList>& params = Teuchos::null); 00975 00985 void 00986 replaceColMap (const Teuchos::RCP<const map_type>& newColMap); 00987 01069 void 01070 reindexColumns (crs_graph_type* const graph, 01071 const Teuchos::RCP<const map_type>& newColMap, 01072 const Teuchos::RCP<const import_type>& newImport = Teuchos::null, 01073 const bool sortEachRow = true); 01074 01087 void 01088 replaceDomainMapAndImporter (const Teuchos::RCP<const map_type>& newDomainMap, 01089 Teuchos::RCP<const import_type>& newImporter); 01090 01104 virtual void 01105 removeEmptyProcessesInPlace (const Teuchos::RCP<const map_type>& newMap); 01106 01108 01109 01110 01112 RCP<const Comm<int> > getComm() const; 01113 01115 RCP<node_type> getNode () const; 01116 01118 RCP<const map_type> getRowMap () const; 01119 01121 RCP<const map_type> getColMap () const; 01122 01124 RCP<const RowGraph<LocalOrdinal,GlobalOrdinal,node_type> > getGraph () const; 01125 01127 RCP<const crs_graph_type> getCrsGraph () const; 01128 01130 k_local_matrix_type getLocalMatrix () const {return k_lclMatrix_; } 01131 01132 01152 global_size_t getGlobalNumRows() const; 01153 01159 global_size_t getGlobalNumCols() const; 01160 01167 size_t getNodeNumRows() const; 01168 01172 size_t getNodeNumCols() const; 01173 01175 GlobalOrdinal getIndexBase() const; 01176 01178 global_size_t getGlobalNumEntries() const; 01179 01181 size_t getNodeNumEntries() const; 01182 01184 01185 size_t getNumEntriesInGlobalRow(GlobalOrdinal globalRow) const; 01186 01188 01189 size_t getNumEntriesInLocalRow(LocalOrdinal localRow) const; 01190 01192 01194 global_size_t getGlobalNumDiags() const; 01195 01197 01199 size_t getNodeNumDiags() const; 01200 01202 01204 size_t getGlobalMaxNumRowEntries() const; 01205 01207 01209 size_t getNodeMaxNumRowEntries() const; 01210 01212 bool hasColMap() const; 01213 01215 01217 bool isLowerTriangular() const; 01218 01220 01222 bool isUpperTriangular() const; 01223 01243 bool isLocallyIndexed() const; 01244 01264 bool isGloballyIndexed() const; 01265 01288 bool isFillComplete() const; 01289 01312 bool isFillActive() const; 01313 01315 01321 bool isStorageOptimized() const; 01322 01324 ProfileType getProfileType() const; 01325 01327 bool isStaticGraph() const; 01328 01340 typename ScalarTraits<Scalar>::magnitudeType getFrobeniusNorm() const; 01341 01343 virtual bool supportsRowViews() const; 01344 01346 01356 void 01357 getGlobalRowCopy (GlobalOrdinal GlobalRow, 01358 const ArrayView<GlobalOrdinal> &Indices, 01359 const ArrayView<Scalar> &Values, 01360 size_t &NumEntries) const; 01361 01363 01375 void 01376 getLocalRowCopy (LocalOrdinal LocalRow, 01377 const ArrayView<LocalOrdinal> &Indices, 01378 const ArrayView<Scalar> &Values, 01379 size_t &NumEntries) const; 01380 01394 void 01395 getGlobalRowView (GlobalOrdinal GlobalRow, 01396 ArrayView<const GlobalOrdinal> &indices, 01397 ArrayView<const Scalar> &values) const; 01398 01412 void 01413 getLocalRowView (LocalOrdinal LocalRow, 01414 ArrayView<const LocalOrdinal> &indices, 01415 ArrayView<const Scalar> &values) const; 01416 01422 void getLocalDiagCopy (Vector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>& diag) const; 01423 01461 void getLocalDiagOffsets (Teuchos::ArrayRCP<size_t>& offsets) const; 01462 01481 void 01482 getLocalDiagCopy (Vector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>& diag, 01483 const Teuchos::ArrayView<const size_t>& offsets) const; 01484 01486 void leftScale(const Vector<Scalar, LocalOrdinal, GlobalOrdinal, node_type>& x); 01487 01489 void rightScale(const Vector<Scalar, LocalOrdinal, GlobalOrdinal, node_type>& x); 01490 01492 01493 01494 01543 template <class DomainScalar, class RangeScalar> 01544 void 01545 localMultiply (const MultiVector<DomainScalar,LocalOrdinal,GlobalOrdinal,node_type>& X, 01546 MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,node_type>& Y, 01547 Teuchos::ETransp trans, 01548 RangeScalar alpha, 01549 RangeScalar beta) const; 01550 01575 template <class DomainScalar, class RangeScalar> 01576 void 01577 localGaussSeidel (const MultiVector<DomainScalar,LocalOrdinal,GlobalOrdinal,node_type> &B, 01578 MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,node_type> &X, 01579 const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type> &D, 01580 const RangeScalar& dampingFactor, 01581 const KokkosClassic::ESweepDirection direction) const; 01582 01609 template <class DomainScalar, class RangeScalar> 01610 void 01611 reorderedLocalGaussSeidel (const MultiVector<DomainScalar,LocalOrdinal,GlobalOrdinal,node_type>& B, 01612 MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,node_type>& X, 01613 const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>& D, 01614 const ArrayView<LocalOrdinal>& rowIndices, 01615 const RangeScalar& dampingFactor, 01616 const KokkosClassic::ESweepDirection direction) const; 01617 01634 template <class DomainScalar, class RangeScalar> 01635 void 01636 localSolve (const MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,node_type>& Y, 01637 MultiVector<DomainScalar,LocalOrdinal,GlobalOrdinal,node_type>& X, 01638 Teuchos::ETransp trans) const; 01639 01641 template <class T> 01642 RCP<CrsMatrix<T,LocalOrdinal,GlobalOrdinal,node_type> > convert() const; 01643 01645 01646 01647 01658 void 01659 apply (const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>& X, 01660 MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>&Y, 01661 Teuchos::ETransp mode = Teuchos::NO_TRANS, 01662 Scalar alpha = ScalarTraits<Scalar>::one(), 01663 Scalar beta = ScalarTraits<Scalar>::zero()) const; 01664 01666 bool hasTransposeApply() const; 01667 01671 RCP<const map_type> getDomainMap () const; 01672 01676 RCP<const map_type> getRangeMap () const; 01677 01679 01680 01681 01746 void 01747 gaussSeidel (const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type> &B, 01748 MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type> &X, 01749 const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type> &D, 01750 const Scalar& dampingFactor, 01751 const ESweepDirection direction, 01752 const int numSweeps) const; 01753 01820 void 01821 reorderedGaussSeidel (const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>& B, 01822 MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>& X, 01823 const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>& D, 01824 const ArrayView<LocalOrdinal>& rowIndices, 01825 const Scalar& dampingFactor, 01826 const ESweepDirection direction, 01827 const int numSweeps) const; 01828 01857 void 01858 gaussSeidelCopy (MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type> &X, 01859 const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type> &B, 01860 const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type> &D, 01861 const Scalar& dampingFactor, 01862 const ESweepDirection direction, 01863 const int numSweeps, 01864 const bool zeroInitialGuess) const; 01865 01895 void 01896 reorderedGaussSeidelCopy (MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>& X, 01897 const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>& B, 01898 const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,node_type>& D, 01899 const ArrayView<LocalOrdinal>& rowIndices, 01900 const Scalar& dampingFactor, 01901 const ESweepDirection direction, 01902 const int numSweeps, 01903 const bool zeroInitialGuess) const; 01904 01915 virtual Teuchos::RCP<RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, node_type> > 01916 add (const Scalar& alpha, 01917 const RowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, node_type>& A, 01918 const Scalar& beta, 01919 const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, node_type> >& domainMap, 01920 const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, node_type> >& rangeMap, 01921 const Teuchos::RCP<Teuchos::ParameterList>& params) const; 01922 01924 01925 01926 01928 std::string description() const; 01929 01931 void 01932 describe (Teuchos::FancyOStream &out, 01933 const Teuchos::EVerbosityLevel verbLevel = 01934 Teuchos::Describable::verbLevel_default) const; 01935 01937 01938 01939 01940 virtual bool 01941 checkSizes (const SrcDistObject& source); 01942 01943 virtual void 01944 copyAndPermute (const SrcDistObject& source, 01945 size_t numSameIDs, 01946 const ArrayView<const LocalOrdinal> &permuteToLIDs, 01947 const ArrayView<const LocalOrdinal> &permuteFromLIDs); 01948 01949 virtual void 01950 packAndPrepare (const SrcDistObject& source, 01951 const Teuchos::ArrayView<const LocalOrdinal>& exportLIDs, 01952 Teuchos::Array<char>& exports, 01953 const Teuchos::ArrayView<size_t>& numPacketsPerLID, 01954 size_t& constantNumPackets, 01955 Distributor& distor); 01956 01966 void 01967 unpackAndCombine (const Teuchos::ArrayView<const LocalOrdinal> &importLIDs, 01968 const Teuchos::ArrayView<const char> &imports, 01969 const Teuchos::ArrayView<size_t> &numPacketsPerLID, 01970 size_t constantNumPackets, 01971 Distributor &distor, 01972 CombineMode combineMode); 01974 01975 01976 01985 virtual void 01986 pack (const Teuchos::ArrayView<const LocalOrdinal>& exportLIDs, 01987 Teuchos::Array<char>& exports, 01988 const Teuchos::ArrayView<size_t>& numPacketsPerLID, 01989 size_t& constantNumPackets, 01990 Distributor& distor) const; 01991 01993 01995 t_ValuesType getLocalValuesView() const { return k_values1D_; } 01996 01997 private: 01998 // Friend declaration for nonmember function. 01999 template<class CrsMatrixType> 02000 friend Teuchos::RCP<CrsMatrixType> 02001 importAndFillCompleteCrsMatrix (const Teuchos::RCP<const CrsMatrixType>& sourceMatrix, 02002 const Import<typename CrsMatrixType::local_ordinal_type, 02003 typename CrsMatrixType::global_ordinal_type, 02004 typename CrsMatrixType::node_type>& importer, 02005 const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type, 02006 typename CrsMatrixType::global_ordinal_type, 02007 typename CrsMatrixType::node_type> >& domainMap, 02008 const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type, 02009 typename CrsMatrixType::global_ordinal_type, 02010 typename CrsMatrixType::node_type> >& rangeMap, 02011 const Teuchos::RCP<Teuchos::ParameterList>& params); 02012 02013 // Friend declaration for nonmember function. 02014 template<class CrsMatrixType> 02015 friend Teuchos::RCP<CrsMatrixType> 02016 exportAndFillCompleteCrsMatrix (const Teuchos::RCP<const CrsMatrixType>& sourceMatrix, 02017 const Export<typename CrsMatrixType::local_ordinal_type, 02018 typename CrsMatrixType::global_ordinal_type, 02019 typename CrsMatrixType::node_type>& exporter, 02020 const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type, 02021 typename CrsMatrixType::global_ordinal_type, 02022 typename CrsMatrixType::node_type> >& domainMap, 02023 const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type, 02024 typename CrsMatrixType::global_ordinal_type, 02025 typename CrsMatrixType::node_type> >& rangeMap, 02026 const Teuchos::RCP<Teuchos::ParameterList>& params); 02027 02028 public: 02044 void 02045 importAndFillComplete (Teuchos::RCP<CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, node_type> >& destMatrix, 02046 const import_type& importer, 02047 const Teuchos::RCP<const map_type>& domainMap, 02048 const Teuchos::RCP<const map_type>& rangeMap, 02049 const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null) const; 02050 02066 void 02067 exportAndFillComplete (Teuchos::RCP<CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, node_type> >& destMatrix, 02068 const export_type& exporter, 02069 const Teuchos::RCP<const map_type>& domainMap = Teuchos::null, 02070 const Teuchos::RCP<const map_type>& rangeMap = Teuchos::null, 02071 const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null) const; 02072 02073 private: 02094 void 02095 transferAndFillComplete (Teuchos::RCP<CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, node_type> >& destMatrix, 02096 const ::Tpetra::Details::Transfer<LocalOrdinal, GlobalOrdinal, node_type>& rowTransfer, 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 02101 // We forbid copy construction by declaring this method private 02102 // and not implementing it. 02103 CrsMatrix (const CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, node_type> &rhs); 02104 02105 // We forbid assignment (operator=) by declaring this method 02106 // private and not implementing it. 02107 CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, node_type>& 02108 operator= (const CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, node_type>& rhs); 02109 02115 void 02116 insertGlobalValuesFiltered (const GlobalOrdinal globalRow, 02117 const ArrayView<const GlobalOrdinal> &indices, 02118 const ArrayView<const Scalar> &values); 02119 02125 void 02126 insertLocalValuesFiltered (const LocalOrdinal localRow, 02127 const ArrayView<const LocalOrdinal> &indices, 02128 const ArrayView<const Scalar> &values); 02129 02137 void 02138 combineGlobalValues (const GlobalOrdinal globalRowIndex, 02139 const Teuchos::ArrayView<const GlobalOrdinal> columnIndices, 02140 const Teuchos::ArrayView<const Scalar> values, 02141 const Tpetra::CombineMode combineMode); 02142 02172 template<class BinaryFunction> 02173 LocalOrdinal 02174 transformGlobalValues (GlobalOrdinal globalRow, 02175 const Teuchos::ArrayView<const GlobalOrdinal>& indices, 02176 const Teuchos::ArrayView<const Scalar> & values, 02177 BinaryFunction f) 02178 { 02179 typedef LocalOrdinal LO; 02180 typedef GlobalOrdinal GO; 02181 using Teuchos::Array; 02182 using Teuchos::ArrayView; 02183 typedef typename ArrayView<const GO>::size_type size_type; 02184 02185 if (! isFillActive ()) { 02186 // Fill must be active in order to call this method. 02187 return Teuchos::OrdinalTraits<LO>::invalid (); 02188 } 02189 else if (values.size () != indices.size ()) { 02190 // The sizes of values and indices must match. 02191 return Teuchos::OrdinalTraits<LO>::invalid (); 02192 } 02193 02194 const LO lrow = this->getRowMap ()->getLocalElement (globalRow); 02195 if (lrow == Teuchos::OrdinalTraits<LO>::invalid ()) { 02196 // We don't own the row, so we're not allowed to modify its values. 02197 return Teuchos::OrdinalTraits<LO>::invalid (); 02198 } 02199 02200 if (staticGraph_.is_null ()) { 02201 return Teuchos::OrdinalTraits<LO>::invalid (); 02202 } 02203 const crs_graph_type& graph = *staticGraph_; 02204 RowInfo rowInfo = graph.getRowInfo (lrow); 02205 if (indices.size () == 0) { 02206 return static_cast<LO> (0); 02207 } 02208 else { 02209 ArrayView<Scalar> curVals = this->getViewNonConst (rowInfo); 02210 if (isLocallyIndexed ()) { 02211 // Convert the given global indices to local indices. 02212 // 02213 // FIXME (mfh 08 Jul 2014) Why can't we ask the graph to do 02214 // that? It could do the conversions in place, so that we 02215 // wouldn't need temporary storage. 02216 const map_type& colMap = * (this->getColMap ()); 02217 const size_type numInds = indices.size (); 02218 Array<LO> lclInds (numInds); 02219 for (size_type k = 0; k < numInds; ++k) { 02220 // There is no need to filter out indices not in the 02221 // column Map. Those that aren't will be mapped to 02222 // invalid(), which the graph's transformGlobalValues() 02223 // will filter out (but not count in its return value). 02224 lclInds[k] = colMap.getLocalElement (indices[k]); 02225 } 02226 return graph.template transformLocalValues<Scalar, BinaryFunction> (rowInfo, curVals, 02227 lclInds (), values, f); 02228 } 02229 else if (isGloballyIndexed ()) { 02230 return graph.template transformGlobalValues<Scalar, BinaryFunction> (rowInfo, curVals, 02231 indices, values, f); 02232 } 02233 else { 02234 // If the graph is neither locally nor globally indexed on 02235 // the calling process, that means that the calling process 02236 // can't possibly have any entries in the owned row. Thus, 02237 // there are no entries to transform, so we return zero. 02238 return static_cast<LO> (0); 02239 } 02240 } 02241 } 02242 02243 private: 02246 void 02247 insertNonownedGlobalValues (const GlobalOrdinal globalRow, 02248 const Teuchos::ArrayView<const GlobalOrdinal>& indices, 02249 const Teuchos::ArrayView<const Scalar>& values); 02250 02251 protected: 02252 // useful typedefs 02253 typedef OrdinalTraits<LocalOrdinal> OTL; 02254 typedef ScalarTraits<Scalar> STS; 02255 typedef typename STS::magnitudeType Magnitude; 02256 typedef ScalarTraits<Magnitude> STM; 02257 typedef MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, node_type> MV; 02258 typedef Vector<Scalar, LocalOrdinal, GlobalOrdinal, node_type> V; 02259 typedef crs_graph_type Graph; 02260 02261 // Enums 02262 enum GraphAllocationStatus { 02263 GraphAlreadyAllocated, 02264 GraphNotYetAllocated 02265 }; 02266 02283 void allocateValues (ELocalGlobal lg, GraphAllocationStatus gas); 02284 02290 void sortEntries(); 02291 02297 void mergeRedundantEntries(); 02298 02306 void clearGlobalConstants(); 02307 02316 void computeGlobalConstants(); 02317 02330 mutable Teuchos::RCP<MV> importMV_; 02331 02344 mutable Teuchos::RCP<MV> exportMV_; 02345 02365 Teuchos::RCP<MV> 02366 getColumnMapMultiVector (const MV& X_domainMap, 02367 const bool force = false) const; 02368 02390 Teuchos::RCP<MV> 02391 getRowMapMultiVector (const MV& Y_rangeMap, 02392 const bool force = false) const; 02393 02395 void 02396 applyNonTranspose (const MV& X_in, 02397 MV& Y_in, 02398 Scalar alpha, 02399 Scalar beta) const; 02400 02402 void 02403 applyTranspose (const MV& X_in, 02404 MV& Y_in, 02405 const Teuchos::ETransp mode, 02406 Scalar alpha, 02407 Scalar beta) const; 02408 02409 // matrix data accessors 02410 02412 Teuchos::ArrayView<const Scalar> getView (RowInfo rowinfo) const; 02413 02415 Teuchos::ArrayView<Scalar> getViewNonConst (RowInfo rowinfo); 02416 02417 // local Kokkos objects 02418 02424 void fillLocalMatrix (const Teuchos::RCP<Teuchos::ParameterList>& params); 02425 02431 void fillLocalGraphAndMatrix (const Teuchos::RCP<Teuchos::ParameterList>& params); 02432 02434 void checkInternalState() const; 02435 02447 02448 RCP<const Graph> staticGraph_; 02449 RCP< Graph> myGraph_; 02451 02453 k_local_matrix_type k_lclMatrix_; 02454 02467 02468 ArrayRCP<Scalar> values1D_; 02469 t_ValuesType k_values1D_; 02470 ArrayRCP<Array<Scalar> > values2D_; 02472 02482 Details::EStorageStatus storageStatus_; 02483 02485 bool fillComplete_; 02486 02514 std::map<GlobalOrdinal, std::pair<Teuchos::Array<GlobalOrdinal>, 02515 Teuchos::Array<Scalar> > > nonlocals_; 02516 02523 mutable Magnitude frobNorm_; 02524 02525 public: 02526 // FIXME (mfh 24 Feb 2014) Is it _really_ necessary to make this a 02527 // public inner class of CrsMatrix? It looks like it doesn't 02528 // depend on any implementation details of CrsMatrix at all. It 02529 // should really be declared and defined outside of CrsMatrix. 02530 template<class ViewType, class OffsetViewType> 02531 struct pack_functor { 02532 typedef typename ViewType::device_type device_type; 02533 ViewType src_; 02534 ViewType dst_; 02535 OffsetViewType src_offset_; 02536 OffsetViewType dst_offset_; 02537 typedef typename OffsetViewType::non_const_value_type scalar_index_type; 02538 02539 pack_functor (ViewType dst, ViewType src, 02540 OffsetViewType dst_offset, OffsetViewType src_offset) : 02541 src_ (src), 02542 dst_ (dst), 02543 src_offset_ (src_offset), 02544 dst_offset_ (dst_offset) 02545 {} 02546 02547 KOKKOS_INLINE_FUNCTION 02548 void operator () (const LocalOrdinal row) const { 02549 scalar_index_type srcPos = src_offset_(row); 02550 const scalar_index_type dstEnd = dst_offset_(row+1); 02551 scalar_index_type dstPos = dst_offset_(row); 02552 for ( ; dstPos < dstEnd; ++dstPos, ++srcPos) { 02553 dst_(dstPos) = src_(srcPos); 02554 } 02555 } 02556 }; 02557 }; // class CrsMatrix 02558 02559 } // namespace Tpetra 02560 02561 #endif // TPETRA_KOKKOSREFACTOR_CRSMATRIX_DECL_HPP
1.7.6.1
