MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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"
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 
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 
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(),
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(),
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(),
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(),
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(),
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(),
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_ */
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Specifies the data that this class needs, and the factories that generate that data.
Exception indicating invalid cast attempted.
GlobalOrdinal GO
iterator begin() const
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)
T * get() 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)
iterator end() const
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
#define SET_VALID_ENTRY(name)
Exception throws to report errors in the internal logical of the program.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
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()