All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups 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"
59 
60 
65 namespace Xpetra {
66 
67  typedef std::string viewLabel_t;
68 
69  template <class Scalar,
70  class LocalOrdinal,
71  class GlobalOrdinal,
72  class Node>
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 
103  : Xpetra::BlockedMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(rangeMap,bvec->getNumVectors(), false) {
104  brm_ = brm;
105  fullVec_ = bvec;
106  }
107 
108  //protected:
109 
112  brm_ = Teuchos::null;
113  fullVec_ = Teuchos::null;
114  }
115 
117 
118  private:
120  RCP<const BlockedMap> bmap = fullVec_->getBlockedMap();
121 
122  // number of sub blocks
123  size_t numBlocks = brm->GetNumBlocks();
124 
125  Teuchos::RCP<const Map> map = Teuchos::null;
126 
127  if(numBlocks == 0) {
128  // it is a leaf node
129  Teuchos::RCP<const Xpetra::BlockReorderLeaf> leaf = Teuchos::rcp_dynamic_cast<const Xpetra::BlockReorderLeaf>(brm);
130 
131  // never extract Thyra style maps (since we have to merge them)
132  map = bmap->getMap(Teuchos::as<size_t>(leaf->GetIndex()), false);
133  } else {
134  // initialize vector for sub maps
135  std::vector<Teuchos::RCP<const Map> > subMaps (numBlocks, Teuchos::null);
136 
137  for(size_t i = 0; i < numBlocks; i++) {
138  Teuchos::RCP<const Xpetra::BlockReorderManager> blkMgr = brm->GetBlock(Teuchos::as<int>(i));
139  subMaps[i] = mergeSubBlockMaps(blkMgr);
140  TEUCHOS_ASSERT(subMaps[i].is_null()==false);
141  }
142 
143  map = MapUtils::concatenateMaps(subMaps);
144  }
145  TEUCHOS_ASSERT(map.is_null()==false);
146  return map;
147  }
148 
149  public:
151 
152 
154  std::string description() const { return "ReorderedBlockedMultiVector"; }
155 
158  TEUCHOS_ASSERT(brm_ != Teuchos::null);
159  out << description() << ": " << brm_->toString() << std::endl;
160  fullVec_->describe(out,verbLevel);
161  }
162 
164 
165  private:
168 
169 
170 };
171 
172 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
175 
176  // TODO distinguish between range and domain map extractor! provide MapExtractor as parameter!
178 
179  // number of sub blocks
180  size_t numBlocks = brm->GetNumBlocks();
181 
183 
184  if(numBlocks == 0) {
185  // it is a leaf node
186  Teuchos::RCP<const Xpetra::BlockReorderLeaf> leaf = Teuchos::rcp_dynamic_cast<const Xpetra::BlockReorderLeaf>(brm);
187 
188  map = bmap->getMap(Teuchos::as<size_t>(leaf->GetIndex()), bThyraMode);
189  } else {
190  // initialize vector for sub maps
191  std::vector<Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > > subMaps (numBlocks, Teuchos::null);
192 
193  for(size_t i = 0; i < numBlocks; i++) {
194  Teuchos::RCP<const Xpetra::BlockReorderManager> blkMgr = brm->GetBlock(Teuchos::as<int>(i));
195  subMaps[i] = mergeSubBlockMaps(blkMgr,bvec,bThyraMode);
196  TEUCHOS_ASSERT(subMaps[i].is_null()==false);
197  }
198 #if 1
199  // concatenate submaps
200  // for Thyra mode this map isn't important
202 
203  // create new BlockedMap (either in Thyra Mode or Xpetra mode)
204  map = Teuchos::rcp(new Xpetra::BlockedMap<LocalOrdinal,GlobalOrdinal,Node>(fullMap, subMaps, bThyraMode));
205 #else
206  // TAW: 11/27/16 we just concatenate the submaps to one monolithic Map object.
207  // Alternatively, we could create a new BlockedMap using the concatenated map and the submaps
208  // However, the block smoothers only need the concatenated map for creating MultiVectors...
209  // But for the Thyra mode version concatenating would not be ok for the whole map
210  map = MapUtils::concatenateMaps(subMaps);
211 #endif
212  }
213  TEUCHOS_ASSERT(map.is_null()==false);
214  return map;
215 }
216 
217 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
219 
226 
227  // number of sub blocks
228  size_t rowSz = rowMgr->GetNumBlocks();
229 
230  Teuchos::RCP<BlockedMultiVector> rbvec = Teuchos::null;
231 
232  if(rowSz == 0) {
233 
234  // it is a leaf node
235  Teuchos::RCP<const Xpetra::BlockReorderLeaf> rowleaf = Teuchos::rcp_dynamic_cast<const Xpetra::BlockReorderLeaf>(rowMgr);
236 
237  // extract leaf node
238  Teuchos::RCP<MultiVector> vec = bvec->getMultiVector(rowleaf->GetIndex(),false);
239 
240  TEUCHOS_ASSERT(vec != Teuchos::null);
241 
242  // check, whether leaf node is of type Xpetra::CrsMatrixWrap
243  Teuchos::RCP<BlockedMultiVector> subBVec = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(vec);
244  if(subBVec == Teuchos::null) {
245 
246  // DEBUG
247  /*{
248  std::cout << "MultiVector:" << std::endl;
249  Teuchos::ArrayRCP<const Scalar> vData = vec->getData(0);
250  for(size_t j=0; j< vec->getMap()->getNodeNumElements(); j++) {
251  std::cout << j << ": " << vec->getMap()->getGlobalElement(j) << ": " << vData[j] << std::endl;
252  }
253  }*/
254  // END DEBUG
255 
256  // If the leaf node is of type Xpetra::MultiVector. Wrap it into a ReorderedBlockMultiVector
257  // with the corresponding MapExtractors for translating Thyra to Xpetra GIDs if necessary
258  RCP<const BlockedMap> fullBlockedMap = bvec->getBlockedMap();
259  Teuchos::RCP<const Map> submap = fullBlockedMap->getMap(rowleaf->GetIndex(),false);
260  std::vector<Teuchos::RCP<const Map> > rowSubMaps (1, submap);
261  Teuchos::RCP<const BlockedMap> bbmap = Teuchos::rcp(new BlockedMap(submap, rowSubMaps, false));
262 
263  rbvec = Teuchos::rcp(new ReorderedBlockedMultiVector(bbmap, rowMgr, bvec));
264  rbvec->setMultiVector(0,Teuchos::rcp_const_cast<MultiVector>(vec),false);
265 
266  } else {
267  // If leaf node is already wrapped into a blocked matrix do not wrap it again.
268  rbvec = subBVec;
269  TEUCHOS_ASSERT(rbvec != Teuchos::null);
270  }
271  } else {
272  // create the map extractors
273  // we cannot create block matrix in thyra mode since merged maps might not start with 0 GID
274  Teuchos::RCP<const BlockedMap> rgBlockedMap = Teuchos::null;
275  std::vector<Teuchos::RCP<const Map> > rowSubMaps (rowSz, Teuchos::null);
276  for(size_t i = 0; i < rowSz; i++) {
277  Teuchos::RCP<const Xpetra::BlockReorderManager> rowSubMgr = rowMgr->GetBlock(Teuchos::as<int>(i));
278  rowSubMaps[i] = mergeSubBlockMaps(rowSubMgr,bvec,false /*xpetra*/);
279  TEUCHOS_ASSERT(rowSubMaps[i].is_null()==false);
280  }
281  Teuchos::RCP<const Map> rgMergedSubMaps = MapUtils::concatenateMaps(rowSubMaps);
282  rgBlockedMap = Teuchos::rcp(new BlockedMap(rgMergedSubMaps, rowSubMaps, false));
283  rbvec = Teuchos::rcp(new ReorderedBlockedMultiVector(rgBlockedMap, rowMgr,bvec));
284 
285  for(size_t i = 0; i < rowSz; i++) {
286  Teuchos::RCP<const Xpetra::BlockReorderManager> rowSubMgr = rowMgr->GetBlock(Teuchos::as<int>(i));
287  Teuchos::RCP<const MultiVector> subvec = mergeSubBlocks(rowSubMgr, bvec);
288  rbvec->setMultiVector(i,Teuchos::rcp_const_cast<MultiVector>(subvec),false);
289  }
290 
291  }
292  return rbvec;
293 }
294 
295 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
302 
303  TEUCHOS_ASSERT(bvec->getBlockedMap()->getThyraMode() == true);
304 
305  // number of sub blocks
306  size_t rowSz = rowMgr->GetNumBlocks();
307 
308  Teuchos::RCP<BlockedMultiVector> rbvec = Teuchos::null;
309 
310  if(rowSz == 0) {
311  // it is a leaf node
312  Teuchos::RCP<const Xpetra::BlockReorderLeaf> rowleaf = Teuchos::rcp_dynamic_cast<const Xpetra::BlockReorderLeaf>(rowMgr);
313 
314  // this MultiVector uses Thyra style GIDs as global row indices
315  Teuchos::RCP<MultiVector> vec = bvec->getMultiVector(rowleaf->GetIndex(),true);
316 
317  TEUCHOS_ASSERT(vec.is_null() == false);
318 
319  // check, whether leaf node is of type Xpetra::CrsMatrixWrap
320  Teuchos::RCP<BlockedMultiVector> bbvec = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(vec);
321  if(bbvec == Teuchos::null) {
323  // build blocked map
324  RCP<const BlockedMap> fullBlockedRangeMap = bvec->getBlockedMap();
325  // extract Xpetra and Thyra based GIDs
326  Teuchos::RCP<const Map> xpsubmap = fullBlockedRangeMap->getMap(rowleaf->GetIndex(),false);
327  Teuchos::RCP<const Map> thysubmap = fullBlockedRangeMap->getMap(rowleaf->GetIndex(),true);
328  std::vector<Teuchos::RCP<const Map> > rowXpSubMaps (1, xpsubmap);
329  std::vector<Teuchos::RCP<const Map> > rowTySubMaps (1, thysubmap);
330  // use expert constructor
331  Teuchos::RCP<const BlockedMap> rgBlockedMap = Teuchos::rcp(new BlockedMap(rowXpSubMaps, rowTySubMaps));
332 
334  // build reordered blocked multi vector
335  rbvec = Teuchos::rcp(new ReorderedBlockedMultiVector(rgBlockedMap, rowMgr, bvec));
336  rbvec->setMultiVector(0,vec,true);
337  } else {
338  // If leaf node is already wrapped into a blocked matrix do not wrap it again.
339  rbvec = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(vec);
340  }
341  } else {
342  // create the blocked map
343  // we cannot create block multivector in thyra mode since merged maps might not start with 0 GID
344 
345  std::vector<Teuchos::RCP<const Map> > rowXpSubMaps (rowSz, Teuchos::null);
346  std::vector<Teuchos::RCP<const Map> > rowTySubMaps (rowSz, Teuchos::null);
347  for(size_t i = 0; i < rowSz; i++) {
348  Teuchos::RCP<const Xpetra::BlockReorderManager> rowSubMgr = rowMgr->GetBlock(Teuchos::as<int>(i));
349  // extract Xpetra and Thyra based merged GIDs
350  rowXpSubMaps[i] = mergeSubBlockMaps(rowSubMgr,bvec,false);
351  rowTySubMaps[i] = mergeSubBlockMaps(rowSubMgr,bvec,true);
352  TEUCHOS_ASSERT(rowXpSubMaps[i].is_null()==false);
353  TEUCHOS_ASSERT(rowTySubMaps[i].is_null()==false);
354  }
355  // use expert constructor
356  Teuchos::RCP<const BlockedMap> rgBlockedMap = Teuchos::rcp(new BlockedMap(rowXpSubMaps, rowTySubMaps));
357 
358  rbvec = Teuchos::rcp(new ReorderedBlockedMultiVector(rgBlockedMap, rowMgr,bvec));
359 
360  for(size_t i = 0; i < rowSz; i++) {
361  Teuchos::RCP<const Xpetra::BlockReorderManager> rowSubMgr = rowMgr->GetBlock(Teuchos::as<int>(i));
362  Teuchos::RCP<const MultiVector> subvec = mergeSubBlocksThyra(rowSubMgr, bvec);
363  rbvec->setMultiVector(i,Teuchos::rcp_const_cast<MultiVector>(subvec),true);
364  }
365 
366  }
367  return rbvec;
368 }
369 
370 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
372 
374  if(bvec->getBlockedMap()->getThyraMode() == false) {
375  rbvec = mergeSubBlocks(brm,bvec);
376  } else {
377  rbvec = mergeSubBlocksThyra(brm,bvec);
378  }
379  TEUCHOS_ASSERT(rbvec.is_null() == false);
380  return rbvec;
381 }
382 
383 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
385 
388  Teuchos::RCP<const MultiVector> rbvec = Teuchos::null;
389  if(bvec->getBlockedMap()->getThyraMode() == false) {
390  rbvec = mergeSubBlocks(brm,Teuchos::rcp_const_cast<const BlockedMultiVector>(bvec));
391  } else {
392  rbvec = mergeSubBlocksThyra(brm,Teuchos::rcp_const_cast<const BlockedMultiVector>(bvec));
393  }
394  TEUCHOS_ASSERT(rbvec.is_null() == false);
395  return Teuchos::rcp_const_cast<MultiVector>(rbvec);
396 }
397 
398 
399 } //namespace Xpetra
400 
401 #define XPETRA_REORDEREDBLOCKEDMULTIVECTOR_SHORT
402 #endif /* XPETRA_REORDEREDBLOCKEDMULTIVECTOR_HPP */
bool is_null(const boost::shared_ptr< T > &p)
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.
int GetIndex() const
Get the index that is stored in this block/leaf.
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 size_t GetNumBlocks() const
Returns number of subblocks.
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)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
virtual const Teuchos::RCP< BlockReorderManager > GetBlock(int blockIndex)
Get a particular block. If there is no block at this index location return a new one.
ReorderedBlockedMultiVector(Teuchos::RCP< const BlockedMap > &rangeMap, Teuchos::RCP< const Xpetra::BlockReorderManager > brm, Teuchos::RCP< const Xpetra::BlockedMultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > bvec)
Constructor.
void setMultiVector(size_t r, Teuchos::RCP< const MultiVector > v, bool bThyraMode)
set partial multivector associated with block row r
std::string viewLabel_t
static const EVerbosityLevel verbLevel_default
virtual RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > getMap() const
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)
#define TEUCHOS_ASSERT(assertion_test)
std::string description() const
Return a simple one-line description of this object.
Teuchos::RCP< const Xpetra::BlockReorderManager > brm_
bool is_null() const