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 
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(),
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(),
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(),
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(),
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(),
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(),
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.
T & get(const std::string &name, T def_value)
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Specifies the data that this class needs, and the factories that generate that data.
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)
Namespace for MueLu classes and methods.
T * get() const
iterator begin() 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.
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.
iterator end() const
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()