MueLu Version of the Day
Loading...
Searching...
No Matches
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"
52#include "Xpetra_MultiVectorFactory.hpp"
53#include "Xpetra_Vector.hpp"
54#include "Xpetra_VectorFactory.hpp"
55#include <Xpetra_Matrix.hpp>
56#include <Xpetra_BlockedCrsMatrix.hpp>
57#include <Xpetra_MapFactory.hpp>
58#include <Xpetra_MapExtractor.hpp>
59#include <Xpetra_MapExtractorFactory.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
73namespace MueLu {
74
75template <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
93template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
95 FactManager_.push_back(FactManager);
96}
97
98template <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
119template <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
130 RCP<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > bOriginalTransferOp =
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()->getLocalElementList();
221 stridedRgMap = StridedMapFactory::Build(
222 originalTransferOp->getRangeMap()->lib(),
223 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
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()->getLocalElementList();
237 stridedDoMap = StridedMapFactory::Build(
238 originalTransferOp->getDomainMap()->lib(),
239 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
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()->getLocalElementList();
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()->getLocalElementList();
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(),
325 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
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(),
336 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
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(),
351 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
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(),
363 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
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_ */
#define SET_VALID_ENTRY(name)
Exception indicating invalid cast attempted.
Exception throws to report errors in the internal logical of the program.
Timer to be used in factories. Similar to Monitor but with additional timers.
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access)....
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
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.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
void AddFactoryManager(RCP< const FactoryManagerBase > FactManager)
Add a factory manager.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
An exception safe way to call the method 'Level::SetFactoryManager()'.
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
Namespace for MueLu classes and methods.
@ Statistics0
Print statistics that do not involve significant additional computation.