MueLu Version of the Day
Loading...
Searching...
No Matches
MueLu_RebalanceTransferFactory_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_REBALANCETRANSFERFACTORY_DEF_HPP
47#define MUELU_REBALANCETRANSFERFACTORY_DEF_HPP
48
49#include <sstream>
50#include <Teuchos_Tuple.hpp>
51
52#include "Xpetra_MultiVector.hpp"
53#include "Xpetra_MultiVectorFactory.hpp"
54#include "Xpetra_Vector.hpp"
55#include "Xpetra_VectorFactory.hpp"
56#include <Xpetra_Matrix.hpp>
57#include <Xpetra_MapFactory.hpp>
58#include <Xpetra_MatrixFactory.hpp>
59#include <Xpetra_Import.hpp>
60#include <Xpetra_ImportFactory.hpp>
61#include <Xpetra_IO.hpp>
62
64
65#include "MueLu_Level.hpp"
66#include "MueLu_MasterList.hpp"
67#include "MueLu_Monitor.hpp"
68#include "MueLu_PerfUtils.hpp"
69
70namespace MueLu {
71
72 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
74 RCP<ParameterList> validParamList = rcp(new ParameterList());
75
76#define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
77 SET_VALID_ENTRY("repartition: rebalance P and R");
78 SET_VALID_ENTRY("repartition: rebalance Nullspace");
79 SET_VALID_ENTRY("transpose: use implicit");
80 SET_VALID_ENTRY("repartition: use subcommunicators");
81#undef SET_VALID_ENTRY
82
83 {
84 typedef Teuchos::StringToIntegralParameterEntryValidator<int> validatorType;
85 RCP<validatorType> typeValidator = rcp (new validatorType(Teuchos::tuple<std::string>("Interpolation", "Restriction"), "type"));
86 validParamList->set("type", "Interpolation", "Type of the transfer operator that need to be rebalanced (Interpolation or Restriction)", typeValidator);
87 }
88
89 validParamList->set< RCP<const FactoryBase> >("P", null, "Factory of the prolongation operator that need to be rebalanced (only used if type=Interpolation)");
90 validParamList->set< RCP<const FactoryBase> >("R", null, "Factory of the restriction operator that need to be rebalanced (only used if type=Restriction)");
91 validParamList->set< RCP<const FactoryBase> >("Nullspace", null, "Factory of the nullspace that need to be rebalanced (only used if type=Interpolation)");
92 validParamList->set< RCP<const FactoryBase> >("Coordinates", null, "Factory of the coordinates that need to be rebalanced (only used if type=Interpolation)");
93 validParamList->set< RCP<const FactoryBase> >("BlockNumber", null, "Factory of the block ids that need to be rebalanced (only used if type=Interpolation)");
94 validParamList->set< RCP<const FactoryBase> >("Importer", null, "Factory of the importer object used for the rebalancing");
95 validParamList->set< int > ("write start", -1, "First level at which coordinates should be written to file");
96 validParamList->set< int > ("write end", -1, "Last level at which coordinates should be written to file");
97
98 // TODO validation: "P" parameter valid only for type="Interpolation" and "R" valid only for type="Restriction". Like so:
99 // if (paramList.isEntry("type") && paramList.get("type) == "Interpolation) {
100 // validParamList->set< RCP<const FactoryBase> >("P", Teuchos::null, "Factory of the prolongation operator that need to be rebalanced (only used if type=Interpolation)");
101
102 return validParamList;
103 }
104
105 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
107 const ParameterList& pL = GetParameterList();
108
109 if (pL.get<std::string>("type") == "Interpolation") {
110 Input(coarseLevel, "P");
111 if (pL.get<bool>("repartition: rebalance Nullspace"))
112 Input(coarseLevel, "Nullspace");
113 if (pL.get< RCP<const FactoryBase> >("Coordinates") != Teuchos::null)
114 Input(coarseLevel, "Coordinates");
115 if (pL.get< RCP<const FactoryBase> >("BlockNumber") != Teuchos::null)
116 Input(coarseLevel, "BlockNumber");
117
118 } else {
119 if (pL.get<bool>("transpose: use implicit") == false)
120 Input(coarseLevel, "R");
121 }
122
123 Input(coarseLevel, "Importer");
124 }
125
126 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
128 FactoryMonitor m(*this, "Build", coarseLevel);
129 typedef Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType, LO, GO, NO> xdMV;
130
131 const ParameterList& pL = GetParameterList();
132
133 RCP<Matrix> originalP = Get< RCP<Matrix> >(coarseLevel, "P");
134 // If we don't have a valid P (e.g., # global aggregates is 0), skip this rebalancing. This level will
135 // ultimately be removed in MueLu_Hierarchy_defs.h via a resize()
136 if (originalP == Teuchos::null) {
137 Set(coarseLevel, "P", originalP);
138 return;
139 }
140 int implicit = !pL.get<bool>("repartition: rebalance P and R");
141 int writeStart = pL.get<int> ("write start");
142 int writeEnd = pL.get<int> ("write end");
143
144 if (writeStart == 0 && fineLevel.GetLevelID() == 0 && writeStart <= writeEnd && IsAvailable(fineLevel, "Coordinates")) {
145 std::string fileName = "coordinates_level_0.m";
146 RCP<xdMV> fineCoords = fineLevel.Get< RCP<xdMV> >("Coordinates");
147 if (fineCoords != Teuchos::null)
148 Xpetra::IO<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO>::Write(fileName, *fineCoords);
149 }
150
151 if (writeStart == 0 && fineLevel.GetLevelID() == 0 && writeStart <= writeEnd && IsAvailable(fineLevel, "BlockNumber")) {
152 std::string fileName = "BlockNumber_level_0.m";
153 RCP<LocalOrdinalVector> fineBlockNumber = fineLevel.Get< RCP<LocalOrdinalVector> >("BlockNumber");
154 if (fineBlockNumber != Teuchos::null)
155 Xpetra::IO<LO,LO,GO,NO>::Write(fileName, *fineBlockNumber);
156 }
157
158 RCP<const Import> importer = Get<RCP<const Import> >(coarseLevel, "Importer");
159 if (implicit) {
160 // Save the importer, we'll need it for solve
161 coarseLevel.Set("Importer", importer, NoFactory::get());
162 }
163
164 RCP<ParameterList> params = rcp(new ParameterList());
165 if (IsPrint(Statistics2)) {
166 params->set("printLoadBalancingInfo", true);
167 params->set("printCommInfo", true);
168 }
169
170 std::string transferType = pL.get<std::string>("type");
171 if (transferType == "Interpolation") {
172 originalP = Get< RCP<Matrix> >(coarseLevel, "P");
173
174 {
175 // This line must be after the Get call
176 SubFactoryMonitor m1(*this, "Rebalancing prolongator", coarseLevel);
177
178 if (implicit || importer.is_null()) {
179 GetOStream(Runtime0) << "Using original prolongator" << std::endl;
180 Set(coarseLevel, "P", originalP);
181
182 } else {
183 // P is the transfer operator from the coarse grid to the fine grid.
184 // P must transfer the data from the newly reordered coarse A to the
185 // (unchanged) fine A. This means that the domain map (coarse) of P
186 // must be changed according to the new partition. The range map
187 // (fine) is kept unchanged.
188 //
189 // The domain map of P must match the range map of R. See also note
190 // below about domain/range map of R and its implications for P.
191 //
192 // To change the domain map of P, P needs to be fillCompleted again
193 // with the new domain map. To achieve this, P is copied into a new
194 // matrix that is not fill-completed. The doImport() operation is
195 // just used here to make a copy of P: the importer is trivial and
196 // there is no data movement involved. The reordering actually
197 // happens during the fillComplete() with domainMap == importer->getTargetMap().
198 RCP<Matrix> rebalancedP = originalP;
199 RCP<const CrsMatrixWrap> crsOp = rcp_dynamic_cast<const CrsMatrixWrap>(originalP);
200 TEUCHOS_TEST_FOR_EXCEPTION(crsOp == Teuchos::null, Exceptions::BadCast, "Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
201
202 RCP<CrsMatrix> rebalancedP2 = crsOp->getCrsMatrix();
203 TEUCHOS_TEST_FOR_EXCEPTION(rebalancedP2 == Teuchos::null, std::runtime_error, "Xpetra::CrsMatrixWrap doesn't have a CrsMatrix");
204
205 {
206 SubFactoryMonitor subM(*this, "Rebalancing prolongator -- fast map replacement", coarseLevel);
207
208 RCP<const Import> newImporter;
209 {
210 SubFactoryMonitor(*this, "Import construction", coarseLevel);
211 newImporter = ImportFactory::Build(importer->getTargetMap(), rebalancedP->getColMap());
212 }
213 rebalancedP2->replaceDomainMapAndImporter(importer->getTargetMap(), newImporter);
214 }
215
217 // TODO FIXME somehow we have to transfer the striding information of the permuted domain/range maps.
218 // That is probably something for an external permutation factory
219 // if (originalP->IsView("stridedMaps"))
220 // rebalancedP->CreateView("stridedMaps", originalP);
222 if(!rebalancedP.is_null()) {std::ostringstream oss; oss << "P_" << coarseLevel.GetLevelID(); rebalancedP->setObjectLabel(oss.str());}
223 Set(coarseLevel, "P", rebalancedP);
224
225 if (IsPrint(Statistics2))
226 GetOStream(Statistics2) << PerfUtils::PrintMatrixInfo(*rebalancedP, "P (rebalanced)", params);
227 }
228 }
229
230 if (importer.is_null()) {
231 if (IsAvailable(coarseLevel, "Nullspace"))
232 Set(coarseLevel, "Nullspace", Get<RCP<MultiVector> >(coarseLevel, "Nullspace"));
233
234 if (pL.isParameter("Coordinates") && pL.get< RCP<const FactoryBase> >("Coordinates") != Teuchos::null)
235 if (IsAvailable(coarseLevel, "Coordinates"))
236 Set(coarseLevel, "Coordinates", Get< RCP<xdMV> >(coarseLevel, "Coordinates"));
237
238 if (pL.isParameter("BlockNumber") && pL.get< RCP<const FactoryBase> >("BlockNumber") != Teuchos::null)
239 if (IsAvailable(coarseLevel, "BlockNumber"))
240 Set(coarseLevel, "BlockNumber", Get< RCP<LocalOrdinalVector> >(coarseLevel, "BlockNumber"));
241
242 return;
243 }
244
245 if (pL.isParameter("Coordinates") &&
246 pL.get< RCP<const FactoryBase> >("Coordinates") != Teuchos::null &&
247 IsAvailable(coarseLevel, "Coordinates")) {
248 RCP<xdMV> coords = Get<RCP<xdMV> >(coarseLevel, "Coordinates");
249
250 // This line must be after the Get call
251 SubFactoryMonitor subM(*this, "Rebalancing coordinates", coarseLevel);
252
253 LO nodeNumElts = coords->getMap()->getLocalNumElements();
254
255 // If a process has no matrix rows, then we can't calculate blocksize using the formula below.
256 LO myBlkSize = 0, blkSize = 0;
257 if (nodeNumElts > 0)
258 myBlkSize = importer->getSourceMap()->getLocalNumElements() / nodeNumElts;
259 MueLu_maxAll(coords->getMap()->getComm(), myBlkSize, blkSize);
260
261 RCP<const Import> coordImporter;
262 if (blkSize == 1) {
263 coordImporter = importer;
264
265 } else {
266 // NOTE: there is an implicit assumption here: we assume that dof any node are enumerated consequently
267 // Proper fix would require using decomposition similar to how we construct importer in the
268 // RepartitionFactory
269 RCP<const Map> origMap = coords->getMap();
270 GO indexBase = origMap->getIndexBase();
271
272 ArrayView<const GO> OEntries = importer->getTargetMap()->getLocalElementList();
273 LO numEntries = OEntries.size()/blkSize;
274 ArrayRCP<GO> Entries(numEntries);
275 for (LO i = 0; i < numEntries; i++)
276 Entries[i] = (OEntries[i*blkSize]-indexBase)/blkSize + indexBase;
277
278 RCP<const Map> targetMap = MapFactory::Build(origMap->lib(), origMap->getGlobalNumElements(), Entries(), indexBase, origMap->getComm());
279 coordImporter = ImportFactory::Build(origMap, targetMap);
280 }
281
282 RCP<xdMV> permutedCoords = Xpetra::MultiVectorFactory<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO>::Build(coordImporter->getTargetMap(), coords->getNumVectors());
283 permutedCoords->doImport(*coords, *coordImporter, Xpetra::INSERT);
284
285 if (pL.isParameter("repartition: use subcommunicators") == true && pL.get<bool>("repartition: use subcommunicators") == true)
286 permutedCoords->replaceMap(permutedCoords->getMap()->removeEmptyProcesses());
287
288 if (permutedCoords->getMap() == Teuchos::null)
289 permutedCoords = Teuchos::null;
290
291 Set(coarseLevel, "Coordinates", permutedCoords);
292
293 std::string fileName = "rebalanced_coordinates_level_" + toString(coarseLevel.GetLevelID()) + ".m";
294 if (writeStart <= coarseLevel.GetLevelID() && coarseLevel.GetLevelID() <= writeEnd && permutedCoords->getMap() != Teuchos::null)
295 Xpetra::IO<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO>::Write(fileName, *permutedCoords);
296 }
297
298 if (pL.isParameter("BlockNumber") &&
299 pL.get< RCP<const FactoryBase> >("BlockNumber") != Teuchos::null &&
300 IsAvailable(coarseLevel, "BlockNumber")) {
301 RCP<LocalOrdinalVector> BlockNumber = Get<RCP<LocalOrdinalVector> >(coarseLevel, "BlockNumber");
302
303 // This line must be after the Get call
304 SubFactoryMonitor subM(*this, "Rebalancing BlockNumber", coarseLevel);
305
306 RCP<LocalOrdinalVector> permutedBlockNumber = LocalOrdinalVectorFactory::Build(importer->getTargetMap(), false);
307 permutedBlockNumber->doImport(*BlockNumber, *importer, Xpetra::INSERT);
308
309 if (pL.isParameter("repartition: use subcommunicators") == true && pL.get<bool>("repartition: use subcommunicators") == true)
310 permutedBlockNumber->replaceMap(permutedBlockNumber->getMap()->removeEmptyProcesses());
311
312 if (permutedBlockNumber->getMap() == Teuchos::null)
313 permutedBlockNumber = Teuchos::null;
314
315 Set(coarseLevel, "BlockNumber", permutedBlockNumber);
316
317 std::string fileName = "rebalanced_BlockNumber_level_" + toString(coarseLevel.GetLevelID()) + ".m";
318 if (writeStart <= coarseLevel.GetLevelID() && coarseLevel.GetLevelID() <= writeEnd && permutedBlockNumber->getMap() != Teuchos::null)
319 Xpetra::IO<LO,LO,GO,NO>::Write(fileName, *permutedBlockNumber);
320 }
321
322 if (IsAvailable(coarseLevel, "Nullspace")) {
323 RCP<MultiVector> nullspace = Get< RCP<MultiVector> >(coarseLevel, "Nullspace");
324
325 // This line must be after the Get call
326 SubFactoryMonitor subM(*this, "Rebalancing nullspace", coarseLevel);
327
328 RCP<MultiVector> permutedNullspace = MultiVectorFactory::Build(importer->getTargetMap(), nullspace->getNumVectors());
329 permutedNullspace->doImport(*nullspace, *importer, Xpetra::INSERT);
330
331 if (pL.get<bool>("repartition: use subcommunicators") == true)
332 permutedNullspace->replaceMap(permutedNullspace->getMap()->removeEmptyProcesses());
333
334 if (permutedNullspace->getMap() == Teuchos::null)
335 permutedNullspace = Teuchos::null;
336
337 Set(coarseLevel, "Nullspace", permutedNullspace);
338 }
339
340 } else {
341 if (pL.get<bool>("transpose: use implicit") == false) {
342 RCP<Matrix> originalR = Get< RCP<Matrix> >(coarseLevel, "R");
343
344 SubFactoryMonitor m2(*this, "Rebalancing restrictor", coarseLevel);
345
346 if (implicit || importer.is_null()) {
347 GetOStream(Runtime0) << "Using original restrictor" << std::endl;
348 Set(coarseLevel, "R", originalR);
349
350 } else {
351 RCP<Matrix> rebalancedR;
352 {
353 SubFactoryMonitor subM(*this, "Rebalancing restriction -- fusedImport", coarseLevel);
354
355 RCP<Map> dummy; // meaning: use originalR's domain map.
356 Teuchos::ParameterList listLabel;
357 listLabel.set("Timer Label","MueLu::RebalanceR-" + Teuchos::toString(coarseLevel.GetLevelID()));
358 rebalancedR = MatrixFactory::Build(originalR, *importer, dummy, importer->getTargetMap(),Teuchos::rcp(&listLabel,false));
359 }
360 if(!rebalancedR.is_null()) {std::ostringstream oss; oss << "R_" << coarseLevel.GetLevelID(); rebalancedR->setObjectLabel(oss.str());}
361 Set(coarseLevel, "R", rebalancedR);
362
364 // TODO FIXME somehow we have to transfer the striding information of the permuted domain/range maps.
365 // That is probably something for an external permutation factory
366 // if (originalR->IsView("stridedMaps"))
367 // rebalancedR->CreateView("stridedMaps", originalR);
369
370 if (IsPrint(Statistics2))
371 GetOStream(Statistics2) << PerfUtils::PrintMatrixInfo(*rebalancedR, "R (rebalanced)", params);
372 }
373 }
374 }
375 }
376
377} // namespace MueLu
378
379#endif // MUELU_REBALANCETRANSFERFACTORY_DEF_HPP
#define SET_VALID_ENTRY(name)
#define MueLu_maxAll(rcpComm, in, out)
Exception indicating invalid cast attempted.
Timer to be used in factories. Similar to Monitor but with additional timers.
Class that holds all level-specific information.
int GetLevelID() const
Return level number.
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 const NoFactory * get()
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
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.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
Namespace for MueLu classes and methods.
@ Statistics2
Print even more statistics.
@ Runtime0
One-liner description of what is happening.
std::string toString(const T &what)
Little helper function to convert non-string types to strings.