MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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"
53 #include "Xpetra_Vector.hpp"
54 #include "Xpetra_VectorFactory.hpp"
55 #include <Xpetra_Matrix.hpp>
57 #include <Xpetra_MapFactory.hpp>
58 #include <Xpetra_MapExtractor.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);
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 
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()->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<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()->getNodeElementList();
321  stridedRgMap = StridedMapFactory::Build(
322  originalTransferOp->getRangeMap()->lib(),
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(),
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(),
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(),
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(),
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(),
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)
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Specifies the data that this class needs, and the factories that generate that data.
iterator begin() const
GlobalOrdinal GO
T & get(const std::string &name, T def_value)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Timer to be used in factories. Similar to Monitor but with additional timers.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
size_type size() const
LocalOrdinal LO
bool isParameter(const std::string &name) const
Print statistics that do not involve significant additional computation.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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.
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
iterator end() const
Node NO
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()