MueLu  Version of the Day
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 #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"
57 #include <Xpetra_Matrix.hpp>
59 #include <Xpetra_MapFactory.hpp>
60 #include <Xpetra_MapExtractor.hpp>
62 #include <Xpetra_MatrixFactory.hpp>
63 #include <Xpetra_Import.hpp>
64 #include <Xpetra_ImportFactory.hpp>
65 
67 
69 #include "MueLu_HierarchyUtils.hpp"
70 #include "MueLu_Level.hpp"
71 #include "MueLu_Monitor.hpp"
72 #include "MueLu_PerfUtils.hpp"
73 
74 namespace MueLu {
75 
76 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
78  RCP<ParameterList> validParamList = rcp(new ParameterList());
79 
80  validParamList->set< RCP<const FactoryBase> >("P", Teuchos::null, "Factory of the prolongation operator that need to be rebalanced (only used if type=Interpolation)");
81  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.");
82 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
83  // SET_VALID_ENTRY("repartition: use subcommunicators");
84 #undef SET_VALID_ENTRY
85 
86  // TODO validation: "P" parameter valid only for type="Interpolation" and "R" valid only for type="Restriction". Like so:
87  // if (paramList.isEntry("type") && paramList.get("type) == "Interpolation) {
88  // validParamList->set< RCP<const FactoryBase> >("P", Teuchos::null, "Factory of the prolongation operator that need to be rebalanced (only used if type=Interpolation)");
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 
101  Input(coarseLevel, "P");
102  Input(coarseLevel, "A"); // we request the non-rebalanced coarse level A since we have to make sure it is calculated before rebalancing starts!
103 
104  std::vector<Teuchos::RCP<const FactoryManagerBase> >::const_iterator it;
105  for(it = FactManager_.begin(); it!=FactManager_.end(); ++it) {
106  SetFactoryManager fineSFM (rcpFromRef(fineLevel), *it);
107  SetFactoryManager coarseSFM(rcpFromRef(coarseLevel), *it);
108  coarseLevel.DeclareInput("Importer",(*it)->GetFactory("Importer").get(), this);
109  }
110 }
111 
112 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
114  FactoryMonitor m(*this, "Build", coarseLevel);
115 
116  RCP<Teuchos::FancyOStream> out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
117 
118  Teuchos::RCP<Matrix> nonrebCoarseA = Get< RCP<Matrix> >(coarseLevel, "A");
119 
120  Teuchos::RCP<Matrix> originalTransferOp = Teuchos::null;
121  originalTransferOp = Get< RCP<Matrix> >(coarseLevel, "P");
122 
123  RCP<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > bOriginalTransferOp =
124  Teuchos::rcp_dynamic_cast<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > (originalTransferOp);
125  TEUCHOS_TEST_FOR_EXCEPTION(bOriginalTransferOp==Teuchos::null, Exceptions::BadCast, "MueLu::RebalanceBlockTransferFactory::Build: input matrix P or R is not of type BlockedCrsMatrix! error.");
126 
127  RCP<const MapExtractor> rangeMapExtractor = bOriginalTransferOp->getRangeMapExtractor();
128  RCP<const MapExtractor> domainMapExtractor = bOriginalTransferOp->getDomainMapExtractor();
129 
130  // check if GIDs for full maps have to be sorted:
131  // For the Thyra mode ordering they do not have to be sorted since the GIDs are
132  // numbered as 0...n1,0...,n2 (starting with zero for each subblock). The MapExtractor
133  // generates unique GIDs during the construction.
134  // For Xpetra style, the GIDs have to be reordered. Such that one obtains a ordered
135  // list of GIDs in an increasing ordering. In Xpetra, the GIDs are all unique through
136  // out all submaps.
137  bool bThyraRangeGIDs = rangeMapExtractor->getThyraMode();
138  bool bThyraDomainGIDs = domainMapExtractor->getThyraMode();
139 
140  // declare variables for maps of blocked rebalanced prolongation operator
141  std::vector<GO> fullRangeMapVector; // contains all range GIDs on current processor
142  std::vector<GO> fullDomainMapVector; // contains all domain GIDs on current processor
143  std::vector<RCP<const Map> > subBlockPRangeMaps;
144  std::vector<RCP<const Map> > subBlockPDomainMaps;
145  subBlockPRangeMaps.reserve(bOriginalTransferOp->Rows()); // reserve size for block P operators
146  subBlockPDomainMaps.reserve(bOriginalTransferOp->Cols()); // reserve size for block P operators
147 
148  std::vector<Teuchos::RCP<Matrix> > subBlockRebP;
149  subBlockRebP.reserve(bOriginalTransferOp->Rows());
150 
151  int curBlockId = 0;
152  Teuchos::RCP<const Import> rebalanceImporter = Teuchos::null;
153  std::vector<Teuchos::RCP<const FactoryManagerBase> >::const_iterator it;
154  for (it = FactManager_.begin(); it != FactManager_.end(); ++it) {
155  // begin SubFactoryManager environment
156  SetFactoryManager fineSFM (rcpFromRef(fineLevel), *it);
157  SetFactoryManager coarseSFM(rcpFromRef(coarseLevel), *it);
158 
159  rebalanceImporter = coarseLevel.Get<Teuchos::RCP<const Import> >("Importer", (*it)->GetFactory("Importer").get());
160 
161  // extract diagonal matrix block
162  Teuchos::RCP<Matrix> Pii = bOriginalTransferOp->getMatrix(curBlockId, curBlockId);
163  Teuchos::RCP<CrsMatrixWrap> Pwii = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(Pii);
164  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!");
165 
166  // TODO run this only in the debug version
167  TEUCHOS_TEST_FOR_EXCEPTION(bThyraRangeGIDs == true && Pii->getRowMap()->getMinAllGlobalIndex() != 0,
169  "MueLu::RebalanceBlockInterpolationFactory::Build: inconsistent Thyra GIDs. Thyra global ids for block range " << curBlockId << " start with " << Pii->getRowMap()->getMinAllGlobalIndex() << " but should start with 0");
170  TEUCHOS_TEST_FOR_EXCEPTION(bThyraDomainGIDs == true && Pii->getColMap()->getMinAllGlobalIndex() != 0,
172  "MueLu::RebalanceBlockInterpolationFactory::Build: inconsistent Thyra GIDs. Thyra global ids for block domain " << curBlockId << " start with " << Pii->getColMap()->getMinAllGlobalIndex() << " but should start with 0");
173 
174  // rebalance P11
175  if(rebalanceImporter != Teuchos::null) {
176  std::stringstream ss; ss << "Rebalancing prolongator block P(" << curBlockId << "," << curBlockId << ")";
177  SubFactoryMonitor m1(*this, ss.str(), coarseLevel);
178 
179  // P is the transfer operator from the coarse grid to the fine grid.
180  // P must transfer the data from the newly reordered coarse A to the (unchanged) fine A.
181  // This means that the domain map (coarse) of P must be changed according to the new partition. The range map (fine) is kept unchanged.
182  //
183  // The domain map of P must match the range map of R.
184  // See also note below about domain/range map of R and its implications for P.
185  //
186  // To change the domain map of P, P needs to be fillCompleted again with the new domain map.
187  // To achieve this, P is copied into a new matrix that is not fill-completed.
188  // The doImport() operation is just used here to make a copy of P: the importer is trivial and there is no data movement involved.
189  // The reordering actually happens during the fillComplete() with domainMap == rebalanceImporter->getTargetMap().
190  RCP<const Import> newImporter;
191  {
192  SubFactoryMonitor subM(*this, "Rebalancing prolongator -- fast map replacement", coarseLevel);
193  newImporter = ImportFactory::Build(rebalanceImporter->getTargetMap(), Pii->getColMap());
194  Pwii->getCrsMatrix()->replaceDomainMapAndImporter(rebalanceImporter->getTargetMap(), newImporter);
195  }
196 
197  RCP<ParameterList> params = rcp(new ParameterList());
198  params->set("printLoadBalancingInfo", true);
199  std::stringstream ss2; ss2 << "P(" << curBlockId << "," << curBlockId << ") rebalanced:";
200  GetOStream(Statistics0) << PerfUtils::PrintMatrixInfo(*Pii, ss2.str(), params);
201 
202  // store rebalanced P block
203  subBlockRebP.push_back(Pii);
204  } // end rebalance P(1,1)
205  else {
206  RCP<ParameterList> params = rcp(new ParameterList());
207  params->set("printLoadBalancingInfo", true);
208  std::stringstream ss; ss << "P(" << curBlockId << "," << curBlockId << ") not rebalanced:";
209  GetOStream(Statistics0) << PerfUtils::PrintMatrixInfo(*Pii, ss.str(), params);
210  // store rebalanced P block
211  subBlockRebP.push_back(Pii);
212  }
213 
214  // fix striding information for rebalanced diagonal block Pii
215  //RCP<const Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node> > rgPMapExtractor = bOriginalTransferOp->getRangeMapExtractor(); // original map extractor
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 = Pii->getRangeMap()->getNodeElementList();
221  stridedRgMap = StridedMapFactory::Build(
222  originalTransferOp->getRangeMap()->lib(),
223  Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
224  nodeRangeMapii,
225  Pii->getRangeMap()->getIndexBase(),
226  stridingData,
227  originalTransferOp->getRangeMap()->getComm(),
228  orig_stridedRgMap->getStridedBlockId(),
229  orig_stridedRgMap->getOffset());
230  } else stridedRgMap = Pii->getRangeMap();
231 
232  //RCP<const Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node> > doPMapExtractor = bOriginalTransferOp->getDomainMapExtractor(); // original map extractor
233  Teuchos::RCP<const StridedMap> orig_stridedDoMap = Teuchos::rcp_dynamic_cast<const StridedMap>(domainMapExtractor->getMap(Teuchos::as<size_t>(curBlockId),domainMapExtractor->getThyraMode()));
234 
235  Teuchos::RCP<const Map> stridedDoMap = Teuchos::null;
236  if(orig_stridedDoMap != Teuchos::null) {
237  std::vector<size_t> stridingData = orig_stridedDoMap->getStridingData();
238  Teuchos::ArrayView< const GlobalOrdinal > nodeDomainMapii = Pii->getDomainMap()->getNodeElementList();
239  stridedDoMap = StridedMapFactory::Build(originalTransferOp->getDomainMap()->lib(),
240  Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
241  nodeDomainMapii,
242  Pii->getDomainMap()->getIndexBase(),
243  stridingData,
244  originalTransferOp->getDomainMap()->getComm(),
245  orig_stridedDoMap->getStridedBlockId(),
246  orig_stridedDoMap->getOffset());
247  } else stridedDoMap = Pii->getDomainMap();
248 
249  TEUCHOS_TEST_FOR_EXCEPTION(stridedRgMap == Teuchos::null,Exceptions::RuntimeError, "MueLu::RebalanceBlockInterpolationFactory::Build: failed to generate striding information. error.");
250  TEUCHOS_TEST_FOR_EXCEPTION(stridedDoMap == Teuchos::null,Exceptions::RuntimeError, "MueLu::RebalanceBlockInterpolationFactory::Build: failed to generate striding information. error.");
251 
252  // replace stridedMaps view in diagonal sub block
253  if(Pii->IsView("stridedMaps")) Pii->RemoveView("stridedMaps");
254  Pii->CreateView("stridedMaps", stridedRgMap, stridedDoMap);
255 
256  // append strided row map (= range map) to list of range maps.
257  Teuchos::RCP<const Map> rangeMapii = Pii->getRowMap("stridedMaps"); // strided range map
258  subBlockPRangeMaps.push_back(rangeMapii);
259  Teuchos::ArrayView< const GlobalOrdinal > nodeRangeMapii = Pii->getRangeMap()->getNodeElementList();
260  // append the GIDs in the end. Do not sort if we have Thyra style GIDs
261  fullRangeMapVector.insert(fullRangeMapVector.end(), nodeRangeMapii.begin(), nodeRangeMapii.end());
262  if(bThyraRangeGIDs == false)
263  sort(fullRangeMapVector.begin(), fullRangeMapVector.end());
264 
265  // append strided col map (= domain map) to list of range maps.
266  Teuchos::RCP<const Map> domainMapii = Pii->getColMap("stridedMaps"); // strided domain map
267  subBlockPDomainMaps.push_back(domainMapii);
268  Teuchos::ArrayView< const GlobalOrdinal > nodeDomainMapii = Pii->getDomainMap()->getNodeElementList();
269  // append the GIDs in the end. Do not sort if we have Thyra style GIDs
270  fullDomainMapVector.insert(fullDomainMapVector.end(), nodeDomainMapii.begin(), nodeDomainMapii.end());
271  if(bThyraDomainGIDs == false)
272  sort(fullDomainMapVector.begin(), fullDomainMapVector.end());
273 
274  curBlockId++; // increase block id index
275 
276  } // end SubFactoryManager environment
277 
278  // extract map index base from maps of blocked P
279  GO rangeIndexBase = originalTransferOp->getRangeMap()->getIndexBase();
280  GO domainIndexBase= originalTransferOp->getDomainMap()->getIndexBase();
281 
282  RCP<const Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node> > rangePMapExtractor = bOriginalTransferOp->getRangeMapExtractor(); // original map extractor
283  Teuchos::ArrayView<GO> fullRangeMapGIDs(fullRangeMapVector.size() ? &fullRangeMapVector[0] : 0,fullRangeMapVector.size());
284  Teuchos::RCP<const StridedMap> stridedRgFullMap = Teuchos::rcp_dynamic_cast<const StridedMap>(rangePMapExtractor->getFullMap());
285  Teuchos::RCP<const Map > fullRangeMap = Teuchos::null;
286  if(stridedRgFullMap != Teuchos::null) {
287  std::vector<size_t> stridedData = stridedRgFullMap->getStridingData();
288  fullRangeMap =
289  StridedMapFactory::Build(
290  originalTransferOp->getRangeMap()->lib(),
291  Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
292  fullRangeMapGIDs,
293  rangeIndexBase,
294  stridedData,
295  originalTransferOp->getRangeMap()->getComm(),
296  stridedRgFullMap->getStridedBlockId(),
297  stridedRgFullMap->getOffset());
298  } else {
299  fullRangeMap =
300  MapFactory::Build(
301  originalTransferOp->getRangeMap()->lib(),
302  Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
303  fullRangeMapGIDs,
304  rangeIndexBase,
305  originalTransferOp->getRangeMap()->getComm());
306  }
307 
308  RCP<const Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node> > domainAMapExtractor = bOriginalTransferOp->getDomainMapExtractor();
309  Teuchos::ArrayView<GO> fullDomainMapGIDs(fullDomainMapVector.size() ? &fullDomainMapVector[0] : 0,fullDomainMapVector.size());
310  Teuchos::RCP<const StridedMap> stridedDoFullMap = Teuchos::rcp_dynamic_cast<const StridedMap>(domainAMapExtractor->getFullMap());
311  Teuchos::RCP<const Map > fullDomainMap = Teuchos::null;
312  if(stridedDoFullMap != Teuchos::null) {
313  TEUCHOS_TEST_FOR_EXCEPTION(stridedDoFullMap==Teuchos::null, Exceptions::BadCast, "MueLu::BlockedPFactory::Build: full map in domain map extractor has no striding information! error.");
314  std::vector<size_t> stridedData2 = stridedDoFullMap->getStridingData();
315  fullDomainMap =
316  StridedMapFactory::Build(
317  originalTransferOp->getDomainMap()->lib(),
318  Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
319  fullDomainMapGIDs,
320  domainIndexBase,
321  stridedData2,
322  originalTransferOp->getDomainMap()->getComm(),
323  stridedDoFullMap->getStridedBlockId(),
324  stridedDoFullMap->getOffset());
325  } else {
326 
327  fullDomainMap =
328  MapFactory::Build(
329  originalTransferOp->getDomainMap()->lib(),
330  Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
331  fullDomainMapGIDs,
332  domainIndexBase,
333  originalTransferOp->getDomainMap()->getComm());
334  }
335 
336  // build map extractors
337  Teuchos::RCP<const MapExtractor> rebrangeMapExtractor =
338  MapExtractorFactory::Build(fullRangeMap, subBlockPRangeMaps, bThyraRangeGIDs);
339  Teuchos::RCP<const MapExtractor> rebdomainMapExtractor =
340  MapExtractorFactory::Build(fullDomainMap, subBlockPDomainMaps, bThyraDomainGIDs);
341 
342  Teuchos::RCP<BlockedCrsMatrix> bRebP = Teuchos::rcp(new BlockedCrsMatrix(rebrangeMapExtractor,rebdomainMapExtractor,10));
343  for(size_t i = 0; i<subBlockPRangeMaps.size(); i++) {
344  Teuchos::RCP<CrsMatrixWrap> crsOpii = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(subBlockRebP[i]);
345  TEUCHOS_TEST_FOR_EXCEPTION(crsOpii == Teuchos::null,Xpetra::Exceptions::BadCast, "MueLu::RebalanceBlockTransferFactory::Build: block P" << i << " is not of type CrsMatrixWrap.");
346  bRebP->setMatrix(i,i,crsOpii);
347  }
348  bRebP->fillComplete();
349 
350  Set(coarseLevel, "P", Teuchos::rcp_dynamic_cast<Matrix>(bRebP));
351 } // Build
352 
353 } // namespace MueLu
354 
355 #endif /* HAVE_MUELU_EXPERIMENTAL */
356 #endif /* MUELU_REBALANCEBLOCKINTERPOLATIONFACTORY_DEF_HPP_ */
Exception indicating invalid cast attempted.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
GlobalOrdinal GO
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)
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Specifies the data that this class needs, and the factories that generate that data.
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()