Xpetra Version of the Day
Loading...
Searching...
No Matches
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 <KokkosCompat_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
65namespace 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
84
85 public:
86
88
89
91
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:
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
157 TEUCHOS_ASSERT(brm_ != Teuchos::null);
158 out << description() << ": " << brm_->toString() << std::endl;
159 fullVec_->describe(out,verbLevel);
160 }
161
163
164 private:
167
168
169};
170
171template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
174
175 // TODO distinguish between range and domain map extractor! provide MapExtractor as parameter!
177
178 // number of sub blocks
179 size_t numBlocks = brm->GetNumBlocks();
180
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
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
216template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
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()->getLocalNumElements(); 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
294template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
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));
362 rbvec->setMultiVector(i,Teuchos::rcp_const_cast<MultiVector>(subvec),true);
363 }
364
365 }
366 return rbvec;
367}
368
369template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
381
382template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
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 */
static const EVerbosityLevel verbLevel_default
bool is_null() const
virtual size_t getNumVectors() const
Number of columns in the multivector.
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.
std::string description() const
Return a simple one-line description of this object.
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.
ReorderedBlockedMultiVector(Teuchos::RCP< const BlockedMap > &rangeMap, Teuchos::RCP< const Xpetra::BlockReorderManager > brm, Teuchos::RCP< const Xpetra::BlockedMultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > bvec)
Constructor.
Teuchos::RCP< const Xpetra::BlockReorderManager > brm_
Teuchos::RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > mergeSubBlockMaps(Teuchos::RCP< const Xpetra::BlockReorderManager > brm)
#define TEUCHOS_ASSERT(assertion_test)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Xpetra namespace
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)
std::string viewLabel_t
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::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::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > buildReorderedBlockedMultiVector(Teuchos::RCP< const Xpetra::BlockReorderManager > brm, Teuchos::RCP< const Xpetra::BlockedMultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > bvec)