Xpetra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Xpetra_ReorderedBlockedMultiVector.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_REORDEREDBLOCKEDMULTIVECTOR_HPP
47 #define XPETRA_REORDEREDBLOCKEDMULTIVECTOR_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 
57 #include "Xpetra_BlockedMap.hpp"
58 #include "Xpetra_BlockedMultiVector.hpp"
59 
60 
65 namespace Xpetra {
66 
67  typedef std::string viewLabel_t;
68 
69  template <class Scalar,
70  class LocalOrdinal,
71  class GlobalOrdinal,
74  public BlockedMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> {
75  public:
76  typedef Scalar scalar_type;
77  typedef LocalOrdinal local_ordinal_type;
78  typedef GlobalOrdinal global_ordinal_type;
79  typedef Node node_type;
80 
81  private:
82 #undef XPETRA_REORDEREDBLOCKEDMULTIVECTOR_SHORT
83 #include "Xpetra_UseShortNames.hpp"
84 
85  public:
86 
88 
89 
91 
99  (Teuchos::RCP<const BlockedMap>& rangeMap,
100  Teuchos::RCP<const Xpetra::BlockReorderManager> brm,
102  : Xpetra::BlockedMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(rangeMap,bvec->getNumVectors(), false) {
103  brm_ = brm;
104  fullVec_ = bvec;
105  }
106 
107  //protected:
108 
111  brm_ = Teuchos::null;
112  fullVec_ = Teuchos::null;
113  }
114 
116 
117  private:
118  Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > mergeSubBlockMaps(Teuchos::RCP<const Xpetra::BlockReorderManager> brm) {
119  RCP<const BlockedMap> bmap = fullVec_->getBlockedMap();
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 = bmap->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  std::string description() const { return "ReorderedBlockedMultiVector"; }
154 
156  void describe(Teuchos::FancyOStream& out, const Teuchos::EVerbosityLevel verbLevel = Teuchos::Describable::verbLevel_default) const {
157  TEUCHOS_ASSERT(brm_ != Teuchos::null);
158  out << description() << ": " << brm_->toString() << std::endl;
159  fullVec_->describe(out,verbLevel);
160  }
161 
163 
164  private:
165  Teuchos::RCP<const Xpetra::BlockReorderManager > brm_;
166  Teuchos::RCP<const Xpetra::BlockedMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > fullVec_;
167 
168 
169 };
170 
171 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
172 Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > mergeSubBlockMaps(Teuchos::RCP<const Xpetra::BlockReorderManager> brm, Teuchos::RCP<const Xpetra::BlockedMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > bvec, bool bThyraMode) {
174 
175  // TODO distinguish between range and domain map extractor! provide MapExtractor as parameter!
176  RCP<const Xpetra::BlockedMap<LocalOrdinal,GlobalOrdinal,Node> > bmap = bvec->getBlockedMap();
177 
178  // number of sub blocks
179  size_t numBlocks = brm->GetNumBlocks();
180 
181  Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node>> map = Teuchos::null;
182 
183  if(numBlocks == 0) {
184  // it is a leaf node
185  Teuchos::RCP<const Xpetra::BlockReorderLeaf> leaf = Teuchos::rcp_dynamic_cast<const Xpetra::BlockReorderLeaf>(brm);
186 
187  map = bmap->getMap(Teuchos::as<size_t>(leaf->GetIndex()), bThyraMode);
188  } else {
189  // initialize vector for sub maps
190  std::vector<Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > > subMaps (numBlocks, Teuchos::null);
191 
192  for(size_t i = 0; i < numBlocks; i++) {
193  Teuchos::RCP<const Xpetra::BlockReorderManager> blkMgr = brm->GetBlock(Teuchos::as<int>(i));
194  subMaps[i] = mergeSubBlockMaps(blkMgr,bvec,bThyraMode);
195  TEUCHOS_ASSERT(subMaps[i].is_null()==false);
196  }
197 #if 1
198  // concatenate submaps
199  // for Thyra mode this map isn't important
200  Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > fullMap = MapUtils::concatenateMaps(subMaps);
201 
202  // create new BlockedMap (either in Thyra Mode or Xpetra mode)
203  map = Teuchos::rcp(new Xpetra::BlockedMap<LocalOrdinal,GlobalOrdinal,Node>(fullMap, subMaps, bThyraMode));
204 #else
205  // TAW: 11/27/16 we just concatenate the submaps to one monolithic Map object.
206  // Alternatively, we could create a new BlockedMap using the concatenated map and the submaps
207  // However, the block smoothers only need the concatenated map for creating MultiVectors...
208  // But for the Thyra mode version concatenating would not be ok for the whole map
209  map = MapUtils::concatenateMaps(subMaps);
210 #endif
211  }
212  TEUCHOS_ASSERT(map.is_null()==false);
213  return map;
214 }
215 
216 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
217 Teuchos::RCP<const Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > mergeSubBlocks(Teuchos::RCP<const Xpetra::BlockReorderManager> rowMgr, Teuchos::RCP<const Xpetra::BlockedMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > bvec) {
218 
225 
226  // number of sub blocks
227  size_t rowSz = rowMgr->GetNumBlocks();
228 
229  Teuchos::RCP<BlockedMultiVector> rbvec = Teuchos::null;
230 
231  if(rowSz == 0) {
232 
233  // it is a leaf node
234  Teuchos::RCP<const Xpetra::BlockReorderLeaf> rowleaf = Teuchos::rcp_dynamic_cast<const Xpetra::BlockReorderLeaf>(rowMgr);
235 
236  // extract leaf node
237  Teuchos::RCP<MultiVector> vec = bvec->getMultiVector(rowleaf->GetIndex(),false);
238 
239  TEUCHOS_ASSERT(vec != Teuchos::null);
240 
241  // check, whether leaf node is of type Xpetra::CrsMatrixWrap
242  Teuchos::RCP<BlockedMultiVector> subBVec = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(vec);
243  if(subBVec == Teuchos::null) {
244 
245  // DEBUG
246  /*{
247  std::cout << "MultiVector:" << std::endl;
248  Teuchos::ArrayRCP<const Scalar> vData = vec->getData(0);
249  for(size_t j=0; j< vec->getMap()->getNodeNumElements(); j++) {
250  std::cout << j << ": " << vec->getMap()->getGlobalElement(j) << ": " << vData[j] << std::endl;
251  }
252  }*/
253  // END DEBUG
254 
255  // If the leaf node is of type Xpetra::MultiVector. Wrap it into a ReorderedBlockMultiVector
256  // with the corresponding MapExtractors for translating Thyra to Xpetra GIDs if necessary
257  RCP<const BlockedMap> fullBlockedMap = bvec->getBlockedMap();
258  Teuchos::RCP<const Map> submap = fullBlockedMap->getMap(rowleaf->GetIndex(),false);
259  std::vector<Teuchos::RCP<const Map> > rowSubMaps (1, submap);
260  Teuchos::RCP<const BlockedMap> bbmap = Teuchos::rcp(new BlockedMap(submap, rowSubMaps, false));
261 
262  rbvec = Teuchos::rcp(new ReorderedBlockedMultiVector(bbmap, rowMgr, bvec));
263  rbvec->setMultiVector(0,Teuchos::rcp_const_cast<MultiVector>(vec),false);
264 
265  } else {
266  // If leaf node is already wrapped into a blocked matrix do not wrap it again.
267  rbvec = subBVec;
268  TEUCHOS_ASSERT(rbvec != Teuchos::null);
269  }
270  } else {
271  // create the map extractors
272  // we cannot create block matrix in thyra mode since merged maps might not start with 0 GID
273  Teuchos::RCP<const BlockedMap> rgBlockedMap = Teuchos::null;
274  std::vector<Teuchos::RCP<const Map> > rowSubMaps (rowSz, Teuchos::null);
275  for(size_t i = 0; i < rowSz; i++) {
276  Teuchos::RCP<const Xpetra::BlockReorderManager> rowSubMgr = rowMgr->GetBlock(Teuchos::as<int>(i));
277  rowSubMaps[i] = mergeSubBlockMaps(rowSubMgr,bvec,false /*xpetra*/);
278  TEUCHOS_ASSERT(rowSubMaps[i].is_null()==false);
279  }
280  Teuchos::RCP<const Map> rgMergedSubMaps = MapUtils::concatenateMaps(rowSubMaps);
281  rgBlockedMap = Teuchos::rcp(new BlockedMap(rgMergedSubMaps, rowSubMaps, false));
282  rbvec = Teuchos::rcp(new ReorderedBlockedMultiVector(rgBlockedMap, rowMgr,bvec));
283 
284  for(size_t i = 0; i < rowSz; i++) {
285  Teuchos::RCP<const Xpetra::BlockReorderManager> rowSubMgr = rowMgr->GetBlock(Teuchos::as<int>(i));
286  Teuchos::RCP<const MultiVector> subvec = mergeSubBlocks(rowSubMgr, bvec);
287  rbvec->setMultiVector(i,Teuchos::rcp_const_cast<MultiVector>(subvec),false);
288  }
289 
290  }
291  return rbvec;
292 }
293 
294 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
295 Teuchos::RCP<const Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > mergeSubBlocksThyra(Teuchos::RCP<const Xpetra::BlockReorderManager> rowMgr, Teuchos::RCP<const Xpetra::BlockedMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > bvec) {
301 
302  TEUCHOS_ASSERT(bvec->getBlockedMap()->getThyraMode() == true);
303 
304  // number of sub blocks
305  size_t rowSz = rowMgr->GetNumBlocks();
306 
307  Teuchos::RCP<BlockedMultiVector> rbvec = Teuchos::null;
308 
309  if(rowSz == 0) {
310  // it is a leaf node
311  Teuchos::RCP<const Xpetra::BlockReorderLeaf> rowleaf = Teuchos::rcp_dynamic_cast<const Xpetra::BlockReorderLeaf>(rowMgr);
312 
313  // this MultiVector uses Thyra style GIDs as global row indices
314  Teuchos::RCP<MultiVector> vec = bvec->getMultiVector(rowleaf->GetIndex(),true);
315 
316  TEUCHOS_ASSERT(vec.is_null() == false);
317 
318  // check, whether leaf node is of type Xpetra::CrsMatrixWrap
319  Teuchos::RCP<BlockedMultiVector> bbvec = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(vec);
320  if(bbvec == Teuchos::null) {
322  // build blocked map
323  RCP<const BlockedMap> fullBlockedRangeMap = bvec->getBlockedMap();
324  // extract Xpetra and Thyra based GIDs
325  Teuchos::RCP<const Map> xpsubmap = fullBlockedRangeMap->getMap(rowleaf->GetIndex(),false);
326  Teuchos::RCP<const Map> thysubmap = fullBlockedRangeMap->getMap(rowleaf->GetIndex(),true);
327  std::vector<Teuchos::RCP<const Map> > rowXpSubMaps (1, xpsubmap);
328  std::vector<Teuchos::RCP<const Map> > rowTySubMaps (1, thysubmap);
329  // use expert constructor
330  Teuchos::RCP<const BlockedMap> rgBlockedMap = Teuchos::rcp(new BlockedMap(rowXpSubMaps, rowTySubMaps));
331 
333  // build reordered blocked multi vector
334  rbvec = Teuchos::rcp(new ReorderedBlockedMultiVector(rgBlockedMap, rowMgr, bvec));
335  rbvec->setMultiVector(0,vec,true);
336  } else {
337  // If leaf node is already wrapped into a blocked matrix do not wrap it again.
338  rbvec = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(vec);
339  }
340  } else {
341  // create the blocked map
342  // we cannot create block multivector in thyra mode since merged maps might not start with 0 GID
343 
344  std::vector<Teuchos::RCP<const Map> > rowXpSubMaps (rowSz, Teuchos::null);
345  std::vector<Teuchos::RCP<const Map> > rowTySubMaps (rowSz, Teuchos::null);
346  for(size_t i = 0; i < rowSz; i++) {
347  Teuchos::RCP<const Xpetra::BlockReorderManager> rowSubMgr = rowMgr->GetBlock(Teuchos::as<int>(i));
348  // extract Xpetra and Thyra based merged GIDs
349  rowXpSubMaps[i] = mergeSubBlockMaps(rowSubMgr,bvec,false);
350  rowTySubMaps[i] = mergeSubBlockMaps(rowSubMgr,bvec,true);
351  TEUCHOS_ASSERT(rowXpSubMaps[i].is_null()==false);
352  TEUCHOS_ASSERT(rowTySubMaps[i].is_null()==false);
353  }
354  // use expert constructor
355  Teuchos::RCP<const BlockedMap> rgBlockedMap = Teuchos::rcp(new BlockedMap(rowXpSubMaps, rowTySubMaps));
356 
357  rbvec = Teuchos::rcp(new ReorderedBlockedMultiVector(rgBlockedMap, rowMgr,bvec));
358 
359  for(size_t i = 0; i < rowSz; i++) {
360  Teuchos::RCP<const Xpetra::BlockReorderManager> rowSubMgr = rowMgr->GetBlock(Teuchos::as<int>(i));
361  Teuchos::RCP<const MultiVector> subvec = mergeSubBlocksThyra(rowSubMgr, bvec);
362  rbvec->setMultiVector(i,Teuchos::rcp_const_cast<MultiVector>(subvec),true);
363  }
364 
365  }
366  return rbvec;
367 }
368 
369 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
370 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) {
371 
372  Teuchos::RCP<const Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > rbvec = Teuchos::null;
373  if(bvec->getBlockedMap()->getThyraMode() == false) {
374  rbvec = mergeSubBlocks(brm,bvec);
375  } else {
376  rbvec = mergeSubBlocksThyra(brm,bvec);
377  }
378  TEUCHOS_ASSERT(rbvec.is_null() == false);
379  return rbvec;
380 }
381 
382 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
383 Teuchos::RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > buildReorderedBlockedMultiVector(Teuchos::RCP<const Xpetra::BlockReorderManager> brm, Teuchos::RCP<Xpetra::BlockedMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > bvec) {
384 
387  Teuchos::RCP<const MultiVector> rbvec = Teuchos::null;
388  if(bvec->getBlockedMap()->getThyraMode() == false) {
389  rbvec = mergeSubBlocks(brm,Teuchos::rcp_const_cast<const BlockedMultiVector>(bvec));
390  } else {
391  rbvec = mergeSubBlocksThyra(brm,Teuchos::rcp_const_cast<const BlockedMultiVector>(bvec));
392  }
393  TEUCHOS_ASSERT(rbvec.is_null() == false);
394  return Teuchos::rcp_const_cast<MultiVector>(rbvec);
395 }
396 
397 
398 } //namespace Xpetra
399 
400 #define XPETRA_REORDEREDBLOCKEDMULTIVECTOR_SHORT
401 #endif /* XPETRA_REORDEREDBLOCKEDMULTIVECTOR_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)
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::BlockedMultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > fullVec_
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::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< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > mergeSubBlockMaps(Teuchos::RCP< const Xpetra::BlockReorderManager > brm)
ReorderedBlockedMultiVector(Teuchos::RCP< const BlockedMap > &rangeMap, Teuchos::RCP< const Xpetra::BlockReorderManager > brm, Teuchos::RCP< const Xpetra::BlockedMultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > bvec)
Constructor.
std::string viewLabel_t
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)
std::string description() const
Return a simple one-line description of this object.
Teuchos::RCP< const Xpetra::BlockReorderManager > brm_