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 
51 #include <Teuchos_SerialDenseMatrix.hpp>
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"
61 #include "Xpetra_CrsGraph.hpp"
62 #include "Xpetra_CrsMatrix.hpp"
64 
65 #include "Xpetra_MapExtractor.hpp"
66 
67 #include "Xpetra_Matrix.hpp"
68 #include "Xpetra_MatrixFactory.hpp"
69 #include "Xpetra_CrsMatrixWrap.hpp"
70 
71 #ifdef HAVE_XPETRA_THYRA
72 #include <Thyra_ProductVectorSpaceBase.hpp>
73 #include <Thyra_VectorSpaceBase.hpp>
74 #include <Thyra_LinearOpBase.hpp>
75 #include <Thyra_BlockedLinearOpBase.hpp>
76 #include <Thyra_PhysicallyBlockedLinearOpBase.hpp>
77 #include "Xpetra_ThyraUtils.hpp"
78 #endif
79 
80 
85 namespace Xpetra {
86 
87  typedef std::string viewLabel_t;
88 
89  template <class Scalar,
90  class LocalOrdinal, // = typename Matrix<Scalar>::local_ordinal_type,
91  class GlobalOrdinal, // = typename Matrix<Scalar, LocalOrdinal>::global_ordinal_type,
92  class Node> // = 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 
117  BlockedCrsMatrix(Teuchos::RCP<const MapExtractor>& rangeMaps,
118  Teuchos::RCP<const MapExtractor>& domainMaps,
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(MatrixFactory::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.
161  Teuchos::RCP<const Xpetra::CrsMatrix<Scalar,LO,GO,Node> > xop =
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.
187  Teuchos::RCP<const Xpetra::CrsMatrix<Scalar,LO,GO,Node> > xop =
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  // TODO we do not support nested Thyra operators here!
214  Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > const_op = thyraOp->getBlock(r,c); // nonConst access is not allowed.
215  Teuchos::RCP<Thyra::LinearOpBase<Scalar> > op = Teuchos::rcp_const_cast<Thyra::LinearOpBase<Scalar> >(const_op); // cast away const
216  Teuchos::RCP<Xpetra::CrsMatrix<Scalar,LO,GO,Node> > xop =
218  Teuchos::RCP<Xpetra::CrsMatrixWrap<Scalar,LO,GO,Node> > xwrap =
219  Teuchos::rcp(new Xpetra::CrsMatrixWrap<Scalar,LO,GO,Node>(xop));
220  blocks_.push_back(xwrap);
221  } else {
222  // add empty block
224  }
225  }
226  }
227  // Default view
229  }
230 
231  private:
233 
240  Teuchos::RCP<const Map> mergeMaps(std::vector<Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > > & subMaps) {
241  // TODO merging for Thyra mode is missing (similar to what we do in constructor of MapExtractor
242 
243  // merge submaps to global map
244  std::vector<GlobalOrdinal> gids;
245  for(size_t tt = 0; tt<subMaps.size(); ++tt) {
246  Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > subMap = subMaps[tt];
247  for(LocalOrdinal l = 0; l < Teuchos::as<LocalOrdinal>(subMap->getNodeNumElements()); ++l) {
248  GlobalOrdinal gid = subMap->getGlobalElement(l);
249  gids.push_back(gid);
250  }
251  }
252 
253  // we have to sort the matrix entries and get rid of the double entries
254  // since we use this to detect Thyra-style numbering or Xpetra-style
255  // numbering. In Thyra-style numbering mode, the Xpetra::MapExtractor builds
256  // the correct row maps.
257  const GO INVALID = Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid();
258  std::sort(gids.begin(), gids.end());
259  gids.erase(std::unique(gids.begin(), gids.end()), gids.end());
260  Teuchos::ArrayView<GO> gidsView(&gids[0], gids.size());
261  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());
262  return fullMap;
263  }
264 
265  public:
266 #endif
267 
269  virtual ~BlockedCrsMatrix() {}
270 
272 
273 
275 
276 
278 
303  void insertGlobalValues(GlobalOrdinal globalRow, const ArrayView<const GlobalOrdinal>& cols, const ArrayView<const Scalar>& vals) {
304  XPETRA_MONITOR("XpetraBlockedCrsMatrix::insertGlobalValues");
305  if (Rows() == 1 && Cols () == 1) {
306  getMatrix(0,0)->insertGlobalValues(globalRow, cols, vals);
307  return;
308  }
309  throw Xpetra::Exceptions::RuntimeError("insertGlobalValues not supported by BlockedCrsMatrix");
310  }
311 
313 
323  void insertLocalValues(LocalOrdinal localRow, const ArrayView<const LocalOrdinal>& cols, const ArrayView<const Scalar>& vals) {
324  XPETRA_MONITOR("XpetraBlockedCrsMatrix::insertLocalValues");
325  if (Rows() == 1 && Cols () == 1) {
326  getMatrix(0,0)->insertLocalValues(localRow, cols, vals);
327  return;
328  }
329  throw Xpetra::Exceptions::RuntimeError("insertLocalValues not supported by BlockedCrsMatrix");
330  }
331 
332  void removeEmptyProcessesInPlace(const Teuchos::RCP<const Map>& newMap) {
333  XPETRA_MONITOR("XpetraBlockedCrsMatrix::removeEmptyProcessesInPlace");
334  if (Rows() == 1 && Cols () == 1) {
335  getMatrix(0,0)->removeEmptyProcessesInPlace(newMap);
336  return;
337  }
338  throw Xpetra::Exceptions::RuntimeError("removeEmptyProcesses not supported by BlockedCrsMatrix");
339  }
340 
342 
350  void replaceGlobalValues(GlobalOrdinal globalRow,
351  const ArrayView<const GlobalOrdinal> &cols,
352  const ArrayView<const Scalar> &vals) {
353  XPETRA_MONITOR("XpetraBlockedCrsMatrix::replaceGlobalValues");
354  if (Rows() == 1 && Cols () == 1) {
355  getMatrix(0,0)->replaceGlobalValues(globalRow,cols,vals);
356  return;
357  }
358  throw Xpetra::Exceptions::RuntimeError("replaceGlobalValues not supported by BlockedCrsMatrix");
359  }
360 
362 
366  void replaceLocalValues(LocalOrdinal localRow,
367  const ArrayView<const LocalOrdinal> &cols,
368  const ArrayView<const Scalar> &vals) {
369  XPETRA_MONITOR("XpetraBlockedCrsMatrix::replaceLocalValues");
370  if (Rows() == 1 && Cols () == 1) {
371  getMatrix(0,0)->replaceLocalValues(localRow,cols,vals);
372  return;
373  }
374  throw Xpetra::Exceptions::RuntimeError("replaceLocalValues not supported by BlockedCrsMatrix");
375  }
376 
378  virtual void setAllToScalar(const Scalar& alpha) {
379  XPETRA_MONITOR("XpetraBlockedCrsMatrix::setAllToScalar");
380  for (size_t row = 0; row < Rows(); row++) {
381  for (size_t col = 0; col < Cols(); col++) {
382  if (!getMatrix(row,col).is_null()) {
383  getMatrix(row,col)->setAllToScalar(alpha);
384  }
385  }
386  }
387  }
388 
390  void scale(const Scalar& alpha) {
391  XPETRA_MONITOR("XpetraBlockedCrsMatrix::scale");
392  for (size_t row = 0; row < Rows(); row++) {
393  for (size_t col = 0; col < Cols(); col++) {
394  if (!getMatrix(row,col).is_null()) {
395  getMatrix(row,col)->scale(alpha);
396  }
397  }
398  }
399  }
400 
402 
404 
405 
413  void resumeFill(const RCP< ParameterList >& params = null) {
414  XPETRA_MONITOR("XpetraBlockedCrsMatrix::resumeFill");
415  for (size_t row = 0; row < Rows(); row++) {
416  for (size_t col = 0; col < Cols(); col++) {
417  if (!getMatrix(row,col).is_null()) {
418  getMatrix(row,col)->resumeFill(params);
419  }
420  }
421  }
422  }
423 
429  void fillComplete(const RCP<const Map>& domainMap, const RCP<const Map>& rangeMap, const RCP<ParameterList>& params = null) {
430  XPETRA_MONITOR("XpetraBlockedCrsMatrix::fillComplete");
431  if (Rows() == 1 && Cols () == 1) {
432  getMatrix(0,0)->fillComplete(domainMap, rangeMap, params);
433  return;
434  }
435  fillComplete(params);
436  }
437 
452  void fillComplete(const RCP<ParameterList>& params = null) {
453  XPETRA_MONITOR("XpetraBlockedCrsMatrix::fillComplete");
454  TEUCHOS_TEST_FOR_EXCEPTION(rangemaps_==Teuchos::null, Xpetra::Exceptions::RuntimeError,"BlockedCrsMatrix::fillComplete: rangemaps_ is not set. Error.");
455 
456  for (size_t r = 0; r < Rows(); ++r)
457  for (size_t c = 0; c < Cols(); ++c) {
458  if(getMatrix(r,c) != Teuchos::null) {
459  Teuchos::RCP<Matrix> Ablock = getMatrix(r,c);
460 
461  if (Ablock != Teuchos::null && !Ablock->isFillComplete())
462  Ablock->fillComplete(getDomainMap(c, bDomainThyraMode_), getRangeMap(r, bRangeThyraMode_), params);
463  }
464  }
465 
466 #if 0
467  // get full row map
468  RCP<const Map> rangeMap = rangemaps_->getFullMap();
469  fullrowmap_ = MapFactory::Build(rangeMap()->lib(), rangeMap()->getGlobalNumElements(), rangeMap()->getNodeElementList(), rangeMap()->getIndexBase(), rangeMap()->getComm());
470 
471  // build full col map
472  fullcolmap_ = Teuchos::null; // delete old full column map
473 
474  std::vector<GO> colmapentries;
475  for (size_t c = 0; c < Cols(); ++c) {
476  // copy all local column lids of all block rows to colset
477  std::set<GO> colset;
478  for (size_t r = 0; r < Rows(); ++r) {
479  Teuchos::RCP<CrsMatrix> Ablock = getMatrix(r,c);
480 
481  if (Ablock != Teuchos::null) {
482  Teuchos::ArrayView<const GO> colElements = Ablock->getColMap()->getNodeElementList();
483  Teuchos::RCP<const Map> colmap = Ablock->getColMap();
484  copy(colElements.getRawPtr(), colElements.getRawPtr() + colElements.size(), inserter(colset, colset.begin()));
485  }
486  }
487 
488  // remove duplicates (entries which are in column maps of more than one block row)
489  colmapentries.reserve(colmapentries.size() + colset.size());
490  copy(colset.begin(), colset.end(), back_inserter(colmapentries));
491  sort(colmapentries.begin(), colmapentries.end());
492  typename std::vector<GO>::iterator gendLocation;
493  gendLocation = std::unique(colmapentries.begin(), colmapentries.end());
494  colmapentries.erase(gendLocation,colmapentries.end());
495  }
496 
497  // sum up number of local elements
498  size_t numGlobalElements = 0;
499  Teuchos::reduceAll(*(rangeMap->getComm()), Teuchos::REDUCE_SUM, colmapentries.size(), Teuchos::outArg(numGlobalElements));
500 
501  // store global full column map
502  const Teuchos::ArrayView<const GO> aView = Teuchos::ArrayView<const GO>(colmapentries);
503  fullcolmap_ = MapFactory::Build(rangeMap->lib(), numGlobalElements, aView, 0, rangeMap->getComm());
504 #endif
505  }
506 
508 
510 
513  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getGlobalNumRows");
514  global_size_t globalNumRows = 0;
515 
516  for (size_t row = 0; row < Rows(); row++)
517  for (size_t col = 0; col < Cols(); col++)
518  if (!getMatrix(row,col).is_null()) {
519  globalNumRows += getMatrix(row,col)->getGlobalNumRows();
520  break; // we need only one non-null matrix in a row
521  }
522 
523  return globalNumRows;
524  }
525 
527 
530  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getGlobalNumCols");
531  global_size_t globalNumCols = 0;
532 
533  for (size_t col = 0; col < Cols(); col++)
534  for (size_t row = 0; row < Rows(); row++)
535  if (!getMatrix(row,col).is_null()) {
536  globalNumCols += getMatrix(row,col)->getGlobalNumCols();
537  break; // we need only one non-null matrix in a col
538  }
539 
540  return globalNumCols;
541  }
542 
544  size_t getNodeNumRows() const {
545  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getNodeNumRows");
546  global_size_t nodeNumRows = 0;
547 
548  for (size_t row = 0; row < Rows(); ++row)
549  for (size_t col = 0; col < Cols(); col++)
550  if (!getMatrix(row,col).is_null()) {
551  nodeNumRows += getMatrix(row,col)->getNodeNumRows();
552  break; // we need only one non-null matrix in a row
553  }
554 
555  return nodeNumRows;
556  }
557 
560  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getGlobalNumEntries");
561  global_size_t globalNumEntries = 0;
562 
563  for (size_t row = 0; row < Rows(); ++row)
564  for (size_t col = 0; col < Cols(); ++col)
565  if (!getMatrix(row,col).is_null())
566  globalNumEntries += getMatrix(row,col)->getGlobalNumEntries();
567 
568  return globalNumEntries;
569  }
570 
572  size_t getNodeNumEntries() const {
573  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getNodeNumEntries");
574  global_size_t nodeNumEntries = 0;
575 
576  for (size_t row = 0; row < Rows(); ++row)
577  for (size_t col = 0; col < Cols(); ++col)
578  if (!getMatrix(row,col).is_null())
579  nodeNumEntries += getMatrix(row,col)->getNodeNumEntries();
580 
581  return nodeNumEntries;
582  }
583 
585 
586  size_t getNumEntriesInLocalRow(LocalOrdinal localRow) const {
587  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getNumEntriesInLocalRow");
588  if (Rows() == 1 && Cols () == 1) {
589  return getMatrix(0,0)->getNumEntriesInLocalRow(localRow);
590  }
591  throw Xpetra::Exceptions::RuntimeError("getNumEntriesInLocalRow not supported by BlockedCrsMatrix");
592  }
593 
595 
598  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getGlobalNumDiags");
599  if (Rows() == 1 && Cols () == 1) {
600  return getMatrix(0,0)->getGlobalNumDiags();
601  }
602  throw Xpetra::Exceptions::RuntimeError("getGlobalNumDiags() not supported by BlockedCrsMatrix");
603  }
604 
606 
608  size_t getNodeNumDiags() const {
609  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getNodeNumDiags");
610  if (Rows() == 1 && Cols () == 1) {
611  return getMatrix(0,0)->getNodeNumDiags();
612  }
613  throw Xpetra::Exceptions::RuntimeError("getNodeNumDiags() not supported by BlockedCrsMatrix");
614  }
615 
617 
619  size_t getGlobalMaxNumRowEntries() const {
620  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getGlobalMaxNumRowEntries");
621  global_size_t globalMaxEntries = 0;
622 
623  for (size_t row = 0; row < Rows(); row++) {
624  global_size_t globalMaxEntriesBlockRows = 0;
625  for (size_t col = 0; col < Cols(); col++) {
626  if (!getMatrix(row,col).is_null()) {
627  globalMaxEntriesBlockRows += getMatrix(row,col)->getGlobalMaxNumRowEntries();
628  }
629  }
630  if(globalMaxEntriesBlockRows > globalMaxEntries)
631  globalMaxEntries = globalMaxEntriesBlockRows;
632  }
633  return globalMaxEntries;
634  }
635 
637 
639  size_t getNodeMaxNumRowEntries() const {
640  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getNodeMaxNumRowEntries");
641  size_t localMaxEntries = 0;
642 
643  for (size_t row = 0; row < Rows(); row++) {
644  size_t localMaxEntriesBlockRows = 0;
645  for (size_t col = 0; col < Cols(); col++) {
646  if (!getMatrix(row,col).is_null()) {
647  localMaxEntriesBlockRows += getMatrix(row,col)->getNodeMaxNumRowEntries();
648  }
649  }
650  if(localMaxEntriesBlockRows > localMaxEntries)
651  localMaxEntries = localMaxEntriesBlockRows;
652  }
653  return localMaxEntries;
654  }
655 
657 
660  bool isLocallyIndexed() const {
661  XPETRA_MONITOR("XpetraBlockedCrsMatrix::isLocallyIndexed");
662  for (size_t i = 0; i < blocks_.size(); ++i)
663  if (blocks_[i] != Teuchos::null && !blocks_[i]->isLocallyIndexed())
664  return false;
665  return true;
666  }
667 
669 
672  bool isGloballyIndexed() const {
673  XPETRA_MONITOR("XpetraBlockedCrsMatrix::isGloballyIndexed");
674  for (size_t i = 0; i < blocks_.size(); i++)
675  if (blocks_[i] != Teuchos::null && !blocks_[i]->isGloballyIndexed())
676  return false;
677  return true;
678  }
679 
681  bool isFillComplete() const {
682  XPETRA_MONITOR("XpetraBlockedCrsMatrix::isFillComplete");
683  for (size_t i = 0; i < blocks_.size(); i++)
684  if (blocks_[i] != Teuchos::null && !blocks_[i]->isFillComplete())
685  return false;
686  return true;
687  }
688 
690 
704  virtual void getLocalRowCopy(LocalOrdinal LocalRow,
705  const ArrayView<LocalOrdinal>& Indices,
706  const ArrayView<Scalar>& Values,
707  size_t &NumEntries) const {
708  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getLocalRowCopy");
709  if (Rows() == 1 && Cols () == 1) {
710  getMatrix(0,0)->getLocalRowCopy(LocalRow, Indices, Values, NumEntries);
711  return;
712  }
713  throw Xpetra::Exceptions::RuntimeError("getLocalRowCopy not supported by BlockedCrsMatrix" );
714  }
715 
717 
726  void getGlobalRowView(GlobalOrdinal GlobalRow, ArrayView<const GlobalOrdinal>& indices, ArrayView<const Scalar>& values) const {
727  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getGlobalRowView");
728  if (Rows() == 1 && Cols () == 1) {
729  getMatrix(0,0)->getGlobalRowView(GlobalRow, indices, values);
730  return;
731  }
732  throw Xpetra::Exceptions::RuntimeError("getGlobalRowView not supported by BlockedCrsMatrix");
733  }
734 
736 
745  void getLocalRowView(LocalOrdinal LocalRow, ArrayView<const LocalOrdinal>& indices, ArrayView<const Scalar>& values) const {
746  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getLocalRowView");
747  if (Rows() == 1 && Cols () == 1) {
748  getMatrix(0,0)->getLocalRowView(LocalRow, indices, values);
749  return;
750  }
751  throw Xpetra::Exceptions::RuntimeError("getLocalRowView not supported by BlockedCrsMatrix");
752  }
753 
755 
758  void getLocalDiagCopy(Vector& diag) const {
759  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getLocalDiagCopy");
760  XPETRA_TEST_FOR_EXCEPTION(diag.getMap()->isSameAs(*rangemaps_->getFullMap()) == false, Xpetra::Exceptions::RuntimeError,
761  "BlockedCrsMatrix::getLocalDiagCopy(): the map of the vector diag is not compatible with the full map of the blocked operator." );
762 
763  TEUCHOS_TEST_FOR_EXCEPTION(Rows() != Cols(), Xpetra::Exceptions::RuntimeError,
764  "BlockedCrsMatrix::getLocalDiagCopy(): you cannot extract the diagonal of a "<< Rows() << "x"<<Cols()<<" blocked matrix." );
765 
766  for (size_t row = 0; row < Rows(); ++row) {
767  if (!getMatrix(row,row).is_null()) {
768  // if we are in Thyra mode, but the block (row,row) is again a blocked operator, we have to use (pseudo) Xpetra-style GIDs with offset!
769  bool bThyraMode = rangemaps_->getThyraMode() && (Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(getMatrix(row,row)) == Teuchos::null);
770  RCP<Vector> dd = rangemaps_->getVector(row, bThyraMode);
771  getMatrix(row,row)->getLocalDiagCopy(*dd);
772  rangemaps_->InsertVector(*dd,row,diag,bThyraMode);
773  }
774  }
775  }
776 
778  void leftScale (const Vector& x) {
779  XPETRA_MONITOR("XpetraBlockedCrsMatrix::leftScale");
780  XPETRA_TEST_FOR_EXCEPTION(x.getMap()->isSameAs(*rangemaps_->getFullMap()) == false, Xpetra::Exceptions::RuntimeError,
781  "BlockedCrsMatrix::leftScale(): the map of the vector x is not compatible with the full map of the blocked operator." );
782 
783  RCP<const Vector> rcpx = Teuchos::rcpFromRef(x);
784 
785  for (size_t row = 0; row < Rows(); ++row) {
786  for (size_t col = 0; col < Cols(); ++col) {
787  if(getMatrix(row,col)!=Teuchos::null) {
788  // if we are in Thyra mode, but the block (row,row) is again a blocked operator, we have to use (pseudo) Xpetra-style GIDs with offset!
789  bool bThyraMode = rangemaps_->getThyraMode() && (Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(getMatrix(row,col)) == Teuchos::null);
790  RCP<Vector> xx = rangemaps_->ExtractVector(rcpx,row,bThyraMode);
791  getMatrix(row,col)->leftScale(*xx);
792  }
793  }
794  }
795  }
796 
798  void rightScale (const Vector& x) {
799  XPETRA_MONITOR("XpetraBlockedCrsMatrix::rightScale");
800  XPETRA_TEST_FOR_EXCEPTION(x.getMap()->isSameAs(*domainmaps_->getFullMap()) == false, Xpetra::Exceptions::RuntimeError,
801  "BlockedCrsMatrix::rightScale(): the map of the vector x is not compatible with the full map of the blocked operator." );
802 
803  RCP<const Vector> rcpx = Teuchos::rcpFromRef(x);
804 
805  for (size_t col = 0; col < Cols(); ++col) {
806  for (size_t row = 0; row < Rows(); ++row) {
807  if(getMatrix(row,col)!=Teuchos::null) {
808  // if we are in Thyra mode, but the block (row,row) is again a blocked operator, we have to use (pseudo) Xpetra-style GIDs with offset!
809  bool bThyraMode = domainmaps_->getThyraMode() && (Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(getMatrix(row,col)) == Teuchos::null);
810  RCP<Vector> xx = domainmaps_->ExtractVector(rcpx,col,bThyraMode);
811  getMatrix(row,col)->rightScale(*xx);
812  }
813  }
814  }
815  }
816 
817 
819  virtual typename ScalarTraits<Scalar>::magnitudeType getFrobeniusNorm() const {
820  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getFrobeniusNorm");
821  typename ScalarTraits<Scalar>::magnitudeType ret = Teuchos::ScalarTraits<Scalar>::magnitude(Teuchos::ScalarTraits<Scalar>::zero());
822  for (size_t col = 0; col < Cols(); ++col) {
823  for (size_t row = 0; row < Rows(); ++row) {
824  if(getMatrix(row,col)!=Teuchos::null) {
825  typename ScalarTraits<Scalar>::magnitudeType n = getMatrix(row,col)->getFrobeniusNorm();
826  ret += n * n;
827  }
828  }
829  }
830  return Teuchos::ScalarTraits< typename ScalarTraits<Scalar>::magnitudeType >::squareroot(ret);
831  }
832 
834 
836 
837 
839 
859 
861 
862 
864 
867  virtual void apply(const MultiVector& X, MultiVector& Y,
868  Teuchos::ETransp mode = Teuchos::NO_TRANS,
869  Scalar alpha = ScalarTraits<Scalar>::one(),
870  Scalar beta = ScalarTraits<Scalar>::zero()) const
871  {
872  XPETRA_MONITOR("XpetraBlockedCrsMatrix::apply");
873  //using Teuchos::RCP;
874 
875  TEUCHOS_TEST_FOR_EXCEPTION(mode != Teuchos::NO_TRANS && mode != Teuchos::TRANS, Xpetra::Exceptions::RuntimeError,
876  "apply() only supports the following modes: NO_TRANS and TRANS." );
877 
878  // check whether input parameters are blocked or not
879  RCP<const MultiVector> refX = rcpFromRef(X);
880  RCP<const BlockedMultiVector> refbX = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(refX);
881  RCP<MultiVector> tmpY = rcpFromRef(Y);
882  RCP<BlockedMultiVector> tmpbY = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(tmpY);
883 
884  bool bBlockedX = (refbX != Teuchos::null) ? true : false;
885  bool bBlockedY = (tmpbY != Teuchos::null) ? true : false;
886 
887  // create (temporary) vectors for output
888  // In the end we call Y.update(alpha, *tmpY, beta). Therefore we need a new vector storing the temporary results
889  tmpY = MultiVectorFactory::Build(Y.getMap(), Y.getNumVectors(), true);
890  if (bBlockedY == true) tmpbY = Teuchos::rcp(new BlockedMultiVector(rangemaps_,tmpY));
891 
892  SC one = ScalarTraits<SC>::one();
893 
894  if (mode == Teuchos::NO_TRANS) {
895 
896  for (size_t row = 0; row < Rows(); row++) {
897  RCP<MultiVector> Yblock = rangemaps_->getVector(row, Y.getNumVectors(), bRangeThyraMode_, true);
898  for (size_t col = 0; col < Cols(); col++) {
899 
900  // extract matrix block
901  RCP<Matrix> Ablock = getMatrix(row, col);
902 
903  if (Ablock.is_null())
904  continue;
905 
906  // check whether Ablock is itself a blocked operator
907  // If it is a blocked operator we have to provide Xpetra style GIDs, i.e. we have to transform GIDs
908  bool bBlockedSubMatrix = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(Ablock) == Teuchos::null ? false : true;
909 
910  // input/output vectors for local block operation
911  RCP<const MultiVector> Xblock = Teuchos::null; // subpart of X vector to be applied to subblock of A
912 
913  // extract sub part of X using Xpetra or Thyra GIDs
914  // if submatrix is again blocked, we extract it using Xpetra style gids. If it is a single
915  // block matrix we use the Thyra or Xpetra style GIDs that are used to store the matrix
916  if(bBlockedX) Xblock = domainmaps_->ExtractVector(refbX, col, bBlockedSubMatrix == true ? false : bDomainThyraMode_);
917  else Xblock = domainmaps_->ExtractVector(refX, col, bBlockedSubMatrix == true ? false : bDomainThyraMode_);
918  RCP<MultiVector> tmpYblock = rangemaps_->getVector(row, Y.getNumVectors(), bBlockedSubMatrix == true ? false : bRangeThyraMode_, false); // subpart of Y vector containing part of solution of Xblock applied to Ablock
919  Ablock->apply(*Xblock, *tmpYblock);
920 
921  // If Ablock is a blocked operator the local vectors are using (pseudo) Xpetra-style gids
922  // that have to be translated to Thyra based GIDs if bRangeThyraMode is set
923  if(bBlockedSubMatrix == true && bRangeThyraMode_ == true) {
924  tmpYblock->replaceMap(rangemaps_->getMap(row, true)); // switch to Thyra maps (compatible to Yblock)
925  }
926  Yblock->update(one, *tmpYblock, one);
927  }
928  if(bBlockedY) rangemaps_->InsertVector(Yblock, row, tmpbY, bRangeThyraMode_);
929  else rangemaps_->InsertVector(Yblock, row, tmpY, bRangeThyraMode_);
930  }
931 
932  } else if (mode == Teuchos::TRANS) {
933  // TODO: test me!
934  for (size_t col = 0; col < Cols(); col++) {
935  RCP<MultiVector> Yblock = domainmaps_->getVector(col, Y.getNumVectors(), bDomainThyraMode_, true);
936 
937  for (size_t row = 0; row<Rows(); row++) {
938  RCP<Matrix> Ablock = getMatrix(row, col);
939 
940  if (Ablock.is_null())
941  continue;
942 
943  // check whether Ablock is itself a blocked operator
944  // If it is a blocked operator we have to provide Xpetra style GIDs, i.e. we have to transform GIDs
945  bool bBlockedSubMatrix = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(Ablock) == Teuchos::null ? false : true;
946 
947  RCP<const MultiVector> Xblock = Teuchos::null;
948 
949  // extract sub part of X using Xpetra or Thyra GIDs
950  if(bBlockedX) Xblock = rangemaps_->ExtractVector(refbX, row, bBlockedSubMatrix == true ? false : bRangeThyraMode_);
951  else Xblock = rangemaps_->ExtractVector(refX, row, bBlockedSubMatrix == true ? false : bRangeThyraMode_);
952  RCP<MultiVector> tmpYblock = domainmaps_->getVector(col, Y.getNumVectors(), bBlockedSubMatrix == true ? false : bDomainThyraMode_, false);
953  Ablock->apply(*Xblock, *tmpYblock, Teuchos::TRANS);
954 
955  // If Ablock is a blocked operator the local vectors are using (pseudo) Xpetra-style gids
956  // that have to be translated to Thyra based GIDs if bRangeThyraMode is set
957  if(bBlockedSubMatrix == true && bDomainThyraMode_ == true) {
958  tmpYblock->replaceMap(domainmaps_->getMap(col, true)); // switch to Thyra maps (compatible to Yblock)
959  }
960 
961  Yblock->update(one, *tmpYblock, one);
962  }
963  if(bBlockedY) domainmaps_->InsertVector(Yblock, col, tmpbY, bDomainThyraMode_);
964  else domainmaps_->InsertVector(Yblock, col, tmpY, bDomainThyraMode_);
965  }
966  }
967  if(bBlockedY) Y.update(alpha, *tmpbY, beta);
968  else Y.update(alpha, *tmpY, beta);
969  }
970 
972  RCP<const Map > getDomainMap() const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getDomainMap()"); return domainmaps_->getFullMap(); }
973 
975  RCP<const Map > getDomainMap(size_t i) const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getDomainMap(size_t)"); return domainmaps_->getMap(i, bDomainThyraMode_); }
976 
978  RCP<const Map > getDomainMap(size_t i, bool bThyraMode) const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getDomainMap(size_t,bool)"); return domainmaps_->getMap(i, bThyraMode); }
979 
981  RCP<const Map > getRangeMap() const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getRangeMap()"); return rangemaps_->getFullMap(); }
982 
984  RCP<const Map > getRangeMap(size_t i) const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getRangeMap(size_t)"); return rangemaps_->getMap(i, bRangeThyraMode_); }
985 
987  RCP<const Map > getRangeMap(size_t i, bool bThyraMode) const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getRangeMap(size_t,bool)"); return rangemaps_->getMap(i, bThyraMode); }
988 
990  RCP<const MapExtractor> getRangeMapExtractor() const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getRangeMapExtractor()"); return rangemaps_; }
991 
993  RCP<const MapExtractor> getDomainMapExtractor() const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getDomainMapExtractor()"); return domainmaps_; }
994 
996 
998  //{@
999 
1001  const Teuchos::RCP< const Map > getMap() const {
1002  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getMap");
1003  if (Rows() == 1 && Cols () == 1) {
1004  return getMatrix(0,0)->getMap();
1005  }
1006  throw Xpetra::Exceptions::RuntimeError("BlockedCrsMatrix::getMap(): operation not supported.");
1007  }
1008 
1010  void doImport(const Matrix &source, const Import& importer, CombineMode CM) {
1011  XPETRA_MONITOR("XpetraBlockedCrsMatrix::doImport");
1012  if (Rows() == 1 && Cols () == 1) {
1013  getMatrix(0,0)->doImport(source, importer, CM);
1014  return;
1015  }
1016  throw Xpetra::Exceptions::RuntimeError("BlockedCrsMatrix::doImport(): operation not supported.");
1017  }
1018 
1020  void doExport(const Matrix& dest, const Import& importer, CombineMode CM) {
1021  XPETRA_MONITOR("XpetraBlockedCrsMatrix::doExport");
1022  if (Rows() == 1 && Cols () == 1) {
1023  getMatrix(0,0)->doExport(dest, importer, CM);
1024  return;
1025  }
1026  throw Xpetra::Exceptions::RuntimeError("BlockedCrsMatrix::doExport(): operation not supported.");
1027  }
1028 
1030  void doImport(const Matrix& source, const Export& exporter, CombineMode CM) {
1031  XPETRA_MONITOR("XpetraBlockedCrsMatrix::doImport");
1032  if (Rows() == 1 && Cols () == 1) {
1033  getMatrix(0,0)->doImport(source, exporter, CM);
1034  return;
1035  }
1036  throw Xpetra::Exceptions::RuntimeError("BlockedCrsMatrix::doImport(): operation not supported.");
1037  }
1038 
1040  void doExport(const Matrix& dest, const Export& exporter, CombineMode CM) {
1041  XPETRA_MONITOR("XpetraBlockedCrsMatrix::doExport");
1042  if (Rows() == 1 && Cols () == 1) {
1043  getMatrix(0,0)->doExport(dest, exporter, CM);
1044  return;
1045  }
1046  throw Xpetra::Exceptions::RuntimeError("BlockedCrsMatrix::doExport(): operation not supported.");
1047  }
1048 
1049  // @}
1050 
1052 
1053 
1055  std::string description() const { return "Xpetra_BlockedCrsMatrix.description()"; }
1056 
1058  void describe(Teuchos::FancyOStream& out, const Teuchos::EVerbosityLevel verbLevel = Teuchos::Describable::verbLevel_default) const {
1059  out << "Xpetra::BlockedCrsMatrix: " << Rows() << " x " << Cols() << std::endl;
1060 
1061  if (isFillComplete()) {
1062  out << "BlockMatrix is fillComplete" << std::endl;
1063 
1064  /*if(fullrowmap_ != Teuchos::null) {
1065  out << "fullRowMap" << std::endl;
1066  fullrowmap_->describe(out,verbLevel);
1067  } else {
1068  out << "fullRowMap not set. Check whether block matrix is properly fillCompleted!" << std::endl;
1069  }*/
1070 
1071  //out << "fullColMap" << std::endl;
1072  //fullcolmap_->describe(out,verbLevel);
1073 
1074  } else {
1075  out << "BlockMatrix is NOT fillComplete" << std::endl;
1076  }
1077 
1078  for (size_t r = 0; r < Rows(); ++r)
1079  for (size_t c = 0; c < Cols(); ++c) {
1080  if(getMatrix(r,c)!=Teuchos::null) {
1081  out << "Block(" << r << "," << c << ")" << std::endl;
1082  getMatrix(r,c)->describe(out,verbLevel);
1083  } else out << "Block(" << r << "," << c << ") = null" << std::endl;
1084  }
1085  }
1086 
1088  RCP<const CrsGraph> getCrsGraph() const {
1089  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getCrsGraph");
1090  if (Rows() == 1 && Cols () == 1) {
1091  return getMatrix(0,0)->getCrsGraph();
1092  }
1093  throw Xpetra::Exceptions::RuntimeError("getCrsGraph() not supported by BlockedCrsMatrix");
1094  }
1095 
1097 
1099 
1100 
1102  virtual size_t Rows() const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::Rows"); return rangemaps_->NumMaps(); }
1103 
1105  virtual size_t Cols() const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::Cols"); return domainmaps_->NumMaps(); }
1106 
1108  Teuchos::RCP<Matrix> getCrsMatrix() const {
1109  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getCrsMatrix");
1110  TEUCHOS_TEST_FOR_EXCEPTION(Rows()!=1, std::out_of_range, "Can only unwrap a 1x1 blocked matrix. The matrix has " << Rows() << " block rows, though.");
1111  TEUCHOS_TEST_FOR_EXCEPTION(Cols()!=1, std::out_of_range, "Can only unwrap a 1x1 blocked matrix. The matrix has " << Cols() << " block columns, though.");
1112 
1113  RCP<Matrix> mat = getMatrix(0,0);
1114  RCP<BlockedCrsMatrix> bmat = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(mat);
1115  if (bmat == Teuchos::null) return mat;
1116  return bmat->getCrsMatrix();
1117  }
1118 
1119  Teuchos::RCP<Matrix> getInnermostCrsMatrix() {
1120  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getInnermostCrsMatrix");
1121  size_t row = Rows()+1, col = Cols()+1;
1122  for (size_t r = 0; r < Rows(); ++r)
1123  for(size_t c = 0; c < Cols(); ++c)
1124  if (getMatrix(r,c) != Teuchos::null) {
1125  row = r;
1126  col = c;
1127  break;
1128  }
1129  TEUCHOS_TEST_FOR_EXCEPTION(row == Rows()+1 || col == Cols()+1, Xpetra::Exceptions::Incompatible, "Xpetra::BlockedCrsMatrix::getInnermostCrsMatrix: Could not find a non-zero sub-block in blocked operator.")
1130  RCP<Matrix> mm = getMatrix(row,col);
1131  RCP<BlockedCrsMatrix> bmat = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(mm);
1132  if (bmat == Teuchos::null) return mm;
1133  return bmat->getInnermostCrsMatrix();
1134  }
1135 
1137  Teuchos::RCP<Matrix> getMatrix(size_t r, size_t c) const {
1138  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getMatrix");
1139  TEUCHOS_TEST_FOR_EXCEPTION(r > Rows(), std::out_of_range, "Error, r = " << Rows() << " is too big");
1140  TEUCHOS_TEST_FOR_EXCEPTION(c > Cols(), std::out_of_range, "Error, c = " << Cols() << " is too big");
1141 
1142  // transfer strided/blocked map information
1143  if (blocks_[r*Cols()+c] != Teuchos::null &&
1144  blocks_[r*Cols()+c]->IsView("stridedMaps") == false)
1145  blocks_[r*Cols()+c]->CreateView("stridedMaps", getRangeMap(r,bRangeThyraMode_), getDomainMap(c,bDomainThyraMode_));
1146  return blocks_[r*Cols()+c];
1147  }
1148 
1150  //void setMatrix(size_t r, size_t c, Teuchos::RCP<CrsMatrix> mat) {
1151  void setMatrix(size_t r, size_t c, Teuchos::RCP<Matrix> mat) {
1152  XPETRA_MONITOR("XpetraBlockedCrsMatrix::setMatrix");
1153  // TODO: if filled -> return error
1154 
1155  TEUCHOS_TEST_FOR_EXCEPTION(r > Rows(), std::out_of_range, "Error, r = " << Rows() << " is too big");
1156  TEUCHOS_TEST_FOR_EXCEPTION(c > Cols(), std::out_of_range, "Error, c = " << Cols() << " is too big");
1157 
1158  // set matrix
1159  blocks_[r*Cols() + c] = mat;
1160  }
1161 
1163  // NOTE: This is a rather expensive operation, since all blocks are copied
1164  // into a new big CrsMatrix
1165  Teuchos::RCP<Matrix> Merge() const {
1166  XPETRA_MONITOR("XpetraBlockedCrsMatrix::Merge");
1167  using Teuchos::RCP;
1168  using Teuchos::rcp_dynamic_cast;
1169  Scalar one = ScalarTraits<SC>::one();
1170 
1172  "BlockedCrsMatrix::Merge: only implemented for Xpetra-style or Thyra-style numbering. No mixup allowed!" );
1173 
1174  TEUCHOS_TEST_FOR_EXCEPTION(isFillComplete() == false, Xpetra::Exceptions::RuntimeError,
1175  "BlockedCrsMatrix::Merge: BlockMatrix must be fill-completed." );
1176 
1177  RCP<Matrix> sparse = MatrixFactory::Build(getRangeMapExtractor()->getFullMap(), 33);
1178 
1179  if(bRangeThyraMode_ == false) {
1180  // Xpetra mode
1181  for (size_t i = 0; i < Rows(); i++) {
1182  for (size_t j = 0; j < Cols(); j++) {
1183  if (getMatrix(i,j) != Teuchos::null) {
1184  RCP<const Matrix> mat = getMatrix(i,j);
1185 
1186  // recursively call Merge routine
1187  RCP<const BlockedCrsMatrix> bMat = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(mat);
1188  if (bMat != Teuchos::null) mat = bMat->Merge();
1189 
1190  bMat = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(mat);
1191  TEUCHOS_TEST_FOR_EXCEPTION(bMat != Teuchos::null, Xpetra::Exceptions::RuntimeError,
1192  "BlockedCrsMatrix::Merge: Merging of blocked sub-operators failed?!" );
1193 
1194  // jump over empty blocks
1195  if(mat->getNodeNumEntries() == 0) continue;
1196 
1197  this->Add(*mat, one, *sparse, one);
1198  }
1199  }
1200  }
1201  } else {
1202  // Thyra mode
1203  for (size_t i = 0; i < Rows(); i++) {
1204  for (size_t j = 0; j < Cols(); j++) {
1205  if (getMatrix(i,j) != Teuchos::null) {
1206  RCP<const Matrix> mat = getMatrix(i,j);
1207  // recursively call Merge routine
1208  RCP<const BlockedCrsMatrix> bMat = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(mat);
1209  if (bMat != Teuchos::null) mat = bMat->Merge();
1210 
1211  bMat = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(mat);
1212  TEUCHOS_TEST_FOR_EXCEPTION(bMat != Teuchos::null, Xpetra::Exceptions::RuntimeError,
1213  "BlockedCrsMatrix::Merge: Merging of blocked sub-operators failed?!" );
1214 
1215  // check whether we have a CrsMatrix block (no blocked operator)
1216  RCP<const CrsMatrixWrap> crsMat = Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(mat);
1217  TEUCHOS_ASSERT(crsMat != Teuchos::null);
1218 
1219  // these are the thyra style maps of the matrix
1220  RCP<const Map> trowMap = mat->getRowMap();
1221  RCP<const Map> tcolMap = mat->getColMap();
1222  RCP<const Map> tdomMap = mat->getDomainMap();
1223 
1224  // get Xpetra maps
1225  RCP<const Map> xrowMap = getRangeMapExtractor()->getMap(i,false);
1226  RCP<const Map> xdomMap = getDomainMapExtractor()->getMap(j,false);
1227 
1228  // generate column map with Xpetra GIDs
1229  // We have to do this separately for each block since the column
1230  // map of each block might be different in the same block column
1231  Teuchos::RCP<Map> xcolMap = MapUtils::transformThyra2XpetraGIDs(
1232  *tcolMap,
1233  *tdomMap,
1234  *xdomMap);
1235 
1236  // jump over empty blocks
1237  if(mat->getNodeNumEntries() == 0) continue;
1238 
1239  size_t maxNumEntries = mat->getNodeMaxNumRowEntries();
1240 
1241  size_t numEntries;
1242  Array<GO> inds (maxNumEntries);
1243  Array<GO> inds2(maxNumEntries);
1244  Array<SC> vals (maxNumEntries);
1245 
1246  // loop over all rows and add entries
1247  for (size_t k = 0; k < mat->getNodeNumRows(); k++) {
1248  GlobalOrdinal rowTGID = trowMap->getGlobalElement(k);
1249  crsMat->getCrsMatrix()->getGlobalRowCopy(rowTGID, inds(), vals(), numEntries);
1250 
1251  // create new indices array
1252  for (size_t l = 0; l < numEntries; ++l) {
1253  LocalOrdinal lid = tcolMap->getLocalElement(inds[l]);
1254  inds2[l] = xcolMap->getGlobalElement(lid);
1255  }
1256 
1257  GlobalOrdinal rowXGID = xrowMap->getGlobalElement(k);
1258  sparse->insertGlobalValues(
1259  rowXGID, inds2(0, numEntries),
1260  vals(0, numEntries));
1261  }
1262  }
1263  }
1264  }
1265  }
1266 
1267  sparse->fillComplete(getDomainMap(), getRangeMap());
1268 
1269  TEUCHOS_TEST_FOR_EXCEPTION(sparse->getNodeNumEntries() != getNodeNumEntries(), Xpetra::Exceptions::RuntimeError,
1270  "BlockedCrsMatrix::Merge: Local number of entries of merged matrix does not coincide with local number of entries of blocked operator." );
1271 
1272  TEUCHOS_TEST_FOR_EXCEPTION(sparse->getGlobalNumEntries() != getGlobalNumEntries(), Xpetra::Exceptions::RuntimeError,
1273  "BlockedCrsMatrix::Merge: Global number of entries of merged matrix does not coincide with global number of entries of blocked operator." );
1274 
1275  return sparse;
1276  }
1278 
1279 #ifdef HAVE_XPETRA_KOKKOS_REFACTOR
1280  typedef typename CrsMatrix::local_matrix_type local_matrix_type;
1281 
1283  local_matrix_type getLocalMatrix () const {
1284  if (Rows() == 1 && Cols () == 1) {
1285  return getMatrix(0,0)->getLocalMatrix();
1286  }
1287  throw Xpetra::Exceptions::RuntimeError("BlockedCrsMatrix::getLocalMatrix(): operation not supported.");
1288  }
1289 #endif
1290 
1291 #ifdef HAVE_XPETRA_THYRA
1292  Teuchos::RCP<Thyra::BlockedLinearOpBase<Scalar> > getThyraOperator() {
1293  Teuchos::RCP<Thyra::LinearOpBase<Scalar> > thOp =
1294  Xpetra::ThyraUtils<Scalar,LO,GO,Node>::toThyra(Teuchos::rcpFromRef(*this));
1295  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thOp));
1296 
1297  Teuchos::RCP<Thyra::BlockedLinearOpBase<Scalar> > thbOp =
1298  Teuchos::rcp_dynamic_cast<Thyra::BlockedLinearOpBase<Scalar> >(thOp);
1299  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thbOp));
1300  return thbOp;
1301  }
1302 #endif
1303 
1304  private:
1305 
1308 
1310 
1320  void Add(const Matrix& A, const Scalar scalarA, Matrix& B, const Scalar scalarB) const {
1321  XPETRA_MONITOR("XpetraBlockedCrsMatrix::Add");
1322  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Xpetra::Exceptions::RuntimeError,
1323  "Matrix A is not completed");
1324  using Teuchos::Array;
1325  using Teuchos::ArrayView;
1326 
1327  B.scale(scalarB);
1328 
1329  Scalar one = ScalarTraits<SC>::one();
1330  Scalar zero = ScalarTraits<SC>::zero();
1331 
1332  if (scalarA == zero)
1333  return;
1334 
1335  Teuchos::RCP<const Matrix> rcpA = Teuchos::rcpFromRef(A);
1336  Teuchos::RCP<const CrsMatrixWrap> rcpAwrap = Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(rcpA);
1337  TEUCHOS_TEST_FOR_EXCEPTION(rcpAwrap == Teuchos::null, Xpetra::Exceptions::BadCast,
1338  "BlockedCrsMatrix::Add: matrix A must be of type CrsMatrixWrap.");
1339  Teuchos::RCP<const CrsMatrix> crsA = rcpAwrap->getCrsMatrix();
1340 
1341  size_t maxNumEntries = crsA->getNodeMaxNumRowEntries();
1342 
1343  size_t numEntries;
1344  Array<GO> inds(maxNumEntries);
1345  Array<SC> vals(maxNumEntries);
1346 
1347  RCP<const Map> rowMap = crsA->getRowMap();
1348  RCP<const Map> colMap = crsA->getColMap();
1349 
1350  ArrayView<const GO> rowGIDs = crsA->getRowMap()->getNodeElementList();
1351  for (size_t i = 0; i < crsA->getNodeNumRows(); i++) {
1352  GO row = rowGIDs[i];
1353  crsA->getGlobalRowCopy(row, inds(), vals(), numEntries);
1354 
1355  if (scalarA != one)
1356  for (size_t j = 0; j < numEntries; ++j)
1357  vals[j] *= scalarA;
1358 
1359  B.insertGlobalValues(row, inds(0, numEntries), vals(0, numEntries)); // insert should be ok, since blocks in BlockedCrsOpeartor do not overlap!
1360  }
1361  }
1362 
1364 
1365  // Default view is created after fillComplete()
1366  // Because ColMap might not be available before fillComplete().
1368  XPETRA_MONITOR("XpetraBlockedCrsMatrix::CreateDefaultView");
1369 
1370  // Create default view
1371  this->defaultViewLabel_ = "point";
1373 
1374  // Set current view
1375  this->currentViewLabel_ = this->GetDefaultViewLabel();
1376  }
1377 
1378  private:
1379  Teuchos::RCP<const MapExtractor> domainmaps_; // full domain map together with all partial domain maps
1380  Teuchos::RCP<const MapExtractor> rangemaps_; // full range map together with all partial domain maps
1381 
1382  std::vector<Teuchos::RCP<Matrix> > blocks_; // row major matrix block storage
1383 #ifdef HAVE_XPETRA_THYRA
1384  Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar> > thyraOp_;
1385 #endif
1388 
1389 };
1390 
1391 } //namespace Xpetra
1392 
1393 #define XPETRA_BLOCKEDCRSMATRIX_SHORT
1394 #endif /* XPETRA_BLOCKEDCRSMATRIX_HPP */
Teuchos::RCP< Matrix > Merge() const
merge BlockedCrsMatrix blocks in a CrsMatrix
size_t getNodeMaxNumRowEntries() const
Returns the maximum number of entries across all rows/columns on this node.
RCP< const CrsGraph > getCrsGraph() const
Returns the CrsGraph associated with this matrix.
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.
const Teuchos::RCP< const Map > getMap() const
Implements DistObject interface.
bool isLocallyIndexed() const
If matrix indices of all matrix blocks are in the local range, this function returns true...
bool isFillComplete() const
Returns true if fillComplete() has been called and the matrix is in compute mode. ...
RCP< const Map > getRangeMap(size_t i, bool bThyraMode) const
Returns the Map associated with the i&#39;th block range of this operator.
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.
virtual void setAllToScalar(const Scalar &alpha)
Set all matrix entries equal to scalar.
GlobalOrdinal GO
RCP< const Map > getDomainMap(size_t i) const
Returns the Map associated with the i&#39;th block domain of this operator.
global_size_t getGlobalNumDiags() const
Returns the number of global diagonal entries, based on global row/column index comparisons.
Teuchos::RCP< const MapExtractor > rangemaps_
Xpetra namespace
void doImport(const Matrix &source, const Import &importer, CombineMode CM)
Import.
Exception throws to report errors in the internal logical of the program.
virtual ScalarTraits< Scalar >::magnitudeType getFrobeniusNorm() const
Get Frobenius norm of the matrix.
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 ~BlockedCrsMatrix()
Destructor.
Teuchos::RCP< Matrix > getCrsMatrix() const
return unwrap 1x1 blocked operators
void doExport(const Matrix &dest, const Import &importer, CombineMode CM)
Export.
size_t getGlobalMaxNumRowEntries() const
Returns the maximum number of entries across all rows/columns on all nodes.
global_size_t getGlobalNumCols() const
Returns the number of global columns in the matrix.
void resumeFill(const RCP< ParameterList > &params=null)
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.
size_t getNodeNumDiags() const
Returns the number of local diagonal entries, based on global row/column index comparisons.
void rightScale(const Vector &x)
Right scale matrix using the given vector entries.
RCP< const MapExtractor > getRangeMapExtractor() const
Returns map extractor class for range map.
Exception indicating invalid cast attempted.
void insertGlobalValues(GlobalOrdinal globalRow, const ArrayView< const GlobalOrdinal > &cols, const ArrayView< const Scalar > &vals)
Insert matrix entries, using global IDs.
viewLabel_t defaultViewLabel_
virtual void insertGlobalValues(GlobalOrdinal globalRow, const ArrayView< const GlobalOrdinal > &cols, const ArrayView< const Scalar > &vals)=0
Insert matrix entries, using global IDs.
size_t getNodeNumRows() const
Returns the number of matrix rows owned on the calling node.
void doExport(const Matrix &dest, const Export &exporter, CombineMode CM)
Export (using an Importer).
void leftScale(const Vector &x)
Left scale matrix using the given vector entries.
Teuchos::RCP< Matrix > getInnermostCrsMatrix()
void insertLocalValues(LocalOrdinal localRow, const ArrayView< const LocalOrdinal > &cols, const ArrayView< const Scalar > &vals)
Insert matrix entries, using local IDs.
bool isGloballyIndexed() const
If matrix indices are in the global range, this function returns true. Otherwise, this function retur...
RCP< const Map > getDomainMap(size_t i, bool bThyraMode) const
Returns the Map associated with the i&#39;th block domain of this operator.
virtual bool isFillComplete() const =0
Returns true if fillComplete() has been called and the matrix is in compute mode. ...
const viewLabel_t & GetDefaultViewLabel() const
void getLocalDiagCopy(Vector &diag) const
Get a copy of the diagonal entries owned by this node, with local row indices.
void Add(const Matrix &A, const Scalar scalarA, Matrix &B, const Scalar scalarB) const
Add a Xpetra::CrsMatrix to another: B = B*scalarB + A*scalarA.
static Teuchos::RCP< Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > transformThyra2XpetraGIDs(const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &input, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &nonOvlInput, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &nonOvlReferenceInput)
replace set of global ids by new global ids
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_
void doImport(const Matrix &source, const Export &exporter, CombineMode CM)
Import (using an Exporter).
RCP< const Map > getRangeMap(size_t i) const
Returns the Map associated with the i&#39;th block range of this operator.
std::string viewLabel_t
size_t global_size_t
Global size_t object.
global_size_t getGlobalNumEntries() const
Returns the global number of entries in this matrix.
std::vector< Teuchos::RCP< Matrix > > blocks_
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...
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.
RCP< const CrsGraph< int, GlobalOrdinal, Node > > toXpetra(const Epetra_CrsGraph &g)
#define XPETRA_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void setMatrix(size_t r, size_t c, Teuchos::RCP< Matrix > mat)
set matrix block
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.
Exception throws to report incompatible objects (like maps).
size_t getNodeNumEntries() const
Returns the local number of entries in this matrix.
RCP< const Map > getRangeMap() const
Returns the Map associated with the full range of this operator.
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.
Concrete implementation of Xpetra::Matrix.
Teuchos::RCP< Matrix > getMatrix(size_t r, size_t c) const
return block (r,c)
std::string description() const
Return a simple one-line description of this object.
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.
bool IsView(const viewLabel_t viewLabel) const
CombineMode
Xpetra::Combine Mode enumerable type.
global_size_t getGlobalNumRows() const
Returns the number of global rows.
#define XPETRA_MONITOR(funcName)
RCP< const MapExtractor > getDomainMapExtractor() const
Returns map extractor for domain map.
static RCP< Matrix > Build(const RCP< const Map > &rowMap, size_t maxNumEntriesPerRow, Xpetra::ProfileType pftype=Xpetra::DynamicProfile)
Constructor specifying the number of non-zeros for all rows.
virtual void scale(const Scalar &alpha)=0
Scale the current values of a matrix, this = alpha*this.
void fillComplete(const RCP< ParameterList > &params=null)
Signal that data entry is complete.
virtual size_t Rows() const
number of row blocks
virtual size_t Cols() const
number of column blocks
virtual size_t getNumVectors() const =0
Number of columns in the multivector.
RCP< const Map > getDomainMap() const
Returns the Map associated with the full domain of this operator.
void removeEmptyProcessesInPlace(const Teuchos::RCP< const Map > &newMap)
size_t getNumEntriesInLocalRow(LocalOrdinal localRow) const
Returns the current number of entries on this node in the specified local row.
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_
bool bRangeThyraMode_
boolean flag, which is true, if BlockedCrsMatrix has been created using Thyra-style numbering for sub...
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.