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