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_CrsMatrixWrap.hpp"
58 
59 
64 namespace Xpetra {
65 
66  typedef std::string viewLabel_t;
67 
68  template <class Scalar,
69  class LocalOrdinal,
70  class GlobalOrdinal,
71  class Node>
72  class ReorderedBlockedCrsMatrix :
73  public BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> {
74  public:
75  typedef Scalar scalar_type;
76  typedef LocalOrdinal local_ordinal_type;
77  typedef GlobalOrdinal global_ordinal_type;
78  typedef Node node_type;
79 
80  private:
81 #undef XPETRA_REORDEREDBLOCKEDCRSMATRIX_SHORT
82 #include "Xpetra_UseShortNames.hpp"
83 
84  public:
85 
87 
88 
90 
99  (Teuchos::RCP<const MapExtractor>& rangeMaps,
100  Teuchos::RCP<const MapExtractor>& domainMaps,
101  size_t npr,
102  Teuchos::RCP<const Xpetra::BlockReorderManager> brm,
105  : Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>(rangeMaps, domainMaps, npr,pftype) {
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  {
162  }
163 
164  // @}
165 
167 
168 
170  std::string description() const { return "ReorderedBlockedCrsMatrix"; }
171 
173  void describe(Teuchos::FancyOStream& out, const Teuchos::EVerbosityLevel verbLevel = Teuchos::Describable::verbLevel_default) const {
175 
176  out << "Xpetra::ReorderedBlockedCrsMatrix: " << BlockedCrsMatrix::Rows() << " x " << BlockedCrsMatrix::Cols() << std::endl;
177 
179  out << "ReorderedBlockMatrix is fillComplete" << std::endl;
180 
181  out << "fullRowMap" << std::endl;
182  BlockedCrsMatrix::getRangeMap(0,false)->describe(out,verbLevel);
183 
184  //out << "fullColMap" << std::endl;
185  //fullcolmap_->describe(out,verbLevel);
186 
187  } else {
188  out << "Xpetra::ReorderedBlockedCrsMatrix is NOT fillComplete" << std::endl;
189  }
190 
191  for (size_t r = 0; r < BlockedCrsMatrix::Rows(); ++r)
192  for (size_t c = 0; c < BlockedCrsMatrix::Cols(); ++c) {
193  out << "Block(" << r << "," << c << ")" << std::endl;
194  BlockedCrsMatrix::getMatrix(r,c)->describe(out,verbLevel);
195  }
196  }
197 
199 
200  private:
201  Teuchos::RCP<const Xpetra::BlockReorderManager > brm_;
202  Teuchos::RCP<const Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > fullOp_;
203 
204 
205 };
206 
207 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
208 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) {
210 
211  // TODO distinguish between range and domain map extractor! provide MapExtractor as parameter!
212  RCP<const Xpetra::MapExtractor<Scalar,LocalOrdinal,GlobalOrdinal,Node> > fullRangeMapExtractor = bmat->getRangeMapExtractor();
213 
214  // number of sub blocks
215  size_t numBlocks = brm->GetNumBlocks();
216 
217  Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node>> map = Teuchos::null;
218 
219  if(numBlocks == 0) {
220  // it is a leaf node
221  Teuchos::RCP<const Xpetra::BlockReorderLeaf> leaf = Teuchos::rcp_dynamic_cast<const Xpetra::BlockReorderLeaf>(brm);
222 
223  map = fullRangeMapExtractor->getMap(Teuchos::as<size_t>(leaf->GetIndex()), bThyraMode);
224  } else {
225  // initialize vector for sub maps
226  std::vector<Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > > subMaps (numBlocks, Teuchos::null);
227 
228  for(size_t i = 0; i < numBlocks; i++) {
229  Teuchos::RCP<const Xpetra::BlockReorderManager> blkMgr = brm->GetBlock(Teuchos::as<int>(i));
230  subMaps[i] = mergeSubBlockMaps(blkMgr,bmat,bThyraMode);
231  TEUCHOS_ASSERT(subMaps[i].is_null()==false);
232  }
233 
234  map = MapUtils::concatenateMaps(subMaps);
235  }
236  TEUCHOS_ASSERT(map.is_null()==false);
237  return map;
238 }
239 
240 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
241 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) {
242 
249 
250  // number of sub blocks
251  size_t rowSz = rowMgr->GetNumBlocks();
252  size_t colSz = colMgr->GetNumBlocks();
253 
254  Teuchos::RCP<BlockedCrsMatrix> rbmat = Teuchos::null;
255 
256  if(rowSz == 0 && colSz == 0) {
257  // it is a leaf node
258  Teuchos::RCP<const Xpetra::BlockReorderLeaf> rowleaf = Teuchos::rcp_dynamic_cast<const Xpetra::BlockReorderLeaf>(rowMgr);
259  Teuchos::RCP<const Xpetra::BlockReorderLeaf> colleaf = Teuchos::rcp_dynamic_cast<const Xpetra::BlockReorderLeaf>(colMgr);
260 
261  // extract leaf node
262  Teuchos::RCP<Matrix> mat = bmat->getMatrix(rowleaf->GetIndex(), colleaf->GetIndex());
263 
264  if (mat == Teuchos::null) return Teuchos::null;
265 
266  // check, whether leaf node is of type Xpetra::CrsMatrixWrap
267  Teuchos::RCP<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node> > matwrap = Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(mat);
268  if(matwrap != Teuchos::null) {
269  // If the leaf node is of type Xpetra::CrsMatrixWrap wrap it into a 1x1 ReorderedBlockMatrix
270  // with the corresponding MapExtractors for translating Thyra to Xpetra GIDs if necessary
271  RCP<const MapExtractor> fullRangeMapExtractor = bmat->getRangeMapExtractor();
272  Teuchos::RCP<const Map> submap = fullRangeMapExtractor->getMap(rowleaf->GetIndex(),false);
273  std::vector<Teuchos::RCP<const Map> > rowSubMaps (1, submap);
274  Teuchos::RCP<const MapExtractor> rgMapExtractor = Teuchos::rcp(new MapExtractor(submap, rowSubMaps, false));
275 
276  RCP<const MapExtractor> fullDomainMapExtractor = bmat->getDomainMapExtractor();
277  Teuchos::RCP<const Map> submap2 = fullDomainMapExtractor->getMap(colleaf->GetIndex(),false);
278  std::vector<Teuchos::RCP<const Map> > colSubMaps (1, submap2);
279  Teuchos::RCP<const MapExtractor> doMapExtractor = Teuchos::rcp(new MapExtractor(submap2, colSubMaps, false));
280 
281  rbmat = Teuchos::rcp(new ReorderedBlockedCrsMatrix(rgMapExtractor,doMapExtractor, 33, rowMgr,bmat));
282  rbmat->setMatrix(0,0,mat);
283  } else {
284  // If leaf node is already wrapped into a blocked matrix do not wrap it again.
285  rbmat = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(mat);
286  TEUCHOS_ASSERT(rbmat != Teuchos::null);
287  }
288  TEUCHOS_ASSERT(mat->getNodeNumEntries() == rbmat->getNodeNumEntries());
289  } else {
290  // create the map extractors
291  // we cannot create block matrix in thyra mode since merged maps might not start with 0 GID
292  Teuchos::RCP<const MapExtractor> rgMapExtractor = Teuchos::null;
293  if(rowSz > 0) {
294  std::vector<Teuchos::RCP<const Map> > rowSubMaps (rowSz, Teuchos::null);
295  for(size_t i = 0; i < rowSz; i++) {
296  Teuchos::RCP<const Xpetra::BlockReorderManager> rowSubMgr = rowMgr->GetBlock(Teuchos::as<int>(i));
297  rowSubMaps[i] = mergeSubBlockMaps(rowSubMgr,bmat,false /*xpetra*/);
298  TEUCHOS_ASSERT(rowSubMaps[i].is_null()==false);
299  }
300  Teuchos::RCP<const Map> rgMergedSubMaps = MapUtils::concatenateMaps(rowSubMaps);
301  rgMapExtractor = Teuchos::rcp(new MapExtractor(rgMergedSubMaps, rowSubMaps, false));
302  } else {
303  Teuchos::RCP<const Xpetra::BlockReorderLeaf> rowleaf = Teuchos::rcp_dynamic_cast<const Xpetra::BlockReorderLeaf>(rowMgr);
304  RCP<const MapExtractor> fullRangeMapExtractor = bmat->getRangeMapExtractor();
305  // TODO think about Thyra style maps: we cannot use thyra style maps when recombining several blocks!!!
306  // The GIDs might not start with 0 and may not be consecutive!
307  Teuchos::RCP<const Map> submap = fullRangeMapExtractor->getMap(rowleaf->GetIndex(),false);
308  std::vector<Teuchos::RCP<const Map> > rowSubMaps (1, submap);
309  rgMapExtractor = Teuchos::rcp(new MapExtractor(submap, rowSubMaps, false));
310  }
311 
312  Teuchos::RCP<const MapExtractor> doMapExtractor = Teuchos::null;
313  if(colSz > 0) {
314  std::vector<Teuchos::RCP<const Map> > colSubMaps (colSz, Teuchos::null);
315  for(size_t j = 0; j < colSz; j++) {
316  Teuchos::RCP<const Xpetra::BlockReorderManager> colSubMgr = colMgr->GetBlock(Teuchos::as<int>(j));
317  colSubMaps[j] = mergeSubBlockMaps(colSubMgr,bmat,false/*xpetra*/);
318  TEUCHOS_ASSERT(colSubMaps[j].is_null()==false);
319  }
320  Teuchos::RCP<const Map> doMergedSubMaps = MapUtils::concatenateMaps(colSubMaps);
321  doMapExtractor = Teuchos::rcp(new MapExtractor(doMergedSubMaps, colSubMaps, false));
322  } else {
323  Teuchos::RCP<const Xpetra::BlockReorderLeaf> colleaf = Teuchos::rcp_dynamic_cast<const Xpetra::BlockReorderLeaf>(colMgr);
324  RCP<const MapExtractor> fullDomainMapExtractor = bmat->getDomainMapExtractor();
325  // TODO think about Thyra style maps: we cannot use thyra style maps when recombining several blocks!!!
326  // The GIDs might not start with 0 and may not be consecutive!
327  Teuchos::RCP<const Map> submap = fullDomainMapExtractor->getMap(colleaf->GetIndex(),false);
328  std::vector<Teuchos::RCP<const Map> > colSubMaps (1, submap);
329  doMapExtractor = Teuchos::rcp(new MapExtractor(submap, colSubMaps, false));
330  }
331 
332  rbmat = Teuchos::rcp(new ReorderedBlockedCrsMatrix(rgMapExtractor,doMapExtractor, 33, rowMgr,bmat));
333 
334  size_t cntNNZ = 0;
335 
336  if (rowSz == 0 && colSz > 0) {
337  for(size_t j = 0; j < colSz; j++) {
338  Teuchos::RCP<const Xpetra::BlockReorderManager> colSubMgr = colMgr->GetBlock(Teuchos::as<int>(j));
339  Teuchos::RCP<const Matrix> submat = mergeSubBlocks(rowMgr, colSubMgr, bmat);
340  rbmat->setMatrix(0,j,Teuchos::rcp_const_cast<Matrix>(submat));
341  if(submat != Teuchos::null) cntNNZ += submat->getNodeNumEntries();
342  }
343  } else if (rowSz > 0 && colSz == 0) {
344  for(size_t i = 0; i < rowSz; i++) {
345  Teuchos::RCP<const Xpetra::BlockReorderManager> rowSubMgr = rowMgr->GetBlock(Teuchos::as<int>(i));
346  Teuchos::RCP<const Matrix> submat = mergeSubBlocks(rowSubMgr, colMgr, bmat);
347  rbmat->setMatrix(i,0,Teuchos::rcp_const_cast<Matrix>(submat));
348  if(submat != Teuchos::null) cntNNZ += submat->getNodeNumEntries();
349  }
350  } else {
351  for(size_t i = 0; i < rowSz; i++) {
352  Teuchos::RCP<const Xpetra::BlockReorderManager> rowSubMgr = rowMgr->GetBlock(Teuchos::as<int>(i));
353  for(size_t j = 0; j < colSz; j++) {
354  Teuchos::RCP<const Xpetra::BlockReorderManager> colSubMgr = colMgr->GetBlock(Teuchos::as<int>(j));
355  Teuchos::RCP<const Matrix> submat = mergeSubBlocks(rowSubMgr, colSubMgr, bmat);
356  rbmat->setMatrix(i,j,Teuchos::rcp_const_cast<Matrix>(submat));
357  if(submat != Teuchos::null) cntNNZ += submat->getNodeNumEntries();
358  }
359  }
360  }
361  TEUCHOS_ASSERT(rbmat->getNodeNumEntries() == cntNNZ);
362  }
363  rbmat->fillComplete();
364  return rbmat;
365 }
366 
367  //MapExtractor(const std::vector<RCP<const Map> >& maps, const std::vector<RCP<const Map> >& thyramaps);
368 
369 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
370 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) {
371 
377 
378  TEUCHOS_ASSERT(bmat->getRangeMapExtractor()->getThyraMode() == true);
379  TEUCHOS_ASSERT(bmat->getDomainMapExtractor()->getThyraMode() == true);
380 
381  // number of sub blocks
382  size_t rowSz = rowMgr->GetNumBlocks();
383  size_t colSz = colMgr->GetNumBlocks();
384 
385  Teuchos::RCP<BlockedCrsMatrix> rbmat = Teuchos::null;
386 
387  if(rowSz == 0 && colSz == 0) {
388  // it is a leaf node
389  Teuchos::RCP<const Xpetra::BlockReorderLeaf> rowleaf = Teuchos::rcp_dynamic_cast<const Xpetra::BlockReorderLeaf>(rowMgr);
390  Teuchos::RCP<const Xpetra::BlockReorderLeaf> colleaf = Teuchos::rcp_dynamic_cast<const Xpetra::BlockReorderLeaf>(colMgr);
391 
392  // this matrix uses Thyra style GIDs as global row, range, domain and column indices
393  Teuchos::RCP<Matrix> mat = bmat->getMatrix(rowleaf->GetIndex(), colleaf->GetIndex());
394 
395  if(mat == Teuchos::null) return Teuchos::null; //std::cout << "Block " << rowleaf->GetIndex() << "," << colleaf->GetIndex() << " is zero" << std::endl;
396 
397  // check, whether leaf node is of type Xpetra::CrsMatrixWrap
398  Teuchos::RCP<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node> > matwrap = Teuchos::rcp_dynamic_cast<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(mat);
399  if(matwrap != Teuchos::null) {
401  // build map extractors
402  RCP<const MapExtractor> fullRangeMapExtractor = bmat->getRangeMapExtractor();
403  // extract Xpetra and Thyra based GIDs
404  Teuchos::RCP<const Map> xpsubmap = fullRangeMapExtractor->getMap(rowleaf->GetIndex(),false);
405  Teuchos::RCP<const Map> thysubmap = fullRangeMapExtractor->getMap(rowleaf->GetIndex(),true);
406  std::vector<Teuchos::RCP<const Map> > rowXpSubMaps (1, xpsubmap);
407  std::vector<Teuchos::RCP<const Map> > rowTySubMaps (1, thysubmap);
408  // use expert constructor
409  Teuchos::RCP<const MapExtractor> rgMapExtractor = Teuchos::rcp(new MapExtractor(rowXpSubMaps, rowTySubMaps));
410 
411  RCP<const MapExtractor> fullDomainMapExtractor = bmat->getDomainMapExtractor();
412  // extract Xpetra and Thyra based GIDs
413  Teuchos::RCP<const Map> xpsubmap2 = fullDomainMapExtractor->getMap(colleaf->GetIndex(),false);
414  Teuchos::RCP<const Map> tysubmap2 = fullDomainMapExtractor->getMap(colleaf->GetIndex(),true);
415  std::vector<Teuchos::RCP<const Map> > colXpSubMaps (1, xpsubmap2);
416  std::vector<Teuchos::RCP<const Map> > colTySubMaps (1, tysubmap2);
417  // use expert constructor
418  Teuchos::RCP<const MapExtractor> doMapExtractor = Teuchos::rcp(new MapExtractor(colXpSubMaps, colTySubMaps));
419 
421  // build reordered block operator
422  rbmat = Teuchos::rcp(new ReorderedBlockedCrsMatrix(rgMapExtractor,doMapExtractor, 33, rowMgr,bmat));
423  rbmat->setMatrix(0,0,mat);
424  } else {
425  // If leaf node is already wrapped into a blocked matrix do not wrap it again.
426  rbmat = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(mat);
427  TEUCHOS_ASSERT(rbmat != Teuchos::null);
428  }
429  TEUCHOS_ASSERT(mat->getNodeNumEntries() == rbmat->getNodeNumEntries());
430  } else {
431  // create the map extractors
432  // we cannot create block matrix in thyra mode since merged maps might not start with 0 GID
433  Teuchos::RCP<const MapExtractor> rgMapExtractor = Teuchos::null;
434  if(rowSz > 0) {
435  std::vector<Teuchos::RCP<const Map> > rowXpSubMaps (rowSz, Teuchos::null);
436  std::vector<Teuchos::RCP<const Map> > rowTySubMaps (rowSz, Teuchos::null);
437  for(size_t i = 0; i < rowSz; i++) {
438  Teuchos::RCP<const Xpetra::BlockReorderManager> rowSubMgr = rowMgr->GetBlock(Teuchos::as<int>(i));
439  // extract Xpetra and Thyra based merged GIDs
440  rowXpSubMaps[i] = mergeSubBlockMaps(rowSubMgr,bmat,false);
441  rowTySubMaps[i] = mergeSubBlockMaps(rowSubMgr,bmat,true);
442  TEUCHOS_ASSERT(rowXpSubMaps[i].is_null()==false);
443  TEUCHOS_ASSERT(rowTySubMaps[i].is_null()==false);
444  }
445  // use expert constructor
446  rgMapExtractor = Teuchos::rcp(new MapExtractor(rowXpSubMaps, rowTySubMaps));
447  } else {
448  Teuchos::RCP<const Xpetra::BlockReorderLeaf> rowleaf = Teuchos::rcp_dynamic_cast<const Xpetra::BlockReorderLeaf>(rowMgr);
449  RCP<const MapExtractor> fullRangeMapExtractor = bmat->getRangeMapExtractor();
450  // extract Xpetra and Thyra based GIDs
451  Teuchos::RCP<const Map> xpsubmap = fullRangeMapExtractor->getMap(rowleaf->GetIndex(),false);
452  Teuchos::RCP<const Map> thysubmap = fullRangeMapExtractor->getMap(rowleaf->GetIndex(),true);
453  std::vector<Teuchos::RCP<const Map> > rowXpSubMaps (1, xpsubmap);
454  std::vector<Teuchos::RCP<const Map> > rowTySubMaps (1, thysubmap);
455  // use expert constructor
456  rgMapExtractor = Teuchos::rcp(new MapExtractor(rowXpSubMaps, rowTySubMaps));
457  }
458 
459  Teuchos::RCP<const MapExtractor> doMapExtractor = Teuchos::null;
460  if(colSz > 0) {
461  std::vector<Teuchos::RCP<const Map> > colXpSubMaps (colSz, Teuchos::null);
462  std::vector<Teuchos::RCP<const Map> > colTySubMaps (colSz, Teuchos::null);
463  for(size_t j = 0; j < colSz; j++) {
464  Teuchos::RCP<const Xpetra::BlockReorderManager> colSubMgr = colMgr->GetBlock(Teuchos::as<int>(j));
465  // extract Xpetra and Thyra based merged GIDs
466  colXpSubMaps[j] = mergeSubBlockMaps(colSubMgr,bmat,false);
467  colTySubMaps[j] = mergeSubBlockMaps(colSubMgr,bmat,true);
468  TEUCHOS_ASSERT(colXpSubMaps[j].is_null()==false);
469  TEUCHOS_ASSERT(colTySubMaps[j].is_null()==false);
470  }
471  // use expert constructor
472  doMapExtractor = Teuchos::rcp(new MapExtractor(colXpSubMaps,colTySubMaps));
473  } else {
474  Teuchos::RCP<const Xpetra::BlockReorderLeaf> colleaf = Teuchos::rcp_dynamic_cast<const Xpetra::BlockReorderLeaf>(colMgr);
475  RCP<const MapExtractor> fullDomainMapExtractor = bmat->getDomainMapExtractor();
476  // extract Xpetra and Thyra based GIDs
477  Teuchos::RCP<const Map> xpsubmap = fullDomainMapExtractor->getMap(colleaf->GetIndex(),false);
478  Teuchos::RCP<const Map> tysubmap = fullDomainMapExtractor->getMap(colleaf->GetIndex(),true);
479  std::vector<Teuchos::RCP<const Map> > colXpSubMaps (1, xpsubmap);
480  std::vector<Teuchos::RCP<const Map> > colTySubMaps (1, tysubmap);
481  // use expert constructor
482  doMapExtractor = Teuchos::rcp(new MapExtractor(colXpSubMaps, colTySubMaps));
483  }
484 
485  // TODO matrix should have both rowMgr and colMgr??
486  rbmat = Teuchos::rcp(new ReorderedBlockedCrsMatrix(rgMapExtractor,doMapExtractor, 33, rowMgr,bmat));
487 
488  size_t cntNNZ = 0;
489 
490  if (rowSz == 0 && colSz > 0) {
491  for(size_t j = 0; j < colSz; j++) {
492  Teuchos::RCP<const Xpetra::BlockReorderManager> colSubMgr = colMgr->GetBlock(Teuchos::as<int>(j));
493  Teuchos::RCP<const Matrix> submat = mergeSubBlocksThyra(rowMgr, colSubMgr, bmat);
494  rbmat->setMatrix(0,j,Teuchos::rcp_const_cast<Matrix>(submat));
495  if(submat != Teuchos::null) cntNNZ += submat->getNodeNumEntries();
496  }
497  } else if (rowSz > 0 && colSz == 0) {
498  for(size_t i = 0; i < rowSz; i++) {
499  Teuchos::RCP<const Xpetra::BlockReorderManager> rowSubMgr = rowMgr->GetBlock(Teuchos::as<int>(i));
500  Teuchos::RCP<const Matrix> submat = mergeSubBlocksThyra(rowSubMgr, colMgr, bmat);
501  rbmat->setMatrix(i,0,Teuchos::rcp_const_cast<Matrix>(submat));
502  if(submat != Teuchos::null) cntNNZ += submat->getNodeNumEntries();
503  }
504  } else {
505  for(size_t i = 0; i < rowSz; i++) {
506  Teuchos::RCP<const Xpetra::BlockReorderManager> rowSubMgr = rowMgr->GetBlock(Teuchos::as<int>(i));
507  for(size_t j = 0; j < colSz; j++) {
508  Teuchos::RCP<const Xpetra::BlockReorderManager> colSubMgr = colMgr->GetBlock(Teuchos::as<int>(j));
509  Teuchos::RCP<const Matrix> submat = mergeSubBlocksThyra(rowSubMgr, colSubMgr, bmat);
510  rbmat->setMatrix(i,j,Teuchos::rcp_const_cast<Matrix>(submat));
511  if(submat != Teuchos::null) cntNNZ += submat->getNodeNumEntries();
512  }
513  }
514  }
515  TEUCHOS_ASSERT(rbmat->getNodeNumEntries() == cntNNZ);
516  }
517 
518  rbmat->fillComplete();
519  return rbmat;
520 }
521 
522 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
523 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) {
524  TEUCHOS_ASSERT(bmat->getRangeMapExtractor()->getThyraMode() == bmat->getDomainMapExtractor()->getThyraMode());
525  Teuchos::RCP<const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > rbmat = Teuchos::null;
526  if(bmat->getRangeMapExtractor()->getThyraMode() == false) {
527  rbmat = mergeSubBlocks(brm, brm, bmat);
528  } else {
529  rbmat = mergeSubBlocksThyra(brm, brm, bmat);
530  }
531 
532  // TAW, 6/7/2016: rbmat might be Teuchos::null for empty blocks!
533  return rbmat;
534 }
535 
536 } //namespace Xpetra
537 
538 #define XPETRA_REORDEREDBLOCKEDCRSMATRIX_SHORT
539 #endif /* XPETRA_REORDEREDBLOCKEDCRSMATRIX_HPP */
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::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. ...
Teuchos::RCP< const Xpetra::BlockedCrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > fullOp_
Xpetra namespace
const Teuchos::RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > getMap() const
Implements DistObject interface.
Xpetra utility class for common map-related routines.
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)
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)
std::string description() const
Return a simple one-line description of this object.
std::string viewLabel_t
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)
BlockedCrsMatrix(Teuchos::RCP< const MapExtractor > &rangeMaps, Teuchos::RCP< const MapExtractor > &domainMaps, size_t npr, Xpetra::ProfileType pftype=Xpetra::DynamicProfile)
Constructor.
RCP< const Map > getRangeMap() const
Returns the Map associated with the full range of this operator.
Concrete implementation of Xpetra::Matrix.
Teuchos::RCP< Matrix > getMatrix(size_t r, size_t c) const
return block (r,c)
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, Xpetra::ProfileType pftype=Xpetra::DynamicProfile)
Constructor.
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.
virtual size_t Rows() const
number of row blocks
virtual size_t Cols() const
number of column blocks
Xpetra-specific matrix class.
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.