Xpetra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Xpetra_ReorderedBlockedCrsMatrix.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_REORDEREDBLOCKEDCRSMATRIX_HPP
47 #define XPETRA_REORDEREDBLOCKEDCRSMATRIX_HPP
48 
49 #include <Kokkos_DefaultNode.hpp>
50 
51 #include "Xpetra_ConfigDefs.hpp"
52 #include "Xpetra_Exceptions.hpp"
53 
54 #include "Xpetra_MapUtils.hpp"
55 
56 #include "Xpetra_BlockedMultiVector.hpp"
58 #include "Xpetra_CrsMatrixWrap.hpp"
60 
61 
66 namespace Xpetra {
67 
68  typedef std::string viewLabel_t;
69 
70  template <class Scalar,
71  class LocalOrdinal,
72  class GlobalOrdinal,
75  public BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> {
76  public:
77  typedef Scalar scalar_type;
78  typedef LocalOrdinal local_ordinal_type;
79  typedef GlobalOrdinal global_ordinal_type;
80  typedef Node node_type;
81 
82  private:
83 #undef XPETRA_REORDEREDBLOCKEDCRSMATRIX_SHORT
84 #include "Xpetra_UseShortNames.hpp"
85 
86  public:
87 
89 
90 
92 
100  (Teuchos::RCP<const MapExtractor>& rangeMaps,
101  Teuchos::RCP<const MapExtractor>& domainMaps,
102  size_t npr,
103  Teuchos::RCP<const Xpetra::BlockReorderManager> brm,
106  brm_ = brm;
107  fullOp_ = bmat;
108  }
109 
110  //protected:
111 
114 
116 
117  private:
118  Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > mergeSubBlockMaps(Teuchos::RCP<const Xpetra::BlockReorderManager> brm) {
119  RCP<const MapExtractor> fullRangeMapExtractor = fullOp_->getRangeMapExtractor();
120 
121  // number of sub blocks
122  size_t numBlocks = brm->GetNumBlocks();
123 
124  Teuchos::RCP<const Map> map = Teuchos::null;
125 
126  if(numBlocks == 0) {
127  // it is a leaf node
128  Teuchos::RCP<const Xpetra::BlockReorderLeaf> leaf = Teuchos::rcp_dynamic_cast<const Xpetra::BlockReorderLeaf>(brm);
129 
130  // never extract Thyra style maps (since we have to merge them)
131  map = fullRangeMapExtractor->getMap(Teuchos::as<size_t>(leaf->GetIndex()), false);
132  } else {
133  // initialize vector for sub maps
134  std::vector<Teuchos::RCP<const Map> > subMaps (numBlocks, Teuchos::null);
135 
136  for(size_t i = 0; i < numBlocks; i++) {
137  Teuchos::RCP<const Xpetra::BlockReorderManager> blkMgr = brm->GetBlock(Teuchos::as<int>(i));
138  subMaps[i] = mergeSubBlockMaps(blkMgr);
139  TEUCHOS_ASSERT(subMaps[i].is_null()==false);
140  }
141 
142  map = MapUtils::concatenateMaps(subMaps);
143  }
144  TEUCHOS_ASSERT(map.is_null()==false);
145  return map;
146  }
147 
148  public:
150 
151 
153 
156  virtual void apply(const MultiVector& X, MultiVector& Y,
157  Teuchos::ETransp mode = Teuchos::NO_TRANS,
158  Scalar alpha = ScalarTraits<Scalar>::one(),
159  Scalar beta = ScalarTraits<Scalar>::zero()) const
160  {
161 
162  // Nested sub-operators should just use the provided X and B vectors
163  if(fullOp_->getNodeNumRows() != this->getNodeNumRows()) {
165  return;
166  }
167 
168  // Special handling for the top level of the nested operator
169 
170  // check whether input parameters are blocked or not
171  RCP<const MultiVector> refX = rcpFromRef(X);
172  RCP<const BlockedMultiVector> refbX = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(refX);
173  RCP<MultiVector> tmpY = rcpFromRef(Y);
174  RCP<BlockedMultiVector> tmpbY = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(tmpY);
175 
176 
177  bool bCopyResultX = false;
178  bool bCopyResultY = false;
179 
180  // TODO create a nested ReorderedBlockedMultiVector out of a nested BlockedMultiVector!
181  // probably not necessary, is it?
182  // check whether X and B are blocked but not ReorderedBlocked operators
183  /*if (refbX != Teuchos::null && fullOp_->getNodeNumRows() == this->getNodeNumRows()) {
184  RCP<const ReorderedBlockedMultiVector> rbCheck = Teuchos::rcp_dynamic_cast<const ReorderedBlockedMultiVector>(refbX);
185  if(rbCheck == Teuchos::null) {
186  RCP<const BlockedMultiVector> bX =
187  Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(Xpetra::buildReorderedBlockedMultiVector(brm_, refbX));
188  TEUCHOS_ASSERT(bX.is_null()==false);
189  refbX.swap(bX);
190  }
191  }
192  if (tmpbY != Teuchos::null && fullOp_->getNodeNumRows() == this->getNodeNumRows()) {
193  RCP<ReorderedBlockedMultiVector> rbCheck = Teuchos::rcp_dynamic_cast<ReorderedBlockedMultiVector>(tmpbY);
194  if(rbCheck == Teuchos::null) {
195  RCP<BlockedMultiVector> bY =
196  Teuchos::rcp_dynamic_cast<BlockedMultiVector>(Xpetra::buildReorderedBlockedMultiVector(brm_, tmpbY));
197  TEUCHOS_ASSERT(bY.is_null()==false);
198  tmpbY.swap(bY);
199  }
200  }*/
201 
202  // if X and B are not blocked, create a blocked version here and use the blocked vectors
203  // for the internal (nested) apply call.
204 
205  // Check whether "this" operator is the reordered variant of the underlying fullOp_.
206  // Note, that nested ReorderedBlockedCrsMatrices always have the same full operator "fullOp_"
207  // stored underneath for being able to "translate" the block ids.
208  if (refbX == Teuchos::null && fullOp_->getNodeNumRows() == this->getNodeNumRows())
209  {
210  // create a new (non-nested) blocked multi vector (using the blocked range map of fullOp_)
211  RCP<const BlockedMap> blkRgMap = Teuchos::rcp_dynamic_cast<const BlockedMap>(fullOp_->getRangeMap());
212  TEUCHOS_ASSERT(blkRgMap.is_null()==false);
213  RCP<const BlockedMultiVector> bXtemp = Teuchos::rcp(new BlockedMultiVector(blkRgMap, refX));
214  TEUCHOS_ASSERT(bXtemp.is_null()==false);
215  RCP<const BlockedMultiVector> bX =
216  Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(Xpetra::buildReorderedBlockedMultiVector(brm_, bXtemp));
217  TEUCHOS_ASSERT(bX.is_null()==false);
218  refbX.swap(bX);
219  bCopyResultX = true;
220  }
221 
222  if (tmpbY == Teuchos::null && fullOp_->getNodeNumRows() == this->getNodeNumRows()) {
223  // create a new (non-nested) blocked multi vector (using the blocked range map of fullOp_)
224  RCP<const BlockedMap> blkRgMap = Teuchos::rcp_dynamic_cast<const BlockedMap>(fullOp_->getRangeMap());
225  TEUCHOS_ASSERT(blkRgMap.is_null()==false);
226  RCP<BlockedMultiVector> tmpbYtemp = Teuchos::rcp(new BlockedMultiVector(blkRgMap, tmpY));
227  TEUCHOS_ASSERT(tmpbYtemp.is_null()==false);
228  RCP<BlockedMultiVector> bY =
229  Teuchos::rcp_dynamic_cast<BlockedMultiVector>(Xpetra::buildReorderedBlockedMultiVector(brm_, tmpbYtemp));
230  TEUCHOS_ASSERT(bY.is_null()==false);
231  tmpbY.swap(bY);
232  bCopyResultY = true;
233  }
234 
235  TEUCHOS_ASSERT(refbX.is_null()==false);
236  TEUCHOS_ASSERT(tmpbY.is_null()==false);
237 
239 
240  if (bCopyResultX == true) {
241  RCP<const MultiVector> Xmerged = refbX->Merge();
242  RCP<MultiVector> nonconstX = Teuchos::rcp_const_cast<MultiVector>(refX);
243  nonconstX->update(Teuchos::ScalarTraits<Scalar>::one(), *Xmerged, Teuchos::ScalarTraits<Scalar>::zero());
244  }
245  if (bCopyResultY == true) {
246  RCP< MultiVector> Ymerged = tmpbY->Merge();
247  Y.update(Teuchos::ScalarTraits<Scalar>::one(), *Ymerged, Teuchos::ScalarTraits<Scalar>::zero());
248  }
249 
250  }
251 
252  // @}
253 
254 
255 
257 
258 
260  Teuchos::RCP<const Xpetra::BlockReorderManager > getBlockReorderManager() { return brm_; }
261 
262  // @}
263 
265 
266 
268  std::string description() const { return "ReorderedBlockedCrsMatrix"; }
269 
271  void describe(Teuchos::FancyOStream& out, const Teuchos::EVerbosityLevel verbLevel = Teuchos::Describable::verbLevel_default) const {
272 
273  out << "Xpetra::ReorderedBlockedCrsMatrix: " << BlockedCrsMatrix::Rows() << " x " << BlockedCrsMatrix::Cols() << std::endl;
274 
276  out << "ReorderedBlockMatrix is fillComplete" << std::endl;
277 
278  out << "fullRowMap" << std::endl;
279  BlockedCrsMatrix::getRangeMap(0,false)->describe(out,verbLevel);
280 
281  //out << "fullColMap" << std::endl;
282  //fullcolmap_->describe(out,verbLevel);
283 
284  } else {
285  out << "Xpetra::ReorderedBlockedCrsMatrix is NOT fillComplete" << std::endl;
286  }
287 
288  for (size_t r = 0; r < BlockedCrsMatrix::Rows(); ++r)
289  for (size_t c = 0; c < BlockedCrsMatrix::Cols(); ++c) {
290  out << "Block(" << r << "," << c << ")" << std::endl;
291  BlockedCrsMatrix::getMatrix(r,c)->describe(out,verbLevel);
292  }
293  }
294 
296 
297  private:
298  Teuchos::RCP<const Xpetra::BlockReorderManager > brm_;
299  Teuchos::RCP<const Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > fullOp_;
300 
301 
302 };
303 
304 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
305 Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > mergeSubBlockMaps(Teuchos::RCP<const Xpetra::BlockReorderManager> brm, Teuchos::RCP<const Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > bmat, bool bThyraMode) {
307 
308  // TODO distinguish between range and domain map extractor! provide MapExtractor as parameter!
309  RCP<const Xpetra::MapExtractor<Scalar,LocalOrdinal,GlobalOrdinal,Node> > fullRangeMapExtractor = bmat->getRangeMapExtractor();
310 
311  // number of sub blocks
312  size_t numBlocks = brm->GetNumBlocks();
313 
314  Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node>> map = Teuchos::null;
315 
316  if(numBlocks == 0) {
317  // it is a leaf node
318  Teuchos::RCP<const Xpetra::BlockReorderLeaf> leaf = Teuchos::rcp_dynamic_cast<const Xpetra::BlockReorderLeaf>(brm);
319 
320  map = fullRangeMapExtractor->getMap(Teuchos::as<size_t>(leaf->GetIndex()), bThyraMode);
321  } else {
322  // initialize vector for sub maps
323  std::vector<Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > > subMaps (numBlocks, Teuchos::null);
324 
325  for(size_t i = 0; i < numBlocks; i++) {
326  Teuchos::RCP<const Xpetra::BlockReorderManager> blkMgr = brm->GetBlock(Teuchos::as<int>(i));
327  subMaps[i] = mergeSubBlockMaps(blkMgr,bmat,bThyraMode);
328  TEUCHOS_ASSERT(subMaps[i].is_null()==false);
329  }
330 
331 #if 1
332  // concatenate submaps
333  // for Thyra mode this map isn't important
334  Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > fullMap = MapUtils::concatenateMaps(subMaps);
335 
336  // create new BlockedMap (either in Thyra Mode or Xpetra mode)
337  map = Teuchos::rcp(new Xpetra::BlockedMap<LocalOrdinal,GlobalOrdinal,Node>(fullMap, subMaps, bThyraMode));
338 #else
339  // TAW: 11/27/16 we just concatenate the submaps to one monolithic Map object.
340  // Alternatively, we could create a new BlockedMap using the concatenated map and the submaps
341  // However, the block smoothers only need the concatenated map for creating MultiVectors...
342  // But for the Thyra mode version concatenating would not be ok for the whole map
343  map = MapUtils::concatenateMaps(subMaps);
344 #endif
345  }
346  TEUCHOS_ASSERT(map.is_null()==false);
347  return map;
348 }
349 
350 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
351 Teuchos::RCP<const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > mergeSubBlocks(Teuchos::RCP<const Xpetra::BlockReorderManager> rowMgr, Teuchos::RCP<const Xpetra::BlockReorderManager> colMgr, Teuchos::RCP<const Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > bmat) {
352 
359 
360  // number of sub blocks
361  size_t rowSz = rowMgr->GetNumBlocks();
362  size_t colSz = colMgr->GetNumBlocks();
363 
364  Teuchos::RCP<BlockedCrsMatrix> rbmat = Teuchos::null;
365 
366  if(rowSz == 0 && colSz == 0) {
367  // it is a leaf node
368  Teuchos::RCP<const Xpetra::BlockReorderLeaf> rowleaf = Teuchos::rcp_dynamic_cast<const Xpetra::BlockReorderLeaf>(rowMgr);
369  Teuchos::RCP<const Xpetra::BlockReorderLeaf> colleaf = Teuchos::rcp_dynamic_cast<const Xpetra::BlockReorderLeaf>(colMgr);
370 
371  // extract leaf node
372  Teuchos::RCP<Matrix> mat = bmat->getMatrix(rowleaf->GetIndex(), colleaf->GetIndex());
373 
374  if (mat == Teuchos::null) return Teuchos::null;
375 
376  // check, whether leaf node is of type Xpetra::CrsMatrixWrap
377  Teuchos::RCP<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node> > matwrap = Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(mat);
378  if(matwrap != Teuchos::null) {
379  // If the leaf node is of type Xpetra::CrsMatrixWrap wrap it into a 1x1 ReorderedBlockMatrix
380  // with the corresponding MapExtractors for translating Thyra to Xpetra GIDs if necessary
381  RCP<const MapExtractor> fullRangeMapExtractor = bmat->getRangeMapExtractor();
382  Teuchos::RCP<const Map> submap = fullRangeMapExtractor->getMap(rowleaf->GetIndex(),false);
383  std::vector<Teuchos::RCP<const Map> > rowSubMaps (1, submap);
384  Teuchos::RCP<const MapExtractor> rgMapExtractor = Teuchos::rcp(new MapExtractor(submap, rowSubMaps, false));
385 
386  RCP<const MapExtractor> fullDomainMapExtractor = bmat->getDomainMapExtractor();
387  Teuchos::RCP<const Map> submap2 = fullDomainMapExtractor->getMap(colleaf->GetIndex(),false);
388  std::vector<Teuchos::RCP<const Map> > colSubMaps (1, submap2);
389  Teuchos::RCP<const MapExtractor> doMapExtractor = Teuchos::rcp(new MapExtractor(submap2, colSubMaps, false));
390 
391  rbmat = Teuchos::rcp(new ReorderedBlockedCrsMatrix(rgMapExtractor,doMapExtractor, 33, rowMgr,bmat));
392  rbmat->setMatrix(0,0,mat);
393  } else {
394  // If leaf node is already wrapped into a blocked matrix do not wrap it again.
395  rbmat = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(mat);
396  TEUCHOS_ASSERT(rbmat != Teuchos::null);
397  }
398  TEUCHOS_ASSERT(mat->getNodeNumEntries() == rbmat->getNodeNumEntries());
399  } else {
400  // create the map extractors
401  // we cannot create block matrix in thyra mode since merged maps might not start with 0 GID
402  Teuchos::RCP<const MapExtractor> rgMapExtractor = Teuchos::null;
403  if(rowSz > 0) {
404  std::vector<Teuchos::RCP<const Map> > rowSubMaps (rowSz, Teuchos::null);
405  for(size_t i = 0; i < rowSz; i++) {
406  Teuchos::RCP<const Xpetra::BlockReorderManager> rowSubMgr = rowMgr->GetBlock(Teuchos::as<int>(i));
407  rowSubMaps[i] = mergeSubBlockMaps(rowSubMgr,bmat,false /*xpetra*/);
408  TEUCHOS_ASSERT(rowSubMaps[i].is_null()==false);
409  }
410  Teuchos::RCP<const Map> rgMergedSubMaps = MapUtils::concatenateMaps(rowSubMaps);
411  rgMapExtractor = Teuchos::rcp(new MapExtractor(rgMergedSubMaps, rowSubMaps, false));
412  } else {
413  Teuchos::RCP<const Xpetra::BlockReorderLeaf> rowleaf = Teuchos::rcp_dynamic_cast<const Xpetra::BlockReorderLeaf>(rowMgr);
414  RCP<const MapExtractor> fullRangeMapExtractor = bmat->getRangeMapExtractor();
415  // TODO think about Thyra style maps: we cannot use thyra style maps when recombining several blocks!!!
416  // The GIDs might not start with 0 and may not be consecutive!
417  Teuchos::RCP<const Map> submap = fullRangeMapExtractor->getMap(rowleaf->GetIndex(),false);
418  std::vector<Teuchos::RCP<const Map> > rowSubMaps (1, submap);
419  rgMapExtractor = Teuchos::rcp(new MapExtractor(submap, rowSubMaps, false));
420  }
421 
422  Teuchos::RCP<const MapExtractor> doMapExtractor = Teuchos::null;
423  if(colSz > 0) {
424  std::vector<Teuchos::RCP<const Map> > colSubMaps (colSz, Teuchos::null);
425  for(size_t j = 0; j < colSz; j++) {
426  Teuchos::RCP<const Xpetra::BlockReorderManager> colSubMgr = colMgr->GetBlock(Teuchos::as<int>(j));
427  colSubMaps[j] = mergeSubBlockMaps(colSubMgr,bmat,false/*xpetra*/);
428  TEUCHOS_ASSERT(colSubMaps[j].is_null()==false);
429  }
430  Teuchos::RCP<const Map> doMergedSubMaps = MapUtils::concatenateMaps(colSubMaps);
431  doMapExtractor = Teuchos::rcp(new MapExtractor(doMergedSubMaps, colSubMaps, false));
432  } else {
433  Teuchos::RCP<const Xpetra::BlockReorderLeaf> colleaf = Teuchos::rcp_dynamic_cast<const Xpetra::BlockReorderLeaf>(colMgr);
434  RCP<const MapExtractor> fullDomainMapExtractor = bmat->getDomainMapExtractor();
435  // TODO think about Thyra style maps: we cannot use thyra style maps when recombining several blocks!!!
436  // The GIDs might not start with 0 and may not be consecutive!
437  Teuchos::RCP<const Map> submap = fullDomainMapExtractor->getMap(colleaf->GetIndex(),false);
438  std::vector<Teuchos::RCP<const Map> > colSubMaps (1, submap);
439  doMapExtractor = Teuchos::rcp(new MapExtractor(submap, colSubMaps, false));
440  }
441 
442  rbmat = Teuchos::rcp(new ReorderedBlockedCrsMatrix(rgMapExtractor,doMapExtractor, 33, rowMgr,bmat));
443 
444  size_t cntNNZ = 0;
445 
446  if (rowSz == 0 && colSz > 0) {
447  for(size_t j = 0; j < colSz; j++) {
448  Teuchos::RCP<const Xpetra::BlockReorderManager> colSubMgr = colMgr->GetBlock(Teuchos::as<int>(j));
449  Teuchos::RCP<const Matrix> submat = mergeSubBlocks(rowMgr, colSubMgr, bmat);
450  rbmat->setMatrix(0,j,Teuchos::rcp_const_cast<Matrix>(submat));
451  if(submat != Teuchos::null) cntNNZ += submat->getNodeNumEntries();
452  }
453  } else if (rowSz > 0 && colSz == 0) {
454  for(size_t i = 0; i < rowSz; i++) {
455  Teuchos::RCP<const Xpetra::BlockReorderManager> rowSubMgr = rowMgr->GetBlock(Teuchos::as<int>(i));
456  Teuchos::RCP<const Matrix> submat = mergeSubBlocks(rowSubMgr, colMgr, bmat);
457  rbmat->setMatrix(i,0,Teuchos::rcp_const_cast<Matrix>(submat));
458  if(submat != Teuchos::null) cntNNZ += submat->getNodeNumEntries();
459  }
460  } else {
461  for(size_t i = 0; i < rowSz; i++) {
462  Teuchos::RCP<const Xpetra::BlockReorderManager> rowSubMgr = rowMgr->GetBlock(Teuchos::as<int>(i));
463  for(size_t j = 0; j < colSz; j++) {
464  Teuchos::RCP<const Xpetra::BlockReorderManager> colSubMgr = colMgr->GetBlock(Teuchos::as<int>(j));
465  Teuchos::RCP<const Matrix> submat = mergeSubBlocks(rowSubMgr, colSubMgr, bmat);
466  rbmat->setMatrix(i,j,Teuchos::rcp_const_cast<Matrix>(submat));
467  if(submat != Teuchos::null) cntNNZ += submat->getNodeNumEntries();
468  }
469  }
470  }
471  TEUCHOS_ASSERT(rbmat->getNodeNumEntries() == cntNNZ);
472  }
473  rbmat->fillComplete();
474  return rbmat;
475 }
476 
477  //MapExtractor(const std::vector<RCP<const Map> >& maps, const std::vector<RCP<const Map> >& thyramaps);
478 
479 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
480 Teuchos::RCP<const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > mergeSubBlocksThyra(Teuchos::RCP<const Xpetra::BlockReorderManager> rowMgr, Teuchos::RCP<const Xpetra::BlockReorderManager> colMgr, Teuchos::RCP<const Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > bmat) {
481 
487 
488  TEUCHOS_ASSERT(bmat->getRangeMapExtractor()->getThyraMode() == true);
489  TEUCHOS_ASSERT(bmat->getDomainMapExtractor()->getThyraMode() == true);
490 
491  // number of sub blocks
492  size_t rowSz = rowMgr->GetNumBlocks();
493  size_t colSz = colMgr->GetNumBlocks();
494 
495  Teuchos::RCP<BlockedCrsMatrix> rbmat = Teuchos::null;
496 
497  if(rowSz == 0 && colSz == 0) {
498  // it is a leaf node
499  Teuchos::RCP<const Xpetra::BlockReorderLeaf> rowleaf = Teuchos::rcp_dynamic_cast<const Xpetra::BlockReorderLeaf>(rowMgr);
500  Teuchos::RCP<const Xpetra::BlockReorderLeaf> colleaf = Teuchos::rcp_dynamic_cast<const Xpetra::BlockReorderLeaf>(colMgr);
501 
502  // this matrix uses Thyra style GIDs as global row, range, domain and column indices
503  Teuchos::RCP<Matrix> mat = bmat->getMatrix(rowleaf->GetIndex(), colleaf->GetIndex());
504 
505  if(mat == Teuchos::null) return Teuchos::null; //std::cout << "Block " << rowleaf->GetIndex() << "," << colleaf->GetIndex() << " is zero" << std::endl;
506 
507  // check, whether leaf node is of type Xpetra::CrsMatrixWrap
508  Teuchos::RCP<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node> > matwrap = Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(mat);
509  if(matwrap != Teuchos::null) {
511  // build map extractors
512  RCP<const MapExtractor> fullRangeMapExtractor = bmat->getRangeMapExtractor();
513  // extract Xpetra and Thyra based GIDs
514  Teuchos::RCP<const Map> xpsubmap = fullRangeMapExtractor->getMap(rowleaf->GetIndex(),false);
515  Teuchos::RCP<const Map> thysubmap = fullRangeMapExtractor->getMap(rowleaf->GetIndex(),true);
516  std::vector<Teuchos::RCP<const Map> > rowXpSubMaps (1, xpsubmap);
517  std::vector<Teuchos::RCP<const Map> > rowTySubMaps (1, thysubmap);
518  // use expert constructor
519  Teuchos::RCP<const MapExtractor> rgMapExtractor = Teuchos::rcp(new MapExtractor(rowXpSubMaps, rowTySubMaps));
520 
521  RCP<const MapExtractor> fullDomainMapExtractor = bmat->getDomainMapExtractor();
522  // extract Xpetra and Thyra based GIDs
523  Teuchos::RCP<const Map> xpsubmap2 = fullDomainMapExtractor->getMap(colleaf->GetIndex(),false);
524  Teuchos::RCP<const Map> tysubmap2 = fullDomainMapExtractor->getMap(colleaf->GetIndex(),true);
525  std::vector<Teuchos::RCP<const Map> > colXpSubMaps (1, xpsubmap2);
526  std::vector<Teuchos::RCP<const Map> > colTySubMaps (1, tysubmap2);
527  // use expert constructor
528  Teuchos::RCP<const MapExtractor> doMapExtractor = Teuchos::rcp(new MapExtractor(colXpSubMaps, colTySubMaps));
529 
531  // build reordered block operator
532  rbmat = Teuchos::rcp(new ReorderedBlockedCrsMatrix(rgMapExtractor,doMapExtractor, 33, rowMgr,bmat));
533  rbmat->setMatrix(0,0,mat);
534  } else {
535  // If leaf node is already wrapped into a blocked matrix do not wrap it again.
536  rbmat = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(mat);
537  TEUCHOS_ASSERT(rbmat != Teuchos::null);
538  }
539  TEUCHOS_ASSERT(mat->getNodeNumEntries() == rbmat->getNodeNumEntries());
540  } else {
541  // create the map extractors
542  // we cannot create block matrix in thyra mode since merged maps might not start with 0 GID
543  Teuchos::RCP<const MapExtractor> rgMapExtractor = Teuchos::null;
544  if(rowSz > 0) {
545  std::vector<Teuchos::RCP<const Map> > rowXpSubMaps (rowSz, Teuchos::null);
546  std::vector<Teuchos::RCP<const Map> > rowTySubMaps (rowSz, Teuchos::null);
547  for(size_t i = 0; i < rowSz; i++) {
548  Teuchos::RCP<const Xpetra::BlockReorderManager> rowSubMgr = rowMgr->GetBlock(Teuchos::as<int>(i));
549  // extract Xpetra and Thyra based merged GIDs
550  rowXpSubMaps[i] = mergeSubBlockMaps(rowSubMgr,bmat,false);
551  rowTySubMaps[i] = mergeSubBlockMaps(rowSubMgr,bmat,true);
552  TEUCHOS_ASSERT(rowXpSubMaps[i].is_null()==false);
553  TEUCHOS_ASSERT(rowTySubMaps[i].is_null()==false);
554  }
555  // use expert constructor
556  rgMapExtractor = Teuchos::rcp(new MapExtractor(rowXpSubMaps, rowTySubMaps));
557  } else {
558  Teuchos::RCP<const Xpetra::BlockReorderLeaf> rowleaf = Teuchos::rcp_dynamic_cast<const Xpetra::BlockReorderLeaf>(rowMgr);
559  RCP<const MapExtractor> fullRangeMapExtractor = bmat->getRangeMapExtractor();
560  // extract Xpetra and Thyra based GIDs
561  Teuchos::RCP<const Map> xpsubmap = fullRangeMapExtractor->getMap(rowleaf->GetIndex(),false);
562  Teuchos::RCP<const Map> thysubmap = fullRangeMapExtractor->getMap(rowleaf->GetIndex(),true);
563  std::vector<Teuchos::RCP<const Map> > rowXpSubMaps (1, xpsubmap);
564  std::vector<Teuchos::RCP<const Map> > rowTySubMaps (1, thysubmap);
565  // use expert constructor
566  rgMapExtractor = Teuchos::rcp(new MapExtractor(rowXpSubMaps, rowTySubMaps));
567  }
568 
569  Teuchos::RCP<const MapExtractor> doMapExtractor = Teuchos::null;
570  if(colSz > 0) {
571  std::vector<Teuchos::RCP<const Map> > colXpSubMaps (colSz, Teuchos::null);
572  std::vector<Teuchos::RCP<const Map> > colTySubMaps (colSz, Teuchos::null);
573  for(size_t j = 0; j < colSz; j++) {
574  Teuchos::RCP<const Xpetra::BlockReorderManager> colSubMgr = colMgr->GetBlock(Teuchos::as<int>(j));
575  // extract Xpetra and Thyra based merged GIDs
576  colXpSubMaps[j] = mergeSubBlockMaps(colSubMgr,bmat,false);
577  colTySubMaps[j] = mergeSubBlockMaps(colSubMgr,bmat,true);
578  TEUCHOS_ASSERT(colXpSubMaps[j].is_null()==false);
579  TEUCHOS_ASSERT(colTySubMaps[j].is_null()==false);
580  }
581  // use expert constructor
582  doMapExtractor = Teuchos::rcp(new MapExtractor(colXpSubMaps,colTySubMaps));
583  } else {
584  Teuchos::RCP<const Xpetra::BlockReorderLeaf> colleaf = Teuchos::rcp_dynamic_cast<const Xpetra::BlockReorderLeaf>(colMgr);
585  RCP<const MapExtractor> fullDomainMapExtractor = bmat->getDomainMapExtractor();
586  // extract Xpetra and Thyra based GIDs
587  Teuchos::RCP<const Map> xpsubmap = fullDomainMapExtractor->getMap(colleaf->GetIndex(),false);
588  Teuchos::RCP<const Map> tysubmap = fullDomainMapExtractor->getMap(colleaf->GetIndex(),true);
589  std::vector<Teuchos::RCP<const Map> > colXpSubMaps (1, xpsubmap);
590  std::vector<Teuchos::RCP<const Map> > colTySubMaps (1, tysubmap);
591  // use expert constructor
592  doMapExtractor = Teuchos::rcp(new MapExtractor(colXpSubMaps, colTySubMaps));
593  }
594 
595  // TODO matrix should have both rowMgr and colMgr??
596  rbmat = Teuchos::rcp(new ReorderedBlockedCrsMatrix(rgMapExtractor,doMapExtractor, 33, rowMgr,bmat));
597 
598  size_t cntNNZ = 0;
599 
600  if (rowSz == 0 && colSz > 0) {
601  for(size_t j = 0; j < colSz; j++) {
602  Teuchos::RCP<const Xpetra::BlockReorderManager> colSubMgr = colMgr->GetBlock(Teuchos::as<int>(j));
603  Teuchos::RCP<const Matrix> submat = mergeSubBlocksThyra(rowMgr, colSubMgr, bmat);
604  rbmat->setMatrix(0,j,Teuchos::rcp_const_cast<Matrix>(submat));
605  if(submat != Teuchos::null) cntNNZ += submat->getNodeNumEntries();
606  }
607  } else if (rowSz > 0 && colSz == 0) {
608  for(size_t i = 0; i < rowSz; i++) {
609  Teuchos::RCP<const Xpetra::BlockReorderManager> rowSubMgr = rowMgr->GetBlock(Teuchos::as<int>(i));
610  Teuchos::RCP<const Matrix> submat = mergeSubBlocksThyra(rowSubMgr, colMgr, bmat);
611  rbmat->setMatrix(i,0,Teuchos::rcp_const_cast<Matrix>(submat));
612  if(submat != Teuchos::null) cntNNZ += submat->getNodeNumEntries();
613  }
614  } else {
615  for(size_t i = 0; i < rowSz; i++) {
616  Teuchos::RCP<const Xpetra::BlockReorderManager> rowSubMgr = rowMgr->GetBlock(Teuchos::as<int>(i));
617  for(size_t j = 0; j < colSz; j++) {
618  Teuchos::RCP<const Xpetra::BlockReorderManager> colSubMgr = colMgr->GetBlock(Teuchos::as<int>(j));
619  Teuchos::RCP<const Matrix> submat = mergeSubBlocksThyra(rowSubMgr, colSubMgr, bmat);
620  rbmat->setMatrix(i,j,Teuchos::rcp_const_cast<Matrix>(submat));
621  if(submat != Teuchos::null) cntNNZ += submat->getNodeNumEntries();
622  }
623  }
624  }
625  TEUCHOS_ASSERT(rbmat->getNodeNumEntries() == cntNNZ);
626  }
627 
628  rbmat->fillComplete();
629  return rbmat;
630 }
631 
632 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
633 Teuchos::RCP<const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > buildReorderedBlockedCrsMatrix(Teuchos::RCP<const Xpetra::BlockReorderManager> brm, Teuchos::RCP<const Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > bmat) {
634  TEUCHOS_ASSERT(bmat->getRangeMapExtractor()->getThyraMode() == bmat->getDomainMapExtractor()->getThyraMode());
635  Teuchos::RCP<const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > rbmat = Teuchos::null;
636  if(bmat->getRangeMapExtractor()->getThyraMode() == false) {
637  rbmat = mergeSubBlocks(brm, brm, bmat);
638  } else {
639  rbmat = mergeSubBlocksThyra(brm, brm, bmat);
640  }
641 
642  // TAW, 6/7/2016: rbmat might be Teuchos::null for empty blocks!
643  return rbmat;
644 }
645 
646 } //namespace Xpetra
647 
648 #define XPETRA_REORDEREDBLOCKEDCRSMATRIX_SHORT
649 #endif /* XPETRA_REORDEREDBLOCKEDCRSMATRIX_HPP */
Teuchos::RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > mergeSubBlockMaps(Teuchos::RCP< const Xpetra::BlockReorderManager > brm, Teuchos::RCP< const Xpetra::BlockedCrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > bmat, bool bThyraMode)
bool isFillComplete() const
Returns true if fillComplete() has been called and the matrix is in compute mode. ...
RCP< const Map > getRangeMap() const
Returns the Map associated with the range of this operator.
const Teuchos::RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > getMap() const
Implements DistObject interface.
Teuchos::RCP< const Xpetra::BlockedCrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > fullOp_
Xpetra utility class for common map-related routines.
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.
Teuchos::RCP< const Xpetra::BlockReorderManager > getBlockReorderManager()
Returns internal BlockReorderManager object.
static Teuchos::RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > concatenateMaps(const std::vector< Teuchos::RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > > &subMaps)
Helper function to concatenate several maps.
Teuchos::RCP< const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > mergeSubBlocksThyra(Teuchos::RCP< const Xpetra::BlockReorderManager > rowMgr, Teuchos::RCP< const Xpetra::BlockReorderManager > colMgr, Teuchos::RCP< const Xpetra::BlockedCrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > bmat)
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.
Teuchos::RCP< const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > mergeSubBlocks(Teuchos::RCP< const Xpetra::BlockReorderManager > rowMgr, Teuchos::RCP< const Xpetra::BlockReorderManager > colMgr, Teuchos::RCP< const Xpetra::BlockedCrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > bmat)
Teuchos::RCP< Matrix > getMatrix(size_t r, size_t c) const
return block (r,c)
virtual size_t Cols() const
number of column blocks
std::string viewLabel_t
std::string description() const
Return a simple one-line description of this object.
Teuchos::RCP< const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > buildReorderedBlockedCrsMatrix(Teuchos::RCP< const Xpetra::BlockReorderManager > brm, Teuchos::RCP< const Xpetra::BlockedCrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > bmat)
Teuchos::RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > mergeSubBlockMaps(Teuchos::RCP< const Xpetra::BlockReorderManager > brm)
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.
Concrete implementation of Xpetra::Matrix.
virtual size_t Rows() const
number of row blocks
Teuchos::RCP< const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > buildReorderedBlockedMultiVector(Teuchos::RCP< const Xpetra::BlockReorderManager > brm, Teuchos::RCP< const Xpetra::BlockedMultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > bvec)
size_t getNodeNumRows() const
Returns the number of matrix rows owned on the calling node.
ReorderedBlockedCrsMatrix(Teuchos::RCP< const MapExtractor > &rangeMaps, Teuchos::RCP< const MapExtractor > &domainMaps, size_t npr, Teuchos::RCP< const Xpetra::BlockReorderManager > brm, Teuchos::RCP< const Xpetra::BlockedCrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > bmat)
Constructor.
Xpetra-specific matrix class.
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.
Teuchos::RCP< const Xpetra::BlockReorderManager > brm_