MueLu  Version of the Day
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 
72 namespace MueLu {
73 
74 template <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 
96 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
98  FactManager_.push_back(FactManager);
99 }
100 
101 template <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 
128 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
130  FactoryMonitor m(*this, "Build", coarseLevel);
131  typedef Xpetra::MultiVector<double, LO, GO, NO> xdMV;
132  typedef Xpetra::BlockedMultiVector<double,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<double,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()->getNodeNumElements();
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()->getNodeNumElements() % nodeNumElts != 0,
262  "MueLu::RebalanceBlockInterpolationFactory::Build: block size. " << rebalanceImporter->getSourceMap()->getNodeNumElements() << " not divisable by " << nodeNumElts);
263  myBlkSize = rebalanceImporter->getSourceMap()->getNodeNumElements() / 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()->getNodeElementList();
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<double,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()->getNodeElementList();
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()->getNodeElementList();
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()->getNodeElementList();
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()->getNodeElementList();
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_ */
Exception indicating invalid cast attempted.
#define MueLu_maxAll(rcpComm, in, out)
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
Timer to be used in factories. Similar to Monitor but with additional timers.
Namespace for MueLu classes and methods.
Print statistics that do not involve significant additional computation.
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
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.
Exception throws to report errors in the internal logical of the program.
#define MUELU_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
An exception safe way to call the method &#39;Level::SetFactoryManager()&#39;.
void AddFactoryManager(RCP< const FactoryManagerBase > FactManager)
Add a factory manager.
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()