Xpetra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
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"
59 #include "Xpetra_BlockedMultiVector.hpp"
60 #include "Xpetra_MultiVectorFactory.hpp"
61 #include "Xpetra_BlockedVector.hpp"
62 #include "Xpetra_CrsGraph.hpp"
63 #include "Xpetra_CrsMatrix.hpp"
65 
66 #include "Xpetra_MapExtractor.hpp"
67 
68 #include "Xpetra_Matrix.hpp"
69 #include "Xpetra_MatrixFactory.hpp"
70 #include "Xpetra_CrsMatrixWrap.hpp"
71 
72 #ifdef HAVE_XPETRA_THYRA
73 #include <Thyra_ProductVectorSpaceBase.hpp>
74 #include <Thyra_VectorSpaceBase.hpp>
75 #include <Thyra_LinearOpBase.hpp>
76 #include <Thyra_BlockedLinearOpBase.hpp>
77 #include <Thyra_PhysicallyBlockedLinearOpBase.hpp>
78 #include "Xpetra_ThyraUtils.hpp"
79 #endif
80 
81 #include "Xpetra_VectorFactory.hpp"
82 
83 
88 namespace Xpetra {
89 
90  typedef std::string viewLabel_t;
91 
92  template <class Scalar,
93  class LocalOrdinal,
94  class GlobalOrdinal,
97  public Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> {
98  public:
99  typedef Scalar scalar_type;
100  typedef LocalOrdinal local_ordinal_type;
101  typedef GlobalOrdinal global_ordinal_type;
102  typedef Node node_type;
103 
104  private:
105 #undef XPETRA_BLOCKEDCRSMATRIX_SHORT
106 #include "Xpetra_UseShortNames.hpp"
107 
108  public:
109 
111 
112 
114 
119  BlockedCrsMatrix(const Teuchos::RCP<const BlockedMap>& rangeMaps,
120  const Teuchos::RCP<const BlockedMap>& domainMaps,
121  size_t npr) {
122  is_diagonal_ = true;
123  domainmaps_ = Teuchos::rcp(new MapExtractor(domainMaps));
124  rangemaps_ = Teuchos::rcp(new MapExtractor(rangeMaps));
125  bRangeThyraMode_ = rangeMaps->getThyraMode();
126  bDomainThyraMode_ = domainMaps->getThyraMode();
127 
128  blocks_.reserve(Rows()*Cols());
129 
130  // add CrsMatrix objects in row,column order
131  for (size_t r = 0; r < Rows(); ++r)
132  for (size_t c = 0; c < Cols(); ++c)
134 
135  // Default view
137  }
138 
140 
147  BlockedCrsMatrix(Teuchos::RCP<const MapExtractor>& rangeMaps,
148  Teuchos::RCP<const MapExtractor>& domainMaps,
149  size_t npr)
150  : is_diagonal_(true), domainmaps_(domainMaps), rangemaps_(rangeMaps)
151  {
152  bRangeThyraMode_ = rangeMaps->getThyraMode();
153  bDomainThyraMode_ = domainMaps->getThyraMode();
154 
155  blocks_.reserve(Rows()*Cols());
156 
157  // add CrsMatrix objects in row,column order
158  for (size_t r = 0; r < Rows(); ++r)
159  for (size_t c = 0; c < Cols(); ++c)
161 
162  // Default view
164  }
165 
166 #ifdef HAVE_XPETRA_THYRA
167 
173  BlockedCrsMatrix(const Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar> >& thyraOp, const Teuchos::RCP<const Teuchos::Comm<int> >& /* comm */)
174  : is_diagonal_(true), thyraOp_(thyraOp)
175  {
176  // extract information from Thyra blocked operator and rebuilt information
177  const Teuchos::RCP<const Thyra::ProductVectorSpaceBase<Scalar> > productRangeSpace = thyraOp->productRange();
178  const Teuchos::RCP<const Thyra::ProductVectorSpaceBase<Scalar> > productDomainSpace = thyraOp->productDomain();
179 
180  int numRangeBlocks = productRangeSpace->numBlocks();
181  int numDomainBlocks = productDomainSpace->numBlocks();
182 
183  // build range map extractor from Thyra::BlockedLinearOpBase object
184  std::vector<Teuchos::RCP<const Map> > subRangeMaps(numRangeBlocks);
185  for (size_t r=0; r<Teuchos::as<size_t>(numRangeBlocks); ++r) {
186  for (size_t c=0; c<Teuchos::as<size_t>(numDomainBlocks); ++c) {
187  if (thyraOp->blockExists(r,c)) {
188  // we only need at least one block in each block row to extract the range map
189  Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > const_op = thyraOp->getBlock(r,c); // nonConst access is not allowed.
190  Teuchos::RCP<const Xpetra::CrsMatrix<Scalar,LO,GO,Node> > xop =
192  subRangeMaps[r] = xop->getRangeMap();
193  if(r!=c) is_diagonal_ = false;
194  break;
195  }
196  }
197  }
198  Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > fullRangeMap = mergeMaps(subRangeMaps);
199 
200  // check whether the underlying Thyra operator uses Thyra-style numbering (default for most applications) or
201  // Xpetra style numbering
202  bool bRangeUseThyraStyleNumbering = false;
203  size_t numAllElements = 0;
204  for(size_t v = 0; v < subRangeMaps.size(); ++v) {
205  numAllElements += subRangeMaps[v]->getGlobalNumElements();
206  }
207  if ( fullRangeMap->getGlobalNumElements() != numAllElements) bRangeUseThyraStyleNumbering = true;
208  rangemaps_ = Teuchos::rcp(new Xpetra::MapExtractor<Scalar,LocalOrdinal,GlobalOrdinal,Node>(fullRangeMap, subRangeMaps, bRangeUseThyraStyleNumbering));
209 
210  // build domain map extractor from Thyra::BlockedLinearOpBase object
211  std::vector<Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > > subDomainMaps(numDomainBlocks);
212  for (size_t c=0; c<Teuchos::as<size_t>(numDomainBlocks); ++c) {
213  for (size_t r=0; r<Teuchos::as<size_t>(numRangeBlocks); ++r) {
214  if (thyraOp->blockExists(r,c)) {
215  // we only need at least one block in each block row to extract the range map
216  Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > const_op = thyraOp->getBlock(r,c); // nonConst access is not allowed.
217  Teuchos::RCP<const Xpetra::CrsMatrix<Scalar,LO,GO,Node> > xop =
219  subDomainMaps[c] = xop->getDomainMap();
220  break;
221  }
222  }
223  }
224  Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > fullDomainMap = mergeMaps(subDomainMaps);
225  // plausibility check for numbering style (Xpetra or Thyra)
226  bool bDomainUseThyraStyleNumbering = false;
227  numAllElements = 0;
228  for(size_t v = 0; v < subDomainMaps.size(); ++v) {
229  numAllElements += subDomainMaps[v]->getGlobalNumElements();
230  }
231  if (fullDomainMap->getGlobalNumElements() != numAllElements) bDomainUseThyraStyleNumbering = true;
232  domainmaps_ = Teuchos::rcp(new Xpetra::MapExtractor<Scalar,LocalOrdinal,GlobalOrdinal,Node>(fullDomainMap, subDomainMaps, bDomainUseThyraStyleNumbering));
233 
234  // store numbering mode
235  bRangeThyraMode_ = bRangeUseThyraStyleNumbering;
236  bDomainThyraMode_ = bDomainUseThyraStyleNumbering;
237 
238  // add CrsMatrix objects in row,column order
239  blocks_.reserve(Rows()*Cols());
240  for (size_t r = 0; r < Rows(); ++r) {
241  for (size_t c = 0; c < Cols(); ++c) {
242  if(thyraOp->blockExists(r,c)) {
243  // TODO we do not support nested Thyra operators here!
244  Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > const_op = thyraOp->getBlock(r,c); // nonConst access is not allowed.
245  Teuchos::RCP<Thyra::LinearOpBase<Scalar> > op = Teuchos::rcp_const_cast<Thyra::LinearOpBase<Scalar> >(const_op); // cast away const
246  Teuchos::RCP<Xpetra::CrsMatrix<Scalar,LO,GO,Node> > xop =
248  Teuchos::RCP<Xpetra::CrsMatrixWrap<Scalar,LO,GO,Node> > xwrap =
249  Teuchos::rcp(new Xpetra::CrsMatrixWrap<Scalar,LO,GO,Node>(xop));
250  blocks_.push_back(xwrap);
251  } else {
252  // add empty block
254  }
255  }
256  }
257  // Default view
259  }
260 
261  private:
263 
270  Teuchos::RCP<const Map> mergeMaps(std::vector<Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > > & subMaps) {
271  // TODO merging for Thyra mode is missing (similar to what we do in constructor of MapExtractor
272 
273  // merge submaps to global map
274  std::vector<GlobalOrdinal> gids;
275  for(size_t tt = 0; tt<subMaps.size(); ++tt) {
276  Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > subMap = subMaps[tt];
277 #if 1
278  Teuchos::ArrayView< const GlobalOrdinal > subMapGids = subMap->getNodeElementList();
279  gids.insert(gids.end(), subMapGids.begin(), subMapGids.end());
280 #else
281  size_t myNumElements = subMap->getNodeNumElements();
282  for(LocalOrdinal l = 0; l < Teuchos::as<LocalOrdinal>(myNumElements); ++l) {
283  GlobalOrdinal gid = subMap->getGlobalElement(l);
284  gids.push_back(gid);
285  }
286 #endif
287  }
288 
289  // we have to sort the matrix entries and get rid of the double entries
290  // since we use this to detect Thyra-style numbering or Xpetra-style
291  // numbering. In Thyra-style numbering mode, the Xpetra::MapExtractor builds
292  // the correct row maps.
293  const GO INVALID = Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid();
294  std::sort(gids.begin(), gids.end());
295  gids.erase(std::unique(gids.begin(), gids.end()), gids.end());
296  Teuchos::ArrayView<GO> gidsView(&gids[0], gids.size());
297  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());
298  return fullMap;
299  }
300 
301  public:
302 #endif
303 
305  virtual ~BlockedCrsMatrix() {}
306 
308 
309 
311 
312 
314 
339  void insertGlobalValues(GlobalOrdinal globalRow, const ArrayView<const GlobalOrdinal>& cols, const ArrayView<const Scalar>& vals) {
340  XPETRA_MONITOR("XpetraBlockedCrsMatrix::insertGlobalValues");
341  if (Rows() == 1 && Cols () == 1) {
342  getMatrix(0,0)->insertGlobalValues(globalRow, cols, vals);
343  return;
344  }
345  throw Xpetra::Exceptions::RuntimeError("insertGlobalValues not supported by BlockedCrsMatrix");
346  }
347 
349 
359  void insertLocalValues(LocalOrdinal localRow, const ArrayView<const LocalOrdinal>& cols, const ArrayView<const Scalar>& vals) {
360  XPETRA_MONITOR("XpetraBlockedCrsMatrix::insertLocalValues");
361  if (Rows() == 1 && Cols () == 1) {
362  getMatrix(0,0)->insertLocalValues(localRow, cols, vals);
363  return;
364  }
365  throw Xpetra::Exceptions::RuntimeError("insertLocalValues not supported by BlockedCrsMatrix");
366  }
367 
368  void removeEmptyProcessesInPlace(const Teuchos::RCP<const Map>& newMap) {
369  XPETRA_MONITOR("XpetraBlockedCrsMatrix::removeEmptyProcessesInPlace");
370  if (Rows() == 1 && Cols () == 1) {
371  getMatrix(0,0)->removeEmptyProcessesInPlace(newMap);
372  return;
373  }
374  throw Xpetra::Exceptions::RuntimeError("removeEmptyProcesses not supported by BlockedCrsMatrix");
375  }
376 
378 
386  void replaceGlobalValues(GlobalOrdinal globalRow,
387  const ArrayView<const GlobalOrdinal> &cols,
388  const ArrayView<const Scalar> &vals) {
389  XPETRA_MONITOR("XpetraBlockedCrsMatrix::replaceGlobalValues");
390  if (Rows() == 1 && Cols () == 1) {
391  getMatrix(0,0)->replaceGlobalValues(globalRow,cols,vals);
392  return;
393  }
394  throw Xpetra::Exceptions::RuntimeError("replaceGlobalValues not supported by BlockedCrsMatrix");
395  }
396 
398 
402  void replaceLocalValues(LocalOrdinal localRow,
403  const ArrayView<const LocalOrdinal> &cols,
404  const ArrayView<const Scalar> &vals) {
405  XPETRA_MONITOR("XpetraBlockedCrsMatrix::replaceLocalValues");
406  if (Rows() == 1 && Cols () == 1) {
407  getMatrix(0,0)->replaceLocalValues(localRow,cols,vals);
408  return;
409  }
410  throw Xpetra::Exceptions::RuntimeError("replaceLocalValues not supported by BlockedCrsMatrix");
411  }
412 
414  virtual void setAllToScalar(const Scalar& alpha) {
415  XPETRA_MONITOR("XpetraBlockedCrsMatrix::setAllToScalar");
416  for (size_t row = 0; row < Rows(); row++) {
417  for (size_t col = 0; col < Cols(); col++) {
418  if (!getMatrix(row,col).is_null()) {
419  getMatrix(row,col)->setAllToScalar(alpha);
420  }
421  }
422  }
423  }
424 
426  void scale(const Scalar& alpha) {
427  XPETRA_MONITOR("XpetraBlockedCrsMatrix::scale");
428  for (size_t row = 0; row < Rows(); row++) {
429  for (size_t col = 0; col < Cols(); col++) {
430  if (!getMatrix(row,col).is_null()) {
431  getMatrix(row,col)->scale(alpha);
432  }
433  }
434  }
435  }
436 
438 
440 
441 
449  void resumeFill(const RCP< ParameterList >& params = null) {
450  XPETRA_MONITOR("XpetraBlockedCrsMatrix::resumeFill");
451  for (size_t row = 0; row < Rows(); row++) {
452  for (size_t col = 0; col < Cols(); col++) {
453  if (!getMatrix(row,col).is_null()) {
454  getMatrix(row,col)->resumeFill(params);
455  }
456  }
457  }
458  }
459 
465  void fillComplete(const RCP<const Map>& domainMap, const RCP<const Map>& rangeMap, const RCP<ParameterList>& params = null) {
466  XPETRA_MONITOR("XpetraBlockedCrsMatrix::fillComplete");
467  if (Rows() == 1 && Cols () == 1) {
468  getMatrix(0,0)->fillComplete(domainMap, rangeMap, params);
469  return;
470  }
471  fillComplete(params);
472  }
473 
488  void fillComplete(const RCP<ParameterList>& params = null) {
489  XPETRA_MONITOR("XpetraBlockedCrsMatrix::fillComplete");
490  TEUCHOS_TEST_FOR_EXCEPTION(rangemaps_==Teuchos::null, Xpetra::Exceptions::RuntimeError,"BlockedCrsMatrix::fillComplete: rangemaps_ is not set. Error.");
491 
492  for (size_t r = 0; r < Rows(); ++r)
493  for (size_t c = 0; c < Cols(); ++c) {
494  if(getMatrix(r,c) != Teuchos::null) {
495  Teuchos::RCP<Matrix> Ablock = getMatrix(r,c);
496  if(r!=c) is_diagonal_ = false;
497  if (Ablock != Teuchos::null && !Ablock->isFillComplete())
498  Ablock->fillComplete(getDomainMap(c, bDomainThyraMode_), getRangeMap(r, bRangeThyraMode_), params);
499  }
500  }
501 
502 #if 0
503  // get full row map
504  RCP<const Map> rangeMap = rangemaps_->getFullMap();
505  fullrowmap_ = MapFactory::Build(rangeMap()->lib(), rangeMap()->getGlobalNumElements(), rangeMap()->getNodeElementList(), rangeMap()->getIndexBase(), rangeMap()->getComm());
506 
507  // build full col map
508  fullcolmap_ = Teuchos::null; // delete old full column map
509 
510  std::vector<GO> colmapentries;
511  for (size_t c = 0; c < Cols(); ++c) {
512  // copy all local column lids of all block rows to colset
513  std::set<GO> colset;
514  for (size_t r = 0; r < Rows(); ++r) {
515  Teuchos::RCP<CrsMatrix> Ablock = getMatrix(r,c);
516 
517  if (Ablock != Teuchos::null) {
518  Teuchos::ArrayView<const GO> colElements = Ablock->getColMap()->getNodeElementList();
519  Teuchos::RCP<const Map> colmap = Ablock->getColMap();
520  copy(colElements.getRawPtr(), colElements.getRawPtr() + colElements.size(), inserter(colset, colset.begin()));
521  }
522  }
523 
524  // remove duplicates (entries which are in column maps of more than one block row)
525  colmapentries.reserve(colmapentries.size() + colset.size());
526  copy(colset.begin(), colset.end(), back_inserter(colmapentries));
527  sort(colmapentries.begin(), colmapentries.end());
528  typename std::vector<GO>::iterator gendLocation;
529  gendLocation = std::unique(colmapentries.begin(), colmapentries.end());
530  colmapentries.erase(gendLocation,colmapentries.end());
531  }
532 
533  // sum up number of local elements
534  size_t numGlobalElements = 0;
535  Teuchos::reduceAll(*(rangeMap->getComm()), Teuchos::REDUCE_SUM, colmapentries.size(), Teuchos::outArg(numGlobalElements));
536 
537  // store global full column map
538  const Teuchos::ArrayView<const GO> aView = Teuchos::ArrayView<const GO>(colmapentries);
539  fullcolmap_ = MapFactory::Build(rangeMap->lib(), numGlobalElements, aView, 0, rangeMap->getComm());
540 #endif
541  }
542 
544 
546 
549  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getGlobalNumRows");
550  global_size_t globalNumRows = 0;
551 
552  for (size_t row = 0; row < Rows(); row++)
553  for (size_t col = 0; col < Cols(); col++)
554  if (!getMatrix(row,col).is_null()) {
555  globalNumRows += getMatrix(row,col)->getGlobalNumRows();
556  break; // we need only one non-null matrix in a row
557  }
558 
559  return globalNumRows;
560  }
561 
563 
566  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getGlobalNumCols");
567  global_size_t globalNumCols = 0;
568 
569  for (size_t col = 0; col < Cols(); col++)
570  for (size_t row = 0; row < Rows(); row++)
571  if (!getMatrix(row,col).is_null()) {
572  globalNumCols += getMatrix(row,col)->getGlobalNumCols();
573  break; // we need only one non-null matrix in a col
574  }
575 
576  return globalNumCols;
577  }
578 
580  size_t getNodeNumRows() const {
581  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getNodeNumRows");
582  global_size_t nodeNumRows = 0;
583 
584  for (size_t row = 0; row < Rows(); ++row)
585  for (size_t col = 0; col < Cols(); col++)
586  if (!getMatrix(row,col).is_null()) {
587  nodeNumRows += getMatrix(row,col)->getNodeNumRows();
588  break; // we need only one non-null matrix in a row
589  }
590 
591  return nodeNumRows;
592  }
593 
596  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getGlobalNumEntries");
597  global_size_t globalNumEntries = 0;
598 
599  for (size_t row = 0; row < Rows(); ++row)
600  for (size_t col = 0; col < Cols(); ++col)
601  if (!getMatrix(row,col).is_null())
602  globalNumEntries += getMatrix(row,col)->getGlobalNumEntries();
603 
604  return globalNumEntries;
605  }
606 
608  size_t getNodeNumEntries() const {
609  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getNodeNumEntries");
610  global_size_t nodeNumEntries = 0;
611 
612  for (size_t row = 0; row < Rows(); ++row)
613  for (size_t col = 0; col < Cols(); ++col)
614  if (!getMatrix(row,col).is_null())
615  nodeNumEntries += getMatrix(row,col)->getNodeNumEntries();
616 
617  return nodeNumEntries;
618  }
619 
621 
622  size_t getNumEntriesInLocalRow(LocalOrdinal localRow) const {
623  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getNumEntriesInLocalRow");
624  GlobalOrdinal gid = this->getRowMap()->getGlobalElement(localRow);
625  size_t row = getBlockedRangeMap()->getMapIndexForGID(gid);
626  LocalOrdinal lid = getBlockedRangeMap()->getMap(row)->getLocalElement(gid);
627  size_t numEntriesInLocalRow = 0;
628  for (size_t col = 0; col < Cols(); ++col)
629  if (!getMatrix(row,col).is_null())
630  numEntriesInLocalRow += getMatrix(row,col)->getNumEntriesInLocalRow(lid);
631  return numEntriesInLocalRow;
632  }
633 
635 
636  size_t getNumEntriesInGlobalRow(GlobalOrdinal globalRow) const {
637  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getNumEntriesInGlobalRow");
638  size_t row = getBlockedRangeMap()->getMapIndexForGID(globalRow);
639  size_t numEntriesInGlobalRow = 0;
640  for (size_t col = 0; col < Cols(); ++col)
641  if (!getMatrix(row,col).is_null())
642  numEntriesInGlobalRow += getMatrix(row,col)->getNumEntriesInGlobalRow(globalRow);
643  return numEntriesInGlobalRow;
644  }
645 
647 
649  size_t getGlobalMaxNumRowEntries() const {
650  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getGlobalMaxNumRowEntries");
651  global_size_t globalMaxEntries = 0;
652 
653  for (size_t row = 0; row < Rows(); row++) {
654  global_size_t globalMaxEntriesBlockRows = 0;
655  for (size_t col = 0; col < Cols(); col++) {
656  if (!getMatrix(row,col).is_null()) {
657  globalMaxEntriesBlockRows += getMatrix(row,col)->getGlobalMaxNumRowEntries();
658  }
659  }
660  if(globalMaxEntriesBlockRows > globalMaxEntries)
661  globalMaxEntries = globalMaxEntriesBlockRows;
662  }
663  return globalMaxEntries;
664  }
665 
667 
669  size_t getNodeMaxNumRowEntries() const {
670  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getNodeMaxNumRowEntries");
671  size_t localMaxEntries = 0;
672 
673  for (size_t row = 0; row < Rows(); row++) {
674  size_t localMaxEntriesBlockRows = 0;
675  for (size_t col = 0; col < Cols(); col++) {
676  if (!getMatrix(row,col).is_null()) {
677  localMaxEntriesBlockRows += getMatrix(row,col)->getNodeMaxNumRowEntries();
678  }
679  }
680  if(localMaxEntriesBlockRows > localMaxEntries)
681  localMaxEntries = localMaxEntriesBlockRows;
682  }
683  return localMaxEntries;
684  }
685 
687 
690  bool isLocallyIndexed() const {
691  XPETRA_MONITOR("XpetraBlockedCrsMatrix::isLocallyIndexed");
692  for (size_t i = 0; i < blocks_.size(); ++i)
693  if (blocks_[i] != Teuchos::null && !blocks_[i]->isLocallyIndexed())
694  return false;
695  return true;
696  }
697 
699 
702  bool isGloballyIndexed() const {
703  XPETRA_MONITOR("XpetraBlockedCrsMatrix::isGloballyIndexed");
704  for (size_t i = 0; i < blocks_.size(); i++)
705  if (blocks_[i] != Teuchos::null && !blocks_[i]->isGloballyIndexed())
706  return false;
707  return true;
708  }
709 
711  bool isFillComplete() const {
712  XPETRA_MONITOR("XpetraBlockedCrsMatrix::isFillComplete");
713  for (size_t i = 0; i < blocks_.size(); i++)
714  if (blocks_[i] != Teuchos::null && !blocks_[i]->isFillComplete())
715  return false;
716  return true;
717  }
718 
720 
734  virtual void getLocalRowCopy(LocalOrdinal LocalRow,
735  const ArrayView<LocalOrdinal>& Indices,
736  const ArrayView<Scalar>& Values,
737  size_t &NumEntries) const {
738  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getLocalRowCopy");
739  if (Rows() == 1 && Cols () == 1) {
740  getMatrix(0,0)->getLocalRowCopy(LocalRow, Indices, Values, NumEntries);
741  return;
742  }
743  throw Xpetra::Exceptions::RuntimeError("getLocalRowCopy not supported by BlockedCrsMatrix" );
744  }
745 
747 
756  void getGlobalRowView(GlobalOrdinal GlobalRow, ArrayView<const GlobalOrdinal>& indices, ArrayView<const Scalar>& values) const {
757  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getGlobalRowView");
758  if (Rows() == 1 && Cols () == 1) {
759  getMatrix(0,0)->getGlobalRowView(GlobalRow, indices, values);
760  return;
761  }
762  throw Xpetra::Exceptions::RuntimeError("getGlobalRowView not supported by BlockedCrsMatrix");
763  }
764 
766 
775  void getLocalRowView(LocalOrdinal LocalRow, ArrayView<const LocalOrdinal>& indices, ArrayView<const Scalar>& values) const {
776  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getLocalRowView");
777  if (Rows() == 1 && Cols () == 1) {
778  getMatrix(0,0)->getLocalRowView(LocalRow, indices, values);
779  return;
780  }
781  else if(is_diagonal_) {
782  GlobalOrdinal gid = this->getRowMap()->getGlobalElement(LocalRow);
783  size_t row = getBlockedRangeMap()->getMapIndexForGID(gid);
784  getMatrix(row,row)->getLocalRowView(getMatrix(row,row)->getRowMap()->getLocalElement(gid),indices,values);
785  return;
786  }
787  throw Xpetra::Exceptions::RuntimeError("getLocalRowView not supported by BlockedCrsMatrix");
788  }
789 
791 
794  void getLocalDiagCopy(Vector& diag) const {
795  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getLocalDiagCopy");
796 
797  //RCP<Teuchos::FancyOStream> fancy = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
798 
799  Teuchos::RCP<Vector> rcpdiag = Teuchos::rcpFromRef(diag);
800  Teuchos::RCP<BlockedVector> bdiag = Teuchos::rcp_dynamic_cast<BlockedVector>(rcpdiag);
801 
802  // special treatment for 1x1 block matrices
803  // ReorderedBlockedCrsMatrix object encapsulate single blocks in ReorderedBlocks
804  // BlockedVectors have Vector objects as Leaf objects.
805  if(Rows() == 1 && Cols() == 1 && bdiag.is_null() == true) {
806  Teuchos::RCP<const Matrix> rm = getMatrix(0,0);
807  rm->getLocalDiagCopy(diag);
808  return;
809  }
810 
811  TEUCHOS_TEST_FOR_EXCEPTION(bdiag.is_null() == true, Xpetra::Exceptions::RuntimeError, "BlockedCrsMatrix::getLocalDiagCopy: diag must be a Blocked(Multi)Vector.");
812  TEUCHOS_TEST_FOR_EXCEPTION(bdiag->getNumVectors() != 1, Xpetra::Exceptions::RuntimeError, "BlockedCrsMatrix::getLocalDiagCopy: diag must be a Blocked(Multi)Vector with exactly one vector. However, the number of stored vectors is " << bdiag->getNumVectors());
813  TEUCHOS_TEST_FOR_EXCEPTION(bdiag->getBlockedMap()->getNumMaps() != this->Rows(), Xpetra::Exceptions::RuntimeError,
814  "BlockedCrsMatrix::getLocalDiagCopy(): the number of blocks in diag differ from the number of blocks in this operator." );
815  //XPETRA_TEST_FOR_EXCEPTION(bdiag->getMap()->isSameAs(*(getMap())) == false, Xpetra::Exceptions::RuntimeError,
816  // "BlockedCrsMatrix::getLocalDiagCopy(): the map of the vector diag is not compatible with the map of the blocked operator." );
817 
818  for (size_t row = 0; row < Rows(); row++) {
819  Teuchos::RCP<const Matrix> rm = getMatrix(row,row);
820  if (!rm.is_null()) {
821  Teuchos::RCP<Vector> rv = VectorFactory::Build(bdiag->getBlockedMap()->getMap(row,bdiag->getBlockedMap()->getThyraMode()));
822  rm->getLocalDiagCopy(*rv);
823  bdiag->setMultiVector(row,rv,bdiag->getBlockedMap()->getThyraMode());
824  }
825  }
826  }
827 
829  void leftScale (const Vector& x) {
830  XPETRA_MONITOR("XpetraBlockedCrsMatrix::leftScale");
831 
832  Teuchos::RCP<const Vector> rcpx = Teuchos::rcpFromRef(x);
833  Teuchos::RCP<const BlockedVector> bx = Teuchos::rcp_dynamic_cast<const BlockedVector>(rcpx);
834 
835  //RCP<Teuchos::FancyOStream> fancy = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
836 
837  // special treatment for 1xn block matrices
838  // ReorderedBlockedCrsMatrix object encapsulate single blocks in ReorderedBlocks
839  // BlockedVectors have Vector objects as Leaf objects.
840  if(Rows() == 1 && bx.is_null() == true) {
841  for (size_t col = 0; col < Cols(); ++col) {
842  Teuchos::RCP<Matrix> rm = getMatrix(0,col);
843  rm->leftScale(x);
844  }
845  return;
846  }
847 
848  TEUCHOS_TEST_FOR_EXCEPTION(bx.is_null() == true, Xpetra::Exceptions::RuntimeError, "BlockedCrsMatrix::leftScale: x must be a Blocked(Multi)Vector.");
849  TEUCHOS_TEST_FOR_EXCEPTION(bx->getNumVectors() != 1, Xpetra::Exceptions::RuntimeError, "BlockedCrsMatrix::leftScale: x must be a Blocked(Multi)Vector with exactly one vector. However, the number of stored vectors is " << bx->getNumVectors());
850  TEUCHOS_TEST_FOR_EXCEPTION(bx->getBlockedMap()->getNumMaps() != this->Rows(), Xpetra::Exceptions::RuntimeError,
851  "BlockedCrsMatrix::leftScale(): the number of blocks in diag differ from the number of blocks in this operator." );
852 
853  for (size_t row = 0; row < Rows(); row++) {
854  Teuchos::RCP<const MultiVector> rmv = bx->getMultiVector(row);
855  Teuchos::RCP<const Vector> rscale = rmv->getVector(0);
856  XPETRA_TEST_FOR_EXCEPTION(rscale.is_null()==true, Xpetra::Exceptions::RuntimeError, "BlockedCrsMatrix::leftScale: x must be a Vector.");
857  for (size_t col = 0; col < Cols(); ++col) {
858  Teuchos::RCP<Matrix> rm = getMatrix(row,col);
859  if (!rm.is_null()) {
860  rm->leftScale(*rscale);
861  }
862  }
863  }
864  }
865 
867  void rightScale (const Vector& x) {
868  XPETRA_MONITOR("XpetraBlockedCrsMatrix::rightScale");
869 
870  Teuchos::RCP<const Vector> rcpx = Teuchos::rcpFromRef(x);
871  Teuchos::RCP<const BlockedVector> bx = Teuchos::rcp_dynamic_cast<const BlockedVector>(rcpx);
872 
873  //RCP<Teuchos::FancyOStream> fancy = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
874 
875  // special treatment for nx1 block matrices
876  // ReorderedBlockedCrsMatrix object encapsulate single blocks in ReorderedBlocks
877  // BlockedVectors have Vector objects as Leaf objects.
878  if(Cols() == 1 && bx.is_null() == true) {
879  for (size_t row = 0; row < Rows(); ++row) {
880  Teuchos::RCP<Matrix> rm = getMatrix(row,0);
881  rm->rightScale(x);
882  }
883  return;
884  }
885 
886  TEUCHOS_TEST_FOR_EXCEPTION(bx.is_null() == true, Xpetra::Exceptions::RuntimeError, "BlockedCrsMatrix::rightScale: x must be a Blocked(Multi)Vector.");
887  TEUCHOS_TEST_FOR_EXCEPTION(bx->getNumVectors() != 1, Xpetra::Exceptions::RuntimeError, "BlockedCrsMatrix::rightScale: x must be a Blocked(Multi)Vector with exactly one vector. However, the number of stored vectors is " << bx->getNumVectors());
888  TEUCHOS_TEST_FOR_EXCEPTION(bx->getBlockedMap()->getNumMaps() != this->Cols(), Xpetra::Exceptions::RuntimeError,
889  "BlockedCrsMatrix::rightScale(): the number of blocks in diag differ from the number of blocks in this operator." );
890 
891  for (size_t col = 0; col < Cols(); ++col) {
892  Teuchos::RCP<const MultiVector> rmv = bx->getMultiVector(col);
893  Teuchos::RCP<const Vector> rscale = rmv->getVector(0);
894  XPETRA_TEST_FOR_EXCEPTION(rscale.is_null()==true, Xpetra::Exceptions::RuntimeError, "BlockedCrsMatrix::leftScale: x must be a Vector.");
895  for (size_t row = 0; row < Rows(); row++) {
896  Teuchos::RCP<Matrix> rm = getMatrix(row,col);
897  if (!rm.is_null()) {
898  rm->rightScale(*rscale);
899  }
900  }
901  }
902  }
903 
904 
906  virtual typename ScalarTraits<Scalar>::magnitudeType getFrobeniusNorm() const {
907  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getFrobeniusNorm");
908  typename ScalarTraits<Scalar>::magnitudeType ret = Teuchos::ScalarTraits<Scalar>::magnitude(Teuchos::ScalarTraits<Scalar>::zero());
909  for (size_t col = 0; col < Cols(); ++col) {
910  for (size_t row = 0; row < Rows(); ++row) {
911  if(getMatrix(row,col)!=Teuchos::null) {
912  typename ScalarTraits<Scalar>::magnitudeType n = getMatrix(row,col)->getFrobeniusNorm();
913  ret += n * n;
914  }
915  }
916  }
917  return Teuchos::ScalarTraits< typename ScalarTraits<Scalar>::magnitudeType >::squareroot(ret);
918  }
919 
920 
922  virtual bool haveGlobalConstants() const {return true;}
923 
924 
926 
928 
929 
931 
951 
953 
954 
956 
959  virtual void apply(const MultiVector& X, MultiVector& Y,
960  Teuchos::ETransp mode = Teuchos::NO_TRANS,
961  Scalar alpha = ScalarTraits<Scalar>::one(),
962  Scalar beta = ScalarTraits<Scalar>::zero()) const
963  {
964  XPETRA_MONITOR("XpetraBlockedCrsMatrix::apply");
965  //using Teuchos::RCP;
966 
967  TEUCHOS_TEST_FOR_EXCEPTION(mode != Teuchos::NO_TRANS && mode != Teuchos::TRANS, Xpetra::Exceptions::RuntimeError,
968  "apply() only supports the following modes: NO_TRANS and TRANS." );
969 
970  // check whether input parameters are blocked or not
971  RCP<const MultiVector> refX = rcpFromRef(X);
972  RCP<const BlockedMultiVector> refbX = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(refX);
973  //RCP<MultiVector> tmpY = rcpFromRef(Y);
974  //RCP<BlockedMultiVector> tmpbY = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(tmpY);
975 
976  // TODO get rid of me: adapt MapExtractor
977  bool bBlockedX = (refbX != Teuchos::null) ? true : false;
978 
979  // create (temporary) vectors for output
980  // In the end we call Y.update(alpha, *tmpY, beta). Therefore we need a new vector storing the temporary results
981  RCP<MultiVector> tmpY = MultiVectorFactory::Build(Y.getMap(), Y.getNumVectors(), true);
982 
983  //RCP<Teuchos::FancyOStream> out = rcp(new Teuchos::FancyOStream(rcp(&std::cout,false)));
984 
985  SC one = ScalarTraits<SC>::one();
986 
987  if (mode == Teuchos::NO_TRANS) {
988 
989  for (size_t row = 0; row < Rows(); row++) {
990  RCP<MultiVector> Yblock = rangemaps_->getVector(row, Y.getNumVectors(), bRangeThyraMode_, true);
991  for (size_t col = 0; col < Cols(); col++) {
992 
993  // extract matrix block
994  RCP<Matrix> Ablock = getMatrix(row, col);
995 
996  if (Ablock.is_null())
997  continue;
998 
999  // check whether Ablock is itself a blocked operator
1000  // If it is a blocked operator we have to provide Xpetra style GIDs, i.e. we have to transform GIDs
1001  bool bBlockedSubMatrix = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(Ablock) == Teuchos::null ? false : true;
1002 
1003  // input/output vectors for local block operation
1004  RCP<const MultiVector> Xblock = Teuchos::null; // subpart of X vector to be applied to subblock of A
1005 
1006  // extract sub part of X using Xpetra or Thyra GIDs
1007  // if submatrix is again blocked, we extract it using Xpetra style gids. If it is a single
1008  // block matrix we use the Thyra or Xpetra style GIDs that are used to store the matrix
1009  if(bBlockedX) Xblock = domainmaps_->ExtractVector(refbX, col, bDomainThyraMode_);
1010  else Xblock = domainmaps_->ExtractVector(refX, col, bBlockedSubMatrix == true ? false : bDomainThyraMode_);
1011 
1012  RCP<MultiVector> tmpYblock = rangemaps_->getVector(row, Y.getNumVectors(), bRangeThyraMode_, false); // subpart of Y vector containing part of solution of Xblock applied to Ablock
1013  Ablock->apply(*Xblock, *tmpYblock);
1014 
1015  Yblock->update(one, *tmpYblock, one);
1016  }
1017  rangemaps_->InsertVector(Yblock, row, tmpY, bRangeThyraMode_);
1018  }
1019 
1020  } else if (mode == Teuchos::TRANS) {
1021  // TODO: test me!
1022  for (size_t col = 0; col < Cols(); col++) {
1023  RCP<MultiVector> Yblock = domainmaps_->getVector(col, Y.getNumVectors(), bDomainThyraMode_, true);
1024 
1025  for (size_t row = 0; row<Rows(); row++) {
1026  RCP<Matrix> Ablock = getMatrix(row, col);
1027 
1028  if (Ablock.is_null())
1029  continue;
1030 
1031  // check whether Ablock is itself a blocked operator
1032  // If it is a blocked operator we have to provide Xpetra style GIDs, i.e. we have to transform GIDs
1033  bool bBlockedSubMatrix = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(Ablock) == Teuchos::null ? false : true;
1034 
1035  RCP<const MultiVector> Xblock = Teuchos::null;
1036 
1037  // extract sub part of X using Xpetra or Thyra GIDs
1038  if(bBlockedX) Xblock = rangemaps_->ExtractVector(refbX, row, bRangeThyraMode_);
1039  else Xblock = rangemaps_->ExtractVector(refX, row, bBlockedSubMatrix == true ? false : bRangeThyraMode_);
1040  RCP<MultiVector> tmpYblock = domainmaps_->getVector(col, Y.getNumVectors(), bDomainThyraMode_, false);
1041  Ablock->apply(*Xblock, *tmpYblock, Teuchos::TRANS);
1042 
1043  Yblock->update(one, *tmpYblock, one);
1044  }
1045  domainmaps_->InsertVector(Yblock, col, tmpY, bDomainThyraMode_);
1046  }
1047  }
1048  Y.update(alpha, *tmpY, beta);
1049  }
1050 
1052  RCP<const Map > getFullDomainMap() const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getFullDomainMap()"); return domainmaps_->getFullMap(); }
1053 
1055  RCP<const BlockedMap > getBlockedDomainMap() const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getBlockedDomainMap()"); return domainmaps_->getBlockedMap(); }
1056 
1058  RCP<const Map > getDomainMap() const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getDomainMap()"); return domainmaps_->getMap(); /*domainmaps_->getFullMap();*/ }
1059 
1061  RCP<const Map > getDomainMap(size_t i) const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getDomainMap(size_t)"); return domainmaps_->getMap(i, bDomainThyraMode_); }
1062 
1064  RCP<const Map > getDomainMap(size_t i, bool bThyraMode) const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getDomainMap(size_t,bool)"); return domainmaps_->getMap(i, bThyraMode); }
1065 
1067  RCP<const Map > getFullRangeMap() const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getRangeMap()"); return rangemaps_->getFullMap(); }
1068 
1070  RCP<const BlockedMap > getBlockedRangeMap() const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getBlockedRangeMap()"); return rangemaps_->getBlockedMap(); }
1071 
1073  RCP<const Map > getRangeMap() const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getRangeMap()"); return rangemaps_->getMap(); /*rangemaps_->getFullMap();*/ }
1074 
1076  RCP<const Map > getRangeMap(size_t i) const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getRangeMap(size_t)"); return rangemaps_->getMap(i, bRangeThyraMode_); }
1077 
1079  RCP<const Map > getRangeMap(size_t i, bool bThyraMode) const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getRangeMap(size_t,bool)"); return rangemaps_->getMap(i, bThyraMode); }
1080 
1082  RCP<const MapExtractor> getRangeMapExtractor() const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getRangeMapExtractor()"); return rangemaps_; }
1083 
1085  RCP<const MapExtractor> getDomainMapExtractor() const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getDomainMapExtractor()"); return domainmaps_; }
1086 
1088 
1090  //{@
1091 
1100  virtual void bgs_apply(
1101  const MultiVector& X,
1102  MultiVector& Y,
1103  size_t row,
1104  Teuchos::ETransp mode = Teuchos::NO_TRANS,
1105  Scalar alpha = ScalarTraits<Scalar>::one(),
1106  Scalar beta = ScalarTraits<Scalar>::zero()
1107  ) const
1108  {
1109  XPETRA_MONITOR("XpetraBlockedCrsMatrix::bgs_apply");
1110  //using Teuchos::RCP;
1111 
1112  TEUCHOS_TEST_FOR_EXCEPTION(mode != Teuchos::NO_TRANS && mode != Teuchos::TRANS, Xpetra::Exceptions::RuntimeError,
1113  "apply() only supports the following modes: NO_TRANS and TRANS." );
1114 
1115  // check whether input parameters are blocked or not
1116  RCP<const MultiVector> refX = rcpFromRef(X);
1117  RCP<const BlockedMultiVector> refbX = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(refX);
1118  //RCP<MultiVector> tmpY = rcpFromRef(Y);
1119  //RCP<BlockedMultiVector> tmpbY = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(tmpY);
1120 
1121  bool bBlockedX = (refbX != Teuchos::null) ? true : false;
1122 
1123  // create (temporary) vectors for output
1124  // In the end we call Y.update(alpha, *tmpY, beta). Therefore we need a new vector storing the temporary results
1125  RCP<MultiVector> tmpY = MultiVectorFactory::Build(Y.getMap(), Y.getNumVectors(), true);
1126 
1127  SC one = ScalarTraits<SC>::one();
1128 
1129  if (mode == Teuchos::NO_TRANS) {
1130  RCP<MultiVector> Yblock = rangemaps_->getVector(row, Y.getNumVectors(), bRangeThyraMode_, true);
1131  for (size_t col = 0; col < Cols(); col++) {
1132 
1133  // extract matrix block
1134  RCP<Matrix> Ablock = getMatrix(row, col);
1135 
1136  if (Ablock.is_null())
1137  continue;
1138 
1139  // check whether Ablock is itself a blocked operator
1140  // If it is a blocked operator we have to provide Xpetra style GIDs, i.e. we have to transform GIDs
1141  bool bBlockedSubMatrix = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(Ablock) == Teuchos::null ? false : true;
1142 
1143  // input/output vectors for local block operation
1144  RCP<const MultiVector> Xblock = Teuchos::null; // subpart of X vector to be applied to subblock of A
1145 
1146  // extract sub part of X using Xpetra or Thyra GIDs
1147  // if submatrix is again blocked, we extract it using Xpetra style gids. If it is a single
1148  // block matrix we use the Thyra or Xpetra style GIDs that are used to store the matrix
1149  if(bBlockedX) Xblock = domainmaps_->ExtractVector(refbX, col, bDomainThyraMode_);
1150  else Xblock = domainmaps_->ExtractVector(refX, col, bBlockedSubMatrix == true ? false : bDomainThyraMode_);
1151 
1152  RCP<MultiVector> tmpYblock = rangemaps_->getVector(row, Y.getNumVectors(), bRangeThyraMode_, false); // subpart of Y vector containing part of solution of Xblock applied to Ablock
1153  Ablock->apply(*Xblock, *tmpYblock);
1154 
1155  Yblock->update(one, *tmpYblock, one);
1156  }
1157  rangemaps_->InsertVector(Yblock, row, tmpY, bRangeThyraMode_);
1158  } else {
1159  TEUCHOS_TEST_FOR_EXCEPTION(true,Xpetra::Exceptions::NotImplemented,"Xpetar::BlockedCrsMatrix::bgs_apply: not implemented for transpose case.");
1160  }
1161  Y.update(alpha, *tmpY, beta);
1162  }
1163 
1164 
1166 
1167 
1169  //{@
1170 
1172  const Teuchos::RCP< const Map > getMap() const {
1173  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getMap");
1174  if (Rows() == 1 && Cols () == 1) {
1175  return getMatrix(0,0)->getMap();
1176  }
1177  throw Xpetra::Exceptions::RuntimeError("BlockedCrsMatrix::getMap(): operation not supported.");
1178  }
1179 
1181  void doImport(const Matrix &source, const Import& importer, CombineMode CM) {
1182  XPETRA_MONITOR("XpetraBlockedCrsMatrix::doImport");
1183  if (Rows() == 1 && Cols () == 1) {
1184  getMatrix(0,0)->doImport(source, importer, CM);
1185  return;
1186  }
1187  throw Xpetra::Exceptions::RuntimeError("BlockedCrsMatrix::doImport(): operation not supported.");
1188  }
1189 
1191  void doExport(const Matrix& dest, const Import& importer, CombineMode CM) {
1192  XPETRA_MONITOR("XpetraBlockedCrsMatrix::doExport");
1193  if (Rows() == 1 && Cols () == 1) {
1194  getMatrix(0,0)->doExport(dest, importer, CM);
1195  return;
1196  }
1197  throw Xpetra::Exceptions::RuntimeError("BlockedCrsMatrix::doExport(): operation not supported.");
1198  }
1199 
1201  void doImport(const Matrix& source, const Export& exporter, CombineMode CM) {
1202  XPETRA_MONITOR("XpetraBlockedCrsMatrix::doImport");
1203  if (Rows() == 1 && Cols () == 1) {
1204  getMatrix(0,0)->doImport(source, exporter, CM);
1205  return;
1206  }
1207  throw Xpetra::Exceptions::RuntimeError("BlockedCrsMatrix::doImport(): operation not supported.");
1208  }
1209 
1211  void doExport(const Matrix& dest, const Export& exporter, CombineMode CM) {
1212  XPETRA_MONITOR("XpetraBlockedCrsMatrix::doExport");
1213  if (Rows() == 1 && Cols () == 1) {
1214  getMatrix(0,0)->doExport(dest, exporter, CM);
1215  return;
1216  }
1217  throw Xpetra::Exceptions::RuntimeError("BlockedCrsMatrix::doExport(): operation not supported.");
1218  }
1219 
1220  // @}
1221 
1223 
1224 
1226  std::string description() const { return "Xpetra_BlockedCrsMatrix.description()"; }
1227 
1229  void describe(Teuchos::FancyOStream& out, const Teuchos::EVerbosityLevel verbLevel = Teuchos::Describable::verbLevel_default) const {
1230  out << "Xpetra::BlockedCrsMatrix: " << Rows() << " x " << Cols() << std::endl;
1231 
1232  if (isFillComplete()) {
1233  out << "BlockMatrix is fillComplete" << std::endl;
1234 
1235  /*if(fullrowmap_ != Teuchos::null) {
1236  out << "fullRowMap" << std::endl;
1237  fullrowmap_->describe(out,verbLevel);
1238  } else {
1239  out << "fullRowMap not set. Check whether block matrix is properly fillCompleted!" << std::endl;
1240  }*/
1241 
1242  //out << "fullColMap" << std::endl;
1243  //fullcolmap_->describe(out,verbLevel);
1244 
1245  } else {
1246  out << "BlockMatrix is NOT fillComplete" << std::endl;
1247  }
1248 
1249  for (size_t r = 0; r < Rows(); ++r)
1250  for (size_t c = 0; c < Cols(); ++c) {
1251  if(getMatrix(r,c)!=Teuchos::null) {
1252  out << "Block(" << r << "," << c << ")" << std::endl;
1253  getMatrix(r,c)->describe(out,verbLevel);
1254  } else out << "Block(" << r << "," << c << ") = null" << std::endl;
1255  }
1256  }
1257 
1259 
1260  void setObjectLabel( const std::string &objectLabel ) {
1261  XPETRA_MONITOR("TpetraBlockedCrsMatrix::setObjectLabel");
1262  for (size_t r = 0; r < Rows(); ++r)
1263  for (size_t c = 0; c < Cols(); ++c) {
1264  if(getMatrix(r,c)!=Teuchos::null) {
1265  std::ostringstream oss; oss<< objectLabel << "(" << r << "," << c << ")";
1266  getMatrix(r,c)->setObjectLabel(oss.str());
1267  }
1268  }
1269  }
1271 
1272 
1274  bool hasCrsGraph() const {
1275  if (Rows() == 1 && Cols () == 1) return true;
1276  else return false;
1277  }
1278 
1280  RCP<const CrsGraph> getCrsGraph() const {
1281  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getCrsGraph");
1282  if (Rows() == 1 && Cols () == 1) {
1283  return getMatrix(0,0)->getCrsGraph();
1284  }
1285  throw Xpetra::Exceptions::RuntimeError("getCrsGraph() not supported by BlockedCrsMatrix");
1286  }
1287 
1289 
1291 
1292 
1293  virtual bool isDiagonal() const {return is_diagonal_;}
1294 
1296  virtual size_t Rows() const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::Rows"); return rangemaps_->NumMaps(); }
1297 
1299  virtual size_t Cols() const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::Cols"); return domainmaps_->NumMaps(); }
1300 
1302  Teuchos::RCP<Matrix> getCrsMatrix() const {
1303  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getCrsMatrix");
1304  TEUCHOS_TEST_FOR_EXCEPTION(Rows()!=1, std::out_of_range, "Can only unwrap a 1x1 blocked matrix. The matrix has " << Rows() << " block rows, though.");
1305  TEUCHOS_TEST_FOR_EXCEPTION(Cols()!=1, std::out_of_range, "Can only unwrap a 1x1 blocked matrix. The matrix has " << Cols() << " block columns, though.");
1306 
1307  RCP<Matrix> mat = getMatrix(0,0);
1308  RCP<BlockedCrsMatrix> bmat = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(mat);
1309  if (bmat == Teuchos::null) return mat;
1310  return bmat->getCrsMatrix();
1311  }
1312 
1314  Teuchos::RCP<Matrix> getInnermostCrsMatrix() {
1315  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getInnermostCrsMatrix");
1316  size_t row = Rows()+1, col = Cols()+1;
1317  for (size_t r = 0; r < Rows(); ++r)
1318  for(size_t c = 0; c < Cols(); ++c)
1319  if (getMatrix(r,c) != Teuchos::null) {
1320  row = r;
1321  col = c;
1322  break;
1323  }
1324  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.")
1325  RCP<Matrix> mm = getMatrix(row,col);
1326  RCP<BlockedCrsMatrix> bmat = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(mm);
1327  if (bmat == Teuchos::null) return mm;
1328  return bmat->getInnermostCrsMatrix();
1329  }
1330 
1332  Teuchos::RCP<Matrix> getMatrix(size_t r, size_t c) const {
1333  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getMatrix");
1334  TEUCHOS_TEST_FOR_EXCEPTION(r > Rows(), std::out_of_range, "Error, r = " << Rows() << " is too big");
1335  TEUCHOS_TEST_FOR_EXCEPTION(c > Cols(), std::out_of_range, "Error, c = " << Cols() << " is too big");
1336 
1337  // transfer strided/blocked map information
1338  /* if (blocks_[r*Cols()+c] != Teuchos::null &&
1339  blocks_[r*Cols()+c]->IsView("stridedMaps") == false)
1340  blocks_[r*Cols()+c]->CreateView("stridedMaps", getRangeMap(r,bRangeThyraMode_), getDomainMap(c,bDomainThyraMode_));*/
1341  return blocks_[r*Cols()+c];
1342  }
1343 
1345  //void setMatrix(size_t r, size_t c, Teuchos::RCP<CrsMatrix> mat) {
1346  void setMatrix(size_t r, size_t c, Teuchos::RCP<Matrix> mat) {
1347  XPETRA_MONITOR("XpetraBlockedCrsMatrix::setMatrix");
1348  // TODO: if filled -> return error
1349 
1350  TEUCHOS_TEST_FOR_EXCEPTION(r > Rows(), std::out_of_range, "Error, r = " << Rows() << " is too big");
1351  TEUCHOS_TEST_FOR_EXCEPTION(c > Cols(), std::out_of_range, "Error, c = " << Cols() << " is too big");
1352  if(!mat.is_null() && r!=c) is_diagonal_ = false;
1353  // set matrix
1354  blocks_[r*Cols() + c] = mat;
1355  }
1356 
1358  // NOTE: This is a rather expensive operation, since all blocks are copied
1359  // into a new big CrsMatrix
1360  Teuchos::RCP<Matrix> Merge() const {
1361  XPETRA_MONITOR("XpetraBlockedCrsMatrix::Merge");
1362  using Teuchos::RCP;
1363  using Teuchos::rcp_dynamic_cast;
1364  Scalar one = ScalarTraits<SC>::one();
1365 
1367  "BlockedCrsMatrix::Merge: only implemented for Xpetra-style or Thyra-style numbering. No mixup allowed!" );
1368 
1369  TEUCHOS_TEST_FOR_EXCEPTION(isFillComplete() == false, Xpetra::Exceptions::RuntimeError,
1370  "BlockedCrsMatrix::Merge: BlockMatrix must be fill-completed." );
1371 
1372  LocalOrdinal lclNumRows = getFullRangeMap()->getNodeNumElements();
1373  Teuchos::ArrayRCP<size_t> numEntPerRow (lclNumRows);
1374  for (LocalOrdinal lclRow = 0; lclRow < lclNumRows; ++lclRow)
1375  numEntPerRow[lclRow] = getNumEntriesInLocalRow(lclRow);
1376 
1377  RCP<Matrix> sparse = MatrixFactory::Build(getFullRangeMap(), numEntPerRow);
1378 
1379  if(bRangeThyraMode_ == false) {
1380  // Xpetra mode
1381  for (size_t i = 0; i < Rows(); i++) {
1382  for (size_t j = 0; j < Cols(); j++) {
1383  if (getMatrix(i,j) != Teuchos::null) {
1384  RCP<const Matrix> mat = getMatrix(i,j);
1385 
1386  // recursively call Merge routine
1387  RCP<const BlockedCrsMatrix> bMat = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(mat);
1388  if (bMat != Teuchos::null) mat = bMat->Merge();
1389 
1390  bMat = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(mat);
1391  TEUCHOS_TEST_FOR_EXCEPTION(bMat != Teuchos::null, Xpetra::Exceptions::RuntimeError,
1392  "BlockedCrsMatrix::Merge: Merging of blocked sub-operators failed?!" );
1393 
1394  // jump over empty blocks
1395  if(mat->getNodeNumEntries() == 0) continue;
1396 
1397  this->Add(*mat, one, *sparse, one);
1398  }
1399  }
1400  }
1401  } else {
1402  // Thyra mode
1403  for (size_t i = 0; i < Rows(); i++) {
1404  for (size_t j = 0; j < Cols(); j++) {
1405  if (getMatrix(i,j) != Teuchos::null) {
1406  RCP<const Matrix> mat = getMatrix(i,j);
1407  // recursively call Merge routine
1408  RCP<const BlockedCrsMatrix> bMat = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(mat);
1409  if (bMat != Teuchos::null) mat = bMat->Merge();
1410 
1411  bMat = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(mat);
1412  TEUCHOS_TEST_FOR_EXCEPTION(bMat != Teuchos::null, Xpetra::Exceptions::RuntimeError,
1413  "BlockedCrsMatrix::Merge: Merging of blocked sub-operators failed?!" );
1414 
1415  // check whether we have a CrsMatrix block (no blocked operator)
1416  RCP<const CrsMatrixWrap> crsMat = Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(mat);
1417  TEUCHOS_ASSERT(crsMat != Teuchos::null);
1418 
1419  // these are the thyra style maps of the matrix
1420  RCP<const Map> trowMap = mat->getRowMap();
1421  RCP<const Map> tcolMap = mat->getColMap();
1422  RCP<const Map> tdomMap = mat->getDomainMap();
1423 
1424  // get Xpetra maps
1425  RCP<const Map> xrowMap = getRangeMapExtractor()->getMap(i,false);
1426  RCP<const Map> xdomMap = getDomainMapExtractor()->getMap(j,false);
1427 
1428  // generate column map with Xpetra GIDs
1429  // We have to do this separately for each block since the column
1430  // map of each block might be different in the same block column
1431  Teuchos::RCP<Map> xcolMap = MapUtils::transformThyra2XpetraGIDs(
1432  *tcolMap,
1433  *tdomMap,
1434  *xdomMap);
1435 
1436  // jump over empty blocks
1437  if(mat->getNodeNumEntries() == 0) continue;
1438 
1439  size_t maxNumEntries = mat->getNodeMaxNumRowEntries();
1440 
1441  size_t numEntries;
1442  Array<GO> inds (maxNumEntries);
1443  Array<GO> inds2(maxNumEntries);
1444  Array<SC> vals (maxNumEntries);
1445 
1446  // loop over all rows and add entries
1447  for (size_t k = 0; k < mat->getNodeNumRows(); k++) {
1448  GlobalOrdinal rowTGID = trowMap->getGlobalElement(k);
1449  crsMat->getCrsMatrix()->getGlobalRowCopy(rowTGID, inds(), vals(), numEntries);
1450 
1451  // create new indices array
1452  for (size_t l = 0; l < numEntries; ++l) {
1453  LocalOrdinal lid = tcolMap->getLocalElement(inds[l]);
1454  inds2[l] = xcolMap->getGlobalElement(lid);
1455  }
1456 
1457  GlobalOrdinal rowXGID = xrowMap->getGlobalElement(k);
1458  sparse->insertGlobalValues(
1459  rowXGID, inds2(0, numEntries),
1460  vals(0, numEntries));
1461  }
1462  }
1463  }
1464  }
1465  }
1466 
1467  sparse->fillComplete(getFullDomainMap(), getFullRangeMap());
1468 
1469  TEUCHOS_TEST_FOR_EXCEPTION(sparse->getNodeNumEntries() != getNodeNumEntries(), Xpetra::Exceptions::RuntimeError,
1470  "BlockedCrsMatrix::Merge: Local number of entries of merged matrix does not coincide with local number of entries of blocked operator." );
1471 
1472  TEUCHOS_TEST_FOR_EXCEPTION(sparse->getGlobalNumEntries() != getGlobalNumEntries(), Xpetra::Exceptions::RuntimeError,
1473  "BlockedCrsMatrix::Merge: Global number of entries of merged matrix does not coincide with global number of entries of blocked operator." );
1474 
1475  return sparse;
1476  }
1478 
1479 #ifdef HAVE_XPETRA_KOKKOS_REFACTOR
1480  typedef typename CrsMatrix::local_matrix_type local_matrix_type;
1481 
1483  local_matrix_type getLocalMatrix () const {
1484  if (Rows() == 1 && Cols () == 1) {
1485  return getMatrix(0,0)->getLocalMatrix();
1486  }
1487  throw Xpetra::Exceptions::RuntimeError("BlockedCrsMatrix::getLocalMatrix(): operation not supported.");
1488  }
1489 #endif
1490 
1491 #ifdef HAVE_XPETRA_THYRA
1492  Teuchos::RCP<Thyra::BlockedLinearOpBase<Scalar> > getThyraOperator() {
1493  Teuchos::RCP<Thyra::LinearOpBase<Scalar> > thOp =
1494  Xpetra::ThyraUtils<Scalar,LO,GO,Node>::toThyra(Teuchos::rcpFromRef(*this));
1495  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thOp));
1496 
1497  Teuchos::RCP<Thyra::BlockedLinearOpBase<Scalar> > thbOp =
1498  Teuchos::rcp_dynamic_cast<Thyra::BlockedLinearOpBase<Scalar> >(thOp);
1499  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thbOp));
1500  return thbOp;
1501  }
1502 #endif
1503 
1505  void residual(const MultiVector & X,
1506  const MultiVector & B,
1507  MultiVector & R) const {
1508  using STS = Teuchos::ScalarTraits<Scalar>;
1509  R.update(STS::one(),B,STS::zero());
1510  this->apply (X, R, Teuchos::NO_TRANS, -STS::one(), STS::one());
1511  }
1512 
1513  private:
1514 
1517 
1519 
1529  void Add(const Matrix& A, const Scalar scalarA, Matrix& B, const Scalar scalarB) const {
1530  XPETRA_MONITOR("XpetraBlockedCrsMatrix::Add");
1531  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete(), Xpetra::Exceptions::RuntimeError,
1532  "Matrix A is not completed");
1533  using Teuchos::Array;
1534  using Teuchos::ArrayView;
1535 
1536  B.scale(scalarB);
1537 
1538  Scalar one = ScalarTraits<SC>::one();
1539  Scalar zero = ScalarTraits<SC>::zero();
1540 
1541  if (scalarA == zero)
1542  return;
1543 
1544  Teuchos::RCP<const Matrix> rcpA = Teuchos::rcpFromRef(A);
1545  Teuchos::RCP<const CrsMatrixWrap> rcpAwrap = Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(rcpA);
1546  TEUCHOS_TEST_FOR_EXCEPTION(rcpAwrap == Teuchos::null, Xpetra::Exceptions::BadCast,
1547  "BlockedCrsMatrix::Add: matrix A must be of type CrsMatrixWrap.");
1548  Teuchos::RCP<const CrsMatrix> crsA = rcpAwrap->getCrsMatrix();
1549 
1550  size_t maxNumEntries = crsA->getNodeMaxNumRowEntries();
1551 
1552  size_t numEntries;
1553  Array<GO> inds(maxNumEntries);
1554  Array<SC> vals(maxNumEntries);
1555 
1556  RCP<const Map> rowMap = crsA->getRowMap();
1557  RCP<const Map> colMap = crsA->getColMap();
1558 
1559  ArrayView<const GO> rowGIDs = crsA->getRowMap()->getNodeElementList();
1560  for (size_t i = 0; i < crsA->getNodeNumRows(); i++) {
1561  GO row = rowGIDs[i];
1562  crsA->getGlobalRowCopy(row, inds(), vals(), numEntries);
1563 
1564  if (scalarA != one)
1565  for (size_t j = 0; j < numEntries; ++j)
1566  vals[j] *= scalarA;
1567 
1568  B.insertGlobalValues(row, inds(0, numEntries), vals(0, numEntries)); // insert should be ok, since blocks in BlockedCrsOpeartor do not overlap!
1569  }
1570  }
1571 
1573 
1574  // Default view is created after fillComplete()
1575  // Because ColMap might not be available before fillComplete().
1577  XPETRA_MONITOR("XpetraBlockedCrsMatrix::CreateDefaultView");
1578 
1579  // Create default view
1580  this->defaultViewLabel_ = "point";
1582 
1583  // Set current view
1584  this->currentViewLabel_ = this->GetDefaultViewLabel();
1585  }
1586 
1587  private:
1588  bool is_diagonal_; // If we're diagonal a bunch of the extraction stuff should work
1589  Teuchos::RCP<const MapExtractor> domainmaps_; // full domain map together with all partial domain maps
1590  Teuchos::RCP<const MapExtractor> rangemaps_; // full range map together with all partial domain maps
1591 
1592  std::vector<Teuchos::RCP<Matrix> > blocks_; // row major matrix block storage
1593 #ifdef HAVE_XPETRA_THYRA
1594  Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar> > thyraOp_;
1595 #endif
1598 
1599 };
1600 
1601 } //namespace Xpetra
1602 
1603 #define XPETRA_BLOCKEDCRSMATRIX_SHORT
1604 #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 Map > getDomainMap(size_t i, bool bThyraMode) const
Returns the Map associated with the i&#39;th block domain of this operator.
BlockedCrsMatrix(Teuchos::RCP< const MapExtractor > &rangeMaps, Teuchos::RCP< const MapExtractor > &domainMaps, size_t npr)
Constructor.
RCP< const Map > getFullDomainMap() const
Returns the Map associated with the full domain of this operator.
static Teuchos::RCP< Xpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Build(const Teuchos::RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node >> &map, bool zeroOut=true)
Constructor specifying the number of non-zeros for all rows.
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 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.
virtual void setAllToScalar(const Scalar &alpha)
Set all matrix entries equal to scalar.
GlobalOrdinal GO
RCP< const Map > getRangeMap() const
Returns the Map associated with the range of this operator.
RCP< const Map > getRangeMap(size_t i, bool bThyraMode) const
Returns the Map associated with the i&#39;th block range of this operator.
Teuchos::RCP< const MapExtractor > rangemaps_
void doImport(const Matrix &source, const Import &importer, CombineMode CM)
Import.
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.
static RCP< Matrix > Build(const RCP< const Map > &rowMap)
virtual void bgs_apply(const MultiVector &X, MultiVector &Y, size_t row, Teuchos::ETransp mode=Teuchos::NO_TRANS, Scalar alpha=ScalarTraits< Scalar >::one(), Scalar beta=ScalarTraits< Scalar >::zero()) const
Special multiplication routine (for BGS/Jacobi smoother)
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 > getRangeMap(size_t i) const
Returns the Map associated with the i&#39;th block range of this operator.
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.
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.
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.
void rightScale(const Vector &x)
Right scale matrix using the given vector entries.
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.
Exception indicating invalid cast attempted.
virtual bool haveGlobalConstants() const
Returns true if globalConstants have been computed; false otherwise.
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 indices.
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.
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.
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)
Map constructor with Xpetra-defined contiguous uniform distribution.
void setObjectLabel(const std::string &objectLabel)
RCP< const BlockedMap > getBlockedDomainMap() const
Returns the BlockedMap associated with the domain of this operator.
Teuchos::RCP< Matrix > getInnermostCrsMatrix()
helper routine recursively returns the first inner-most non-null matrix block from a (nested) blocked...
void insertLocalValues(LocalOrdinal localRow, const ArrayView< const LocalOrdinal > &cols, const ArrayView< const Scalar > &vals)
Insert matrix entries, using local IDs.
Teuchos::RCP< Matrix > getMatrix(size_t r, size_t c) const
return block (r,c)
virtual bool isFillComplete() const =0
Returns true if fillComplete() has been called and the matrix is in compute mode. ...
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
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
RCP< const Map > getFullRangeMap() const
Returns the Map associated with the full range of this operator.
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)
void residual(const MultiVector &X, const MultiVector &B, MultiVector &R) const
Compute a residual R = B - (*this) * X.
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).
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.
Exception throws when you call an unimplemented method of Xpetra.
bool hasCrsGraph() const
Supports the getCrsGraph() call.
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.
std::vector< Teuchos::RCP< Matrix > > blocks_
Teuchos::RCP< Matrix > Merge() const
merge BlockedCrsMatrix blocks in a CrsMatrix
RCP< const MapExtractor > getDomainMapExtractor() const
Returns map extractor for domain map.
global_size_t getGlobalNumRows() const
Returns the number of global rows.
RCP< const CrsGraph< int, GlobalOrdinal, Node > > toXpetra(const Epetra_CrsGraph &g)
#define XPETRA_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
size_t getNumEntriesInGlobalRow(GlobalOrdinal globalRow) const
Returns the current number of entries in the specified (locally owned) global row.
void setMatrix(size_t r, size_t c, Teuchos::RCP< Matrix > mat)
set matrix block
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...
bool isGloballyIndexed() const
If matrix indices are in the global range, this function returns true. Otherwise, this function retur...
Exception throws to report incompatible objects (like maps).
RCP< const Map > getDomainMap() const
Returns the Map associated with the domain of this operator.
Concrete implementation of Xpetra::Matrix.
virtual size_t Rows() const
number of row blocks
BlockedCrsMatrix(const Teuchos::RCP< const BlockedMap > &rangeMaps, const Teuchos::RCP< const BlockedMap > &domainMaps, size_t npr)
Constructor.
CombineMode
Xpetra::Combine Mode enumerable type.
#define XPETRA_MONITOR(funcName)
virtual const RCP< const Map > & getRowMap() const
Returns the Map that describes the row distribution in this matrix.
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.
const viewLabel_t & GetDefaultViewLabel() const
RCP< const BlockedMap > getBlockedRangeMap() const
Returns the BlockedMap associated with the range of this operator.
RCP< const Map > getDomainMap(size_t i) const
Returns the Map associated with the i&#39;th block domain of this operator.
virtual size_t getNumVectors() const =0
Number of columns in the multivector.
size_t getNodeNumRows() const
Returns the number of matrix rows owned on the calling node.
virtual ScalarTraits< Scalar >::magnitudeType getFrobeniusNorm() const
Get Frobenius norm of the matrix.
void removeEmptyProcessesInPlace(const Teuchos::RCP< const Map > &newMap)
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...
RCP< const MapExtractor > getRangeMapExtractor() const
Returns map extractor class for range map.