|
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 00045 // TODO: row-wise insertion of entries in globalAssemble() may be more efficient 00046 00047 // TODO: add typeglobs: CrsMatrix<Scalar,typeglob> 00048 // TODO: add template (template) parameter for nonlocal container (this will be part of typeglob) 00049 00050 #include <Kokkos_DefaultNode.hpp> 00051 #include <Kokkos_DefaultKernels.hpp> 00052 00053 #include "Tpetra_ConfigDefs.hpp" 00054 #include "Tpetra_RowMatrix.hpp" 00055 #include "Tpetra_DistObject.hpp" 00056 #include "Tpetra_CrsGraph.hpp" 00057 #include "Tpetra_Vector.hpp" 00058 #include "Tpetra_CrsMatrixMultiplyOp_decl.hpp" 00059 00060 #ifndef DOXYGEN_SHOULD_SKIP_THIS 00061 namespace Tpetra { 00062 // struct for i,j,v triplets 00063 template <class Ordinal, class Scalar> 00064 struct CrsIJV { 00065 CrsIJV(); 00066 CrsIJV(Ordinal row, Ordinal col, const Scalar &val); 00067 Ordinal i,j; 00068 Scalar v; 00069 }; 00070 } 00071 00072 namespace Teuchos { 00073 // SerializationTraits specialization for CrsIJV, using DirectSerialization 00074 template <typename Ordinal, typename Scalar> 00075 class SerializationTraits<int,Tpetra::CrsIJV<Ordinal,Scalar> > 00076 : public DirectSerializationTraits<int,Tpetra::CrsIJV<Ordinal,Scalar> > 00077 {}; 00078 } 00079 00080 namespace std { 00081 template <class Ordinal, class Scalar> 00082 bool operator<(const Tpetra::CrsIJV<Ordinal,Scalar> &ijv1, const Tpetra::CrsIJV<Ordinal,Scalar> &ijv2); 00083 } 00084 #endif 00085 00086 namespace Tpetra { 00087 00089 00180 template <class Scalar, 00181 class LocalOrdinal = int, 00182 class GlobalOrdinal = LocalOrdinal, 00183 class Node = Kokkos::DefaultNode::DefaultNodeType, 00184 class LocalMatOps = typename Kokkos::DefaultKernels<Scalar,LocalOrdinal,Node>::SparseOps > 00185 class CrsMatrix : public RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>, 00186 public DistObject<char, LocalOrdinal,GlobalOrdinal,Node> { 00187 public: 00188 typedef Scalar scalar_type; 00189 typedef LocalOrdinal local_ordinal_type; 00190 typedef GlobalOrdinal global_ordinal_type; 00191 typedef Node node_type; 00192 typedef Map<LocalOrdinal,GlobalOrdinal,Node> map_type; 00193 // backwards compatibility defines both of these 00194 typedef LocalMatOps mat_vec_type; 00195 typedef LocalMatOps mat_solve_type; 00196 00198 00199 00217 CrsMatrix (const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& rowMap, 00218 size_t maxNumEntriesPerRow, 00219 ProfileType pftype = DynamicProfile, 00220 const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null); 00221 00239 CrsMatrix (const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& rowMap, 00240 const ArrayRCP<const size_t>& NumEntriesPerRowToAlloc, 00241 ProfileType pftype = DynamicProfile, 00242 const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null); 00243 00266 CrsMatrix (const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& rowMap, 00267 const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& colMap, 00268 size_t maxNumEntriesPerRow, 00269 ProfileType pftype = DynamicProfile, 00270 const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null); 00271 00294 CrsMatrix (const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& rowMap, 00295 const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> >& colMap, 00296 const ArrayRCP<const size_t>& NumEntriesPerRowToAlloc, 00297 ProfileType pftype = DynamicProfile, 00298 const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null); 00299 00320 explicit CrsMatrix (const Teuchos::RCP<const CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps> >& graph, 00321 const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null); 00322 00324 virtual ~CrsMatrix(); 00325 00327 00328 00329 00331 00340 void insertGlobalValues(GlobalOrdinal globalRow, const ArrayView<const GlobalOrdinal> &cols, const ArrayView<const Scalar> &vals); 00341 00343 00354 void insertLocalValues(LocalOrdinal localRow, 00355 const ArrayView<const LocalOrdinal> &cols, 00356 const ArrayView<const Scalar> &vals); 00357 00359 00364 void replaceGlobalValues(GlobalOrdinal globalRow, 00365 const ArrayView<const GlobalOrdinal> &cols, 00366 const ArrayView<const Scalar> &vals); 00367 00369 00371 void replaceLocalValues(LocalOrdinal localRow, 00372 const ArrayView<const LocalOrdinal> &cols, 00373 const ArrayView<const Scalar> &vals); 00374 00376 00381 void sumIntoGlobalValues(GlobalOrdinal globalRow, 00382 const ArrayView<const GlobalOrdinal> &cols, 00383 const ArrayView<const Scalar> &vals); 00384 00385 00387 00392 void sumIntoLocalValues(LocalOrdinal globalRow, 00393 const ArrayView<const LocalOrdinal> &cols, 00394 const ArrayView<const Scalar> &vals); 00395 00397 void setAllToScalar(const Scalar &alpha); 00398 00400 void scale(const Scalar &alpha); 00401 00403 00404 00405 00412 void globalAssemble(); 00413 00422 void resumeFill(const RCP<ParameterList> ¶ms = null); 00423 00434 void fillComplete(const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &domainMap, 00435 const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &rangeMap, 00436 const RCP<ParameterList> ¶ms = null); 00437 00450 void fillComplete(const RCP<ParameterList> ¶ms = null); 00451 00453 00454 00455 00457 const RCP<const Comm<int> > & getComm() const; 00458 00460 RCP<Node> getNode() const; 00461 00463 const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & getRowMap() const; 00464 00466 const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & getColMap() const; 00467 00469 RCP<const RowGraph<LocalOrdinal,GlobalOrdinal,Node> > getGraph() const; 00470 00472 RCP<const CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps> > getCrsGraph() const; 00473 00493 global_size_t getGlobalNumRows() const; 00494 00499 global_size_t getGlobalNumCols() const; 00500 00502 size_t getNodeNumRows() const; 00503 00505 00507 size_t getNodeNumCols() const; 00508 00510 GlobalOrdinal getIndexBase() const; 00511 00513 global_size_t getGlobalNumEntries() const; 00514 00516 size_t getNodeNumEntries() const; 00517 00519 00520 size_t getNumEntriesInGlobalRow(GlobalOrdinal globalRow) const; 00521 00523 00524 size_t getNumEntriesInLocalRow(LocalOrdinal localRow) const; 00525 00527 00529 global_size_t getGlobalNumDiags() const; 00530 00532 00534 size_t getNodeNumDiags() const; 00535 00537 00539 size_t getGlobalMaxNumRowEntries() const; 00540 00542 00544 size_t getNodeMaxNumRowEntries() const; 00545 00547 bool hasColMap() const; 00548 00550 00552 bool isLowerTriangular() const; 00553 00555 00557 bool isUpperTriangular() const; 00558 00560 bool isLocallyIndexed() const; 00561 00563 bool isGloballyIndexed() const; 00564 00566 bool isFillComplete() const; 00567 00569 bool isFillActive() const; 00570 00572 00578 bool isStorageOptimized() const; 00579 00581 ProfileType getProfileType() const; 00582 00584 bool isStaticGraph() const; 00585 00587 00593 typename ScalarTraits<Scalar>::magnitudeType getFrobeniusNorm() const; 00594 00596 00606 void getGlobalRowCopy(GlobalOrdinal GlobalRow, 00607 const ArrayView<GlobalOrdinal> &Indices, 00608 const ArrayView<Scalar> &Values, 00609 size_t &NumEntries 00610 ) const; 00611 00613 00625 void getLocalRowCopy(LocalOrdinal LocalRow, 00626 const ArrayView<LocalOrdinal> &Indices, 00627 const ArrayView<Scalar> &Values, 00628 size_t &NumEntries 00629 ) const; 00630 00632 00641 void getGlobalRowView(GlobalOrdinal GlobalRow, ArrayView<const GlobalOrdinal> &indices, ArrayView<const Scalar> &values) const; 00642 00644 00653 void getLocalRowView(LocalOrdinal LocalRow, ArrayView<const LocalOrdinal> &indices, ArrayView<const Scalar> &values) const; 00654 00656 00658 void getLocalDiagCopy(Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &diag) const; 00659 00661 void leftScale(const Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& x); 00662 00664 void rightScale(const Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node>& x); 00665 00667 00668 00669 00671 00681 template <class DomainScalar, class RangeScalar> 00682 void 00683 localMultiply (const MultiVector<DomainScalar,LocalOrdinal,GlobalOrdinal,Node>& X, 00684 MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,Node>& Y, 00685 Teuchos::ETransp trans, 00686 RangeScalar alpha, 00687 RangeScalar beta) const; 00688 00705 template <class DomainScalar, class RangeScalar> 00706 void 00707 localSolve (const MultiVector<RangeScalar,LocalOrdinal,GlobalOrdinal,Node>& Y, 00708 MultiVector<DomainScalar,LocalOrdinal,GlobalOrdinal,Node>& X, 00709 Teuchos::ETransp trans) const; 00710 00712 template <class T> 00713 RCP<CrsMatrix<T,LocalOrdinal,GlobalOrdinal,Node> > convert() const; 00714 00716 00717 00718 00720 00723 void apply(const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> & X, MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &Y, 00724 Teuchos::ETransp mode = Teuchos::NO_TRANS, 00725 Scalar alpha = ScalarTraits<Scalar>::one(), 00726 Scalar beta = ScalarTraits<Scalar>::zero()) const; 00727 00729 bool hasTransposeApply() const; 00730 00733 const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & getDomainMap() const; 00734 00737 const RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > & getRangeMap() const; 00738 00740 00741 00742 00744 std::string description() const; 00745 00747 void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const; 00748 00750 00751 00752 00753 bool checkSizes(const DistObject<char, LocalOrdinal,GlobalOrdinal,Node>& source); 00754 00755 void copyAndPermute(const DistObject<char, LocalOrdinal,GlobalOrdinal,Node>& source, 00756 size_t numSameIDs, 00757 const ArrayView<const LocalOrdinal> &permuteToLIDs, 00758 const ArrayView<const LocalOrdinal> &permuteFromLIDs); 00759 00760 void packAndPrepare(const DistObject<char, LocalOrdinal,GlobalOrdinal,Node>& source, 00761 const ArrayView<const LocalOrdinal> &exportLIDs, 00762 Array<char> &exports, 00763 const ArrayView<size_t> & numPacketsPerLID, 00764 size_t& constantNumPackets, 00765 Distributor &distor); 00766 00776 void 00777 unpackAndCombine (const Teuchos::ArrayView<const LocalOrdinal> &importLIDs, 00778 const Teuchos::ArrayView<const char> &imports, 00779 const Teuchos::ArrayView<size_t> &numPacketsPerLID, 00780 size_t constantNumPackets, 00781 Distributor &distor, 00782 CombineMode combineMode); 00784 00785 private: 00786 // We forbid copy construction by declaring this method private 00787 // and not implementing it. 00788 CrsMatrix (const CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps> &rhs); 00789 00790 // We forbid assignment (operator=) by declaring this method 00791 // private and not implementing it. 00792 CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps>& 00793 operator= (const CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps> &rhs); 00794 00802 void 00803 combineGlobalValues (const GlobalOrdinal globalRowIndex, 00804 const Teuchos::ArrayView<const GlobalOrdinal> columnIndices, 00805 const Teuchos::ArrayView<const Scalar> values, 00806 const Tpetra::CombineMode combineMode); 00807 00822 template<class BinaryFunction> 00823 void 00824 transformGlobalValues (GlobalOrdinal globalRow, 00825 const Teuchos::ArrayView<const GlobalOrdinal>& indices, 00826 const Teuchos::ArrayView<const Scalar> & values, 00827 BinaryFunction f) 00828 { 00829 typedef LocalOrdinal LO; 00830 typedef GlobalOrdinal GO; 00831 typedef Node NT; 00832 using Teuchos::Array; 00833 using Teuchos::ArrayView; 00834 00835 TEUCHOS_TEST_FOR_EXCEPTION(values.size() != indices.size(), 00836 std::invalid_argument, "transformGlobalValues: values.size() = " 00837 << values.size() << " != indices.size() = " << indices.size() << "."); 00838 00839 const LO lrow = this->getRowMap()->getLocalElement(globalRow); 00840 00841 TEUCHOS_TEST_FOR_EXCEPTION(lrow == LOT::invalid(), std::invalid_argument, 00842 "transformGlobalValues: The given global row index " << globalRow 00843 << " is not owned by the calling process (rank " 00844 << this->getRowMap()->getComm()->getRank() << ")."); 00845 00846 RowInfo rowInfo = staticGraph_->getRowInfo(lrow); 00847 if (indices.size() > 0) { 00848 if (isLocallyIndexed()) { 00849 // Convert global indices to local indices. 00850 const Map<LO, GO, NT> &colMap = *(this->getColMap()); 00851 Array<LO> lindices (indices.size()); 00852 typename ArrayView<const GO>::iterator gindit = indices.begin(); 00853 typename Array<LO>::iterator lindit = lindices.begin(); 00854 while (gindit != indices.end()) { 00855 // No need to filter before asking the column Map to 00856 // convert GID->LID. If the GID doesn't exist in the 00857 // column Map, the GID will be mapped to invalid(), which 00858 // will not be found in the graph. 00859 *lindit++ = colMap.getLocalElement(*gindit++); 00860 } 00861 typename Graph::SLocalGlobalViews inds_view; 00862 inds_view.linds = lindices(); 00863 staticGraph_->template transformValues<LocalIndices>(rowInfo, inds_view, this->getViewNonConst(rowInfo).begin(), values.begin(), f); 00864 } 00865 else if (isGloballyIndexed()) { 00866 typename Graph::SLocalGlobalViews inds_view; 00867 inds_view.ginds = indices; 00868 staticGraph_->template transformValues<GlobalIndices>(rowInfo, inds_view, this->getViewNonConst(rowInfo).begin(), values.begin(), f); 00869 } 00870 } 00871 } 00872 00873 protected: 00874 // useful typedefs 00875 typedef OrdinalTraits<LocalOrdinal> LOT; 00876 typedef OrdinalTraits<GlobalOrdinal> GOT; 00877 typedef ScalarTraits<Scalar> ST; 00878 typedef typename ST::magnitudeType Magnitude; 00879 typedef MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV; 00880 typedef Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> V; 00881 typedef CrsGraph<LocalOrdinal,GlobalOrdinal,Node,LocalMatOps> Graph; 00882 typedef typename LocalMatOps::template bind_scalar<Scalar>::other_type sparse_ops_type; 00883 typedef typename sparse_ops_type::template graph<LocalOrdinal,Node>::graph_type local_graph_type; 00884 typedef typename sparse_ops_type::template matrix<Scalar,LocalOrdinal,Node>::matrix_type local_matrix_type; 00885 // Enums 00886 enum GraphAllocationStatus { 00887 GraphAlreadyAllocated, 00888 GraphNotYetAllocated 00889 }; 00890 00907 void allocateValues (ELocalGlobal lg, GraphAllocationStatus gas); 00908 00914 void sortEntries(); 00915 00921 void mergeRedundantEntries(); 00922 00930 void clearGlobalConstants(); 00931 00940 void computeGlobalConstants(); 00941 00942 // matrix data accessors 00943 ArrayView<const Scalar> getView(RowInfo rowinfo) const; 00944 ArrayView< Scalar> getViewNonConst(RowInfo rowinfo); 00945 // local Kokkos objects 00946 void fillLocalMatrix(const RCP<ParameterList> ¶ms); 00947 void fillLocalGraphAndMatrix(const RCP<ParameterList> ¶ms); 00948 // debugging 00949 void checkInternalState() const; 00950 00962 00963 RCP<const Graph> staticGraph_; 00964 RCP< Graph> myGraph_; 00966 00967 RCP<sparse_ops_type> lclMatOps_; 00968 RCP<local_matrix_type> lclMatrix_; 00969 00970 // matrix values. before allocation, both are null. 00971 // after allocation, one is null. 00972 // 1D == StaticAllocation, 2D == DynamicAllocation 00973 // The allocation always matches that of graph_, as the graph does the allocation for the matrix. 00974 ArrayRCP<Scalar> values1D_; 00975 ArrayRCP<ArrayRCP<Scalar> > values2D_; 00976 // TODO: these could be allocated at resumeFill() and de-allocated at fillComplete() to make for very fast getView()/getViewNonConst() 00977 // ArrayRCP< typedef ArrayRCP<const Scalar>::iterator > rowPtrs_; 00978 // ArrayRCP< typedef ArrayRCP< Scalar>::iterator > rowPtrsNC_; 00979 00981 bool fillComplete_; 00982 00995 std::map<GlobalOrdinal, Array<std::pair<GlobalOrdinal,Scalar> > > nonlocals_; 00996 00997 // a wrapper around multiply, for use in apply; it contains a non-owning RCP to *this, therefore, it is not allowed 00998 // to persist past the destruction of *this. therefore, WE MAY NOT SHARE THIS POINTER. 00999 RCP< const CrsMatrixMultiplyOp<Scalar,Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps> > sameScalarMultiplyOp_; 01000 01001 // cached frobenius norm: -ST::one() means invalid 01002 mutable Magnitude frobNorm_; 01003 01005 bool insertGlobalValuesWarnedEfficiency_; 01007 bool insertLocalValuesWarnedEfficiency_; 01008 }; // class CrsMatrix 01009 01016 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node> 01017 Teuchos::RCP<CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > 01018 createCrsMatrix (const Teuchos::RCP<const Map<LocalOrdinal,GlobalOrdinal,Node> > &map, 01019 size_t maxNumEntriesPerRow = 0, 01020 const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null) 01021 { 01022 using Teuchos::rcp; 01023 typedef CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> matrix_type; 01024 01025 return rcp (new matrix_type (map, maxNumEntriesPerRow, DynamicProfile, params)); 01026 } 01027 01077 template<class CrsMatrixType> 01078 Teuchos::RCP<CrsMatrixType> 01079 importAndFillCompleteCrsMatrix (const Teuchos::RCP<const CrsMatrixType>& sourceMatrix, 01080 const Import<typename CrsMatrixType::local_ordinal_type, 01081 typename CrsMatrixType::global_ordinal_type, 01082 typename CrsMatrixType::node_type>& importer, 01083 const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type, 01084 typename CrsMatrixType::global_ordinal_type, 01085 typename CrsMatrixType::node_type> >& domainMap = Teuchos::null, 01086 const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type, 01087 typename CrsMatrixType::global_ordinal_type, 01088 typename CrsMatrixType::node_type> >& rangeMap = Teuchos::null, 01089 const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null) 01090 { 01091 using Teuchos::as; 01092 using Teuchos::RCP; 01093 using Teuchos::rcp; 01094 typedef Map<typename CrsMatrixType::local_ordinal_type, 01095 typename CrsMatrixType::global_ordinal_type, 01096 typename CrsMatrixType::node_type> map_type; 01097 01098 // FIXME (mfh 11 Apr 2012) The current implementation of this 01099 // method doesn't actually fuse the Import with fillComplete(). 01100 // This will change in the future. 01101 RCP<CrsMatrixType> destMat = 01102 rcp (new CrsMatrixType (importer.getTargetMap (), 01103 as<size_t> (0), 01104 DynamicProfile, 01105 params)); 01106 destMat->doImport (*sourceMatrix, importer, INSERT); 01107 01108 // Use the source matrix's domain Map as the default. 01109 RCP<const map_type> theDomainMap = 01110 domainMap.is_null () ? sourceMatrix->getDomainMap () : domainMap; 01111 // Use the source matrix's range Map as the default. 01112 RCP<const map_type> theRangeMap = 01113 rangeMap.is_null () ? sourceMatrix->getRangeMap () : rangeMap; 01114 01115 destMat->fillComplete (theDomainMap, theRangeMap); 01116 return destMat; 01117 } 01118 01152 template<class CrsMatrixType> 01153 Teuchos::RCP<CrsMatrixType> 01154 exportAndFillCompleteCrsMatrix (const Teuchos::RCP<const CrsMatrixType>& sourceMatrix, 01155 const Export<typename CrsMatrixType::local_ordinal_type, 01156 typename CrsMatrixType::global_ordinal_type, 01157 typename CrsMatrixType::node_type>& exporter, 01158 const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type, 01159 typename CrsMatrixType::global_ordinal_type, 01160 typename CrsMatrixType::node_type> >& domainMap = Teuchos::null, 01161 const Teuchos::RCP<const Map<typename CrsMatrixType::local_ordinal_type, 01162 typename CrsMatrixType::global_ordinal_type, 01163 typename CrsMatrixType::node_type> >& rangeMap = Teuchos::null, 01164 const Teuchos::RCP<Teuchos::ParameterList>& params = Teuchos::null) 01165 { 01166 using Teuchos::as; 01167 using Teuchos::RCP; 01168 using Teuchos::rcp; 01169 typedef Map<typename CrsMatrixType::local_ordinal_type, 01170 typename CrsMatrixType::global_ordinal_type, 01171 typename CrsMatrixType::node_type> map_type; 01172 01173 // FIXME (mfh 11 Apr 2012) The current implementation of this 01174 // method doesn't actually fuse the Export with fillComplete(). 01175 // This will change in the future. 01176 RCP<CrsMatrixType> destMat = 01177 rcp (new CrsMatrixType (exporter.getTargetMap (), 01178 as<size_t> (0), 01179 DynamicProfile, 01180 params)); 01181 destMat->doExport (*sourceMatrix, exporter, INSERT); 01182 01183 // Use the source matrix's domain Map as the default. 01184 RCP<const map_type> theDomainMap = 01185 domainMap.is_null () ? sourceMatrix->getDomainMap () : domainMap; 01186 // Use the source matrix's range Map as the default. 01187 RCP<const map_type> theRangeMap = 01188 rangeMap.is_null () ? sourceMatrix->getRangeMap () : rangeMap; 01189 01190 destMat->fillComplete (theDomainMap, theRangeMap); 01191 return destMat; 01192 } 01193 } // namespace Tpetra 01194 01205 #endif
1.7.6.1
