MueLu Version of the Day
Loading...
Searching...
No Matches
MueLu_RebalanceBlockInterpolationFactory_def.hpp
Go to the documentation of this file.
1// @HEADER
2//
3// ***********************************************************************
4//
5// MueLu: A package for multigrid based preconditioning
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 MUELU_REBALANCEBLOCKINTERPOLATIONFACTORY_DEF_HPP_
47#define MUELU_REBALANCEBLOCKINTERPOLATIONFACTORY_DEF_HPP_
48
49#include <Teuchos_Tuple.hpp>
50
51#include "Xpetra_MultiVector.hpp"
52#include "Xpetra_MultiVectorFactory.hpp"
53#include "Xpetra_Vector.hpp"
54#include "Xpetra_VectorFactory.hpp"
55#include <Xpetra_Matrix.hpp>
56#include <Xpetra_BlockedCrsMatrix.hpp>
57#include <Xpetra_MapFactory.hpp>
58#include <Xpetra_MapExtractor.hpp>
59#include <Xpetra_MapExtractorFactory.hpp>
60#include <Xpetra_MatrixFactory.hpp>
61#include <Xpetra_Import.hpp>
62#include <Xpetra_ImportFactory.hpp>
63
65
67#include "MueLu_HierarchyUtils.hpp"
68#include "MueLu_Level.hpp"
69#include "MueLu_Monitor.hpp"
70#include "MueLu_PerfUtils.hpp"
71
72namespace MueLu {
73
74template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
76 RCP<ParameterList> validParamList = rcp(new ParameterList());
77
78 validParamList->set< RCP<const FactoryBase> >("P", Teuchos::null, "Factory of the prolongation operator that need to be rebalanced (only used if type=Interpolation)");
79 validParamList->set< RCP<const FactoryBase> >("A", Teuchos::null, "Factory for generating the non-rebalanced coarse level A. We need this to make sure the non-rebalanced coarse A is calculated first before rebalancing takes place.");
80
81 validParamList->set< RCP<const FactoryBase> >("Coordinates", Teuchos::null, "Factory for generating the non-rebalanced Coordinates.");
82 validParamList->set<RCP<const FactoryBase> >("Importer", Teuchos::null, "Generating factory of the matrix Importer for rebalancing");
83 validParamList->set<RCP<const FactoryBase> >("SubImporters", Teuchos::null, "Generating factory of the matrix sub-Importers for rebalancing");
84
85#define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
86 // SET_VALID_ENTRY("repartition: use subcommunicators");
87#undef SET_VALID_ENTRY
88
89 // TODO validation: "P" parameter valid only for type="Interpolation" and "R" valid only for type="Restriction". Like so:
90 // if (paramList.isEntry("type") && paramList.get("type) == "Interpolation) {
91 // validParamList->set< RCP<const FactoryBase> >("P", Teuchos::null, "Factory of the prolongation operator that need to be rebalanced (only used if type=Interpolation)");
92
93 return validParamList;
94}
95
96template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
98 FactManager_.push_back(FactManager);
99}
100
101template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
103 Input(coarseLevel, "P");
104 Input(coarseLevel, "A"); // we request the non-rebalanced coarse level A since we have to make sure it is calculated before rebalancing starts!
105
106 std::vector<Teuchos::RCP<const FactoryManagerBase> >::const_iterator it;
107 for(it = FactManager_.begin(); it!=FactManager_.end(); ++it) {
108 SetFactoryManager fineSFM (rcpFromRef(fineLevel), *it);
109 SetFactoryManager coarseSFM(rcpFromRef(coarseLevel), *it);
110
111 // Request Importer and Coordinates (if defined in xml file)
112 // Note, that we have to use the Level::DeclareInput routine in order to use the FactoryManager *it (rather than the main factory manager)
113 coarseLevel.DeclareInput("Importer",(*it)->GetFactory("Importer").get(), this);
114 if((*it)->hasFactory("Coordinates") == true)
115 coarseLevel.DeclareInput("Coordinates",(*it)->GetFactory("Coordinates").get(), this);
116 }
117
118 // Use the non-manager path if the maps / importers are generated in one place
119 if(FactManager_.size() == 0) {
120 Input(coarseLevel,"Importer");
121 Input(coarseLevel,"SubImporters");
122 Input(coarseLevel,"Coordinates");
123 }
124
125
126}
127
128template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
130 FactoryMonitor m(*this, "Build", coarseLevel);
131 typedef Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType, LO, GO, NO> xdMV;
132 typedef Xpetra::BlockedMultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO> xdBV;
133
134 bool UseSingleSource = FactManager_.size() == 0;
135 //RCP<Teuchos::FancyOStream> out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
136
137 Teuchos::RCP<Matrix> nonrebCoarseA = Get< RCP<Matrix> >(coarseLevel, "A");
138
139 Teuchos::RCP<Matrix> originalTransferOp = Teuchos::null;
140 originalTransferOp = Get< RCP<Matrix> >(coarseLevel, "P");
141
142 RCP<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > bOriginalTransferOp =
143 Teuchos::rcp_dynamic_cast<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > (originalTransferOp);
144 TEUCHOS_TEST_FOR_EXCEPTION(bOriginalTransferOp==Teuchos::null, Exceptions::BadCast, "MueLu::RebalanceBlockTransferFactory::Build: input matrix P or R is not of type BlockedCrsMatrix! error.");
145
146 RCP<const MapExtractor> rangeMapExtractor = bOriginalTransferOp->getRangeMapExtractor();
147 RCP<const MapExtractor> domainMapExtractor = bOriginalTransferOp->getDomainMapExtractor();
148
149 // check if GIDs for full maps have to be sorted:
150 // For the Thyra mode ordering they do not have to be sorted since the GIDs are
151 // numbered as 0...n1,0...,n2 (starting with zero for each subblock). The MapExtractor
152 // generates unique GIDs during the construction.
153 // For Xpetra style, the GIDs have to be reordered. Such that one obtains a ordered
154 // list of GIDs in an increasing ordering. In Xpetra, the GIDs are all unique through
155 // out all submaps.
156 bool bThyraRangeGIDs = rangeMapExtractor->getThyraMode();
157 bool bThyraDomainGIDs = domainMapExtractor->getThyraMode();
158
159 // declare variables for maps of blocked rebalanced prolongation operator
160 std::vector<GO> fullRangeMapVector; // contains all range GIDs on current processor
161 std::vector<GO> fullDomainMapVector; // contains all domain GIDs on current processor
162 std::vector<RCP<const Map> > subBlockPRangeMaps;
163 std::vector<RCP<const Map> > subBlockPDomainMaps;
164 subBlockPRangeMaps.reserve(bOriginalTransferOp->Rows()); // reserve size for block P operators
165 subBlockPDomainMaps.reserve(bOriginalTransferOp->Cols()); // reserve size for block P operators
166
167 std::vector<Teuchos::RCP<Matrix> > subBlockRebP;
168 subBlockRebP.reserve(bOriginalTransferOp->Rows());
169
170 // For use in single-source mode only
171 std::vector<RCP<const Import> > importers = std::vector<RCP<const Import> >(bOriginalTransferOp->Rows(), Teuchos::null);
172 std::vector<RCP<xdMV> > newCoordinates(bOriginalTransferOp->Rows());
173 RCP<xdBV> oldCoordinates;
174 if(UseSingleSource) {
175 importers = Get<std::vector<RCP<const Import> > >(coarseLevel,"SubImporters");
176 oldCoordinates = Get<RCP<xdBV> >(coarseLevel,"Coordinates");
177 }
178
179 int curBlockId = 0;
180 Teuchos::RCP<const Import> rebalanceImporter = Teuchos::null;
181 std::vector<Teuchos::RCP<const FactoryManagerBase> >::const_iterator it;
182 for (it = FactManager_.begin(); it != FactManager_.end(); ++it) {
183 // begin SubFactoryManager environment
184 SetFactoryManager fineSFM (rcpFromRef(fineLevel), *it);
185 SetFactoryManager coarseSFM(rcpFromRef(coarseLevel), *it);
186
187 // TAW: use the Level::Get routine in order to access the data declared in (*it) factory manager (rather than the main factory manager)
188 if(UseSingleSource) rebalanceImporter = importers[curBlockId];
189 else rebalanceImporter = coarseLevel.Get<Teuchos::RCP<const Import> >("Importer", (*it)->GetFactory("Importer").get());
190
191 // extract diagonal matrix block
192 Teuchos::RCP<Matrix> Pii = bOriginalTransferOp->getMatrix(curBlockId, curBlockId);
193 Teuchos::RCP<CrsMatrixWrap> Pwii = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(Pii);
194 TEUCHOS_TEST_FOR_EXCEPTION(Pwii == Teuchos::null,Xpetra::Exceptions::BadCast, "MueLu::RebalanceBlockTransferFactory::Build: block " << curBlockId << " is not of type CrsMatrixWrap. We need an underlying CsrMatrix to replace domain map and importer!");
195
196 MUELU_TEST_FOR_EXCEPTION(bThyraRangeGIDs == true && Pii->getRowMap()->getMinAllGlobalIndex() != 0,
198 "MueLu::RebalanceBlockInterpolationFactory::Build: inconsistent Thyra GIDs. Thyra global ids for block range " << curBlockId << " start with " << Pii->getRowMap()->getMinAllGlobalIndex() << " but should start with 0");
199 MUELU_TEST_FOR_EXCEPTION(bThyraDomainGIDs == true && Pii->getColMap()->getMinAllGlobalIndex() != 0,
201 "MueLu::RebalanceBlockInterpolationFactory::Build: inconsistent Thyra GIDs. Thyra global ids for block domain " << curBlockId << " start with " << Pii->getColMap()->getMinAllGlobalIndex() << " but should start with 0");
202
203 // rebalance P11
204 if(rebalanceImporter != Teuchos::null) {
205 std::stringstream ss; ss << "Rebalancing prolongator block P(" << curBlockId << "," << curBlockId << ")";
206 SubFactoryMonitor m1(*this, ss.str(), coarseLevel);
207
208 // P is the transfer operator from the coarse grid to the fine grid.
209 // P must transfer the data from the newly reordered coarse A to the (unchanged) fine A.
210 // This means that the domain map (coarse) of P must be changed according to the new partition. The range map (fine) is kept unchanged.
211 //
212 // The domain map of P must match the range map of R.
213 // See also note below about domain/range map of R and its implications for P.
214 //
215 // To change the domain map of P, P needs to be fillCompleted again with the new domain map.
216 // To achieve this, P is copied into a new matrix that is not fill-completed.
217 // The doImport() operation is just used here to make a copy of P: the importer is trivial and there is no data movement involved.
218 // The reordering actually happens during the fillComplete() with domainMap == rebalanceImporter->getTargetMap().
219 RCP<const Import> newImporter;
220 {
221 SubFactoryMonitor subM(*this, "Rebalancing prolongator -- fast map replacement", coarseLevel);
222 newImporter = ImportFactory::Build(rebalanceImporter->getTargetMap(), Pii->getColMap());
223 Pwii->getCrsMatrix()->replaceDomainMapAndImporter(rebalanceImporter->getTargetMap(), newImporter);
224 }
225
226 RCP<ParameterList> params = rcp(new ParameterList());
227 params->set("printLoadBalancingInfo", true);
228 std::stringstream ss2; ss2 << "P(" << curBlockId << "," << curBlockId << ") rebalanced:";
229 GetOStream(Statistics0) << PerfUtils::PrintMatrixInfo(*Pii, ss2.str(), params);
230
231 // store rebalanced P block
232 subBlockRebP.push_back(Pii);
233
234 // rebalance coordinates
235 // TAW: Note, that each sub-block manager overwrites the Coordinates. So far we only support one set of Coordinates
236 // for a multiphysics problem (i.e., we only support volume coupled problems with the same mesh)
237 if(UseSingleSource) {
238 RCP<xdMV> localCoords = oldCoordinates->getMultiVector(curBlockId);
239
240 // FIXME: This should be extended to work with blocking
241 RCP<const Import> coordImporter = rebalanceImporter;
242
243 RCP<xdMV> permutedLocalCoords = Xpetra::MultiVectorFactory<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO>::Build(coordImporter->getTargetMap(), localCoords->getNumVectors());
244 permutedLocalCoords->doImport(*localCoords, *coordImporter, Xpetra::INSERT);
245
246 newCoordinates[curBlockId] = permutedLocalCoords;
247 }
248 else if ( (*it)->hasFactory("Coordinates") == true && coarseLevel.IsAvailable("Coordinates",(*it)->GetFactory("Coordinates").get()) == true) {
249 RCP<xdMV> coords = coarseLevel.Get< RCP<xdMV> >("Coordinates",(*it)->GetFactory("Coordinates").get());
250
251 // This line must be after the Get call
252 SubFactoryMonitor subM(*this, "Rebalancing coordinates", coarseLevel);
253
254 LO nodeNumElts = coords->getMap()->getLocalNumElements();
255
256 // If a process has no matrix rows, then we can't calculate blocksize using the formula below.
257 LO myBlkSize = 0, blkSize = 0;
258
259 if (nodeNumElts > 0) {
260 MUELU_TEST_FOR_EXCEPTION(rebalanceImporter->getSourceMap()->getLocalNumElements() % nodeNumElts != 0,
262 "MueLu::RebalanceBlockInterpolationFactory::Build: block size. " << rebalanceImporter->getSourceMap()->getLocalNumElements() << " not divisable by " << nodeNumElts);
263 myBlkSize = rebalanceImporter->getSourceMap()->getLocalNumElements() / nodeNumElts;
264 }
265
266 MueLu_maxAll(coords->getMap()->getComm(), myBlkSize, blkSize);
267
268 RCP<const Import> coordImporter = Teuchos::null;
269 if (blkSize == 1) {
270 coordImporter = rebalanceImporter;
271 } else {
272 // NOTE: there is an implicit assumption here: we assume that dof any node are enumerated consequently
273 // Proper fix would require using decomposition similar to how we construct importer in the
274 // RepartitionFactory
275 RCP<const Map> origMap = coords->getMap();
276 GO indexBase = origMap->getIndexBase();
277
278 ArrayView<const GO> OEntries = rebalanceImporter->getTargetMap()->getLocalElementList();
279 LO numEntries = OEntries.size()/blkSize;
280 ArrayRCP<GO> Entries(numEntries);
281 for (LO i = 0; i < numEntries; i++)
282 Entries[i] = (OEntries[i*blkSize]-indexBase)/blkSize + indexBase;
283
284 RCP<const Map> targetMap = MapFactory::Build(origMap->lib(), origMap->getGlobalNumElements(), Entries(), indexBase, origMap->getComm());
285 coordImporter = ImportFactory::Build(origMap, targetMap);
286 }
287
288 RCP<xdMV> permutedCoords = Xpetra::MultiVectorFactory<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO>::Build(coordImporter->getTargetMap(), coords->getNumVectors());
289 permutedCoords->doImport(*coords, *coordImporter, Xpetra::INSERT);
290
291 const ParameterList& pL = GetParameterList();
292 if (pL.isParameter("repartition: use subcommunicators") == true && pL.get<bool>("repartition: use subcommunicators") == true)
293 permutedCoords->replaceMap(permutedCoords->getMap()->removeEmptyProcesses());
294
295 Set(coarseLevel, "Coordinates", permutedCoords);
296 }
297 } // end rebalance P(1,1)
298 else {
299 RCP<ParameterList> params = rcp(new ParameterList());
300 params->set("printLoadBalancingInfo", true);
301 std::stringstream ss; ss << "P(" << curBlockId << "," << curBlockId << ") not rebalanced:";
302 GetOStream(Statistics0) << PerfUtils::PrintMatrixInfo(*Pii, ss.str(), params);
303 // store rebalanced P block
304 subBlockRebP.push_back(Pii);
305
306 // Store Coordinates on coarse level (generated by this)
307 // TAW: Note, that each sub-block manager overwrites the Coordinates. So far we only support one set of Coordinates
308 // for a multiphysics problem (i.e., we only support volume coupled problems with the same mesh)
309 if((*it)->hasFactory("Coordinates") == true && coarseLevel.IsAvailable("Coordinates",(*it)->GetFactory("Coordinates").get()) == true) {
310 coarseLevel.Set("Coordinates", coarseLevel.Get< RCP<xdMV> >("Coordinates",(*it)->GetFactory("Coordinates").get()),this);
311 }
312 }
313
314 // fix striding information for rebalanced diagonal block Pii
315 //RCP<const Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node> > rgPMapExtractor = bOriginalTransferOp->getRangeMapExtractor(); // original map extractor
316 Teuchos::RCP<const StridedMap> orig_stridedRgMap = Teuchos::rcp_dynamic_cast<const StridedMap>(rangeMapExtractor->getMap(Teuchos::as<size_t>(curBlockId),rangeMapExtractor->getThyraMode()));
317 Teuchos::RCP<const Map> stridedRgMap = Teuchos::null;
318 if(orig_stridedRgMap != Teuchos::null) {
319 std::vector<size_t> stridingData = orig_stridedRgMap->getStridingData();
320 Teuchos::ArrayView< const GlobalOrdinal > nodeRangeMapii = Pii->getRangeMap()->getLocalElementList();
321 stridedRgMap = StridedMapFactory::Build(
322 originalTransferOp->getRangeMap()->lib(),
323 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
324 nodeRangeMapii,
325 Pii->getRangeMap()->getIndexBase(),
326 stridingData,
327 originalTransferOp->getRangeMap()->getComm(),
328 orig_stridedRgMap->getStridedBlockId(),
329 orig_stridedRgMap->getOffset());
330 } else stridedRgMap = Pii->getRangeMap();
331
332 //RCP<const Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node> > doPMapExtractor = bOriginalTransferOp->getDomainMapExtractor(); // original map extractor
333 Teuchos::RCP<const StridedMap> orig_stridedDoMap = Teuchos::rcp_dynamic_cast<const StridedMap>(domainMapExtractor->getMap(Teuchos::as<size_t>(curBlockId),domainMapExtractor->getThyraMode()));
334
335 Teuchos::RCP<const Map> stridedDoMap = Teuchos::null;
336 if(orig_stridedDoMap != Teuchos::null) {
337 std::vector<size_t> stridingData = orig_stridedDoMap->getStridingData();
338 Teuchos::ArrayView< const GlobalOrdinal > nodeDomainMapii = Pii->getDomainMap()->getLocalElementList();
339 stridedDoMap = StridedMapFactory::Build(originalTransferOp->getDomainMap()->lib(),
340 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
341 nodeDomainMapii,
342 Pii->getDomainMap()->getIndexBase(),
343 stridingData,
344 originalTransferOp->getDomainMap()->getComm(),
345 orig_stridedDoMap->getStridedBlockId(),
346 orig_stridedDoMap->getOffset());
347 } else stridedDoMap = Pii->getDomainMap();
348
349 TEUCHOS_TEST_FOR_EXCEPTION(stridedRgMap == Teuchos::null,Exceptions::RuntimeError, "MueLu::RebalanceBlockInterpolationFactory::Build: failed to generate striding information. error.");
350 TEUCHOS_TEST_FOR_EXCEPTION(stridedDoMap == Teuchos::null,Exceptions::RuntimeError, "MueLu::RebalanceBlockInterpolationFactory::Build: failed to generate striding information. error.");
351
352 // replace stridedMaps view in diagonal sub block
353 if(Pii->IsView("stridedMaps")) Pii->RemoveView("stridedMaps");
354 Pii->CreateView("stridedMaps", stridedRgMap, stridedDoMap);
355
356 // append strided row map (= range map) to list of range maps.
357 Teuchos::RCP<const Map> rangeMapii = Pii->getRowMap("stridedMaps"); // strided range map
358 subBlockPRangeMaps.push_back(rangeMapii);
359 Teuchos::ArrayView< const GlobalOrdinal > nodeRangeMapii = Pii->getRangeMap()->getLocalElementList();
360 // append the GIDs in the end. Do not sort if we have Thyra style GIDs
361 fullRangeMapVector.insert(fullRangeMapVector.end(), nodeRangeMapii.begin(), nodeRangeMapii.end());
362 if(bThyraRangeGIDs == false)
363 sort(fullRangeMapVector.begin(), fullRangeMapVector.end());
364
365 // append strided col map (= domain map) to list of range maps.
366 Teuchos::RCP<const Map> domainMapii = Pii->getColMap("stridedMaps"); // strided domain map
367 subBlockPDomainMaps.push_back(domainMapii);
368 Teuchos::ArrayView< const GlobalOrdinal > nodeDomainMapii = Pii->getDomainMap()->getLocalElementList();
369 // append the GIDs in the end. Do not sort if we have Thyra style GIDs
370 fullDomainMapVector.insert(fullDomainMapVector.end(), nodeDomainMapii.begin(), nodeDomainMapii.end());
371 if(bThyraDomainGIDs == false)
372 sort(fullDomainMapVector.begin(), fullDomainMapVector.end());
373
374 curBlockId++; // increase block id index
375
376 } // end SubFactoryManager environment
377
378 // extract map index base from maps of blocked P
379 GO rangeIndexBase = originalTransferOp->getRangeMap()->getIndexBase();
380 GO domainIndexBase= originalTransferOp->getDomainMap()->getIndexBase();
381
382 RCP<const Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node> > rangePMapExtractor = bOriginalTransferOp->getRangeMapExtractor(); // original map extractor
383 Teuchos::ArrayView<GO> fullRangeMapGIDs(fullRangeMapVector.size() ? &fullRangeMapVector[0] : 0,fullRangeMapVector.size());
384 Teuchos::RCP<const StridedMap> stridedRgFullMap = Teuchos::rcp_dynamic_cast<const StridedMap>(rangePMapExtractor->getFullMap());
385 Teuchos::RCP<const Map > fullRangeMap = Teuchos::null;
386 if(stridedRgFullMap != Teuchos::null) {
387 std::vector<size_t> stridedData = stridedRgFullMap->getStridingData();
388 fullRangeMap =
389 StridedMapFactory::Build(
390 originalTransferOp->getRangeMap()->lib(),
391 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
392 fullRangeMapGIDs,
393 rangeIndexBase,
394 stridedData,
395 originalTransferOp->getRangeMap()->getComm(),
396 stridedRgFullMap->getStridedBlockId(),
397 stridedRgFullMap->getOffset());
398 } else {
399 fullRangeMap =
400 MapFactory::Build(
401 originalTransferOp->getRangeMap()->lib(),
402 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
403 fullRangeMapGIDs,
404 rangeIndexBase,
405 originalTransferOp->getRangeMap()->getComm());
406 }
407
408 RCP<const Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node> > domainAMapExtractor = bOriginalTransferOp->getDomainMapExtractor();
409 Teuchos::ArrayView<GO> fullDomainMapGIDs(fullDomainMapVector.size() ? &fullDomainMapVector[0] : 0,fullDomainMapVector.size());
410 Teuchos::RCP<const StridedMap> stridedDoFullMap = Teuchos::rcp_dynamic_cast<const StridedMap>(domainAMapExtractor->getFullMap());
411 Teuchos::RCP<const Map > fullDomainMap = Teuchos::null;
412 if(stridedDoFullMap != Teuchos::null) {
413 TEUCHOS_TEST_FOR_EXCEPTION(stridedDoFullMap==Teuchos::null, Exceptions::BadCast, "MueLu::BlockedPFactory::Build: full map in domain map extractor has no striding information! error.");
414 std::vector<size_t> stridedData2 = stridedDoFullMap->getStridingData();
415 fullDomainMap =
416 StridedMapFactory::Build(
417 originalTransferOp->getDomainMap()->lib(),
418 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
419 fullDomainMapGIDs,
420 domainIndexBase,
421 stridedData2,
422 originalTransferOp->getDomainMap()->getComm(),
423 stridedDoFullMap->getStridedBlockId(),
424 stridedDoFullMap->getOffset());
425 } else {
426
427 fullDomainMap =
428 MapFactory::Build(
429 originalTransferOp->getDomainMap()->lib(),
430 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
431 fullDomainMapGIDs,
432 domainIndexBase,
433 originalTransferOp->getDomainMap()->getComm());
434 }
435
436 // build map extractors
437 Teuchos::RCP<const MapExtractor> rebrangeMapExtractor =
438 MapExtractorFactory::Build(fullRangeMap, subBlockPRangeMaps, bThyraRangeGIDs);
439 Teuchos::RCP<const MapExtractor> rebdomainMapExtractor =
440 MapExtractorFactory::Build(fullDomainMap, subBlockPDomainMaps, bThyraDomainGIDs);
441
442 Teuchos::RCP<BlockedCrsMatrix> bRebP = Teuchos::rcp(new BlockedCrsMatrix(rebrangeMapExtractor,rebdomainMapExtractor,10));
443 for(size_t i = 0; i<subBlockPRangeMaps.size(); i++) {
444 Teuchos::RCP<CrsMatrixWrap> crsOpii = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(subBlockRebP[i]);
445 TEUCHOS_TEST_FOR_EXCEPTION(crsOpii == Teuchos::null,Xpetra::Exceptions::BadCast, "MueLu::RebalanceBlockTransferFactory::Build: block P" << i << " is not of type CrsMatrixWrap.");
446 bRebP->setMatrix(i,i,crsOpii);
447 }
448 bRebP->fillComplete();
449
450 Set(coarseLevel, "P", Teuchos::rcp_dynamic_cast<Matrix>(bRebP));
451
452
453 // Finish up the coordinates (single source only)
454 if(UseSingleSource) {
455 RCP<xdBV> bcoarseCoords = rcp(new xdBV(rebrangeMapExtractor->getBlockedMap(),newCoordinates));
456 Set(coarseLevel,"Coordinates",bcoarseCoords);
457 }
458
459
460} // Build
461
462} // namespace MueLu
463
464#endif /* MUELU_REBALANCEBLOCKINTERPOLATIONFACTORY_DEF_HPP_ */
#define MUELU_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
#define MueLu_maxAll(rcpComm, in, out)
Exception indicating invalid cast attempted.
Exception throws to report errors in the internal logical of the program.
Timer to be used in factories. Similar to Monitor but with additional timers.
Class that holds all level-specific information.
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need's value has been saved.
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access)....
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Specifies the data that this class needs, and the factories that generate that data.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
void AddFactoryManager(RCP< const FactoryManagerBase > FactManager)
Add a factory manager.
An exception safe way to call the method 'Level::SetFactoryManager()'.
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
Namespace for MueLu classes and methods.
@ Statistics0
Print statistics that do not involve significant additional computation.