Xpetra_BlockedCrsMatrix.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // Xpetra: A linear algebra interface package
6 // Copyright 2012 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact
39 // Jonathan Hu (jhu@sandia.gov)
40 // Andrey Prokopenko (aprokop@sandia.gov)
41 // Ray Tuminaro (rstumin@sandia.gov)
42 //
43 // ***********************************************************************
44 //
45 // @HEADER
46 #ifndef XPETRA_BLOCKEDCRSMATRIX_HPP
47 #define XPETRA_BLOCKEDCRSMATRIX_HPP
48 
49 #include <Kokkos_DefaultNode.hpp>
50 
52 #include <Teuchos_Hashtable.hpp>
53 
54 #include "Xpetra_ConfigDefs.hpp"
55 #include "Xpetra_Exceptions.hpp"
56 
57 #include "Xpetra_MapFactory.hpp"
58 #include "Xpetra_MultiVector.hpp"
60 #include "Xpetra_CrsGraph.hpp"
61 #include "Xpetra_CrsMatrix.hpp"
63 
64 #include "Xpetra_MapExtractor.hpp"
65 
66 #include "Xpetra_Matrix.hpp"
67 
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>
74 #include "Xpetra_ThyraUtils.hpp"
75 #endif
76 
77 
82 namespace Xpetra {
83 
84  typedef std::string viewLabel_t;
85 
86  template <class Scalar /*= Matrix<>::scalar_type*/,
87  class LocalOrdinal /*=
88  typename Matrix<Scalar>::local_ordinal_type*/,
89  class GlobalOrdinal /*=
90  typename Matrix<Scalar, LocalOrdinal>::global_ordinal_type*/,
91  class Node /*=
92  typename Matrix<Scalar, LocalOrdinal, GlobalOrdinal>::node_type*/>
93  class BlockedCrsMatrix :
94  public Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> {
95  public:
96  typedef Scalar scalar_type;
97  typedef LocalOrdinal local_ordinal_type;
98  typedef GlobalOrdinal global_ordinal_type;
99  typedef Node node_type;
100 
101  private:
102 #undef XPETRA_BLOCKEDCRSMATRIX_SHORT
103 #include "Xpetra_UseShortNames.hpp"
104 
105  public:
106 
108 
109 
111 
119  size_t npr, Xpetra::ProfileType pftype = Xpetra::DynamicProfile)
120  : domainmaps_(domainMaps), rangemaps_(rangeMaps)
121  {
122  bRangeThyraMode_ = rangeMaps->getThyraMode();
123  bDomainThyraMode_ = domainMaps->getThyraMode();
124 
125  blocks_.reserve(Rows()*Cols());
126 
127  // add CrsMatrix objects in row,column order
128  for (size_t r = 0; r < Rows(); ++r)
129  for (size_t c = 0; c < Cols(); ++c)
130  blocks_.push_back(CrsMatrixFactory::Build(getRangeMap(r,bRangeThyraMode_), npr, pftype));
131 
132  // Default view
134  }
135 
136 #ifdef HAVE_XPETRA_THYRA
137 
144  BlockedCrsMatrix(const Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar> >& thyraOp, const Teuchos::RCP<const Teuchos::Comm<int> >& comm)
145  : thyraOp_(thyraOp)
146  {
147  // extract information from Thyra blocked operator and rebuilt information
148  const Teuchos::RCP<const Thyra::ProductVectorSpaceBase<Scalar> > productRangeSpace = thyraOp->productRange();
149  const Teuchos::RCP<const Thyra::ProductVectorSpaceBase<Scalar> > productDomainSpace = thyraOp->productDomain();
150 
151  int numRangeBlocks = productRangeSpace->numBlocks();
152  int numDomainBlocks = productDomainSpace->numBlocks();
153 
154  // build range map extractor from Thyra::BlockedLinearOpBase object
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)) {
159  // we only need at least one block in each block row to extract the range map
160  Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > const_op = thyraOp->getBlock(r,c); // nonConst access is not allowed.
163  subRangeMaps[r] = xop->getRangeMap();
164  break;
165  }
166  }
167  }
168  Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > fullRangeMap = mergeMaps(subRangeMaps);
169 
170  // check whether the underlying Thyra operator uses Thyra-style numbering (default for most applications) or
171  // Xpetra style numbering
172  bool bRangeUseThyraStyleNumbering = false;
173  size_t numAllElements = 0;
174  for(size_t v = 0; v < subRangeMaps.size(); ++v) {
175  numAllElements += subRangeMaps[v]->getGlobalNumElements();
176  }
177  if ( fullRangeMap->getGlobalNumElements() != numAllElements) bRangeUseThyraStyleNumbering = true;
178  rangemaps_ = Teuchos::rcp(new Xpetra::MapExtractor<Scalar,LocalOrdinal,GlobalOrdinal,Node>(fullRangeMap, subRangeMaps, bRangeUseThyraStyleNumbering));
179 
180  // build domain map extractor from Thyra::BlockedLinearOpBase object
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)) {
185  // we only need at least one block in each block row to extract the range map
186  Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > const_op = thyraOp->getBlock(r,c); // nonConst access is not allowed.
189  subDomainMaps[c] = xop->getDomainMap();
190  break;
191  }
192  }
193  }
194  Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > fullDomainMap = mergeMaps(subDomainMaps);
195  // plausibility check for numbering style (Xpetra or Thyra)
196  bool bDomainUseThyraStyleNumbering = false;
197  numAllElements = 0;
198  for(size_t v = 0; v < subDomainMaps.size(); ++v) {
199  numAllElements += subDomainMaps[v]->getGlobalNumElements();
200  }
201  if (fullDomainMap->getGlobalNumElements() != numAllElements) bDomainUseThyraStyleNumbering = true;
202  domainmaps_ = Teuchos::rcp(new Xpetra::MapExtractor<Scalar,LocalOrdinal,GlobalOrdinal,Node>(fullDomainMap, subDomainMaps, bDomainUseThyraStyleNumbering));
203 
204  // store numbering mode
205  bRangeThyraMode_ = bRangeUseThyraStyleNumbering;
206  bDomainThyraMode_ = bDomainUseThyraStyleNumbering;
207 
208  // add CrsMatrix objects in row,column order
209  blocks_.reserve(Rows()*Cols());
210  for (size_t r = 0; r < Rows(); ++r) {
211  for (size_t c = 0; c < Cols(); ++c) {
212  if(thyraOp->blockExists(r,c)) {
213  Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > const_op = thyraOp->getBlock(r,c); // nonConst access is not allowed.
214  Teuchos::RCP<Thyra::LinearOpBase<Scalar> > op = Teuchos::rcp_const_cast<Thyra::LinearOpBase<Scalar> >(const_op); // cast away const
217  blocks_.push_back(xop);
218  } else {
219  // add empty block
221  }
222  }
223  }
224  // Default view
226  }
227 
228  private:
230 
236  Teuchos::RCP<const Map> mergeMaps(std::vector<Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > > & subMaps) {
237  // TODO merging for Thyra mode is missing (similar to what we do in constructor of MapExtractor
238 
239  // merge submaps to global map
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);
245  gids.push_back(gid);
246  }
247  }
248 
249  // we have to sort the matrix entries and get rid of the double entries
250  // since we use this to detect Thyra-style numbering or Xpetra-style
251  // numbering. In Thyra-style numbering mode, the Xpetra::MapExtractor builds
252  // the correct row maps.
254  std::sort(gids.begin(), gids.end());
255  gids.erase(std::unique(gids.begin(), gids.end()), gids.end());
256  Teuchos::ArrayView<GO> gidsView(&gids[0], gids.size());
257  Teuchos::RCP<Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > fullMap = Xpetra::MapFactory<LocalOrdinal,GlobalOrdinal,Node>::Build(subMaps[0]->lib(), INVALID, gidsView, subMaps[0]->getIndexBase(), subMaps[0]->getComm());
258  return fullMap;
259  }
260 
261  public:
262 #endif
263 
265  virtual ~BlockedCrsMatrix() {}
266 
268 
269 
271 
272 
274 
296  void insertGlobalValues(GlobalOrdinal globalRow, const ArrayView<const GlobalOrdinal>& cols, const ArrayView<const Scalar>& vals) {
297  throw Xpetra::Exceptions::RuntimeError("insertGlobalValues not supported by BlockedCrsMatrix");
298  }
299 
301 
308  void insertLocalValues(LocalOrdinal localRow, const ArrayView<const LocalOrdinal>& cols, const ArrayView<const Scalar>& vals) {
309  throw Xpetra::Exceptions::RuntimeError("insertLocalValues not supported by BlockedCrsMatrix");
310  }
311 
313  throw Xpetra::Exceptions::RuntimeError("removeEmptyProcesses not supported by BlockedCrsMatrix");
314  }
315 
317 
325  void replaceGlobalValues(GlobalOrdinal globalRow,
326  const ArrayView<const GlobalOrdinal> &cols,
327  const ArrayView<const Scalar> &vals) {
328  throw Xpetra::Exceptions::RuntimeError("replaceGlobalValues not supported by BlockedCrsMatrix");
329  }
330 
332 
336  void replaceLocalValues(LocalOrdinal localRow,
337  const ArrayView<const LocalOrdinal> &cols,
338  const ArrayView<const Scalar> &vals) {
339  throw Xpetra::Exceptions::RuntimeError("replaceLocalValues not supported by BlockedCrsMatrix");
340  }
341 
343  virtual void setAllToScalar(const Scalar& alpha) {
344  throw Xpetra::Exceptions::RuntimeError("setAllToScalar not supported by BlockedCrsMatrix");
345  }
346 
348  void scale(const Scalar& alpha) {
349  throw Xpetra::Exceptions::RuntimeError("scale not supported by BlockedCrsMatrix");
350  }
351 
353 
355 
356 
365  void resumeFill(const RCP< ParameterList >& params = null) {
366  throw Xpetra::Exceptions::RuntimeError("resumeFill not supported for block matrices");
367  }
368 
382  void fillComplete(const RCP<const Map>& domainMap, const RCP<const Map>& rangeMap, const RCP<ParameterList>& params = null) {
383  throw Xpetra::Exceptions::RuntimeError("fillComplete with arguments not supported for block matrices");
384  }
385 
400  void fillComplete(const RCP<ParameterList>& params = null) {
401  for (size_t r = 0; r < Rows(); ++r)
402  for (size_t c = 0; c < Cols(); ++c) {
403  Teuchos::RCP<CrsMatrix> Ablock = getMatrix(r,c);
404 
405  if (Ablock != Teuchos::null && !Ablock->isFillComplete())
407  }
408 
409  // get full row map
410  RCP<const Map> rangeMap = rangemaps_->getFullMap();
411  fullrowmap_ = MapFactory::Build(rangeMap()->lib(), rangeMap()->getGlobalNumElements(), rangeMap()->getNodeElementList(), rangeMap()->getIndexBase(), rangeMap()->getComm());
412 
413 #if 0
414  // build full col map
415  fullcolmap_ = Teuchos::null; // delete old full column map
416 
417  std::vector<GO> colmapentries;
418  for (size_t c = 0; c < Cols(); ++c) {
419  // copy all local column lids of all block rows to colset
420  std::set<GO> colset;
421  for (size_t r = 0; r < Rows(); ++r) {
422  Teuchos::RCP<CrsMatrix> Ablock = getMatrix(r,c);
423 
424  if (Ablock != Teuchos::null) {
425  Teuchos::ArrayView<const GO> colElements = Ablock->getColMap()->getNodeElementList();
426  Teuchos::RCP<const Map> colmap = Ablock->getColMap();
427  copy(colElements.getRawPtr(), colElements.getRawPtr() + colElements.size(), inserter(colset, colset.begin()));
428  }
429  }
430 
431  // remove duplicates (entries which are in column maps of more than one block row)
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());
438  }
439 
440  // sum up number of local elements
441  size_t numGlobalElements = 0;
442  Teuchos::reduceAll(*(rangeMap->getComm()), Teuchos::REDUCE_SUM, colmapentries.size(), Teuchos::outArg(numGlobalElements));
443 
444  // store global full column map
445  const Teuchos::ArrayView<const GO> aView = Teuchos::ArrayView<const GO>(colmapentries);
446  fullcolmap_ = MapFactory::Build(rangeMap->lib(), numGlobalElements, aView, 0, rangeMap->getComm());
447 #endif
448  }
449 
451 
453 
456  global_size_t globalNumRows = 0;
457 
458  for (size_t row = 0; row < Rows(); row++)
459  for (size_t col = 0; col < Cols(); col++)
460  if (!getMatrix(row,col).is_null()) {
461  globalNumRows += getMatrix(row,col)->getGlobalNumRows();
462  break; // we need only one non-null matrix in a row
463  }
464 
465  return globalNumRows;
466  }
467 
469 
472  global_size_t globalNumCols = 0;
473 
474  for (size_t col = 0; col < Cols(); col++)
475  for (size_t row = 0; row < Rows(); row++)
476  if (!getMatrix(row,col).is_null()) {
477  globalNumCols += getMatrix(row,col)->getGlobalNumCols();
478  break; // we need only one non-null matrix in a col
479  }
480 
481  return globalNumCols;
482  }
483 
485  size_t getNodeNumRows() const {
486  global_size_t nodeNumRows = 0;
487 
488  for (size_t row = 0; row < Rows(); ++row)
489  for (size_t col = 0; col < Cols(); col++)
490  if (!getMatrix(row,col).is_null()) {
491  nodeNumRows += getMatrix(row,col)->getNodeNumRows();
492  break; // we need only one non-null matrix in a row
493  }
494 
495  return nodeNumRows;
496  }
497 
500  global_size_t globalNumEntries = 0;
501 
502  for (size_t row = 0; row < Rows(); ++row)
503  for (size_t col = 0; col < Cols(); ++col)
504  if (!getMatrix(row,col).is_null())
505  globalNumEntries += getMatrix(row,col)->getGlobalNumEntries();
506 
507  return globalNumEntries;
508  }
509 
511  size_t getNodeNumEntries() const {
512  global_size_t nodeNumEntries = 0;
513 
514  for (size_t row = 0; row < Rows(); ++row)
515  for (size_t col = 0; col < Cols(); ++col)
516  if (!getMatrix(row,col).is_null())
517  nodeNumEntries += getMatrix(row,col)->getNodeNumEntries();
518 
519  return nodeNumEntries;
520  }
521 
523 
524  size_t getNumEntriesInLocalRow(LocalOrdinal localRow) const {
525  throw Xpetra::Exceptions::RuntimeError("getNumEntriesInLocalRow not supported by BlockedCrsMatrix");
526  }
527 
529 
532  throw Xpetra::Exceptions::RuntimeError("getGlobalNumDiags() not supported by BlockedCrsMatrix");
533  }
534 
536 
538  size_t getNodeNumDiags() const {
539  throw Xpetra::Exceptions::RuntimeError("getNodeNumDiags() not supported by BlockedCrsMatrix");
540  }
541 
543 
545  size_t getGlobalMaxNumRowEntries() const {
546  throw Xpetra::Exceptions::RuntimeError("getGlobalMaxNumRowEntries() not supported by BlockedCrsMatrix");
547  }
548 
550 
552  size_t getNodeMaxNumRowEntries() const {
553  throw Xpetra::Exceptions::RuntimeError("getNodeMaxNumRowEntries() not supported by BlockedCrsMatrix");
554  }
555 
557 
560  bool isLocallyIndexed() const {
561  for (size_t i = 0; i < blocks_.size(); ++i)
562  if (blocks_[i] != Teuchos::null && !blocks_[i]->isLocallyIndexed())
563  return false;
564  return true;
565  }
566 
568 
571  bool isGloballyIndexed() const {
572  for (size_t i = 0; i < blocks_.size(); i++)
573  if (blocks_[i] != Teuchos::null && !blocks_[i]->isGloballyIndexed())
574  return false;
575  return true;
576  }
577 
579  bool isFillComplete() const {
580  for (size_t i = 0; i < blocks_.size(); i++)
581  if (blocks_[i] != Teuchos::null && !blocks_[i]->isFillComplete())
582  return false;
583  return true;
584  }
585 
587 
601  virtual void getLocalRowCopy(LocalOrdinal LocalRow,
602  const ArrayView<LocalOrdinal>& Indices,
603  const ArrayView<Scalar>& Values,
604  size_t &NumEntries) const {
605  throw Xpetra::Exceptions::RuntimeError("getLocalRowCopy not supported by BlockedCrsMatrix" );
606  }
607 
609 
618  void getGlobalRowView(GlobalOrdinal GlobalRow, ArrayView<const GlobalOrdinal>& indices, ArrayView<const Scalar>& values) const {
619  throw Xpetra::Exceptions::RuntimeError("getGlobalRowView not supported by BlockedCrsMatrix");
620  }
621 
623 
632  void getLocalRowView(LocalOrdinal LocalRow, ArrayView<const LocalOrdinal>& indices, ArrayView<const Scalar>& values) const {
633  throw Xpetra::Exceptions::RuntimeError("getLocalRowView not supported by BlockedCrsMatrix");
634  }
635 
637 
640  void getLocalDiagCopy(Vector& diag) const {
641  throw Xpetra::Exceptions::RuntimeError("getLocalDiagCopy not supported by BlockedCrsMatrix" );
642  }
643 
646  throw Xpetra::Exceptions::RuntimeError("getFrobeniusNorm() not supported by BlockedCrsMatrix, yet");
647  }
648 
650 
652 
653 
655 
675 
677 
678 
680 
683  virtual void apply(const MultiVector& X, MultiVector& Y,
685  Scalar alpha = ScalarTraits<Scalar>::one(),
686  Scalar beta = ScalarTraits<Scalar>::zero()) const
687  {
688  using Teuchos::RCP;
689 
691  "apply() only supports the following modes: NO_TRANS and TRANS." );
692 
693  RCP<const MultiVector> refX = rcpFromRef(X);
695 
696  SC one = ScalarTraits<SC>::one();
697 
698  if (mode == Teuchos::NO_TRANS) {
699  for (size_t row = 0; row < Rows(); row++) {
700  RCP<MultiVector> Yblock = rangemaps_->getVector(row, Y.getNumVectors(), bRangeThyraMode_);
701  RCP<MultiVector> tmpYblock = rangemaps_->getVector(row, Y.getNumVectors(), bRangeThyraMode_);
702 
703  for (size_t col = 0; col < Cols(); col++) {
704  RCP<const MultiVector> Xblock = domainmaps_->ExtractVector(refX, col, bDomainThyraMode_);
705  RCP<CrsMatrix> Ablock = getMatrix(row, col);
706 
707  if (Ablock.is_null())
708  continue;
709 
710  Ablock->apply(*Xblock, *tmpYblock);
711  Yblock->update(one, *tmpYblock, one);
712  }
713  rangemaps_->InsertVector(Yblock, row, tmpY, bRangeThyraMode_);
714  }
715 
716  } else if (mode == Teuchos::TRANS) {
717  // TODO: test me!
718  for (size_t col = 0; col < Cols(); col++) {
719  RCP<MultiVector> Yblock = domainmaps_->getVector(col, Y.getNumVectors(), bDomainThyraMode_);
720  RCP<MultiVector> tmpYblock = domainmaps_->getVector(col, Y.getNumVectors(), bDomainThyraMode_);
721 
722  for (size_t row = 0; row<Rows(); row++) {
723  RCP<const MultiVector> Xblock = rangemaps_->ExtractVector(refX, row, bRangeThyraMode_);
724  RCP<CrsMatrix> Ablock = getMatrix(row, col);
725 
726  if (Ablock.is_null())
727  continue;
728 
729  Ablock->apply(*Xblock, *tmpYblock, Teuchos::TRANS);
730 
731  Yblock->update(one, *tmpYblock, one);
732  }
733  domainmaps_->InsertVector(Yblock, col, tmpY, bDomainThyraMode_);
734  }
735  }
736 
737  Y.update(alpha, *tmpY, beta);
738  }
739 
742  RCP<const Map > getDomainMap() const { return domainmaps_->getFullMap(); }
743 
746  RCP<const Map > getDomainMap(size_t i, bool bThyraMode = false) const { return domainmaps_->getMap(i, bDomainThyraMode_); }
747 
750  RCP<const Map > getRangeMap() const { return rangemaps_->getFullMap(); }
751 
754  RCP<const Map > getRangeMap(size_t i, bool bThyraMode = false) const { return rangemaps_->getMap(i, bRangeThyraMode_); }
755 
758 
761 
763 
765  //{@
766 
769  throw Xpetra::Exceptions::RuntimeError("BlockedCrsMatrix::getMap(): operation not supported.");
770  }
771 
773  void doImport(const Matrix &source, const Import& importer, CombineMode CM) {
774  throw Xpetra::Exceptions::RuntimeError("BlockedCrsMatrix::doImport(): operation not supported.");
775  }
776 
778  void doExport(const Matrix& dest, const Import& importer, CombineMode CM) {
779  throw Xpetra::Exceptions::RuntimeError("BlockedCrsMatrix::doExport(): operation not supported.");
780  }
781 
783  void doImport(const Matrix& source, const Export& exporter, CombineMode CM) {
784  throw Xpetra::Exceptions::RuntimeError("BlockedCrsMatrix::doImport(): operation not supported.");
785  }
786 
788  void doExport(const Matrix& dest, const Export& exporter, CombineMode CM) {
789  throw Xpetra::Exceptions::RuntimeError("BlockedCrsMatrix::doExport(): operation not supported.");
790  }
791 
792  // @}
793 
795 
796 
798  std::string description() const { return "Xpetra_BlockedCrsMatrix.description()"; }
799 
802  out << "Xpetra::BlockedCrsMatrix: " << Rows() << " x " << Cols() << std::endl;
803 
804  if (isFillComplete()) {
805  out << "BlockMatrix is fillComplete" << std::endl;
806 
807  out << "fullRowMap" << std::endl;
808  fullrowmap_->describe(out,verbLevel);
809 
810  //out << "fullColMap" << std::endl;
811  //fullcolmap_->describe(out,verbLevel);
812 
813  } else {
814  out << "BlockMatrix is NOT fillComplete" << std::endl;
815  }
816 
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;
820  getMatrix(r,c)->describe(out,verbLevel);
821  }
822  }
823 
826  throw Xpetra::Exceptions::RuntimeError("getCrsGraph() not supported by BlockedCrsMatrix");
827  }
828 
830 
832 
833 
835  virtual size_t Rows() const { return rangemaps_->NumMaps(); }
836 
838  virtual size_t Cols() const { return domainmaps_->NumMaps(); }
839 
841  Teuchos::RCP<CrsMatrix> getMatrix(size_t r, size_t c) const { return blocks_[r*Cols()+c]; }
842 
844  void setMatrix(size_t r, size_t c, Teuchos::RCP<CrsMatrix> mat) {
845  // TODO: if filled -> return error
846 
847  TEUCHOS_TEST_FOR_EXCEPTION(r > Rows(), std::out_of_range, "Error, r = " << Rows() << " is too big");
848  TEUCHOS_TEST_FOR_EXCEPTION(c > Cols(), std::out_of_range, "Error, c = " << Cols() << " is too big");
849 
850  // set matrix
851  blocks_[r*Cols() + c] = mat;
852  }
853 
855  // NOTE: This is a rather expensive operation, since all blocks are copied
856  // into a new big CrsMatrix
858  using Teuchos::RCP;
859  using Teuchos::rcp_dynamic_cast;
860  Scalar one = ScalarTraits<SC>::one();
861 
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." );
864 
866 
867  for (size_t i = 0; i < blocks_.size(); i++)
868  if (blocks_[i] != Teuchos::null)
869  this->Add(*blocks_[i], one, *sparse, one);
870 
871  sparse->fillComplete(getDomainMap(), getRangeMap());
872 
873  return sparse;
874  }
876 
877 #ifdef HAVE_XPETRA_KOKKOS_REFACTOR
878  typedef typename CrsMatrix::local_matrix_type local_matrix_type;
879 
881  local_matrix_type getLocalMatrix () const {
882  throw Xpetra::Exceptions::RuntimeError("BlockedCrsMatrix::getLocalMatrix(): operation not supported.");
883  }
884 #endif
885 
886 #ifdef HAVE_XPETRA_THYRA
889  Xpetra::ThyraUtils<Scalar,LO,GO,Node>::toThyra(this);
891 
893  Teuchos::rcp_dynamic_cast<Thyra::BlockedLinearOpBase<Scalar> >(thOp);
895  return thbOp;
896  }
897 #endif
898 
899  private:
900 
903 
905 
915  void Add(const CrsMatrix& A, const Scalar scalarA, CrsMatrix& B, const Scalar scalarB) const {
917  "Matrix A is not completed");
918  using Teuchos::Array;
919  using Teuchos::ArrayView;
920 
921  B.scale(scalarB);
922 
923  Scalar one = ScalarTraits<SC>::one();
924  Scalar zero = ScalarTraits<SC>::zero();
925 
926  if (scalarA == zero)
927  return;
928 
929  size_t maxNumEntries = A.getNodeMaxNumRowEntries();
930 
931  size_t numEntries;
932  Array<GO> inds(maxNumEntries);
933  Array<SC> vals(maxNumEntries);
934 
935  RCP<const Map> rowMap = A.getRowMap();
936  RCP<const Map> colMap = A.getColMap();
937 
938  ArrayView<const GO> rowGIDs = A.getRowMap()->getNodeElementList();
939  for (size_t i = 0; i < A.getNodeNumRows(); i++) {
940  GO row = rowGIDs[i];
941  A.getGlobalRowCopy(row, inds(), vals(), numEntries);
942 
943  if (scalarA != one)
944  for (size_t j = 0; j < numEntries; ++j)
945  vals[j] *= scalarA;
946 
947  B.insertGlobalValues(row, inds(0, numEntries), vals(0, numEntries)); // insert should be ok, since blocks in BlockedCrsOpeartor do not overlap!
948  }
949  }
950 
952 
953  // Default view is created after fillComplete()
954  // Because ColMap might not be available before fillComplete().
956  // Create default view
957  this->defaultViewLabel_ = "point";
959 
960  // Set current view
961  this->currentViewLabel_ = this->GetDefaultViewLabel();
962  }
963 
964  private:
965  Teuchos::RCP<const MapExtractor> domainmaps_; // full domain map together with all partial domain maps
966  Teuchos::RCP<const MapExtractor> rangemaps_; // full range map together with all partial domain maps
967  Teuchos::RCP<Map> fullrowmap_; // full matrix row map
968  //Teuchos::RCP<Map> fullcolmap_; // full matrix column map
969 
970  std::vector<Teuchos::RCP<CrsMatrix> > blocks_; // row major matrix block storage
971 #ifdef HAVE_XPETRA_THYRA
973 #endif
976 
977 };
978 
979 } //namespace Xpetra
980 
981 #define XPETRA_BLOCKEDCRSMATRIX_SHORT
982 #endif /* XPETRA_BLOCKEDCRSMATRIX_HPP */
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&#39;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 > &params=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 > &params=null)=0
Signal that data entry is complete, specifying domain and range maps.
GlobalOrdinal GO
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.
Xpetra namespace
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.
size_type size() const
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.
RCP< const Map > getDomainMap(size_t i, bool bThyraMode=false) const
Returns the Map associated with the i&#39;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.
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 > &params=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_
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.
T * getRawPtr() const
std::string viewLabel_t
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 > &params=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...
bool is_null() const