MueLu  Version of the Day
MueLu_RebalanceBlockRestrictionFactory_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_REBALANCEBLOCKRESTRICTIONFACTORY_DEF_HPP_
47 #define MUELU_REBALANCEBLOCKRESTRICTIONFACTORY_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 
66 #include "MueLu_HierarchyUtils.hpp"
68 #include "MueLu_Level.hpp"
69 #include "MueLu_MasterList.hpp"
70 #include "MueLu_Monitor.hpp"
71 #include "MueLu_PerfUtils.hpp"
72 
73 namespace MueLu {
74 
75 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
77  RCP<ParameterList> validParamList = rcp(new ParameterList());
78 
79 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
80  SET_VALID_ENTRY("repartition: use subcommunicators");
81 #undef SET_VALID_ENTRY
82 
83  //validParamList->set< RCP<const FactoryBase> >("repartition: use subcommunicators", Teuchos::null, "test");
84 
85  validParamList->set< RCP<const FactoryBase> >("R", Teuchos::null, "Factory of the restriction operator that need to be rebalanced (only used if type=Restriction)");
86  validParamList->set<RCP<const FactoryBase> >("Importer", Teuchos::null, "Generating factory of the matrix Importer for rebalancing");
87  validParamList->set<RCP<const FactoryBase> >("SubImporters", Teuchos::null, "Generating factory of the matrix sub-Importers for rebalancing");
88  validParamList->set<RCP<const FactoryBase> >("Nullspace", Teuchos::null, "Generating factory of the Nullspace operator");
89 
90  return validParamList;
91 }
92 
93 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
95  FactManager_.push_back(FactManager);
96 }
97 
98 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
100  Input(coarseLevel, "R");
101 
102  std::vector<Teuchos::RCP<const FactoryManagerBase> >::const_iterator it;
103  for(it = FactManager_.begin(); it!=FactManager_.end(); ++it) {
104  SetFactoryManager fineSFM (rcpFromRef(fineLevel), *it);
105  SetFactoryManager coarseSFM(rcpFromRef(coarseLevel), *it);
106 
107  if(!UseSingleSourceImporters_) coarseLevel.DeclareInput("Importer",(*it)->GetFactory("Importer").get(), this);
108  coarseLevel.DeclareInput("Nullspace",(*it)->GetFactory("Nullspace").get(), this);
109  }
110 
111  // Use the non-manager path if the maps / importers are generated in one place
112  if(UseSingleSourceImporters_) {
113  Input(coarseLevel,"SubImporters");
114  }
115 
116 
117 }
118 
119 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
121  FactoryMonitor m(*this, "Build", coarseLevel);
122 
123 
124 
125  RCP<Teuchos::FancyOStream> out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
126 
127  Teuchos::RCP<Matrix> originalTransferOp = Teuchos::null;
128  originalTransferOp = Get< RCP<Matrix> >(coarseLevel, "R");
129 
130  RCP<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > bOriginalTransferOp =
131  Teuchos::rcp_dynamic_cast<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(originalTransferOp);
132  TEUCHOS_TEST_FOR_EXCEPTION(bOriginalTransferOp==Teuchos::null, Exceptions::BadCast, "MueLu::RebalanceBlockTransferFactory::Build: input matrix P or R is not of type BlockedCrsMatrix! error.");
133 
134  RCP<const MapExtractor> rangeMapExtractor = bOriginalTransferOp->getRangeMapExtractor();
135  RCP<const MapExtractor> domainMapExtractor = bOriginalTransferOp->getDomainMapExtractor();
136 
137  // restrict communicator?
138  bool bRestrictComm = false;
139  const ParameterList& pL = GetParameterList();
140  if (pL.get<bool>("repartition: use subcommunicators") == true)
141  bRestrictComm = true;
142 
143  // check if GIDs for full maps have to be sorted:
144  // For the Thyra mode ordering they do not have to be sorted since the GIDs are
145  // numbered as 0...n1,0...,n2 (starting with zero for each subblock). The MapExtractor
146  // generates unique GIDs during the construction.
147  // For Xpetra style, the GIDs have to be reordered. Such that one obtains a ordered
148  // list of GIDs in an increasing ordering. In Xpetra, the GIDs are all unique through
149  // out all submaps.
150  bool bThyraRangeGIDs = rangeMapExtractor->getThyraMode();
151  bool bThyraDomainGIDs = domainMapExtractor->getThyraMode();
152 
153  // rebuild rebalanced blocked P operator
154  std::vector<GO> fullRangeMapVector;
155  std::vector<GO> fullDomainMapVector;
156  std::vector<RCP<const Map> > subBlockRRangeMaps;
157  std::vector<RCP<const Map> > subBlockRDomainMaps;
158  subBlockRRangeMaps.reserve(bOriginalTransferOp->Rows()); // reserve size for block P operators
159  subBlockRDomainMaps.reserve(bOriginalTransferOp->Cols()); // reserve size for block P operators
160 
161  std::vector<Teuchos::RCP<Matrix> > subBlockRebR;
162  subBlockRebR.reserve(bOriginalTransferOp->Cols());
163 
164  std::vector<RCP<const Import> > importers = std::vector<RCP<const Import> >(bOriginalTransferOp->Rows(), Teuchos::null);
165  if(UseSingleSourceImporters_) {
166  importers = Get<std::vector<RCP<const Import> > >(coarseLevel,"SubImporters");
167  }
168 
169  int curBlockId = 0;
170  Teuchos::RCP<const Import> rebalanceImporter;
171  std::vector<Teuchos::RCP<const FactoryManagerBase> >::const_iterator it;
172  for (it = FactManager_.begin(); it != FactManager_.end(); ++it) {
173  // begin SubFactoryManager environment
174  SetFactoryManager fineSFM (rcpFromRef(fineLevel), *it);
175  SetFactoryManager coarseSFM(rcpFromRef(coarseLevel), *it);
176 
177  if(UseSingleSourceImporters_) rebalanceImporter = importers[curBlockId];
178  else rebalanceImporter = coarseLevel.Get<Teuchos::RCP<const Import> >("Importer", (*it)->GetFactory("Importer").get());
179 
180  // extract matrix block
181  Teuchos::RCP<Matrix> Rii = bOriginalTransferOp->getMatrix(curBlockId, curBlockId);
182 
183  // TODO run this only in the debug version
184  TEUCHOS_TEST_FOR_EXCEPTION(bThyraRangeGIDs == true && Rii->getRowMap()->getMinAllGlobalIndex() != 0,
186  "MueLu::RebalanceBlockRestrictionFactory::Build: inconsistent Thyra GIDs. Thyra global ids for block range " << curBlockId << " start with " << Rii->getRowMap()->getMinAllGlobalIndex() << " but should start with 0");
187  TEUCHOS_TEST_FOR_EXCEPTION(bThyraDomainGIDs == true && Rii->getColMap()->getMinAllGlobalIndex() != 0,
189  "MueLu::RebalanceBlockRestrictionFactory::Build: inconsistent Thyra GIDs. Thyra global ids for block domain " << curBlockId << " start with " << Rii->getColMap()->getMinAllGlobalIndex() << " but should start with 0");
190 
191  Teuchos::RCP<Matrix> rebRii;
192  if(rebalanceImporter != Teuchos::null) {
193  std::stringstream ss; ss << "Rebalancing restriction block R(" << curBlockId << "," << curBlockId << ")";
194  SubFactoryMonitor m1(*this, ss.str(), coarseLevel);
195  {
196  SubFactoryMonitor subM(*this, "Rebalancing restriction -- fusedImport", coarseLevel);
197  // Note: The 3rd argument says to use originalR's domain map.
198 
199  RCP<Map> dummy;
200  rebRii = MatrixFactory::Build(Rii,*rebalanceImporter,dummy,rebalanceImporter->getTargetMap());
201  }
202 
203  RCP<ParameterList> params = rcp(new ParameterList());
204  params->set("printLoadBalancingInfo", true);
205  std::stringstream ss2; ss2 << "R(" << curBlockId << "," << curBlockId << ") rebalanced:";
206  GetOStream(Statistics0) << PerfUtils::PrintMatrixInfo(*rebRii, ss2.str(), params);
207  } else {
208  rebRii = Rii;
209  RCP<ParameterList> params = rcp(new ParameterList());
210  params->set("printLoadBalancingInfo", true);
211  std::stringstream ss2; ss2 << "R(" << curBlockId << "," << curBlockId << ") not rebalanced:";
212  GetOStream(Statistics0) << PerfUtils::PrintMatrixInfo(*rebRii, ss2.str(), params);
213  }
214 
215  // fix striding information for rebalanced diagonal block rebRii
216  Teuchos::RCP<const StridedMap> orig_stridedRgMap = Teuchos::rcp_dynamic_cast<const StridedMap>(rangeMapExtractor->getMap(Teuchos::as<size_t>(curBlockId),rangeMapExtractor->getThyraMode()));
217  Teuchos::RCP<const Map> stridedRgMap = Teuchos::null;
218  if(orig_stridedRgMap != Teuchos::null) {
219  std::vector<size_t> stridingData = orig_stridedRgMap->getStridingData();
220  Teuchos::ArrayView< const GlobalOrdinal > nodeRangeMapii = rebRii->getRangeMap()->getNodeElementList();
221  stridedRgMap = StridedMapFactory::Build(
222  originalTransferOp->getRangeMap()->lib(),
223  Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
224  nodeRangeMapii,
225  rebRii->getRangeMap()->getIndexBase(),
226  stridingData,
227  originalTransferOp->getRangeMap()->getComm(),
228  orig_stridedRgMap->getStridedBlockId(),
229  orig_stridedRgMap->getOffset());
230  } else stridedRgMap = Rii->getRangeMap();
231 
232  Teuchos::RCP<const StridedMap> orig_stridedDoMap = Teuchos::rcp_dynamic_cast<const StridedMap>(domainMapExtractor->getMap(Teuchos::as<size_t>(curBlockId),domainMapExtractor->getThyraMode()));
233  Teuchos::RCP<const Map> stridedDoMap = Teuchos::null;
234  if(orig_stridedDoMap != Teuchos::null) {
235  std::vector<size_t> stridingData = orig_stridedDoMap->getStridingData();
236  Teuchos::ArrayView< const GlobalOrdinal > nodeDomainMapii = rebRii->getDomainMap()->getNodeElementList();
237  stridedDoMap = StridedMapFactory::Build(
238  originalTransferOp->getDomainMap()->lib(),
239  Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
240  nodeDomainMapii,
241  rebRii->getDomainMap()->getIndexBase(),
242  stridingData,
243  originalTransferOp->getDomainMap()->getComm(),
244  orig_stridedDoMap->getStridedBlockId(),
245  orig_stridedDoMap->getOffset());
246  } else stridedDoMap = Rii->getDomainMap();
247 
248  if(bRestrictComm) {
249  stridedRgMap->removeEmptyProcesses();
250  stridedDoMap->removeEmptyProcesses();
251  }
252 
253  TEUCHOS_TEST_FOR_EXCEPTION(stridedRgMap == Teuchos::null,Exceptions::RuntimeError, "MueLu::RebalanceBlockRestrictionFactory::Build: failed to generate striding information. error.");
254  TEUCHOS_TEST_FOR_EXCEPTION(stridedDoMap == Teuchos::null,Exceptions::RuntimeError, "MueLu::RebalanceBlockRestrictionFactory::Build: failed to generate striding information. error.");
255 
256  // replace stridedMaps view in diagonal sub block
257  if(rebRii->IsView("stridedMaps")) rebRii->RemoveView("stridedMaps");
258  rebRii->CreateView("stridedMaps", stridedRgMap, stridedDoMap);
259 
260  // store rebalanced subblock
261  subBlockRebR.push_back(rebRii);
262 
263  // append strided row map (= range map) to list of range maps.
264  Teuchos::RCP<const Map> rangeMapii = rebRii->getRowMap("stridedMaps");
265  subBlockRRangeMaps.push_back(rangeMapii);
266  Teuchos::ArrayView< const GlobalOrdinal > nodeRangeMapii = rebRii->getRangeMap()->getNodeElementList();
267  // append the GIDs in the end. Do not sort if we have Thyra style GIDs
268  fullRangeMapVector.insert(fullRangeMapVector.end(), nodeRangeMapii.begin(), nodeRangeMapii.end());
269  if(bThyraRangeGIDs == false)
270  sort(fullRangeMapVector.begin(), fullRangeMapVector.end());
271 
272  // append strided col map (= domain map) to list of range maps.
273  Teuchos::RCP<const Map> domainMapii = rebRii->getColMap("stridedMaps");
274  subBlockRDomainMaps.push_back(domainMapii);
275  Teuchos::ArrayView< const GlobalOrdinal > nodeDomainMapii = rebRii->getDomainMap()->getNodeElementList();
276  // append the GIDs in the end. Do not sort if we have Thyra style GIDs
277  fullDomainMapVector.insert(fullDomainMapVector.end(), nodeDomainMapii.begin(), nodeDomainMapii.end());
278  if(bThyraDomainGIDs == false)
279  sort(fullDomainMapVector.begin(), fullDomainMapVector.end());
280 
282 
283  // rebalance null space
284  // This rebalances the null space partial vector associated with the current block (generated by the NullspaceFactory
285  // associated with the block)
286  if(rebalanceImporter != Teuchos::null)
287  { // rebalance null space
288  std::stringstream ss2; ss2 << "Rebalancing nullspace block(" << curBlockId << "," << curBlockId << ")";
289  SubFactoryMonitor subM(*this, ss2.str(), coarseLevel);
290 
291  RCP<MultiVector> nullspace = coarseLevel.Get<RCP<MultiVector> >("Nullspace", (*it)->GetFactory("Nullspace").get());
292  RCP<MultiVector> permutedNullspace = MultiVectorFactory::Build(rebalanceImporter->getTargetMap(), nullspace->getNumVectors());
293  permutedNullspace->doImport(*nullspace, *rebalanceImporter, Xpetra::INSERT);
294 
295  // TODO subcomm enabled everywhere or nowhere
296  if (bRestrictComm)
297  permutedNullspace->replaceMap(permutedNullspace->getMap()->removeEmptyProcesses());
298 
299  coarseLevel.Set<RCP<MultiVector> >("Nullspace", permutedNullspace, (*it)->GetFactory("Nullspace").get());
300 
301  } // end rebalance null space
302  else { // do nothing
303  RCP<MultiVector> nullspace = coarseLevel.Get<RCP<MultiVector> >("Nullspace", (*it)->GetFactory("Nullspace").get());
304  coarseLevel.Set<RCP<MultiVector> >("Nullspace", nullspace, (*it)->GetFactory("Nullspace").get());
305  }
306 
308 
309  curBlockId++;
310  } // end for loop
311 
312  // extract map index base from maps of blocked P
313  GO rangeIndexBase = originalTransferOp->getRangeMap()->getIndexBase();
314  GO domainIndexBase= originalTransferOp->getDomainMap()->getIndexBase();
315 
316  // check this
317  Teuchos::ArrayView<GO> fullRangeMapGIDs(fullRangeMapVector.size() ? &fullRangeMapVector[0] : 0,fullRangeMapVector.size());
318  Teuchos::RCP<const StridedMap> stridedRgFullMap = Teuchos::rcp_dynamic_cast<const StridedMap>(rangeMapExtractor->getFullMap());
319  Teuchos::RCP<const Map > fullRangeMap = Teuchos::null;
320  if(stridedRgFullMap != Teuchos::null) {
321  std::vector<size_t> stridedData = stridedRgFullMap->getStridingData();
322  fullRangeMap =
323  StridedMapFactory::Build(
324  originalTransferOp->getRangeMap()->lib(),
325  Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
326  fullRangeMapGIDs,
327  rangeIndexBase,
328  stridedData,
329  originalTransferOp->getRangeMap()->getComm(),
330  stridedRgFullMap->getStridedBlockId(),
331  stridedRgFullMap->getOffset());
332  } else {
333  fullRangeMap =
334  MapFactory::Build(
335  originalTransferOp->getRangeMap()->lib(),
336  Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
337  fullRangeMapGIDs,
338  rangeIndexBase,
339  originalTransferOp->getRangeMap()->getComm());
340  }
341 
342  Teuchos::ArrayView<GO> fullDomainMapGIDs(fullDomainMapVector.size() ? &fullDomainMapVector[0] : 0,fullDomainMapVector.size());
343  Teuchos::RCP<const StridedMap> stridedDoFullMap = Teuchos::rcp_dynamic_cast<const StridedMap>(domainMapExtractor->getFullMap());
344  Teuchos::RCP<const Map > fullDomainMap = Teuchos::null;
345  if(stridedDoFullMap != Teuchos::null) {
346  TEUCHOS_TEST_FOR_EXCEPTION(stridedDoFullMap==Teuchos::null, Exceptions::BadCast, "MueLu::BlockedPFactory::Build: full map in domain map extractor has no striding information! error.");
347  std::vector<size_t> stridedData2 = stridedDoFullMap->getStridingData();
348  fullDomainMap =
349  StridedMapFactory::Build(
350  originalTransferOp->getDomainMap()->lib(),
351  Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
352  fullDomainMapGIDs,
353  domainIndexBase,
354  stridedData2,
355  originalTransferOp->getDomainMap()->getComm(),
356  stridedDoFullMap->getStridedBlockId(),
357  stridedDoFullMap->getOffset());
358  } else {
359 
360  fullDomainMap =
361  MapFactory::Build(
362  originalTransferOp->getDomainMap()->lib(),
363  Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
364  fullDomainMapGIDs,
365  domainIndexBase,
366  originalTransferOp->getDomainMap()->getComm());
367  }
368 
369  if(bRestrictComm) {
370  fullRangeMap->removeEmptyProcesses();
371  fullDomainMap->removeEmptyProcesses();
372  }
373 
374  // build map extractors
375  Teuchos::RCP<const MapExtractor> rebrangeMapExtractor =
376  MapExtractorFactory::Build(fullRangeMap, subBlockRRangeMaps, bThyraRangeGIDs);
377  Teuchos::RCP<const MapExtractor> rebdomainMapExtractor =
378  MapExtractorFactory::Build(fullDomainMap, subBlockRDomainMaps, bThyraDomainGIDs);
379 
380  Teuchos::RCP<BlockedCrsMatrix> bRebR = Teuchos::rcp(new BlockedCrsMatrix(rebrangeMapExtractor,rebdomainMapExtractor,10));
381  for(size_t i = 0; i<subBlockRRangeMaps.size(); i++) {
382  Teuchos::RCP<CrsMatrixWrap> crsOpii = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(subBlockRebR[i]);
383  bRebR->setMatrix(i,i,crsOpii);
384  }
385 
386  bRebR->fillComplete();
387 
388  Set(coarseLevel, "R", Teuchos::rcp_dynamic_cast<Matrix>(bRebR)); // do nothing // TODO remove this!
389 
390 } // Build
391 
392 } // namespace MueLu
393 
394 #endif /* MUELU_REBALANCEBLOCKRESTRICTIONFACTORY_DEF_HPP_ */
Exception indicating invalid cast attempted.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Specifies the data that this class needs, and the factories that generate that data.
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)
#define SET_VALID_ENTRY(name)
Exception throws to report errors in the internal logical of the program.
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()