46 #ifndef XPETRA_BLOCKEDCRSMATRIX_HPP 47 #define XPETRA_BLOCKEDCRSMATRIX_HPP 68 #ifdef HAVE_XPETRA_THYRA 69 #include <Thyra_ProductVectorSpaceBase.hpp> 70 #include <Thyra_VectorSpaceBase.hpp> 71 #include <Thyra_LinearOpBase.hpp> 72 #include <Thyra_BlockedLinearOpBase.hpp> 73 #include <Thyra_PhysicallyBlockedLinearOpBase.hpp> 86 template <
class Scalar ,
94 public Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> {
102 #undef XPETRA_BLOCKEDCRSMATRIX_SHORT 128 for (
size_t r = 0; r <
Rows(); ++r)
129 for (
size_t c = 0; c <
Cols(); ++c)
136 #ifdef HAVE_XPETRA_THYRA 151 int numRangeBlocks = productRangeSpace->numBlocks();
152 int numDomainBlocks = productDomainSpace->numBlocks();
155 std::vector<Teuchos::RCP<const Map> > subRangeMaps(numRangeBlocks);
156 for (
size_t r=0; r<Teuchos::as<size_t>(numRangeBlocks); ++r) {
157 for (
size_t c=0; c<Teuchos::as<size_t>(numDomainBlocks); ++c) {
158 if (thyraOp->blockExists(r,c)) {
163 subRangeMaps[r] = xop->getRangeMap();
172 bool bRangeUseThyraStyleNumbering =
false;
173 size_t numAllElements = 0;
174 for(
size_t v = 0; v < subRangeMaps.size(); ++v) {
175 numAllElements += subRangeMaps[v]->getGlobalNumElements();
177 if ( fullRangeMap->getGlobalNumElements() != numAllElements) bRangeUseThyraStyleNumbering =
true;
181 std::vector<Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > > subDomainMaps(numDomainBlocks);
182 for (
size_t c=0; c<Teuchos::as<size_t>(numDomainBlocks); ++c) {
183 for (
size_t r=0; r<Teuchos::as<size_t>(numRangeBlocks); ++r) {
184 if (thyraOp->blockExists(r,c)) {
189 subDomainMaps[c] = xop->getDomainMap();
196 bool bDomainUseThyraStyleNumbering =
false;
198 for(
size_t v = 0; v < subDomainMaps.size(); ++v) {
199 numAllElements += subDomainMaps[v]->getGlobalNumElements();
201 if (fullDomainMap->getGlobalNumElements() != numAllElements) bDomainUseThyraStyleNumbering =
true;
210 for (
size_t r = 0; r <
Rows(); ++r) {
211 for (
size_t c = 0; c <
Cols(); ++c) {
212 if(thyraOp->blockExists(r,c)) {
240 std::vector<GlobalOrdinal> gids;
241 for(
size_t tt = 0; tt<subMaps.size(); ++tt) {
243 for(LocalOrdinal l = 0; l < Teuchos::as<LocalOrdinal>(subMap->getNodeNumElements()); ++l) {
244 GlobalOrdinal gid = subMap->getGlobalElement(l);
254 std::sort(gids.begin(), gids.end());
255 gids.erase(std::unique(gids.begin(), gids.end()), gids.end());
401 for (
size_t r = 0; r <
Rows(); ++r)
402 for (
size_t c = 0; c <
Cols(); ++c) {
411 fullrowmap_ =
MapFactory::Build(rangeMap()->lib(), rangeMap()->getGlobalNumElements(), rangeMap()->getNodeElementList(), rangeMap()->getIndexBase(), rangeMap()->getComm());
415 fullcolmap_ = Teuchos::null;
417 std::vector<GO> colmapentries;
418 for (
size_t c = 0; c <
Cols(); ++c) {
421 for (
size_t r = 0; r <
Rows(); ++r) {
424 if (Ablock != Teuchos::null) {
427 copy(colElements.
getRawPtr(), colElements.
getRawPtr() + colElements.
size(), inserter(colset, colset.begin()));
432 colmapentries.reserve(colmapentries.size() + colset.size());
433 copy(colset.begin(), colset.end(), back_inserter(colmapentries));
434 sort(colmapentries.begin(), colmapentries.end());
435 typename std::vector<GO>::iterator gendLocation;
436 gendLocation = std::unique(colmapentries.begin(), colmapentries.end());
437 colmapentries.erase(gendLocation,colmapentries.end());
441 size_t numGlobalElements = 0;
442 Teuchos::reduceAll(*(rangeMap->
getComm()), Teuchos::REDUCE_SUM, colmapentries.size(), Teuchos::outArg(numGlobalElements));
458 for (
size_t row = 0; row <
Rows(); row++)
459 for (
size_t col = 0; col <
Cols(); col++)
461 globalNumRows +=
getMatrix(row,col)->getGlobalNumRows();
465 return globalNumRows;
474 for (
size_t col = 0; col <
Cols(); col++)
475 for (
size_t row = 0; row <
Rows(); row++)
477 globalNumCols +=
getMatrix(row,col)->getGlobalNumCols();
481 return globalNumCols;
488 for (
size_t row = 0; row <
Rows(); ++row)
489 for (
size_t col = 0; col <
Cols(); col++)
491 nodeNumRows +=
getMatrix(row,col)->getNodeNumRows();
502 for (
size_t row = 0; row <
Rows(); ++row)
503 for (
size_t col = 0; col <
Cols(); ++col)
505 globalNumEntries +=
getMatrix(row,col)->getGlobalNumEntries();
507 return globalNumEntries;
514 for (
size_t row = 0; row <
Rows(); ++row)
515 for (
size_t col = 0; col <
Cols(); ++col)
517 nodeNumEntries +=
getMatrix(row,col)->getNodeNumEntries();
519 return nodeNumEntries;
561 for (
size_t i = 0; i <
blocks_.size(); ++i)
572 for (
size_t i = 0; i <
blocks_.size(); i++)
580 for (
size_t i = 0; i <
blocks_.size(); i++)
604 size_t &NumEntries)
const {
691 "apply() only supports the following modes: NO_TRANS and TRANS." );
699 for (
size_t row = 0; row <
Rows(); row++) {
703 for (
size_t col = 0; col <
Cols(); col++) {
710 Ablock->
apply(*Xblock, *tmpYblock);
711 Yblock->
update(one, *tmpYblock, one);
718 for (
size_t col = 0; col <
Cols(); col++) {
722 for (
size_t row = 0; row<
Rows(); row++) {
731 Yblock->
update(one, *tmpYblock, one);
737 Y.
update(alpha, *tmpY, beta);
798 std::string
description()
const {
return "Xpetra_BlockedCrsMatrix.description()"; }
802 out <<
"Xpetra::BlockedCrsMatrix: " <<
Rows() <<
" x " <<
Cols() << std::endl;
805 out <<
"BlockMatrix is fillComplete" << std::endl;
807 out <<
"fullRowMap" << std::endl;
814 out <<
"BlockMatrix is NOT fillComplete" << std::endl;
817 for (
size_t r = 0; r <
Rows(); ++r)
818 for (
size_t c = 0; c <
Cols(); ++c) {
819 out <<
"Block(" << r <<
"," << c <<
")" << std::endl;
859 using Teuchos::rcp_dynamic_cast;
863 "BlockedCrsMatrix::Merge: only implemented for Xpetra-style numbering. If you need the merge routine for Thyra-style numbering, report this to the Xpetra developers." );
867 for (
size_t i = 0; i <
blocks_.size(); i++)
868 if (
blocks_[i] != Teuchos::null)
877 #ifdef HAVE_XPETRA_KOKKOS_REFACTOR 878 typedef typename CrsMatrix::local_matrix_type local_matrix_type;
881 local_matrix_type getLocalMatrix ()
const {
886 #ifdef HAVE_XPETRA_THYRA 889 Xpetra::ThyraUtils<Scalar,LO,GO,Node>::toThyra(
this);
893 Teuchos::rcp_dynamic_cast<Thyra::BlockedLinearOpBase<Scalar> >(thOp);
917 "Matrix A is not completed");
944 for (
size_t j = 0; j < numEntries; ++j)
970 std::vector<Teuchos::RCP<CrsMatrix> >
blocks_;
971 #ifdef HAVE_XPETRA_THYRA 981 #define XPETRA_BLOCKEDCRSMATRIX_SHORT virtual void getLocalRowCopy(LocalOrdinal LocalRow, const ArrayView< LocalOrdinal > &Indices, const ArrayView< Scalar > &Values, size_t &NumEntries) const
Extract a list of entries in a specified local row of the matrix. Put into storage allocated by calli...
RCP< const MapExtractor > getRangeMapExtractor()
Returns map extractor class for range map.
virtual void apply(const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &X, MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Y, Teuchos::ETransp mode=Teuchos::NO_TRANS, Scalar alpha=ScalarTraits< Scalar >::one(), Scalar beta=ScalarTraits< Scalar >::zero()) const =0
Computes the sparse matrix-multivector multiplication.
virtual size_t getNodeMaxNumRowEntries() const =0
Returns the maximum number of entries across all rows/columns on this node.
virtual Teuchos::RCP< const Teuchos::Comm< int > > getComm() const =0
Get this Map's Comm object.
static Teuchos::RCP< Map< LocalOrdinal, GlobalOrdinal, Node > > Build(UnderlyingLib lib, global_size_t numGlobalElements, GlobalOrdinal indexBase, const Teuchos::RCP< const Teuchos::Comm< int > > &comm, LocalGlobal lg=Xpetra::GloballyDistributed, const Teuchos::RCP< Node > &node=defaultArgNode())
Map constructor with Xpetra-defined contiguous uniform distribution.
std::vector< Teuchos::RCP< CrsMatrix > > blocks_
virtual const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getColMap() const =0
Returns the Map that describes the column distribution in this matrix.
bool isFillComplete() const
Returns true if fillComplete() has been called and the matrix is in compute mode. ...
void replaceGlobalValues(GlobalOrdinal globalRow, const ArrayView< const GlobalOrdinal > &cols, const ArrayView< const Scalar > &vals)
Replace matrix entries, using global IDs.
void fillComplete(const RCP< const Map > &domainMap, const RCP< const Map > &rangeMap, const RCP< ParameterList > ¶ms=null)
Signal that data entry is complete, specifying domain and range maps.
virtual void setAllToScalar(const Scalar &alpha)
Set all matrix entries equal to scalar.
virtual void scale(const Scalar &alpha)=0
Scale the current values of a matrix, this = alpha*this.
virtual void fillComplete(const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &domainMap, const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &rangeMap, const RCP< ParameterList > ¶ms=null)=0
Signal that data entry is complete, specifying domain and range maps.
virtual void insertGlobalValues(GlobalOrdinal globalRow, const ArrayView< const GlobalOrdinal > &cols, const ArrayView< const Scalar > &vals)=0
Insert matrix entries, using global IDs.
bool is_null(const std::shared_ptr< T > &p)
RCP< const Map > getRangeMap() const
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
static RCP< CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Build(const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &rowMap, size_t maxNumEntriesPerRow, Xpetra::ProfileType pftype=Xpetra::DynamicProfile, const Teuchos::RCP< Teuchos::ParameterList > &plist=Teuchos::null)
Constructor specifying fixed number of entries for each row.
Teuchos::RCP< const MapExtractor > rangemaps_
virtual const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getRowMap() const =0
Returns the Map that describes the row distribution in this matrix.
void doImport(const Matrix &source, const Import &importer, CombineMode CM)
Import.
Teuchos::RCP< CrsMatrix > getMatrix(size_t r, size_t c) const
return block (r,c)
void getLocalRowView(LocalOrdinal LocalRow, ArrayView< const LocalOrdinal > &indices, ArrayView< const Scalar > &values) const
Extract a const, non-persisting view of local indices in a specified row of the matrix.
Exception throws to report errors in the internal logical of the program.
size_t getNodeNumEntries() const
Returns the local number of entries in this matrix.
virtual ~BlockedCrsMatrix()
Destructor.
Teuchos::RCP< Map > fullrowmap_
RCP< const Map > getDomainMap(size_t i, bool bThyraMode=false) const
Returns the Map associated with the i'th block domain of this operator. This will be null until fillC...
size_t getNodeNumDiags() const
Returns the number of local diagonal entries, based on global row/column index comparisons.
void doExport(const Matrix &dest, const Import &importer, CombineMode CM)
Export.
LocalOrdinal local_ordinal_type
void Add(const CrsMatrix &A, const Scalar scalarA, CrsMatrix &B, const Scalar scalarB) const
Add a Xpetra::CrsMatrix to another: B = B*scalarB + A*scalarA.
void setMatrix(size_t r, size_t c, Teuchos::RCP< CrsMatrix > mat)
set matrix block
size_t getGlobalMaxNumRowEntries() const
Returns the maximum number of entries across all rows/columns on all nodes.
void resumeFill(const RCP< ParameterList > ¶ms=null)
global_size_t getGlobalNumCols() const
Returns the number of global columns in the matrix.
virtual void update(const Scalar &alpha, const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Scalar &beta)=0
Update multi-vector values with scaled values of A, this = beta*this + alpha*A.
virtual void getGlobalRowCopy(GlobalOrdinal GlobalRow, const ArrayView< GlobalOrdinal > &indices, const ArrayView< Scalar > &values, size_t &numEntries) const =0
Extract a list of entries in a specified global row of this matrix. Put into pre-allocated storage...
global_size_t getGlobalNumEntries() const
Returns the global number of entries in this matrix.
size_t getNumEntriesInLocalRow(LocalOrdinal localRow) const
Returns the current number of entries on this node in the specified local row.
RCP< const CrsGraph > getCrsGraph() const
Returns the CrsGraph associated with this matrix.
void insertGlobalValues(GlobalOrdinal globalRow, const ArrayView< const GlobalOrdinal > &cols, const ArrayView< const Scalar > &vals)
Insert matrix entries, using global IDs.
void getLocalDiagCopy(Vector &diag) const
Get a copy of the diagonal entries owned by this node, with local row idices.
viewLabel_t defaultViewLabel_
GlobalOrdinal global_ordinal_type
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void doExport(const Matrix &dest, const Export &exporter, CombineMode CM)
Export (using an Importer).
void insertLocalValues(LocalOrdinal localRow, const ArrayView< const LocalOrdinal > &cols, const ArrayView< const Scalar > &vals)
Insert matrix entries, using local IDs.
RCP< const MapExtractor > getDomainMapExtractor()
Returns map extractor for domain map.
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const
Print the object with some verbosity level to an FancyOStream object.
virtual size_t Cols() const
number of column blocks
const Teuchos::RCP< const Map > getMap() const
Implements DistObject interface.
void scale(const Scalar &alpha)
Scale the current values of a matrix, this = alpha*this.
void CreateView(viewLabel_t viewLabel, const RCP< const Map > &rowMap, const RCP< const Map > &colMap)
virtual Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getMap() const =0
The Map describing the parallel distribution of this object.
viewLabel_t currentViewLabel_
TEUCHOS_DEPRECATED void reduceAll(const Comm< Ordinal > &comm, const EReductionType reductType, const Packet &send, Packet *globalReduct)
Teuchos::RCP< CrsMatrix > Merge() const
merge BlockedCrsMatrix blocks in a CrsMatrix
void doImport(const Matrix &source, const Export &exporter, CombineMode CM)
Import (using an Exporter).
void getGlobalRowView(GlobalOrdinal GlobalRow, ArrayView< const GlobalOrdinal > &indices, ArrayView< const Scalar > &values) const
Extract a const, non-persisting view of global indices in a specified row of the matrix.
size_t global_size_t
Global size_t object.
size_t getNodeMaxNumRowEntries() const
Returns the maximum number of entries across all rows/columns on this node.
static const EVerbosityLevel verbLevel_default
global_size_t getGlobalNumDiags() const
Returns the number of global diagonal entries, based on global row/column index comparisons.
global_size_t getGlobalNumRows() const
Returns the number of global rows.
virtual bool isFillComplete() const =0
Returns true if the matrix is in compute mode, i.e. if fillComplete() has been called.
RCP< const CrsGraph< int, GlobalOrdinal, Node > > toXpetra(const Epetra_CrsGraph &g)
std::string description() const
Return a simple one-line description of this object.
bool bDomainThyraMode_
boolean flag, which is true, if BlockedCrsMatrix has been created using Thyra-style numbering for sub...
BlockedCrsMatrix(Teuchos::RCP< const MapExtractor > &rangeMaps, Teuchos::RCP< const MapExtractor > &domainMaps, size_t npr, Xpetra::ProfileType pftype=Xpetra::DynamicProfile)
Constructor.
bool isGloballyIndexed() const
If matrix indices are in the global range, this function returns true. Otherwise, this function retur...
RCP< const Map > getDomainMap() const
Returns the Map associated with the full domain of this operator. This will be null until fillComplet...
virtual size_t Rows() const
number of row blocks
static Teuchos::RCP< MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Build(const Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &map, size_t NumVectors, bool zeroOut=true)
Constructor specifying the number of non-zeros for all rows.
CombineMode
Xpetra::Combine Mode enumerable type.
void fillComplete(const RCP< ParameterList > ¶ms=null)
Signal that data entry is complete.
const viewLabel_t & GetDefaultViewLabel() const
virtual size_t getNumVectors() const =0
Number of columns in the multivector.
virtual size_t getNodeNumRows() const =0
Returns the number of matrix rows owned on the calling node.
size_t getNodeNumRows() const
Returns the number of matrix rows owned on the calling node.
RCP< const Map > getRangeMap(size_t i, bool bThyraMode=false) const
virtual UnderlyingLib lib() const =0
Get the library used by this object (Tpetra or Epetra?)
virtual ScalarTraits< Scalar >::magnitudeType getFrobeniusNorm() const
Get Frobenius norm of the matrix.
void removeEmptyProcessesInPlace(const Teuchos::RCP< const Map > &newMap)
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
bool isLocallyIndexed() const
If matrix indices of all matrix blocks are in the local range, this function returns true...
void replaceLocalValues(LocalOrdinal localRow, const ArrayView< const LocalOrdinal > &cols, const ArrayView< const Scalar > &vals)
Replace matrix entries, using local IDs.
Xpetra-specific matrix class.
Teuchos::RCP< const MapExtractor > domainmaps_
virtual void apply(const MultiVector &X, MultiVector &Y, Teuchos::ETransp mode=Teuchos::NO_TRANS, Scalar alpha=ScalarTraits< Scalar >::one(), Scalar beta=ScalarTraits< Scalar >::zero()) const
Computes the sparse matrix-multivector multiplication.
bool bRangeThyraMode_
boolean flag, which is true, if BlockedCrsMatrix has been created using Thyra-style numbering for sub...