All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups 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 
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_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 
86 namespace Xpetra {
87 
88  typedef std::string viewLabel_t;
89 
90  template <class Scalar,
91  class LocalOrdinal,
92  class GlobalOrdinal,
93  class Node>
95  public Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> {
96  public:
97  typedef Scalar scalar_type;
98  typedef LocalOrdinal local_ordinal_type;
99  typedef GlobalOrdinal global_ordinal_type;
100  typedef Node node_type;
101 
102  private:
103 #undef XPETRA_BLOCKEDCRSMATRIX_SHORT
104 #include "Xpetra_UseShortNames.hpp"
105 
106  public:
107 
109 
110 
112 
119  const Teuchos::RCP<const BlockedMap>& domainMaps,
120  size_t npr, Xpetra::ProfileType pftype = Xpetra::DynamicProfile) {
121  is_diagonal_ = true;
122  domainmaps_ = Teuchos::rcp(new MapExtractor(domainMaps));
123  rangemaps_ = Teuchos::rcp(new MapExtractor(rangeMaps));
124  bRangeThyraMode_ = rangeMaps->getThyraMode();
125  bDomainThyraMode_ = domainMaps->getThyraMode();
126 
127  blocks_.reserve(Rows()*Cols());
128 
129  // add CrsMatrix objects in row,column order
130  for (size_t r = 0; r < Rows(); ++r)
131  for (size_t c = 0; c < Cols(); ++c)
132  blocks_.push_back(MatrixFactory::Build(getRangeMap(r,bRangeThyraMode_), npr, pftype));
133 
134  // Default view
136  }
137 
139 
149  size_t npr, Xpetra::ProfileType pftype = Xpetra::DynamicProfile)
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)
160  blocks_.push_back(MatrixFactory::Build(getRangeMap(r,bRangeThyraMode_), npr, pftype));
161 
162  // Default view
164  }
165 
166 #ifdef HAVE_XPETRA_THYRA
167 
174  BlockedCrsMatrix(const Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar> >& thyraOp, const Teuchos::RCP<const Teuchos::Comm<int> >& /* comm */)
175  : is_diagonal_(true), thyraOp_(thyraOp)
176  {
177  // extract information from Thyra blocked operator and rebuilt information
178  const Teuchos::RCP<const Thyra::ProductVectorSpaceBase<Scalar> > productRangeSpace = thyraOp->productRange();
179  const Teuchos::RCP<const Thyra::ProductVectorSpaceBase<Scalar> > productDomainSpace = thyraOp->productDomain();
180 
181  int numRangeBlocks = productRangeSpace->numBlocks();
182  int numDomainBlocks = productDomainSpace->numBlocks();
183 
184  // build range map extractor from Thyra::BlockedLinearOpBase object
185  std::vector<Teuchos::RCP<const Map> > subRangeMaps(numRangeBlocks);
186  for (size_t r=0; r<Teuchos::as<size_t>(numRangeBlocks); ++r) {
187  for (size_t c=0; c<Teuchos::as<size_t>(numDomainBlocks); ++c) {
188  if (thyraOp->blockExists(r,c)) {
189  // we only need at least one block in each block row to extract the range map
190  Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > const_op = thyraOp->getBlock(r,c); // nonConst access is not allowed.
193  subRangeMaps[r] = xop->getRangeMap();
194  if(r!=c) is_diagonal_ = false;
195  break;
196  }
197  }
198  }
199  Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > fullRangeMap = mergeMaps(subRangeMaps);
200 
201  // check whether the underlying Thyra operator uses Thyra-style numbering (default for most applications) or
202  // Xpetra style numbering
203  bool bRangeUseThyraStyleNumbering = false;
204  size_t numAllElements = 0;
205  for(size_t v = 0; v < subRangeMaps.size(); ++v) {
206  numAllElements += subRangeMaps[v]->getGlobalNumElements();
207  }
208  if ( fullRangeMap->getGlobalNumElements() != numAllElements) bRangeUseThyraStyleNumbering = true;
209  rangemaps_ = Teuchos::rcp(new Xpetra::MapExtractor<Scalar,LocalOrdinal,GlobalOrdinal,Node>(fullRangeMap, subRangeMaps, bRangeUseThyraStyleNumbering));
210 
211  // build domain map extractor from Thyra::BlockedLinearOpBase object
212  std::vector<Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > > subDomainMaps(numDomainBlocks);
213  for (size_t c=0; c<Teuchos::as<size_t>(numDomainBlocks); ++c) {
214  for (size_t r=0; r<Teuchos::as<size_t>(numRangeBlocks); ++r) {
215  if (thyraOp->blockExists(r,c)) {
216  // we only need at least one block in each block row to extract the range map
217  Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > const_op = thyraOp->getBlock(r,c); // nonConst access is not allowed.
220  subDomainMaps[c] = xop->getDomainMap();
221  break;
222  }
223  }
224  }
225  Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > fullDomainMap = mergeMaps(subDomainMaps);
226  // plausibility check for numbering style (Xpetra or Thyra)
227  bool bDomainUseThyraStyleNumbering = false;
228  numAllElements = 0;
229  for(size_t v = 0; v < subDomainMaps.size(); ++v) {
230  numAllElements += subDomainMaps[v]->getGlobalNumElements();
231  }
232  if (fullDomainMap->getGlobalNumElements() != numAllElements) bDomainUseThyraStyleNumbering = true;
233  domainmaps_ = Teuchos::rcp(new Xpetra::MapExtractor<Scalar,LocalOrdinal,GlobalOrdinal,Node>(fullDomainMap, subDomainMaps, bDomainUseThyraStyleNumbering));
234 
235  // store numbering mode
236  bRangeThyraMode_ = bRangeUseThyraStyleNumbering;
237  bDomainThyraMode_ = bDomainUseThyraStyleNumbering;
238 
239  // add CrsMatrix objects in row,column order
240  blocks_.reserve(Rows()*Cols());
241  for (size_t r = 0; r < Rows(); ++r) {
242  for (size_t c = 0; c < Cols(); ++c) {
243  if(thyraOp->blockExists(r,c)) {
244  // TODO we do not support nested Thyra operators here!
245  Teuchos::RCP<const Thyra::LinearOpBase<Scalar> > const_op = thyraOp->getBlock(r,c); // nonConst access is not allowed.
246  Teuchos::RCP<Thyra::LinearOpBase<Scalar> > op = Teuchos::rcp_const_cast<Thyra::LinearOpBase<Scalar> >(const_op); // cast away const
251  blocks_.push_back(xwrap);
252  } else {
253  // add empty block
255  }
256  }
257  }
258  // Default view
260  }
261 
262  private:
264 
271  Teuchos::RCP<const Map> mergeMaps(std::vector<Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > > & subMaps) {
272  // TODO merging for Thyra mode is missing (similar to what we do in constructor of MapExtractor
273 
274  // merge submaps to global map
275  std::vector<GlobalOrdinal> gids;
276  for(size_t tt = 0; tt<subMaps.size(); ++tt) {
278 #if 1
279  Teuchos::ArrayView< const GlobalOrdinal > subMapGids = subMap->getNodeElementList();
280  gids.insert(gids.end(), subMapGids.begin(), subMapGids.end());
281 #else
282  size_t myNumElements = subMap->getNodeNumElements();
283  for(LocalOrdinal l = 0; l < Teuchos::as<LocalOrdinal>(myNumElements); ++l) {
284  GlobalOrdinal gid = subMap->getGlobalElement(l);
285  gids.push_back(gid);
286  }
287 #endif
288  }
289 
290  // we have to sort the matrix entries and get rid of the double entries
291  // since we use this to detect Thyra-style numbering or Xpetra-style
292  // numbering. In Thyra-style numbering mode, the Xpetra::MapExtractor builds
293  // the correct row maps.
295  std::sort(gids.begin(), gids.end());
296  gids.erase(std::unique(gids.begin(), gids.end()), gids.end());
297  Teuchos::ArrayView<GO> gidsView(&gids[0], gids.size());
298  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());
299  return fullMap;
300  }
301 
302  public:
303 #endif
304 
306  virtual ~BlockedCrsMatrix() {}
307 
309 
310 
312 
313 
315 
340  void insertGlobalValues(GlobalOrdinal globalRow, const ArrayView<const GlobalOrdinal>& cols, const ArrayView<const Scalar>& vals) {
341  XPETRA_MONITOR("XpetraBlockedCrsMatrix::insertGlobalValues");
342  if (Rows() == 1 && Cols () == 1) {
343  getMatrix(0,0)->insertGlobalValues(globalRow, cols, vals);
344  return;
345  }
346  throw Xpetra::Exceptions::RuntimeError("insertGlobalValues not supported by BlockedCrsMatrix");
347  }
348 
350 
360  void insertLocalValues(LocalOrdinal localRow, const ArrayView<const LocalOrdinal>& cols, const ArrayView<const Scalar>& vals) {
361  XPETRA_MONITOR("XpetraBlockedCrsMatrix::insertLocalValues");
362  if (Rows() == 1 && Cols () == 1) {
363  getMatrix(0,0)->insertLocalValues(localRow, cols, vals);
364  return;
365  }
366  throw Xpetra::Exceptions::RuntimeError("insertLocalValues not supported by BlockedCrsMatrix");
367  }
368 
370  XPETRA_MONITOR("XpetraBlockedCrsMatrix::removeEmptyProcessesInPlace");
371  if (Rows() == 1 && Cols () == 1) {
372  getMatrix(0,0)->removeEmptyProcessesInPlace(newMap);
373  return;
374  }
375  throw Xpetra::Exceptions::RuntimeError("removeEmptyProcesses not supported by BlockedCrsMatrix");
376  }
377 
379 
387  void replaceGlobalValues(GlobalOrdinal globalRow,
388  const ArrayView<const GlobalOrdinal> &cols,
389  const ArrayView<const Scalar> &vals) {
390  XPETRA_MONITOR("XpetraBlockedCrsMatrix::replaceGlobalValues");
391  if (Rows() == 1 && Cols () == 1) {
392  getMatrix(0,0)->replaceGlobalValues(globalRow,cols,vals);
393  return;
394  }
395  throw Xpetra::Exceptions::RuntimeError("replaceGlobalValues not supported by BlockedCrsMatrix");
396  }
397 
399 
403  void replaceLocalValues(LocalOrdinal localRow,
404  const ArrayView<const LocalOrdinal> &cols,
405  const ArrayView<const Scalar> &vals) {
406  XPETRA_MONITOR("XpetraBlockedCrsMatrix::replaceLocalValues");
407  if (Rows() == 1 && Cols () == 1) {
408  getMatrix(0,0)->replaceLocalValues(localRow,cols,vals);
409  return;
410  }
411  throw Xpetra::Exceptions::RuntimeError("replaceLocalValues not supported by BlockedCrsMatrix");
412  }
413 
415  virtual void setAllToScalar(const Scalar& alpha) {
416  XPETRA_MONITOR("XpetraBlockedCrsMatrix::setAllToScalar");
417  for (size_t row = 0; row < Rows(); row++) {
418  for (size_t col = 0; col < Cols(); col++) {
419  if (!getMatrix(row,col).is_null()) {
420  getMatrix(row,col)->setAllToScalar(alpha);
421  }
422  }
423  }
424  }
425 
427  void scale(const Scalar& alpha) {
428  XPETRA_MONITOR("XpetraBlockedCrsMatrix::scale");
429  for (size_t row = 0; row < Rows(); row++) {
430  for (size_t col = 0; col < Cols(); col++) {
431  if (!getMatrix(row,col).is_null()) {
432  getMatrix(row,col)->scale(alpha);
433  }
434  }
435  }
436  }
437 
439 
441 
442 
450  void resumeFill(const RCP< ParameterList >& params = null) {
451  XPETRA_MONITOR("XpetraBlockedCrsMatrix::resumeFill");
452  for (size_t row = 0; row < Rows(); row++) {
453  for (size_t col = 0; col < Cols(); col++) {
454  if (!getMatrix(row,col).is_null()) {
455  getMatrix(row,col)->resumeFill(params);
456  }
457  }
458  }
459  }
460 
466  void fillComplete(const RCP<const Map>& domainMap, const RCP<const Map>& rangeMap, const RCP<ParameterList>& params = null) {
467  XPETRA_MONITOR("XpetraBlockedCrsMatrix::fillComplete");
468  if (Rows() == 1 && Cols () == 1) {
469  getMatrix(0,0)->fillComplete(domainMap, rangeMap, params);
470  return;
471  }
472  fillComplete(params);
473  }
474 
489  void fillComplete(const RCP<ParameterList>& params = null) {
490  XPETRA_MONITOR("XpetraBlockedCrsMatrix::fillComplete");
491  TEUCHOS_TEST_FOR_EXCEPTION(rangemaps_==Teuchos::null, Xpetra::Exceptions::RuntimeError,"BlockedCrsMatrix::fillComplete: rangemaps_ is not set. Error.");
492 
493  for (size_t r = 0; r < Rows(); ++r)
494  for (size_t c = 0; c < Cols(); ++c) {
495  if(getMatrix(r,c) != Teuchos::null) {
496  Teuchos::RCP<Matrix> Ablock = getMatrix(r,c);
497  if(r!=c) is_diagonal_ = false;
498  if (Ablock != Teuchos::null && !Ablock->isFillComplete())
499  Ablock->fillComplete(getDomainMap(c, bDomainThyraMode_), getRangeMap(r, bRangeThyraMode_), params);
500  }
501  }
502 
503 #if 0
504  // get full row map
505  RCP<const Map> rangeMap = rangemaps_->getFullMap();
506  fullrowmap_ = MapFactory::Build(rangeMap()->lib(), rangeMap()->getGlobalNumElements(), rangeMap()->getNodeElementList(), rangeMap()->getIndexBase(), rangeMap()->getComm());
507 
508  // build full col map
509  fullcolmap_ = Teuchos::null; // delete old full column map
510 
511  std::vector<GO> colmapentries;
512  for (size_t c = 0; c < Cols(); ++c) {
513  // copy all local column lids of all block rows to colset
514  std::set<GO> colset;
515  for (size_t r = 0; r < Rows(); ++r) {
516  Teuchos::RCP<CrsMatrix> Ablock = getMatrix(r,c);
517 
518  if (Ablock != Teuchos::null) {
519  Teuchos::ArrayView<const GO> colElements = Ablock->getColMap()->getNodeElementList();
520  Teuchos::RCP<const Map> colmap = Ablock->getColMap();
521  copy(colElements.getRawPtr(), colElements.getRawPtr() + colElements.size(), inserter(colset, colset.begin()));
522  }
523  }
524 
525  // remove duplicates (entries which are in column maps of more than one block row)
526  colmapentries.reserve(colmapentries.size() + colset.size());
527  copy(colset.begin(), colset.end(), back_inserter(colmapentries));
528  sort(colmapentries.begin(), colmapentries.end());
529  typename std::vector<GO>::iterator gendLocation;
530  gendLocation = std::unique(colmapentries.begin(), colmapentries.end());
531  colmapentries.erase(gendLocation,colmapentries.end());
532  }
533 
534  // sum up number of local elements
535  size_t numGlobalElements = 0;
536  Teuchos::reduceAll(*(rangeMap->getComm()), Teuchos::REDUCE_SUM, colmapentries.size(), Teuchos::outArg(numGlobalElements));
537 
538  // store global full column map
539  const Teuchos::ArrayView<const GO> aView = Teuchos::ArrayView<const GO>(colmapentries);
540  fullcolmap_ = MapFactory::Build(rangeMap->lib(), numGlobalElements, aView, 0, rangeMap->getComm());
541 #endif
542  }
543 
545 
547 
550  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getGlobalNumRows");
551  global_size_t globalNumRows = 0;
552 
553  for (size_t row = 0; row < Rows(); row++)
554  for (size_t col = 0; col < Cols(); col++)
555  if (!getMatrix(row,col).is_null()) {
556  globalNumRows += getMatrix(row,col)->getGlobalNumRows();
557  break; // we need only one non-null matrix in a row
558  }
559 
560  return globalNumRows;
561  }
562 
564 
567  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getGlobalNumCols");
568  global_size_t globalNumCols = 0;
569 
570  for (size_t col = 0; col < Cols(); col++)
571  for (size_t row = 0; row < Rows(); row++)
572  if (!getMatrix(row,col).is_null()) {
573  globalNumCols += getMatrix(row,col)->getGlobalNumCols();
574  break; // we need only one non-null matrix in a col
575  }
576 
577  return globalNumCols;
578  }
579 
581  size_t getNodeNumRows() const {
582  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getNodeNumRows");
583  global_size_t nodeNumRows = 0;
584 
585  for (size_t row = 0; row < Rows(); ++row)
586  for (size_t col = 0; col < Cols(); col++)
587  if (!getMatrix(row,col).is_null()) {
588  nodeNumRows += getMatrix(row,col)->getNodeNumRows();
589  break; // we need only one non-null matrix in a row
590  }
591 
592  return nodeNumRows;
593  }
594 
597  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getGlobalNumEntries");
598  global_size_t globalNumEntries = 0;
599 
600  for (size_t row = 0; row < Rows(); ++row)
601  for (size_t col = 0; col < Cols(); ++col)
602  if (!getMatrix(row,col).is_null())
603  globalNumEntries += getMatrix(row,col)->getGlobalNumEntries();
604 
605  return globalNumEntries;
606  }
607 
609  size_t getNodeNumEntries() const {
610  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getNodeNumEntries");
611  global_size_t nodeNumEntries = 0;
612 
613  for (size_t row = 0; row < Rows(); ++row)
614  for (size_t col = 0; col < Cols(); ++col)
615  if (!getMatrix(row,col).is_null())
616  nodeNumEntries += getMatrix(row,col)->getNodeNumEntries();
617 
618  return nodeNumEntries;
619  }
620 
622 
623  size_t getNumEntriesInLocalRow(LocalOrdinal localRow) const {
624  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getNumEntriesInLocalRow");
625  if (Rows() == 1 && Cols () == 1) {
626  return getMatrix(0,0)->getNumEntriesInLocalRow(localRow);
627  }
628  else if(is_diagonal_){
629  GlobalOrdinal gid = this->getRowMap()->getGlobalElement(localRow);
630  size_t row = getBlockedRangeMap()->getMapIndexForGID(gid);
631  return getMatrix(row,row)->getNumEntriesInLocalRow(getMatrix(row,row)->getRowMap()->getLocalElement(gid));
632  }
633  throw Xpetra::Exceptions::RuntimeError("getNumEntriesInLocalRow() not supported by BlockedCrsMatrix");
634  }
635 
637 
639  size_t getGlobalMaxNumRowEntries() const {
640  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getGlobalMaxNumRowEntries");
641  global_size_t globalMaxEntries = 0;
642 
643  for (size_t row = 0; row < Rows(); row++) {
644  global_size_t globalMaxEntriesBlockRows = 0;
645  for (size_t col = 0; col < Cols(); col++) {
646  if (!getMatrix(row,col).is_null()) {
647  globalMaxEntriesBlockRows += getMatrix(row,col)->getGlobalMaxNumRowEntries();
648  }
649  }
650  if(globalMaxEntriesBlockRows > globalMaxEntries)
651  globalMaxEntries = globalMaxEntriesBlockRows;
652  }
653  return globalMaxEntries;
654  }
655 
657 
659  size_t getNodeMaxNumRowEntries() const {
660  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getNodeMaxNumRowEntries");
661  size_t localMaxEntries = 0;
662 
663  for (size_t row = 0; row < Rows(); row++) {
664  size_t localMaxEntriesBlockRows = 0;
665  for (size_t col = 0; col < Cols(); col++) {
666  if (!getMatrix(row,col).is_null()) {
667  localMaxEntriesBlockRows += getMatrix(row,col)->getNodeMaxNumRowEntries();
668  }
669  }
670  if(localMaxEntriesBlockRows > localMaxEntries)
671  localMaxEntries = localMaxEntriesBlockRows;
672  }
673  return localMaxEntries;
674  }
675 
677 
680  bool isLocallyIndexed() const {
681  XPETRA_MONITOR("XpetraBlockedCrsMatrix::isLocallyIndexed");
682  for (size_t i = 0; i < blocks_.size(); ++i)
683  if (blocks_[i] != Teuchos::null && !blocks_[i]->isLocallyIndexed())
684  return false;
685  return true;
686  }
687 
689 
692  bool isGloballyIndexed() const {
693  XPETRA_MONITOR("XpetraBlockedCrsMatrix::isGloballyIndexed");
694  for (size_t i = 0; i < blocks_.size(); i++)
695  if (blocks_[i] != Teuchos::null && !blocks_[i]->isGloballyIndexed())
696  return false;
697  return true;
698  }
699 
701  bool isFillComplete() const {
702  XPETRA_MONITOR("XpetraBlockedCrsMatrix::isFillComplete");
703  for (size_t i = 0; i < blocks_.size(); i++)
704  if (blocks_[i] != Teuchos::null && !blocks_[i]->isFillComplete())
705  return false;
706  return true;
707  }
708 
710 
724  virtual void getLocalRowCopy(LocalOrdinal LocalRow,
725  const ArrayView<LocalOrdinal>& Indices,
726  const ArrayView<Scalar>& Values,
727  size_t &NumEntries) const {
728  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getLocalRowCopy");
729  if (Rows() == 1 && Cols () == 1) {
730  getMatrix(0,0)->getLocalRowCopy(LocalRow, Indices, Values, NumEntries);
731  return;
732  }
733  throw Xpetra::Exceptions::RuntimeError("getLocalRowCopy not supported by BlockedCrsMatrix" );
734  }
735 
737 
746  void getGlobalRowView(GlobalOrdinal GlobalRow, ArrayView<const GlobalOrdinal>& indices, ArrayView<const Scalar>& values) const {
747  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getGlobalRowView");
748  if (Rows() == 1 && Cols () == 1) {
749  getMatrix(0,0)->getGlobalRowView(GlobalRow, indices, values);
750  return;
751  }
752  throw Xpetra::Exceptions::RuntimeError("getGlobalRowView not supported by BlockedCrsMatrix");
753  }
754 
756 
765  void getLocalRowView(LocalOrdinal LocalRow, ArrayView<const LocalOrdinal>& indices, ArrayView<const Scalar>& values) const {
766  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getLocalRowView");
767  if (Rows() == 1 && Cols () == 1) {
768  getMatrix(0,0)->getLocalRowView(LocalRow, indices, values);
769  return;
770  }
771  else if(is_diagonal_) {
772  //CMS
773  GlobalOrdinal gid = this->getRowMap()->getGlobalElement(LocalRow);
774  size_t row = getBlockedRangeMap()->getMapIndexForGID(gid);
775  getMatrix(row,row)->getLocalRowView(getMatrix(row,row)->getRowMap()->getLocalElement(gid),indices,values);
776  return;
777  }
778  throw Xpetra::Exceptions::RuntimeError("getLocalRowView not supported by BlockedCrsMatrix");
779  }
780 
782 
785  void getLocalDiagCopy(Vector& diag) const {
786  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getLocalDiagCopy");
787 
788  //RCP<Teuchos::FancyOStream> fancy = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
789 
790  Teuchos::RCP<Vector> rcpdiag = Teuchos::rcpFromRef(diag);
791  Teuchos::RCP<BlockedVector> bdiag = Teuchos::rcp_dynamic_cast<BlockedVector>(rcpdiag);
792 
793  // special treatment for 1x1 block matrices
794  // ReorderedBlockedCrsMatrix object encapsulate single blocks in ReorderedBlocks
795  // BlockedVectors have Vector objects as Leaf objects.
796  if(Rows() == 1 && Cols() == 1 && bdiag.is_null() == true) {
798  rm->getLocalDiagCopy(diag);
799  return;
800  }
801 
802  TEUCHOS_TEST_FOR_EXCEPTION(bdiag.is_null() == true, Xpetra::Exceptions::RuntimeError, "BlockedCrsMatrix::getLocalDiagCopy: diag must be a Blocked(Multi)Vector.");
803  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());
804  TEUCHOS_TEST_FOR_EXCEPTION(bdiag->getBlockedMap()->getNumMaps() != this->Rows(), Xpetra::Exceptions::RuntimeError,
805  "BlockedCrsMatrix::getLocalDiagCopy(): the number of blocks in diag differ from the number of blocks in this operator." );
806  //XPETRA_TEST_FOR_EXCEPTION(bdiag->getMap()->isSameAs(*(getMap())) == false, Xpetra::Exceptions::RuntimeError,
807  // "BlockedCrsMatrix::getLocalDiagCopy(): the map of the vector diag is not compatible with the map of the blocked operator." );
808 
809  for (size_t row = 0; row < Rows(); row++) {
811  if (!rm.is_null()) {
812  Teuchos::RCP<Vector> rv = VectorFactory::Build(bdiag->getBlockedMap()->getMap(row,bdiag->getBlockedMap()->getThyraMode()));
813  rm->getLocalDiagCopy(*rv);
814  bdiag->setMultiVector(row,rv,bdiag->getBlockedMap()->getThyraMode());
815  }
816  }
817  }
818 
820  void leftScale (const Vector& x) {
821  XPETRA_MONITOR("XpetraBlockedCrsMatrix::leftScale");
822 
823  Teuchos::RCP<const Vector> rcpx = Teuchos::rcpFromRef(x);
824  Teuchos::RCP<const BlockedVector> bx = Teuchos::rcp_dynamic_cast<const BlockedVector>(rcpx);
825 
826  //RCP<Teuchos::FancyOStream> fancy = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
827 
828  // special treatment for 1xn block matrices
829  // ReorderedBlockedCrsMatrix object encapsulate single blocks in ReorderedBlocks
830  // BlockedVectors have Vector objects as Leaf objects.
831  if(Rows() == 1 && bx.is_null() == true) {
832  for (size_t col = 0; col < Cols(); ++col) {
833  Teuchos::RCP<Matrix> rm = getMatrix(0,col);
834  rm->leftScale(x);
835  }
836  return;
837  }
838 
839  TEUCHOS_TEST_FOR_EXCEPTION(bx.is_null() == true, Xpetra::Exceptions::RuntimeError, "BlockedCrsMatrix::leftScale: x must be a Blocked(Multi)Vector.");
840  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());
841  TEUCHOS_TEST_FOR_EXCEPTION(bx->getBlockedMap()->getNumMaps() != this->Rows(), Xpetra::Exceptions::RuntimeError,
842  "BlockedCrsMatrix::leftScale(): the number of blocks in diag differ from the number of blocks in this operator." );
843 
844  for (size_t row = 0; row < Rows(); row++) {
845  Teuchos::RCP<const MultiVector> rmv = bx->getMultiVector(row);
846  Teuchos::RCP<const Vector> rscale = rmv->getVector(0);
847  XPETRA_TEST_FOR_EXCEPTION(rscale.is_null()==true, Xpetra::Exceptions::RuntimeError, "BlockedCrsMatrix::leftScale: x must be a Vector.");
848  for (size_t col = 0; col < Cols(); ++col) {
849  Teuchos::RCP<Matrix> rm = getMatrix(row,col);
850  if (!rm.is_null()) {
851  rm->leftScale(*rscale);
852  }
853  }
854  }
855  }
856 
858  void rightScale (const Vector& x) {
859  XPETRA_MONITOR("XpetraBlockedCrsMatrix::rightScale");
860 
861  Teuchos::RCP<const Vector> rcpx = Teuchos::rcpFromRef(x);
862  Teuchos::RCP<const BlockedVector> bx = Teuchos::rcp_dynamic_cast<const BlockedVector>(rcpx);
863 
864  //RCP<Teuchos::FancyOStream> fancy = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
865 
866  // special treatment for nx1 block matrices
867  // ReorderedBlockedCrsMatrix object encapsulate single blocks in ReorderedBlocks
868  // BlockedVectors have Vector objects as Leaf objects.
869  if(Cols() == 1 && bx.is_null() == true) {
870  for (size_t row = 0; row < Rows(); ++row) {
871  Teuchos::RCP<Matrix> rm = getMatrix(row,0);
872  rm->rightScale(x);
873  }
874  return;
875  }
876 
877  TEUCHOS_TEST_FOR_EXCEPTION(bx.is_null() == true, Xpetra::Exceptions::RuntimeError, "BlockedCrsMatrix::rightScale: x must be a Blocked(Multi)Vector.");
878  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());
879  TEUCHOS_TEST_FOR_EXCEPTION(bx->getBlockedMap()->getNumMaps() != this->Cols(), Xpetra::Exceptions::RuntimeError,
880  "BlockedCrsMatrix::rightScale(): the number of blocks in diag differ from the number of blocks in this operator." );
881 
882  for (size_t col = 0; col < Cols(); ++col) {
883  Teuchos::RCP<const MultiVector> rmv = bx->getMultiVector(col);
884  Teuchos::RCP<const Vector> rscale = rmv->getVector(0);
885  XPETRA_TEST_FOR_EXCEPTION(rscale.is_null()==true, Xpetra::Exceptions::RuntimeError, "BlockedCrsMatrix::leftScale: x must be a Vector.");
886  for (size_t row = 0; row < Rows(); row++) {
887  Teuchos::RCP<Matrix> rm = getMatrix(row,col);
888  if (!rm.is_null()) {
889  rm->rightScale(*rscale);
890  }
891  }
892  }
893  }
894 
895 
898  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getFrobeniusNorm");
900  for (size_t col = 0; col < Cols(); ++col) {
901  for (size_t row = 0; row < Rows(); ++row) {
902  if(getMatrix(row,col)!=Teuchos::null) {
903  typename ScalarTraits<Scalar>::magnitudeType n = getMatrix(row,col)->getFrobeniusNorm();
904  ret += n * n;
905  }
906  }
907  }
909  }
910 
911 
913  virtual bool haveGlobalConstants() const {return true;}
914 
915 
917 
919 
920 
922 
942 
944 
945 
947 
950  virtual void apply(const MultiVector& X, MultiVector& Y,
952  Scalar alpha = ScalarTraits<Scalar>::one(),
953  Scalar beta = ScalarTraits<Scalar>::zero()) const
954  {
955  XPETRA_MONITOR("XpetraBlockedCrsMatrix::apply");
956  //using Teuchos::RCP;
957 
959  "apply() only supports the following modes: NO_TRANS and TRANS." );
960 
961  // check whether input parameters are blocked or not
962  RCP<const MultiVector> refX = rcpFromRef(X);
963  RCP<const BlockedMultiVector> refbX = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(refX);
964  //RCP<MultiVector> tmpY = rcpFromRef(Y);
965  //RCP<BlockedMultiVector> tmpbY = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(tmpY);
966 
967  // TODO get rid of me: adapt MapExtractor
968  bool bBlockedX = (refbX != Teuchos::null) ? true : false;
969 
970  // create (temporary) vectors for output
971  // In the end we call Y.update(alpha, *tmpY, beta). Therefore we need a new vector storing the temporary results
973 
974  //RCP<Teuchos::FancyOStream> out = rcp(new Teuchos::FancyOStream(rcp(&std::cout,false)));
975 
976  SC one = ScalarTraits<SC>::one();
977 
978  if (mode == Teuchos::NO_TRANS) {
979 
980  for (size_t row = 0; row < Rows(); row++) {
981  RCP<MultiVector> Yblock = rangemaps_->getVector(row, Y.getNumVectors(), bRangeThyraMode_, true);
982  for (size_t col = 0; col < Cols(); col++) {
983 
984  // extract matrix block
985  RCP<Matrix> Ablock = getMatrix(row, col);
986 
987  if (Ablock.is_null())
988  continue;
989 
990  // check whether Ablock is itself a blocked operator
991  // If it is a blocked operator we have to provide Xpetra style GIDs, i.e. we have to transform GIDs
992  bool bBlockedSubMatrix = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(Ablock) == Teuchos::null ? false : true;
993 
994  // input/output vectors for local block operation
995  RCP<const MultiVector> Xblock = Teuchos::null; // subpart of X vector to be applied to subblock of A
996 
997  // extract sub part of X using Xpetra or Thyra GIDs
998  // if submatrix is again blocked, we extract it using Xpetra style gids. If it is a single
999  // block matrix we use the Thyra or Xpetra style GIDs that are used to store the matrix
1000  if(bBlockedX) Xblock = domainmaps_->ExtractVector(refbX, col, bDomainThyraMode_);
1001  else Xblock = domainmaps_->ExtractVector(refX, col, bBlockedSubMatrix == true ? false : bDomainThyraMode_);
1002 
1003  RCP<MultiVector> tmpYblock = rangemaps_->getVector(row, Y.getNumVectors(), bRangeThyraMode_, false); // subpart of Y vector containing part of solution of Xblock applied to Ablock
1004  Ablock->apply(*Xblock, *tmpYblock);
1005 
1006  Yblock->update(one, *tmpYblock, one);
1007  }
1008  rangemaps_->InsertVector(Yblock, row, tmpY, bRangeThyraMode_);
1009  }
1010 
1011  } else if (mode == Teuchos::TRANS) {
1012  // TODO: test me!
1013  for (size_t col = 0; col < Cols(); col++) {
1014  RCP<MultiVector> Yblock = domainmaps_->getVector(col, Y.getNumVectors(), bDomainThyraMode_, true);
1015 
1016  for (size_t row = 0; row<Rows(); row++) {
1017  RCP<Matrix> Ablock = getMatrix(row, col);
1018 
1019  if (Ablock.is_null())
1020  continue;
1021 
1022  // check whether Ablock is itself a blocked operator
1023  // If it is a blocked operator we have to provide Xpetra style GIDs, i.e. we have to transform GIDs
1024  bool bBlockedSubMatrix = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(Ablock) == Teuchos::null ? false : true;
1025 
1026  RCP<const MultiVector> Xblock = Teuchos::null;
1027 
1028  // extract sub part of X using Xpetra or Thyra GIDs
1029  if(bBlockedX) Xblock = rangemaps_->ExtractVector(refbX, row, bRangeThyraMode_);
1030  else Xblock = rangemaps_->ExtractVector(refX, row, bBlockedSubMatrix == true ? false : bRangeThyraMode_);
1031  RCP<MultiVector> tmpYblock = domainmaps_->getVector(col, Y.getNumVectors(), bDomainThyraMode_, false);
1032  Ablock->apply(*Xblock, *tmpYblock, Teuchos::TRANS);
1033 
1034  Yblock->update(one, *tmpYblock, one);
1035  }
1036  domainmaps_->InsertVector(Yblock, col, tmpY, bDomainThyraMode_);
1037  }
1038  }
1039  Y.update(alpha, *tmpY, beta);
1040  }
1041 
1043  RCP<const Map > getFullDomainMap() const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getFullDomainMap()"); return domainmaps_->getFullMap(); }
1044 
1046  RCP<const BlockedMap > getBlockedDomainMap() const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getBlockedDomainMap()"); return domainmaps_->getBlockedMap(); }
1047 
1049  RCP<const Map > getDomainMap() const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getDomainMap()"); return domainmaps_->getMap(); /*domainmaps_->getFullMap();*/ }
1050 
1052  RCP<const Map > getDomainMap(size_t i) const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getDomainMap(size_t)"); return domainmaps_->getMap(i, bDomainThyraMode_); }
1053 
1055  RCP<const Map > getDomainMap(size_t i, bool bThyraMode) const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getDomainMap(size_t,bool)"); return domainmaps_->getMap(i, bThyraMode); }
1056 
1058  RCP<const Map > getFullRangeMap() const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getRangeMap()"); return rangemaps_->getFullMap(); }
1059 
1061  RCP<const BlockedMap > getBlockedRangeMap() const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getBlockedRangeMap()"); return rangemaps_->getBlockedMap(); }
1062 
1064  RCP<const Map > getRangeMap() const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getRangeMap()"); return rangemaps_->getMap(); /*rangemaps_->getFullMap();*/ }
1065 
1067  RCP<const Map > getRangeMap(size_t i) const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getRangeMap(size_t)"); return rangemaps_->getMap(i, bRangeThyraMode_); }
1068 
1070  RCP<const Map > getRangeMap(size_t i, bool bThyraMode) const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getRangeMap(size_t,bool)"); return rangemaps_->getMap(i, bThyraMode); }
1071 
1073  RCP<const MapExtractor> getRangeMapExtractor() const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getRangeMapExtractor()"); return rangemaps_; }
1074 
1076  RCP<const MapExtractor> getDomainMapExtractor() const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::getDomainMapExtractor()"); return domainmaps_; }
1077 
1079 
1081  //{@
1082 
1091  virtual void bgs_apply(
1092  const MultiVector& X,
1093  MultiVector& Y,
1094  size_t row,
1096  Scalar alpha = ScalarTraits<Scalar>::one(),
1097  Scalar beta = ScalarTraits<Scalar>::zero()
1098  ) const
1099  {
1100  XPETRA_MONITOR("XpetraBlockedCrsMatrix::bgs_apply");
1101  //using Teuchos::RCP;
1102 
1104  "apply() only supports the following modes: NO_TRANS and TRANS." );
1105 
1106  // check whether input parameters are blocked or not
1107  RCP<const MultiVector> refX = rcpFromRef(X);
1108  RCP<const BlockedMultiVector> refbX = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(refX);
1109  //RCP<MultiVector> tmpY = rcpFromRef(Y);
1110  //RCP<BlockedMultiVector> tmpbY = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(tmpY);
1111 
1112  bool bBlockedX = (refbX != Teuchos::null) ? true : false;
1113 
1114  // create (temporary) vectors for output
1115  // In the end we call Y.update(alpha, *tmpY, beta). Therefore we need a new vector storing the temporary results
1117 
1118  SC one = ScalarTraits<SC>::one();
1119 
1120  if (mode == Teuchos::NO_TRANS) {
1121  RCP<MultiVector> Yblock = rangemaps_->getVector(row, Y.getNumVectors(), bRangeThyraMode_, true);
1122  for (size_t col = 0; col < Cols(); col++) {
1123 
1124  // extract matrix block
1125  RCP<Matrix> Ablock = getMatrix(row, col);
1126 
1127  if (Ablock.is_null())
1128  continue;
1129 
1130  // check whether Ablock is itself a blocked operator
1131  // If it is a blocked operator we have to provide Xpetra style GIDs, i.e. we have to transform GIDs
1132  bool bBlockedSubMatrix = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(Ablock) == Teuchos::null ? false : true;
1133 
1134  // input/output vectors for local block operation
1135  RCP<const MultiVector> Xblock = Teuchos::null; // subpart of X vector to be applied to subblock of A
1136 
1137  // extract sub part of X using Xpetra or Thyra GIDs
1138  // if submatrix is again blocked, we extract it using Xpetra style gids. If it is a single
1139  // block matrix we use the Thyra or Xpetra style GIDs that are used to store the matrix
1140  if(bBlockedX) Xblock = domainmaps_->ExtractVector(refbX, col, bDomainThyraMode_);
1141  else Xblock = domainmaps_->ExtractVector(refX, col, bBlockedSubMatrix == true ? false : bDomainThyraMode_);
1142 
1143  RCP<MultiVector> tmpYblock = rangemaps_->getVector(row, Y.getNumVectors(), bRangeThyraMode_, false); // subpart of Y vector containing part of solution of Xblock applied to Ablock
1144  Ablock->apply(*Xblock, *tmpYblock);
1145 
1146  Yblock->update(one, *tmpYblock, one);
1147  }
1148  rangemaps_->InsertVector(Yblock, row, tmpY, bRangeThyraMode_);
1149  } else {
1150  TEUCHOS_TEST_FOR_EXCEPTION(true,Xpetra::Exceptions::NotImplemented,"Xpetar::BlockedCrsMatrix::bgs_apply: not implemented for transpose case.");
1151  }
1152  Y.update(alpha, *tmpY, beta);
1153  }
1154 
1155 
1157 
1158 
1160  //{@
1161 
1164  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getMap");
1165  if (Rows() == 1 && Cols () == 1) {
1166  return getMatrix(0,0)->getMap();
1167  }
1168  throw Xpetra::Exceptions::RuntimeError("BlockedCrsMatrix::getMap(): operation not supported.");
1169  }
1170 
1172  void doImport(const Matrix &source, const Import& importer, CombineMode CM) {
1173  XPETRA_MONITOR("XpetraBlockedCrsMatrix::doImport");
1174  if (Rows() == 1 && Cols () == 1) {
1175  getMatrix(0,0)->doImport(source, importer, CM);
1176  return;
1177  }
1178  throw Xpetra::Exceptions::RuntimeError("BlockedCrsMatrix::doImport(): operation not supported.");
1179  }
1180 
1182  void doExport(const Matrix& dest, const Import& importer, CombineMode CM) {
1183  XPETRA_MONITOR("XpetraBlockedCrsMatrix::doExport");
1184  if (Rows() == 1 && Cols () == 1) {
1185  getMatrix(0,0)->doExport(dest, importer, CM);
1186  return;
1187  }
1188  throw Xpetra::Exceptions::RuntimeError("BlockedCrsMatrix::doExport(): operation not supported.");
1189  }
1190 
1192  void doImport(const Matrix& source, const Export& exporter, CombineMode CM) {
1193  XPETRA_MONITOR("XpetraBlockedCrsMatrix::doImport");
1194  if (Rows() == 1 && Cols () == 1) {
1195  getMatrix(0,0)->doImport(source, exporter, CM);
1196  return;
1197  }
1198  throw Xpetra::Exceptions::RuntimeError("BlockedCrsMatrix::doImport(): operation not supported.");
1199  }
1200 
1202  void doExport(const Matrix& dest, const Export& exporter, CombineMode CM) {
1203  XPETRA_MONITOR("XpetraBlockedCrsMatrix::doExport");
1204  if (Rows() == 1 && Cols () == 1) {
1205  getMatrix(0,0)->doExport(dest, exporter, CM);
1206  return;
1207  }
1208  throw Xpetra::Exceptions::RuntimeError("BlockedCrsMatrix::doExport(): operation not supported.");
1209  }
1210 
1211  // @}
1212 
1214 
1215 
1217  std::string description() const { return "Xpetra_BlockedCrsMatrix.description()"; }
1218 
1221  out << "Xpetra::BlockedCrsMatrix: " << Rows() << " x " << Cols() << std::endl;
1222 
1223  if (isFillComplete()) {
1224  out << "BlockMatrix is fillComplete" << std::endl;
1225 
1226  /*if(fullrowmap_ != Teuchos::null) {
1227  out << "fullRowMap" << std::endl;
1228  fullrowmap_->describe(out,verbLevel);
1229  } else {
1230  out << "fullRowMap not set. Check whether block matrix is properly fillCompleted!" << std::endl;
1231  }*/
1232 
1233  //out << "fullColMap" << std::endl;
1234  //fullcolmap_->describe(out,verbLevel);
1235 
1236  } else {
1237  out << "BlockMatrix is NOT fillComplete" << std::endl;
1238  }
1239 
1240  for (size_t r = 0; r < Rows(); ++r)
1241  for (size_t c = 0; c < Cols(); ++c) {
1242  if(getMatrix(r,c)!=Teuchos::null) {
1243  out << "Block(" << r << "," << c << ")" << std::endl;
1244  getMatrix(r,c)->describe(out,verbLevel);
1245  } else out << "Block(" << r << "," << c << ") = null" << std::endl;
1246  }
1247  }
1248 
1250 
1251  void setObjectLabel( const std::string &objectLabel ) {
1252  XPETRA_MONITOR("TpetraBlockedCrsMatrix::setObjectLabel");
1253  for (size_t r = 0; r < Rows(); ++r)
1254  for (size_t c = 0; c < Cols(); ++c) {
1255  if(getMatrix(r,c)!=Teuchos::null) {
1256  std::ostringstream oss; oss<< objectLabel << "(" << r << "," << c << ")";
1257  getMatrix(r,c)->setObjectLabel(oss.str());
1258  }
1259  }
1260  }
1262 
1263 
1265  bool hasCrsGraph() const {
1266  if (Rows() == 1 && Cols () == 1) return true;
1267  else return false;
1268  }
1269 
1272  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getCrsGraph");
1273  if (Rows() == 1 && Cols () == 1) {
1274  return getMatrix(0,0)->getCrsGraph();
1275  }
1276  throw Xpetra::Exceptions::RuntimeError("getCrsGraph() not supported by BlockedCrsMatrix");
1277  }
1278 
1280 
1282 
1283 
1284  virtual bool isDiagonal() const {return is_diagonal_;}
1285 
1287  virtual size_t Rows() const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::Rows"); return rangemaps_->NumMaps(); }
1288 
1290  virtual size_t Cols() const { XPETRA_MONITOR("XpetraBlockedCrsMatrix::Cols"); return domainmaps_->NumMaps(); }
1291 
1294  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getCrsMatrix");
1295  TEUCHOS_TEST_FOR_EXCEPTION(Rows()!=1, std::out_of_range, "Can only unwrap a 1x1 blocked matrix. The matrix has " << Rows() << " block rows, though.");
1296  TEUCHOS_TEST_FOR_EXCEPTION(Cols()!=1, std::out_of_range, "Can only unwrap a 1x1 blocked matrix. The matrix has " << Cols() << " block columns, though.");
1297 
1298  RCP<Matrix> mat = getMatrix(0,0);
1299  RCP<BlockedCrsMatrix> bmat = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(mat);
1300  if (bmat == Teuchos::null) return mat;
1301  return bmat->getCrsMatrix();
1302  }
1303 
1306  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getInnermostCrsMatrix");
1307  size_t row = Rows()+1, col = Cols()+1;
1308  for (size_t r = 0; r < Rows(); ++r)
1309  for(size_t c = 0; c < Cols(); ++c)
1310  if (getMatrix(r,c) != Teuchos::null) {
1311  row = r;
1312  col = c;
1313  break;
1314  }
1315  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.")
1316  RCP<Matrix> mm = getMatrix(row,col);
1317  RCP<BlockedCrsMatrix> bmat = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(mm);
1318  if (bmat == Teuchos::null) return mm;
1319  return bmat->getInnermostCrsMatrix();
1320  }
1321 
1323  Teuchos::RCP<Matrix> getMatrix(size_t r, size_t c) const {
1324  XPETRA_MONITOR("XpetraBlockedCrsMatrix::getMatrix");
1325  TEUCHOS_TEST_FOR_EXCEPTION(r > Rows(), std::out_of_range, "Error, r = " << Rows() << " is too big");
1326  TEUCHOS_TEST_FOR_EXCEPTION(c > Cols(), std::out_of_range, "Error, c = " << Cols() << " is too big");
1327 
1328  // transfer strided/blocked map information
1329  /* if (blocks_[r*Cols()+c] != Teuchos::null &&
1330  blocks_[r*Cols()+c]->IsView("stridedMaps") == false)
1331  blocks_[r*Cols()+c]->CreateView("stridedMaps", getRangeMap(r,bRangeThyraMode_), getDomainMap(c,bDomainThyraMode_));*/
1332  return blocks_[r*Cols()+c];
1333  }
1334 
1336  //void setMatrix(size_t r, size_t c, Teuchos::RCP<CrsMatrix> mat) {
1337  void setMatrix(size_t r, size_t c, Teuchos::RCP<Matrix> mat) {
1338  XPETRA_MONITOR("XpetraBlockedCrsMatrix::setMatrix");
1339  // TODO: if filled -> return error
1340 
1341  TEUCHOS_TEST_FOR_EXCEPTION(r > Rows(), std::out_of_range, "Error, r = " << Rows() << " is too big");
1342  TEUCHOS_TEST_FOR_EXCEPTION(c > Cols(), std::out_of_range, "Error, c = " << Cols() << " is too big");
1343  if(!mat.is_null() && r!=c) is_diagonal_ = false;
1344  // set matrix
1345  blocks_[r*Cols() + c] = mat;
1346  }
1347 
1349  // NOTE: This is a rather expensive operation, since all blocks are copied
1350  // into a new big CrsMatrix
1352  XPETRA_MONITOR("XpetraBlockedCrsMatrix::Merge");
1353  using Teuchos::RCP;
1354  using Teuchos::rcp_dynamic_cast;
1355  Scalar one = ScalarTraits<SC>::one();
1356 
1358  "BlockedCrsMatrix::Merge: only implemented for Xpetra-style or Thyra-style numbering. No mixup allowed!" );
1359 
1361  "BlockedCrsMatrix::Merge: BlockMatrix must be fill-completed." );
1362 
1363  RCP<Matrix> sparse = MatrixFactory::Build(getRangeMapExtractor()->getFullMap(), 33);
1364 
1365  if(bRangeThyraMode_ == false) {
1366  // Xpetra mode
1367  for (size_t i = 0; i < Rows(); i++) {
1368  for (size_t j = 0; j < Cols(); j++) {
1369  if (getMatrix(i,j) != Teuchos::null) {
1370  RCP<const Matrix> mat = getMatrix(i,j);
1371 
1372  // recursively call Merge routine
1373  RCP<const BlockedCrsMatrix> bMat = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(mat);
1374  if (bMat != Teuchos::null) mat = bMat->Merge();
1375 
1376  bMat = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(mat);
1378  "BlockedCrsMatrix::Merge: Merging of blocked sub-operators failed?!" );
1379 
1380  // jump over empty blocks
1381  if(mat->getNodeNumEntries() == 0) continue;
1382 
1383  this->Add(*mat, one, *sparse, one);
1384  }
1385  }
1386  }
1387  } else {
1388  // Thyra mode
1389  for (size_t i = 0; i < Rows(); i++) {
1390  for (size_t j = 0; j < Cols(); j++) {
1391  if (getMatrix(i,j) != Teuchos::null) {
1392  RCP<const Matrix> mat = getMatrix(i,j);
1393  // recursively call Merge routine
1394  RCP<const BlockedCrsMatrix> bMat = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(mat);
1395  if (bMat != Teuchos::null) mat = bMat->Merge();
1396 
1397  bMat = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(mat);
1399  "BlockedCrsMatrix::Merge: Merging of blocked sub-operators failed?!" );
1400 
1401  // check whether we have a CrsMatrix block (no blocked operator)
1402  RCP<const CrsMatrixWrap> crsMat = Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(mat);
1403  TEUCHOS_ASSERT(crsMat != Teuchos::null);
1404 
1405  // these are the thyra style maps of the matrix
1406  RCP<const Map> trowMap = mat->getRowMap();
1407  RCP<const Map> tcolMap = mat->getColMap();
1408  RCP<const Map> tdomMap = mat->getDomainMap();
1409 
1410  // get Xpetra maps
1411  RCP<const Map> xrowMap = getRangeMapExtractor()->getMap(i,false);
1412  RCP<const Map> xdomMap = getDomainMapExtractor()->getMap(j,false);
1413 
1414  // generate column map with Xpetra GIDs
1415  // We have to do this separately for each block since the column
1416  // map of each block might be different in the same block column
1418  *tcolMap,
1419  *tdomMap,
1420  *xdomMap);
1421 
1422  // jump over empty blocks
1423  if(mat->getNodeNumEntries() == 0) continue;
1424 
1425  size_t maxNumEntries = mat->getNodeMaxNumRowEntries();
1426 
1427  size_t numEntries;
1428  Array<GO> inds (maxNumEntries);
1429  Array<GO> inds2(maxNumEntries);
1430  Array<SC> vals (maxNumEntries);
1431 
1432  // loop over all rows and add entries
1433  for (size_t k = 0; k < mat->getNodeNumRows(); k++) {
1434  GlobalOrdinal rowTGID = trowMap->getGlobalElement(k);
1435  crsMat->getCrsMatrix()->getGlobalRowCopy(rowTGID, inds(), vals(), numEntries);
1436 
1437  // create new indices array
1438  for (size_t l = 0; l < numEntries; ++l) {
1439  LocalOrdinal lid = tcolMap->getLocalElement(inds[l]);
1440  inds2[l] = xcolMap->getGlobalElement(lid);
1441  }
1442 
1443  GlobalOrdinal rowXGID = xrowMap->getGlobalElement(k);
1444  sparse->insertGlobalValues(
1445  rowXGID, inds2(0, numEntries),
1446  vals(0, numEntries));
1447  }
1448  }
1449  }
1450  }
1451  }
1452 
1453  sparse->fillComplete(getDomainMap(), getRangeMap());
1454 
1456  "BlockedCrsMatrix::Merge: Local number of entries of merged matrix does not coincide with local number of entries of blocked operator." );
1457 
1459  "BlockedCrsMatrix::Merge: Global number of entries of merged matrix does not coincide with global number of entries of blocked operator." );
1460 
1461  return sparse;
1462  }
1464 
1465 #ifdef HAVE_XPETRA_KOKKOS_REFACTOR
1466  typedef typename CrsMatrix::local_matrix_type local_matrix_type;
1467 
1469  local_matrix_type getLocalMatrix () const {
1470  if (Rows() == 1 && Cols () == 1) {
1471  return getMatrix(0,0)->getLocalMatrix();
1472  }
1473  throw Xpetra::Exceptions::RuntimeError("BlockedCrsMatrix::getLocalMatrix(): operation not supported.");
1474  }
1475 #endif
1476 
1477 #ifdef HAVE_XPETRA_THYRA
1480  Xpetra::ThyraUtils<Scalar,LO,GO,Node>::toThyra(Teuchos::rcpFromRef(*this));
1482 
1484  Teuchos::rcp_dynamic_cast<Thyra::BlockedLinearOpBase<Scalar> >(thOp);
1486  return thbOp;
1487  }
1488 #endif
1489 
1490  private:
1491 
1494 
1496 
1506  void Add(const Matrix& A, const Scalar scalarA, Matrix& B, const Scalar scalarB) const {
1507  XPETRA_MONITOR("XpetraBlockedCrsMatrix::Add");
1509  "Matrix A is not completed");
1510  using Teuchos::Array;
1511  using Teuchos::ArrayView;
1512 
1513  B.scale(scalarB);
1514 
1515  Scalar one = ScalarTraits<SC>::one();
1516  Scalar zero = ScalarTraits<SC>::zero();
1517 
1518  if (scalarA == zero)
1519  return;
1520 
1521  Teuchos::RCP<const Matrix> rcpA = Teuchos::rcpFromRef(A);
1522  Teuchos::RCP<const CrsMatrixWrap> rcpAwrap = Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(rcpA);
1523  TEUCHOS_TEST_FOR_EXCEPTION(rcpAwrap == Teuchos::null, Xpetra::Exceptions::BadCast,
1524  "BlockedCrsMatrix::Add: matrix A must be of type CrsMatrixWrap.");
1525  Teuchos::RCP<const CrsMatrix> crsA = rcpAwrap->getCrsMatrix();
1526 
1527  size_t maxNumEntries = crsA->getNodeMaxNumRowEntries();
1528 
1529  size_t numEntries;
1530  Array<GO> inds(maxNumEntries);
1531  Array<SC> vals(maxNumEntries);
1532 
1533  RCP<const Map> rowMap = crsA->getRowMap();
1534  RCP<const Map> colMap = crsA->getColMap();
1535 
1536  ArrayView<const GO> rowGIDs = crsA->getRowMap()->getNodeElementList();
1537  for (size_t i = 0; i < crsA->getNodeNumRows(); i++) {
1538  GO row = rowGIDs[i];
1539  crsA->getGlobalRowCopy(row, inds(), vals(), numEntries);
1540 
1541  if (scalarA != one)
1542  for (size_t j = 0; j < numEntries; ++j)
1543  vals[j] *= scalarA;
1544 
1545  B.insertGlobalValues(row, inds(0, numEntries), vals(0, numEntries)); // insert should be ok, since blocks in BlockedCrsOpeartor do not overlap!
1546  }
1547  }
1548 
1550 
1551  // Default view is created after fillComplete()
1552  // Because ColMap might not be available before fillComplete().
1554  XPETRA_MONITOR("XpetraBlockedCrsMatrix::CreateDefaultView");
1555 
1556  // Create default view
1557  this->defaultViewLabel_ = "point";
1559 
1560  // Set current view
1561  this->currentViewLabel_ = this->GetDefaultViewLabel();
1562  }
1563 
1564 
1565  private:
1566  bool is_diagonal_; // If we're diagonal a bunch of the extraction stuff should work
1567  Teuchos::RCP<const MapExtractor> domainmaps_; // full domain map together with all partial domain maps
1568  Teuchos::RCP<const MapExtractor> rangemaps_; // full range map together with all partial domain maps
1569 
1570  std::vector<Teuchos::RCP<Matrix> > blocks_; // row major matrix block storage
1571 #ifdef HAVE_XPETRA_THYRA
1573 #endif
1576 
1577 };
1578 
1579 } //namespace Xpetra
1580 
1581 #define XPETRA_BLOCKEDCRSMATRIX_SHORT
1582 #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.
virtual size_t getNodeMaxNumRowEntries() const =0
Returns the maximum number of entries across all rows/columns on this node.
virtual Teuchos::RCP< const Teuchos::Comm< int > > getComm() const =0
Get this Map&#39;s Comm object.
RCP< const Map > getFullDomainMap() const
Returns the Map associated with the full domain of this operator.
virtual const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getColMap() const =0
Returns the Map that describes the column distribution in this matrix.
bool isFillComplete() const
Returns true if fillComplete() has been called and the matrix is in compute mode. ...
void replaceGlobalValues(GlobalOrdinal globalRow, const ArrayView< const GlobalOrdinal > &cols, const ArrayView< const Scalar > &vals)
Replace matrix entries, using global IDs.
void fillComplete(const RCP< const Map > &domainMap, const RCP< const Map > &rangeMap, const RCP< ParameterList > &params=null)
Signal that data entry is complete.
virtual void setAllToScalar(const Scalar &alpha)
Set all matrix entries equal to scalar.
GlobalOrdinal GO
iterator begin() const
bool is_null(const std::shared_ptr< T > &p)
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.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Teuchos::RCP< const MapExtractor > rangemaps_
virtual const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getRowMap() const =0
Returns the Map that describes the row distribution in this matrix.
void doImport(const Matrix &source, const Import &importer, CombineMode CM)
Import.
void getLocalRowView(LocalOrdinal LocalRow, ArrayView< const LocalOrdinal > &indices, ArrayView< const Scalar > &values) const
Extract a const, non-persisting view of local indices in a specified row of the matrix.
size_type size() const
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 LocalOrdinal getLocalElement(GlobalOrdinal globalIndex) const =0
The local index corresponding to the given global index.
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)
BlockedCrsMatrix(const Teuchos::RCP< const BlockedMap > &rangeMaps, const Teuchos::RCP< const BlockedMap > &domainMaps, size_t npr, Xpetra::ProfileType pftype=Xpetra::DynamicProfile)
Constructor.
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.
virtual void getGlobalRowCopy(GlobalOrdinal GlobalRow, const ArrayView< GlobalOrdinal > &indices, const ArrayView< Scalar > &values, size_t &numEntries) const =0
Extract a list of entries in a specified global row of this matrix. Put into pre-allocated storage...
global_size_t getGlobalNumEntries() const
Returns the global number of entries in this matrix.
size_t getNumEntriesInLocalRow(LocalOrdinal localRow) const
Returns the current number of entries on this node in the specified local row.
RCP< const CrsGraph > getCrsGraph() const
Returns the CrsGraph associated with this matrix.
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.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void doExport(const Matrix &dest, const Export &exporter, CombineMode CM)
Export (using an Importer).
void leftScale(const Vector &x)
Left scale matrix using the given vector entries.
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)
virtual Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getMap() const =0
The Map describing the parallel distribution of this object.
viewLabel_t currentViewLabel_
TEUCHOS_DEPRECATED void reduceAll(const Comm< Ordinal > &comm, const EReductionType reductType, const Packet &send, Packet *globalReduct)
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.
virtual GlobalOrdinal getGlobalElement(LocalOrdinal localIndex) const =0
The global index corresponding to the given local index.
Exception throws when you call an unimplemented method of Xpetra.
T * getRawPtr() const
bool hasCrsGraph() const
Supports the getCrsGraph() call.
std::string viewLabel_t
size_t global_size_t
Global size_t object.
static RCP< Vector > Build(const Teuchos::RCP< const Map > &map, bool zeroOut=true)
Constructor specifying the number of non-zeros for all rows.
size_t getNodeMaxNumRowEntries() const
Returns the maximum number of entries across all rows/columns on this node.
static const EVerbosityLevel verbLevel_default
std::vector< Teuchos::RCP< Matrix > > blocks_
Teuchos::RCP< Matrix > Merge() const
merge BlockedCrsMatrix blocks in a CrsMatrix
static magnitudeType magnitude(T a)
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)
void setMatrix(size_t r, size_t c, Teuchos::RCP< Matrix > mat)
set matrix block
iterator end() const
std::string description() const
Return a simple one-line description of this object.
bool bDomainThyraMode_
boolean flag, which is true, if BlockedCrsMatrix has been created using Thyra-style numbering for sub...
BlockedCrsMatrix(Teuchos::RCP< const MapExtractor > &rangeMaps, Teuchos::RCP< const MapExtractor > &domainMaps, size_t npr, Xpetra::ProfileType pftype=Xpetra::DynamicProfile)
Constructor.
bool isGloballyIndexed() const
If matrix indices are in the global range, this function returns true. Otherwise, this function retur...
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
static Teuchos::RCP< MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Build(const Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &map, size_t NumVectors, bool zeroOut=true)
Constructor specifying the number of non-zeros for all rows.
CombineMode
Xpetra::Combine Mode enumerable type.
#define XPETRA_MONITOR(funcName)
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 const RCP< const Map > & getRowMap() const
Returns the Map that describes the row distribution in 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)
Map constructor with Xpetra-defined contiguous uniform distribution.
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.
#define TEUCHOS_ASSERT(assertion_test)
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.
virtual size_t getNodeNumRows() const =0
Returns the number of matrix rows owned on the calling node.
size_t getNodeNumRows() const
Returns the number of matrix rows owned on the calling node.
virtual UnderlyingLib lib() const =0
Get the library used by this object (Tpetra or Epetra?)
virtual ScalarTraits< Scalar >::magnitudeType getFrobeniusNorm() const
Get Frobenius norm of the matrix.
void removeEmptyProcessesInPlace(const Teuchos::RCP< const Map > &newMap)
virtual Teuchos::RCP< const Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getVector(size_t j) const =0
Return a Vector which is a const view of column j.
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
bool isLocallyIndexed() const
If matrix indices of all matrix blocks are in the local range, this function returns true...
void replaceLocalValues(LocalOrdinal localRow, const ArrayView< const LocalOrdinal > &cols, const ArrayView< const Scalar > &vals)
Replace matrix entries, using local IDs.
Xpetra-specific matrix class.
Teuchos::RCP< const MapExtractor > domainmaps_
virtual void apply(const MultiVector &X, MultiVector &Y, Teuchos::ETransp mode=Teuchos::NO_TRANS, Scalar alpha=ScalarTraits< Scalar >::one(), Scalar beta=ScalarTraits< Scalar >::zero()) const
Computes the sparse matrix-multivector multiplication.
bool bRangeThyraMode_
boolean flag, which is true, if BlockedCrsMatrix has been created using Thyra-style numbering for sub...
RCP< const MapExtractor > getRangeMapExtractor() const
Returns map extractor class for range map.
bool is_null() const