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"
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 
69 #include "MueLu_HierarchyHelpers.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> >("R", Teuchos::null, "Factory of the restriction operator that need to be rebalanced (only used if type=Restriction)");
82  // validParamList->set< RCP<const FactoryBase> >("Nullspace", Teuchos::null, "Factory of the nullspace that need to be rebalanced (only used if type=Restriction)");
83  // validParamList->set< RCP<const FactoryBase> >("Coordinates", Teuchos::null, "Factory of the coordinates that need to be rebalanced (only used if type=Restriction)");
84  // validParamList->set< RCP<const FactoryBase> >("Importer", Teuchos::null, "Factory of the importer object used for the rebalancing");
85 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
86  // SET_VALID_ENTRY("repartition: use subcommunicators");
87 #undef SET_VALID_ENTRY
88  // validParamList->set< int > ("write start", -1, "First level at which coordinates should be written to file");
89  // validParamList->set< int > ("write end", -1, "Last level at which coordinates should be written to file");
90 
91  // TODO validation: "P" parameter valid only for type="Interpolation" and "R" valid only for type="Restriction". Like so:
92  // if (paramList.isEntry("type") && paramList.get("type) == "Interpolation) {
93  // validParamList->set< RCP<const FactoryBase> >("P", Teuchos::null, "Factory of the prolongation operator that need to be rebalanced (only used if type=Interpolation)");
94 
95  return validParamList;
96 }
97 
98 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
100  FactManager_.push_back(FactManager);
101 }
102 
103 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
105 
106  Input(coarseLevel, "P");
107 
108  std::vector<Teuchos::RCP<const FactoryManagerBase> >::const_iterator it;
109  for(it = FactManager_.begin(); it!=FactManager_.end(); ++it) {
110  SetFactoryManager fineSFM (rcpFromRef(fineLevel), *it);
111  SetFactoryManager coarseSFM(rcpFromRef(coarseLevel), *it);
112  coarseLevel.DeclareInput("Importer",(*it)->GetFactory("Importer").get(), this);
113  }
114 }
115 
116 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
118  FactoryMonitor m(*this, "Build", coarseLevel);
119 
120  RCP<Teuchos::FancyOStream> out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
121 
122  Teuchos::RCP<Matrix> originalTransferOp = Teuchos::null;
123  originalTransferOp = Get< RCP<Matrix> >(coarseLevel, "P");
124 
125  RCP<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > bOriginalTransferOp =
126  Teuchos::rcp_dynamic_cast<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > (originalTransferOp);
127  TEUCHOS_TEST_FOR_EXCEPTION(bOriginalTransferOp==Teuchos::null, Exceptions::BadCast, "MueLu::RebalanceBlockTransferFactory::Build: input matrix P or R is not of type BlockedCrsMatrix! error.");
128 
129  // plausibility check
130  TEUCHOS_TEST_FOR_EXCEPTION(bOriginalTransferOp->Rows() != 2,Exceptions::RuntimeError, "MueLu::RebalanceBlockTransferFactory::Build: number of block rows of transfer operator is not equal 2. error.");
131  TEUCHOS_TEST_FOR_EXCEPTION(bOriginalTransferOp->Cols() != 2,Exceptions::RuntimeError, "MueLu::RebalanceBlockTransferFactory::Build: number of block columns of transfer operator is not equal 2. error.");
132 
133  // declare variables for maps of blocked rebalanced prolongation operator
134  std::vector<GO> fullRangeMapVector;
135  std::vector<GO> fullDomainMapVector;
136  std::vector<RCP<const Map> > subBlockPRangeMaps;
137  std::vector<RCP<const Map> > subBlockPDomainMaps;
138  subBlockPRangeMaps.reserve(bOriginalTransferOp->Rows()); // reserve size for block P operators
139  subBlockPDomainMaps.reserve(bOriginalTransferOp->Cols()); // reserve size for block P operators
140 
141  std::vector<Teuchos::RCP<Matrix> > subBlockRebP;
142  subBlockRebP.reserve(bOriginalTransferOp->Rows());
143 
144  int curBlockId = 0;
145  Teuchos::RCP<const Import> rebalanceImporter = Teuchos::null;
146  std::vector<Teuchos::RCP<const FactoryManagerBase> >::const_iterator it;
147  for (it = FactManager_.begin(); it != FactManager_.end(); ++it) {
148  // begin SubFactoryManager environment
149  SetFactoryManager fineSFM (rcpFromRef(fineLevel), *it);
150  SetFactoryManager coarseSFM(rcpFromRef(coarseLevel), *it);
151 
152  rebalanceImporter = coarseLevel.Get<Teuchos::RCP<const Import> >("Importer", (*it)->GetFactory("Importer").get());
153 
154  // extract diagonal matrix block
155  Teuchos::RCP<CrsMatrix> Pmii = bOriginalTransferOp->getMatrix(curBlockId, curBlockId);
156  Teuchos::RCP<CrsMatrixWrap> Pwii = Teuchos::rcp(new CrsMatrixWrap(Pmii));
157  Teuchos::RCP<Matrix> Pii = Teuchos::rcp_dynamic_cast<Matrix>(Pwii);
158 
159  // rebalance P11
160  if(rebalanceImporter != Teuchos::null) {
161  std::stringstream ss; ss << "Rebalancing prolongator block P(" << curBlockId << "," << curBlockId << ")";
162  SubFactoryMonitor m1(*this, ss.str(), coarseLevel);
163 
164  // P is the transfer operator from the coarse grid to the fine grid.
165  // P must transfer the data from the newly reordered coarse A to the (unchanged) fine A.
166  // This means that the domain map (coarse) of P must be changed according to the new partition. The range map (fine) is kept unchanged.
167  //
168  // The domain map of P must match the range map of R.
169  // See also note below about domain/range map of R and its implications for P.
170  //
171  // To change the domain map of P, P needs to be fillCompleted again with the new domain map.
172  // To achieve this, P is copied into a new matrix that is not fill-completed.
173  // The doImport() operation is just used here to make a copy of P: the importer is trivial and there is no data movement involved.
174  // The reordering actually happens during the fillComplete() with domainMap == rebalanceImporter->getTargetMap().
175 
176  RCP<const Import> newImporter;
177  {
178  SubFactoryMonitor subM(*this, "Rebalancing prolongator -- fast map replacement", coarseLevel);
179  newImporter = ImportFactory::Build(rebalanceImporter->getTargetMap(), Pii->getColMap());
180  Pmii->replaceDomainMapAndImporter(rebalanceImporter->getTargetMap(), newImporter);
181  }
182 
183  RCP<ParameterList> params = rcp(new ParameterList());
184  params->set("printLoadBalancingInfo", true);
185  std::stringstream ss2; ss2 << "P(" << curBlockId << "," << curBlockId << ") rebalanced:";
186  GetOStream(Statistics0) << PerfUtils::PrintMatrixInfo(*Pii, ss2.str(), params);
187 
188  // store rebalanced P block
189  subBlockRebP.push_back(Pii);
190  } // end rebalance P(1,1)
191  else {
192  RCP<ParameterList> params = rcp(new ParameterList());
193  params->set("printLoadBalancingInfo", true);
194  std::stringstream ss; ss << "P(" << curBlockId << "," << curBlockId << ") not rebalanced:";
195  GetOStream(Statistics0) << PerfUtils::PrintMatrixInfo(*Pii, ss.str(), params);
196  // store rebalanced P block
197  subBlockRebP.push_back(Pii);
198  }
199 
200  // fix striding information for rebalanced diagonal block Pii
201  RCP<const Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node> > rgPMapExtractor = bOriginalTransferOp->getRangeMapExtractor(); // original map extractor
202  Teuchos::RCP<const StridedMap> orig_stridedRgMap = Teuchos::rcp_dynamic_cast<const StridedMap>(rgPMapExtractor->getMap(Teuchos::as<size_t>(curBlockId)));
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 = Pii->getRangeMap()->getNodeElementList();
207  stridedRgMap = StridedMapFactory::Build(
208  originalTransferOp->getRangeMap()->lib(),
209  Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
210  nodeRangeMapii,
211  Pii->getRangeMap()->getIndexBase(),
212  stridingData,
213  originalTransferOp->getRangeMap()->getComm(),
214  orig_stridedRgMap->getStridedBlockId(),
215  orig_stridedRgMap->getOffset());
216  } else stridedRgMap = Pii->getRangeMap();
217  RCP<const Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node> > doPMapExtractor = bOriginalTransferOp->getDomainMapExtractor(); // original map extractor
218  Teuchos::RCP<const StridedMap> orig_stridedDoMap = Teuchos::rcp_dynamic_cast<const StridedMap>(doPMapExtractor->getMap(Teuchos::as<size_t>(curBlockId)));
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 = Pii->getDomainMap()->getNodeElementList();
223  stridedDoMap = StridedMapFactory::Build(originalTransferOp->getDomainMap()->lib(),
224  Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
225  nodeDomainMapii,
226  Pii->getDomainMap()->getIndexBase(),
227  stridingData,
228  originalTransferOp->getDomainMap()->getComm(),
229  orig_stridedDoMap->getStridedBlockId(),
230  orig_stridedDoMap->getOffset());
231  } else stridedDoMap = Pii->getDomainMap();
232 
233  TEUCHOS_TEST_FOR_EXCEPTION(stridedRgMap == Teuchos::null,Exceptions::RuntimeError, "MueLu::RebalanceBlockInterpolationFactory::Build: failed to generate striding information. error.");
234  TEUCHOS_TEST_FOR_EXCEPTION(stridedDoMap == Teuchos::null,Exceptions::RuntimeError, "MueLu::RebalanceBlockInterpolationFactory::Build: failed to generate striding information. error.");
235 
236  // replace stridedMaps view in diagonal sub block
237  if(Pii->IsView("stridedMaps")) Pii->RemoveView("stridedMaps");
238  Pii->CreateView("stridedMaps", stridedRgMap, stridedDoMap);
239 
240  // append strided row map (= range map) to list of range maps.
241  Teuchos::RCP<const Map> rangeMapii = Pii->getRowMap("stridedMaps"); //Pii->getRangeMap();
242  subBlockPRangeMaps.push_back(rangeMapii);
243  Teuchos::ArrayView< const GlobalOrdinal > nodeRangeMapii = Pii->getRangeMap()->getNodeElementList();
244  fullRangeMapVector.insert(fullRangeMapVector.end(), nodeRangeMapii.begin(), nodeRangeMapii.end());
245  sort(fullRangeMapVector.begin(), fullRangeMapVector.end());
246 
247  // append strided col map (= domain map) to list of range maps.
248  Teuchos::RCP<const Map> domainMapii = Pii->getColMap("stridedMaps"); //Pii->getDomainMap();
249  subBlockPDomainMaps.push_back(domainMapii);
250  Teuchos::ArrayView< const GlobalOrdinal > nodeDomainMapii = Pii->getDomainMap()->getNodeElementList();
251  fullDomainMapVector.insert(fullDomainMapVector.end(), nodeDomainMapii.begin(), nodeDomainMapii.end());
252  sort(fullDomainMapVector.begin(), fullDomainMapVector.end());
253 
254  curBlockId++; // increase block id index
255 
256  } // end SubFactoryManager environment
257 
258  // extract map index base from maps of blocked P
259  GO rangeIndexBase = originalTransferOp->getRangeMap()->getIndexBase();
260  GO domainIndexBase= originalTransferOp->getDomainMap()->getIndexBase();
261 
262  RCP<const Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node> > rangePMapExtractor = bOriginalTransferOp->getRangeMapExtractor(); // original map extractor
263  Teuchos::ArrayView<GO> fullRangeMapGIDs(fullRangeMapVector.size() ? &fullRangeMapVector[0] : 0,fullRangeMapVector.size());
264  Teuchos::RCP<const StridedMap> stridedRgFullMap = Teuchos::rcp_dynamic_cast<const StridedMap>(rangePMapExtractor->getFullMap());
265  Teuchos::RCP<const Map > fullRangeMap = Teuchos::null;
266  if(stridedRgFullMap != Teuchos::null) {
267  std::vector<size_t> stridedData = stridedRgFullMap->getStridingData();
268  fullRangeMap =
269  StridedMapFactory::Build(
270  originalTransferOp->getRangeMap()->lib(),
271  Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
272  fullRangeMapGIDs,
273  rangeIndexBase,
274  stridedData,
275  originalTransferOp->getRangeMap()->getComm(),
276  stridedRgFullMap->getStridedBlockId(),
277  stridedRgFullMap->getOffset());
278  } else {
279  fullRangeMap =
280  MapFactory::Build(
281  originalTransferOp->getRangeMap()->lib(),
282  Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
283  fullRangeMapGIDs,
284  rangeIndexBase,
285  originalTransferOp->getRangeMap()->getComm());
286  }
287 
288  RCP<const Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node> > domainAMapExtractor = bOriginalTransferOp->getDomainMapExtractor();
289  Teuchos::ArrayView<GO> fullDomainMapGIDs(fullDomainMapVector.size() ? &fullDomainMapVector[0] : 0,fullDomainMapVector.size());
290  Teuchos::RCP<const StridedMap> stridedDoFullMap = Teuchos::rcp_dynamic_cast<const StridedMap>(domainAMapExtractor->getFullMap());
291  Teuchos::RCP<const Map > fullDomainMap = Teuchos::null;
292  if(stridedDoFullMap != Teuchos::null) {
293  TEUCHOS_TEST_FOR_EXCEPTION(stridedDoFullMap==Teuchos::null, Exceptions::BadCast, "MueLu::BlockedPFactory::Build: full map in domain map extractor has no striding information! error.");
294  std::vector<size_t> stridedData2 = stridedDoFullMap->getStridingData();
295  fullDomainMap =
296  StridedMapFactory::Build(
297  originalTransferOp->getDomainMap()->lib(),
298  Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
299  fullDomainMapGIDs,
300  domainIndexBase,
301  stridedData2,
302  originalTransferOp->getDomainMap()->getComm(),
303  stridedDoFullMap->getStridedBlockId(),
304  stridedDoFullMap->getOffset());
305  } else {
306 
307  fullDomainMap =
308  MapFactory::Build(
309  originalTransferOp->getDomainMap()->lib(),
310  Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
311  fullDomainMapGIDs,
312  domainIndexBase,
313  originalTransferOp->getDomainMap()->getComm());
314  }
315 
316  // build map extractors
317  Teuchos::RCP<const Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node> > rangeMapExtractor =
318  Xpetra::MapExtractorFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(fullRangeMap, subBlockPRangeMaps);
319  Teuchos::RCP<const Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node> > domainMapExtractor =
320  Xpetra::MapExtractorFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(fullDomainMap, subBlockPDomainMaps);
321 
322  Teuchos::RCP<BlockedCrsMatrix> bRebP = Teuchos::rcp(new BlockedCrsMatrix(rangeMapExtractor,domainMapExtractor,10));
323  for(size_t i = 0; i<subBlockPRangeMaps.size(); i++) {
324  Teuchos::RCP<const CrsMatrixWrap> crsOpii = Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(subBlockRebP[i]);
325  Teuchos::RCP<CrsMatrix> crsMatii = crsOpii->getCrsMatrix();
326  bRebP->setMatrix(i,i,crsMatii);
327  }
328  bRebP->fillComplete();
329 
330  Set(coarseLevel, "P", Teuchos::rcp_dynamic_cast<Matrix>(bRebP));
331 
332 
333 } // Build
334 
335 } // namespace MueLu
336 
337 #endif /* HAVE_MUELU_EXPERIMENTAL */
338 #endif /* MUELU_REBALANCEBLOCKINTERPOLATIONFACTORY_DEF_HPP_ */
Exception indicating invalid cast attempted.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
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.
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
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()