MueLu Version of the Day
MueLu_Hierarchy_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
47#ifndef MUELU_HIERARCHY_DEF_HPP
48#define MUELU_HIERARCHY_DEF_HPP
49
50#include <time.h>
51
52#include <algorithm>
53#include <sstream>
54
55#include <Xpetra_Matrix.hpp>
56#include <Xpetra_MultiVectorFactory.hpp>
57#include <Xpetra_Operator.hpp>
58#include <Xpetra_IO.hpp>
59
60#include "MueLu_Hierarchy_decl.hpp"
61
62#include "MueLu_BoostGraphviz.hpp"
63#include "MueLu_FactoryManager.hpp"
64#include "MueLu_HierarchyUtils.hpp"
65#include "MueLu_TopRAPFactory.hpp"
66#include "MueLu_TopSmootherFactory.hpp"
67#include "MueLu_Level.hpp"
68#include "MueLu_Monitor.hpp"
69#include "MueLu_PerfUtils.hpp"
70#include "MueLu_PFactory.hpp"
71#include "MueLu_SmootherFactoryBase.hpp"
72#include "MueLu_SmootherFactory.hpp"
73#include "MueLu_SmootherBase.hpp"
74
75#include "Teuchos_TimeMonitor.hpp"
76
77
78
79namespace MueLu {
80
81 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
82 Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Hierarchy()
83 : maxCoarseSize_(GetDefaultMaxCoarseSize()), implicitTranspose_(GetDefaultImplicitTranspose()),
84 fuseProlongationAndUpdate_(GetDefaultFuseProlongationAndUpdate()),
85 doPRrebalance_(GetDefaultPRrebalance()), isPreconditioner_(true), Cycle_(GetDefaultCycle()), WCycleStartLevel_(0),
86 scalingFactor_(Teuchos::ScalarTraits<double>::one()), lib_(Xpetra::UseTpetra), isDumpingEnabled_(false), dumpLevel_(-2), rate_(-1),
87 sizeOfAllocatedLevelMultiVectors_(0)
88 {
89 AddLevel(rcp(new Level));
90 }
91
92 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
93 Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Hierarchy(const std::string& label)
94 : Hierarchy()
95 {
96 setObjectLabel(label);
97 Levels_[0]->setObjectLabel(label);
98 }
99
100 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
101 Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Hierarchy(const RCP<Matrix>& A)
102 : maxCoarseSize_(GetDefaultMaxCoarseSize()), implicitTranspose_(GetDefaultImplicitTranspose()),
103 fuseProlongationAndUpdate_(GetDefaultFuseProlongationAndUpdate()),
104 doPRrebalance_(GetDefaultPRrebalance()), isPreconditioner_(true), Cycle_(GetDefaultCycle()), WCycleStartLevel_(0),
105 scalingFactor_(Teuchos::ScalarTraits<double>::one()), isDumpingEnabled_(false), dumpLevel_(-2), rate_(-1),
106 sizeOfAllocatedLevelMultiVectors_(0)
107 {
108 lib_ = A->getDomainMap()->lib();
109
110 RCP<Level> Finest = rcp(new Level);
111 AddLevel(Finest);
112
113 Finest->Set("A", A);
114 }
115
116 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
117 Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Hierarchy(const RCP<Matrix>& A, const std::string& label)
118 : Hierarchy(A)
119 {
120 setObjectLabel(label);
121 Levels_[0]->setObjectLabel(label);
122 }
123
124 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
125 void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::AddLevel(const RCP<Level>& level) {
126 int levelID = LastLevelID() + 1; // ID of the inserted level
127
128 if (level->GetLevelID() != -1 && (level->GetLevelID() != levelID))
129 GetOStream(Warnings1) << "Hierarchy::AddLevel(): Level with ID=" << level->GetLevelID() <<
130 " have been added at the end of the hierarchy\n but its ID have been redefined" <<
131 " because last level ID of the hierarchy was " << LastLevelID() << "." << std::endl;
132
133 Levels_.push_back(level);
134 level->SetLevelID(levelID);
135 level->setlib(lib_);
136
137 level->SetPreviousLevel( (levelID == 0) ? Teuchos::null : Levels_[LastLevelID() - 1] );
138 level->setObjectLabel(this->getObjectLabel());
139 }
140
141 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
142 void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::AddNewLevel() {
143 RCP<Level> newLevel = Levels_[LastLevelID()]->Build(); // new coarse level, using copy constructor
144 newLevel->setlib(lib_);
145 this->AddLevel(newLevel); // add to hierarchy
146 }
147
148 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
149 RCP<Level> & Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetLevel(const int levelID) {
150 TEUCHOS_TEST_FOR_EXCEPTION(levelID < 0 || levelID > LastLevelID(), Exceptions::RuntimeError,
151 "MueLu::Hierarchy::GetLevel(): invalid input parameter value: LevelID = " << levelID);
152 return Levels_[levelID];
153 }
154
155 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
156 int Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetNumLevels() const {
157 return Levels_.size();
158 }
159
160 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
161 int Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetGlobalNumLevels() const {
162 RCP<Operator> A = Levels_[0]->template Get<RCP<Operator> >("A");
163 RCP<const Teuchos::Comm<int> > comm = A->getDomainMap()->getComm();
164
165 int numLevels = GetNumLevels();
166 int numGlobalLevels;
167 Teuchos::reduceAll(*comm, Teuchos::REDUCE_MAX, numLevels, Teuchos::ptr(&numGlobalLevels));
168
169 return numGlobalLevels;
170 }
171
172 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
173 double Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetOperatorComplexity() const {
174 double totalNnz = 0, lev0Nnz = 1;
175 for (int i = 0; i < GetNumLevels(); ++i) {
176 TEUCHOS_TEST_FOR_EXCEPTION(!(Levels_[i]->IsAvailable("A")) , Exceptions::RuntimeError,
177 "Operator complexity cannot be calculated because A is unavailable on level " << i);
178 RCP<Operator> A = Levels_[i]->template Get<RCP<Operator> >("A");
179 if (A.is_null())
180 break;
181
182 RCP<Matrix> Am = rcp_dynamic_cast<Matrix>(A);
183 if (Am.is_null()) {
184 GetOStream(Warnings0) << "Some level operators are not matrices, operator complexity calculation aborted" << std::endl;
185 return 0.0;
186 }
187
188 totalNnz += as<double>(Am->getGlobalNumEntries());
189 if (i == 0)
190 lev0Nnz = totalNnz;
191 }
192 return totalNnz / lev0Nnz;
193 }
194
195 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
196 double Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetSmootherComplexity() const {
197 double node_sc = 0, global_sc=0;
198 double a0_nnz =0;
199 const size_t INVALID = Teuchos::OrdinalTraits<size_t>::invalid();
200 // Get cost of fine matvec
201 if (GetNumLevels() <= 0) return -1.0;
202 if (!Levels_[0]->IsAvailable("A")) return -1.0;
203
204 RCP<Operator> A = Levels_[0]->template Get<RCP<Operator> >("A");
205 if (A.is_null()) return -1.0;
206 RCP<Matrix> Am = rcp_dynamic_cast<Matrix>(A);
207 if(Am.is_null()) return -1.0;
208 a0_nnz = as<double>(Am->getGlobalNumEntries());
209
210 // Get smoother complexity at each level
211 for (int i = 0; i < GetNumLevels(); ++i) {
212 size_t level_sc=0;
213 if(!Levels_[i]->IsAvailable("PreSmoother")) continue;
214 RCP<SmootherBase> S = Levels_[i]->template Get<RCP<SmootherBase> >("PreSmoother");
215 if (S.is_null()) continue;
216 level_sc = S->getNodeSmootherComplexity();
217 if(level_sc == INVALID) {global_sc=-1.0;break;}
218
219 node_sc += as<double>(level_sc);
220 }
221
222 double min_sc=0.0;
223 RCP<const Teuchos::Comm<int> > comm =A->getDomainMap()->getComm();
224 Teuchos::reduceAll(*comm,Teuchos::REDUCE_SUM,node_sc,Teuchos::ptr(&global_sc));
225 Teuchos::reduceAll(*comm,Teuchos::REDUCE_MIN,node_sc,Teuchos::ptr(&min_sc));
226
227 if(min_sc < 0.0) return -1.0;
228 else return global_sc / a0_nnz;
229 }
230
231
232
233
234 // Coherence checks todo in Setup() (using an helper function):
235 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
236 void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::CheckLevel(Level& level, int levelID) {
237 TEUCHOS_TEST_FOR_EXCEPTION(level.lib() != lib_, Exceptions::RuntimeError,
238 "MueLu::Hierarchy::CheckLevel(): wrong underlying linear algebra library.");
239 TEUCHOS_TEST_FOR_EXCEPTION(level.GetLevelID() != levelID, Exceptions::RuntimeError,
240 "MueLu::Hierarchy::CheckLevel(): wrong level ID");
241 TEUCHOS_TEST_FOR_EXCEPTION(levelID != 0 && level.GetPreviousLevel() != Levels_[levelID-1], Exceptions::RuntimeError,
242 "MueLu::Hierarchy::Setup(): wrong level parent");
243 }
244
245 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
246 void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::SetMatvecParams(RCP<ParameterList> matvecParams) {
247 for (int i = 0; i < GetNumLevels(); ++i) {
248 RCP<Level> level = Levels_[i];
249 if (level->IsAvailable("A")) {
250 RCP<Operator> Aop = level->Get<RCP<Operator> >("A");
251 RCP<Matrix> A = rcp_dynamic_cast<Matrix>(Aop);
252 if (!A.is_null()) {
253 RCP<const Import> xpImporter = A->getCrsGraph()->getImporter();
254 if (!xpImporter.is_null())
255 xpImporter->setDistributorParameters(matvecParams);
256 RCP<const Export> xpExporter = A->getCrsGraph()->getExporter();
257 if (!xpExporter.is_null())
258 xpExporter->setDistributorParameters(matvecParams);
259 }
260 }
261 if (level->IsAvailable("P")) {
262 RCP<Matrix> P = level->Get<RCP<Matrix> >("P");
263 RCP<const Import> xpImporter = P->getCrsGraph()->getImporter();
264 if (!xpImporter.is_null())
265 xpImporter->setDistributorParameters(matvecParams);
266 RCP<const Export> xpExporter = P->getCrsGraph()->getExporter();
267 if (!xpExporter.is_null())
268 xpExporter->setDistributorParameters(matvecParams);
269 }
270 if (level->IsAvailable("R")) {
271 RCP<Matrix> R = level->Get<RCP<Matrix> >("R");
272 RCP<const Import> xpImporter = R->getCrsGraph()->getImporter();
273 if (!xpImporter.is_null())
274 xpImporter->setDistributorParameters(matvecParams);
275 RCP<const Export> xpExporter = R->getCrsGraph()->getExporter();
276 if (!xpExporter.is_null())
277 xpExporter->setDistributorParameters(matvecParams);
278 }
279 if (level->IsAvailable("Importer")) {
280 RCP<const Import> xpImporter = level->Get< RCP<const Import> >("Importer");
281 if (!xpImporter.is_null())
282 xpImporter->setDistributorParameters(matvecParams);
283 }
284 }
285 }
286
287 // The function uses three managers: fine, coarse and next coarse
288 // We construct the data for the coarse level, and do requests for the next coarse
289 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
290 bool Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Setup(int coarseLevelID,
291 const RCP<const FactoryManagerBase> fineLevelManager,
292 const RCP<const FactoryManagerBase> coarseLevelManager,
293 const RCP<const FactoryManagerBase> nextLevelManager) {
294 // Use PrintMonitor/TimerMonitor instead of just a FactoryMonitor to print "Level 0" instead of Hierarchy(0)
295 // Print is done after the requests for next coarse level
296
297 TEUCHOS_TEST_FOR_EXCEPTION(LastLevelID() < coarseLevelID, Exceptions::RuntimeError,
298 "MueLu::Hierarchy:Setup(): level " << coarseLevelID << " (specified by coarseLevelID argument) "
299 "must be built before calling this function.");
300
301 Level& level = *Levels_[coarseLevelID];
302
303 std::string label = FormattingHelper::getColonLabel(level.getObjectLabel());
304 TimeMonitor m1(*this, label + this->ShortClassName() + ": " + "Setup (total)");
305 TimeMonitor m2(*this, label + this->ShortClassName() + ": " + "Setup" + " (total, level=" + Teuchos::toString(coarseLevelID) + ")");
306
307 // TODO: pass coarseLevelManager by reference
308 TEUCHOS_TEST_FOR_EXCEPTION(coarseLevelManager == Teuchos::null, Exceptions::RuntimeError,
309 "MueLu::Hierarchy::Setup(): argument coarseLevelManager cannot be null");
310
311 typedef MueLu::TopRAPFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node> TopRAPFactory;
312 typedef MueLu::TopSmootherFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node> TopSmootherFactory;
313
314 if (levelManagers_.size() < coarseLevelID+1)
315 levelManagers_.resize(coarseLevelID+1);
316 levelManagers_[coarseLevelID] = coarseLevelManager;
317
318 bool isFinestLevel = (fineLevelManager.is_null());
319 bool isLastLevel = (nextLevelManager.is_null());
320
321 int oldRank = -1;
322 if (isFinestLevel) {
323 RCP<Operator> A = level.Get< RCP<Operator> >("A");
324 RCP<const Map> domainMap = A->getDomainMap();
325 RCP<const Teuchos::Comm<int> > comm = domainMap->getComm();
326
327 // Initialize random seed for reproducibility
328 Utilities::SetRandomSeed(*comm);
329
330 // Record the communicator on the level (used for timers sync)
331 level.SetComm(comm);
332 oldRank = SetProcRankVerbose(comm->getRank());
333
334 // Set the Hierarchy library to match that of the finest level matrix,
335 // even if it was already set
336 lib_ = domainMap->lib();
337 level.setlib(lib_);
338
339 } else {
340 // Permeate library to a coarser level
341 level.setlib(lib_);
342
343 Level& prevLevel = *Levels_[coarseLevelID-1];
344 oldRank = SetProcRankVerbose(prevLevel.GetComm()->getRank());
345 }
346
347 CheckLevel(level, coarseLevelID);
348
349 // Attach FactoryManager to the fine level
350 RCP<SetFactoryManager> SFMFine;
351 if (!isFinestLevel)
352 SFMFine = rcp(new SetFactoryManager(Levels_[coarseLevelID-1], fineLevelManager));
353
354 if (isFinestLevel && Levels_[coarseLevelID]->IsAvailable("Coordinates"))
355 ReplaceCoordinateMap(*Levels_[coarseLevelID]);
356
357 // Attach FactoryManager to the coarse level
358 SetFactoryManager SFMCoarse(Levels_[coarseLevelID], coarseLevelManager);
359
360 if (isDumpingEnabled_ && (dumpLevel_ == 0 || dumpLevel_ == -1) && coarseLevelID == 1)
361 DumpCurrentGraph(0);
362
363 RCP<TopSmootherFactory> coarseFact = rcp(new TopSmootherFactory(coarseLevelManager, "CoarseSolver"));
364 RCP<TopSmootherFactory> smootherFact = rcp(new TopSmootherFactory(coarseLevelManager, "Smoother"));
365
366 int nextLevelID = coarseLevelID + 1;
367
368 RCP<SetFactoryManager> SFMNext;
369 if (isLastLevel == false) {
370 // We are not at the coarsest level, so there is going to be another level ("next coarse") after this one ("coarse")
371 if (nextLevelID > LastLevelID())
372 AddNewLevel();
373 CheckLevel(*Levels_[nextLevelID], nextLevelID);
374
375 // Attach FactoryManager to the next level (level after coarse)
376 SFMNext = rcp(new SetFactoryManager(Levels_[nextLevelID], nextLevelManager));
377 Levels_[nextLevelID]->Request(TopRAPFactory(coarseLevelManager, nextLevelManager));
378
379 // Do smoother requests here. We don't know whether this is going to be
380 // the coarsest level or not, but we need to DeclareInput before we call
381 // coarseRAPFactory.Build(), otherwise some stuff may be erased after
382 // level releases
383 level.Request(*smootherFact);
384
385 } else {
386 // Similar to smoother above, do the coarse solver request here. We don't
387 // know whether this is going to be the coarsest level or not, but we
388 // need to DeclareInput before we call coarseRAPFactory.Build(),
389 // otherwise some stuff may be erased after level releases. This is
390 // actually evident on ProjectorSmoother. It requires both "A" and
391 // "Nullspace". However, "Nullspace" is erased after all releases, so if
392 // we call the coarse factory request after RAP build we would not have
393 // any data, and cannot get it as we don't have previous managers. The
394 // typical trace looks like this:
395 //
396 // MueLu::Level(0)::GetFactory(Aggregates, 0): No FactoryManager
397 // during request for data " Aggregates" on level 0 by factory TentativePFactory
398 // during request for data " P" on level 1 by factory EminPFactory
399 // during request for data " P" on level 1 by factory TransPFactory
400 // during request for data " R" on level 1 by factory RAPFactory
401 // during request for data " A" on level 1 by factory TentativePFactory
402 // during request for data " Nullspace" on level 2 by factory NullspaceFactory
403 // during request for data " Nullspace" on level 2 by factory NullspacePresmoothFactory
404 // during request for data " Nullspace" on level 2 by factory ProjectorSmoother
405 // during request for data " PreSmoother" on level 2 by factory NoFactory
406 level.Request(*coarseFact);
407 }
408
409 PrintMonitor m0(*this, "Level " + Teuchos::toString(coarseLevelID), static_cast<MsgType>(Runtime0 | Test));
410
411 // Build coarse level hierarchy
412 RCP<Operator> Ac = Teuchos::null;
413 TopRAPFactory coarseRAPFactory(fineLevelManager, coarseLevelManager);
414
415 if (level.IsAvailable("A")) {
416 Ac = level.Get<RCP<Operator> >("A");
417 } else if (!isFinestLevel) {
418 // We only build here, the release is done later
419 coarseRAPFactory.Build(*level.GetPreviousLevel(), level);
420 }
421
422 if (level.IsAvailable("A"))
423 Ac = level.Get<RCP<Operator> >("A");
424 RCP<Matrix> Acm = rcp_dynamic_cast<Matrix>(Ac);
425
426 // Record the communicator on the level
427 if (!Ac.is_null())
428 level.SetComm(Ac->getDomainMap()->getComm());
429
430 // Test if we reach the end of the hierarchy
431 bool isOrigLastLevel = isLastLevel;
432 if (isLastLevel) {
433 // Last level as we have achieved the max limit
434 isLastLevel = true;
435
436 } else if (Ac.is_null()) {
437 // Last level for this processor, as it does not belong to the next
438 // subcommunicator. Other processors may continue working on the
439 // hierarchy
440 isLastLevel = true;
441
442 } else {
443 if (!Acm.is_null() && Acm->getGlobalNumRows() <= maxCoarseSize_) {
444 // Last level as the size of the coarse matrix became too small
445 GetOStream(Runtime0) << "Max coarse size (<= " << maxCoarseSize_ << ") achieved" << std::endl;
446 isLastLevel = true;
447 }
448 }
449
450 if (!Ac.is_null() && !isFinestLevel) {
451 RCP<Operator> A = Levels_[coarseLevelID-1]->template Get< RCP<Operator> >("A");
452 RCP<Matrix> Am = rcp_dynamic_cast<Matrix>(A);
453
454 const double maxCoarse2FineRatio = 0.8;
455 if (!Acm.is_null() && !Am.is_null() && Acm->getGlobalNumRows() > maxCoarse2FineRatio * Am->getGlobalNumRows()) {
456 // We could abort here, but for now we simply notify user.
457 // Couple of additional points:
458 // - if repartitioning is delayed until level K, but the aggregation
459 // procedure stagnates between levels K-1 and K. In this case,
460 // repartitioning could enable faster coarsening once again, but the
461 // hierarchy construction will abort due to the stagnation check.
462 // - if the matrix is small enough, we could move it to one processor.
463 GetOStream(Warnings0) << "Aggregation stagnated. Please check your matrix and/or adjust your configuration file."
464 << "Possible fixes:\n"
465 << " - reduce the maximum number of levels\n"
466 << " - enable repartitioning\n"
467 << " - increase the minimum coarse size." << std::endl;
468
469 }
470 }
471
472 if (isLastLevel) {
473 if (!isOrigLastLevel) {
474 // We did not expect to finish this early so we did request a smoother.
475 // We need a coarse solver instead. Do the magic.
476 level.Release(*smootherFact);
477 level.Request(*coarseFact);
478 }
479
480 // Do the actual build, if we have any data.
481 // NOTE: this is not a great check, we may want to call Build() regardless.
482 if (!Ac.is_null())
483 coarseFact->Build(level);
484
485 // Once the dirty deed is done, release stuff. The smoother has already
486 // been released.
487 level.Release(*coarseFact);
488
489 } else {
490 // isLastLevel = false => isOrigLastLevel = false, meaning that we have
491 // requested the smoother. Now we need to build it and to release it.
492 // We don't need to worry about the coarse solver, as we didn't request it.
493 if (!Ac.is_null())
494 smootherFact->Build(level);
495
496 level.Release(*smootherFact);
497 }
498
499 if (isLastLevel == true) {
500 if (isOrigLastLevel == false) {
501 // Earlier in the function, we constructed the next coarse level, and requested data for the that level,
502 // assuming that we are not at the coarsest level. Now, we changed our mind, so we have to release those.
503 Levels_[nextLevelID]->Release(TopRAPFactory(coarseLevelManager, nextLevelManager));
504 }
505 Levels_.resize(nextLevelID);
506 }
507
508 // I think this is the proper place for graph so that it shows every dependence
509 if (isDumpingEnabled_ && ( (dumpLevel_ > 0 && coarseLevelID == dumpLevel_) || dumpLevel_ == -1 ) )
510 DumpCurrentGraph(coarseLevelID);
511
512 if (!isFinestLevel) {
513 // Release the hierarchy data
514 // We release so late to help blocked solvers, as the smoothers for them need A blocks
515 // which we construct in RAPFactory
516 level.Release(coarseRAPFactory);
517 }
518
519 if (oldRank != -1)
520 SetProcRankVerbose(oldRank);
521
522 return isLastLevel;
523 }
524
525 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
526 void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::SetupRe() {
527 int numLevels = Levels_.size();
528 TEUCHOS_TEST_FOR_EXCEPTION(levelManagers_.size() != numLevels, Exceptions::RuntimeError,
529 "Hierarchy::SetupRe: " << Levels_.size() << " levels, but " << levelManagers_.size() << " level factory managers");
530
531 const int startLevel = 0;
532 Clear(startLevel);
533
534#ifdef HAVE_MUELU_DEBUG
535 // Reset factories' data used for debugging
536 for (int i = 0; i < numLevels; i++)
537 levelManagers_[i]->ResetDebugData();
538
539#endif
540
541 int levelID;
542 for (levelID = startLevel; levelID < numLevels;) {
543 bool r = Setup(levelID,
544 (levelID != 0 ? levelManagers_[levelID-1] : Teuchos::null),
545 levelManagers_[levelID],
546 (levelID+1 != numLevels ? levelManagers_[levelID+1] : Teuchos::null));
547 levelID++;
548 if (r) break;
549 }
550 // We may construct fewer levels for some reason, make sure we continue
551 // doing that in the future
552 Levels_ .resize(levelID);
553 levelManagers_.resize(levelID);
554
555 int sizeOfVecs = sizeOfAllocatedLevelMultiVectors_;
556
557 AllocateLevelMultiVectors(sizeOfVecs, true);
558
559 // since the # of levels, etc. may have changed, force re-determination of description during next call to description()
560 ResetDescription();
561
562 describe(GetOStream(Statistics0), GetVerbLevel());
563 }
564
565 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
566 void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Setup(const FactoryManagerBase& manager, int startLevel, int numDesiredLevels) {
567 // Use MueLu::BaseClass::description() to avoid printing "{numLevels = 1}" (numLevels is increasing...)
568 PrintMonitor m0(*this, "Setup (" + this->MueLu::BaseClass::description() + ")", Runtime0);
569
570 Clear(startLevel);
571
572 // Check Levels_[startLevel] exists.
573 TEUCHOS_TEST_FOR_EXCEPTION(Levels_.size() <= startLevel, Exceptions::RuntimeError,
574 "MueLu::Hierarchy::Setup(): fine level (" << startLevel << ") does not exist");
575
576 TEUCHOS_TEST_FOR_EXCEPTION(numDesiredLevels <= 0, Exceptions::RuntimeError,
577 "Constructing non-positive (" << numDesiredLevels << ") number of levels does not make sense.");
578
579 // Check for fine level matrix A
580 TEUCHOS_TEST_FOR_EXCEPTION(!Levels_[startLevel]->IsAvailable("A"), Exceptions::RuntimeError,
581 "MueLu::Hierarchy::Setup(): fine level (" << startLevel << ") has no matrix A! "
582 "Set fine level matrix A using Level.Set()");
583
584 RCP<Operator> A = Levels_[startLevel]->template Get<RCP<Operator> >("A");
585 lib_ = A->getDomainMap()->lib();
586
587 if (IsPrint(Statistics2)) {
588 RCP<Matrix> Amat = rcp_dynamic_cast<Matrix>(A);
589
590 if (!Amat.is_null()) {
591 RCP<ParameterList> params = rcp(new ParameterList());
592 params->set("printLoadBalancingInfo", true);
593 params->set("printCommInfo", true);
594
595 GetOStream(Statistics2) << PerfUtils::PrintMatrixInfo(*Amat, "A0", params);
596 } else {
597 GetOStream(Warnings1) << "Fine level operator is not a matrix, statistics are not available" << std::endl;
598 }
599 }
600
601 RCP<const FactoryManagerBase> rcpmanager = rcpFromRef(manager);
602
603 const int lastLevel = startLevel + numDesiredLevels - 1;
604 GetOStream(Runtime0) << "Setup loop: startLevel = " << startLevel << ", lastLevel = " << lastLevel
605 << " (stop if numLevels = " << numDesiredLevels << " or Ac.size() < " << maxCoarseSize_ << ")" << std::endl;
606
607 // Setup multigrid levels
608 int iLevel = 0;
609 if (numDesiredLevels == 1) {
610 iLevel = 0;
611 Setup(startLevel, Teuchos::null, rcpmanager, Teuchos::null); // setup finest==coarsest level (first and last managers are Teuchos::null)
612
613 } else {
614 bool bIsLastLevel = Setup(startLevel, Teuchos::null, rcpmanager, rcpmanager); // setup finest level (level 0) (first manager is Teuchos::null)
615 if (bIsLastLevel == false) {
616 for (iLevel = startLevel + 1; iLevel < lastLevel; iLevel++) {
617 bIsLastLevel = Setup(iLevel, rcpmanager, rcpmanager, rcpmanager); // setup intermediate levels
618 if (bIsLastLevel == true)
619 break;
620 }
621 if (bIsLastLevel == false)
622 Setup(lastLevel, rcpmanager, rcpmanager, Teuchos::null); // setup coarsest level (last manager is Teuchos::null)
623 }
624 }
625
626 // TODO: some check like this should be done at the beginning of the routine
627 TEUCHOS_TEST_FOR_EXCEPTION(iLevel != Levels_.size() - 1, Exceptions::RuntimeError,
628 "MueLu::Hierarchy::Setup(): number of level");
629
630 // TODO: this is not exception safe: manager will still hold default
631 // factories if you exit this function with an exception
632 manager.Clean();
633
634 describe(GetOStream(Statistics0), GetVerbLevel());
635 }
636
637 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
638 void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Clear(int startLevel) {
639 if (startLevel < GetNumLevels())
640 GetOStream(Runtime0) << "Clearing old data (if any)" << std::endl;
641
642 for (int iLevel = startLevel; iLevel < GetNumLevels(); iLevel++)
643 Levels_[iLevel]->Clear();
644 }
645
646 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
647 void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::ExpertClear() {
648 GetOStream(Runtime0) << "Clearing old data (expert)" << std::endl;
649 for (int iLevel = 0; iLevel < GetNumLevels(); iLevel++)
650 Levels_[iLevel]->ExpertClear();
651 }
652
653#if defined(HAVE_MUELU_EXPERIMENTAL) && defined(HAVE_MUELU_ADDITIVE_VARIANT)
654 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
655 ReturnType Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Iterate(const MultiVector& B, MultiVector& X, ConvData conv,
656 bool InitialGuessIsZero, LO startLevel) {
657 LO nIts = conv.maxIts_;
658 MagnitudeType tol = conv.tol_;
659
660 std::string prefix = this->ShortClassName() + ": ";
661 std::string levelSuffix = " (level=" + toString(startLevel) + ")";
662 std::string levelSuffix1 = " (level=" + toString(startLevel+1) + ")";
663
664 using namespace Teuchos;
665 RCP<Time> CompTime = Teuchos::TimeMonitor::getNewCounter(prefix + "Computational Time (total)");
666 RCP<Time> Concurrent = Teuchos::TimeMonitor::getNewCounter(prefix + "Concurrent portion");
667 RCP<Time> ApplyR = Teuchos::TimeMonitor::getNewCounter(prefix + "R: Computational Time");
668 RCP<Time> ApplyPbar = Teuchos::TimeMonitor::getNewCounter(prefix + "Pbar: Computational Time");
669 RCP<Time> CompFine = Teuchos::TimeMonitor::getNewCounter(prefix + "Fine: Computational Time");
670 RCP<Time> CompCoarse = Teuchos::TimeMonitor::getNewCounter(prefix + "Coarse: Computational Time");
671 RCP<Time> ApplySum = Teuchos::TimeMonitor::getNewCounter(prefix + "Sum: Computational Time");
672 RCP<Time> Synchronize_beginning = Teuchos::TimeMonitor::getNewCounter(prefix + "Synchronize_beginning");
673 RCP<Time> Synchronize_center = Teuchos::TimeMonitor::getNewCounter(prefix + "Synchronize_center");
674 RCP<Time> Synchronize_end = Teuchos::TimeMonitor::getNewCounter(prefix + "Synchronize_end");
675
676 RCP<Level> Fine = Levels_[0];
677 RCP<Level> Coarse;
678
679 RCP<Operator> A = Fine->Get< RCP<Operator> >("A");
680 Teuchos::RCP< const Teuchos::Comm< int > > communicator = A->getDomainMap()->getComm();
681
682
683 //Synchronize_beginning->start();
684 //communicator->barrier();
685 //Synchronize_beginning->stop();
686
687 CompTime->start();
688
689 SC one = STS::one(), zero = STS::zero();
690
691 bool zeroGuess = InitialGuessIsZero;
692
693 // ======= UPFRONT DEFINITION OF COARSE VARIABLES ===========
694
695 //RCP<const Map> origMap;
696 RCP< Operator > P;
697 RCP< Operator > Pbar;
698 RCP< Operator > R;
699 RCP< MultiVector > coarseRhs, coarseX;
700 RCP< Operator > Ac;
701 RCP<SmootherBase> preSmoo_coarse, postSmoo_coarse;
702 bool emptyCoarseSolve = true;
703 RCP<MultiVector> coarseX_prolonged = MultiVectorFactory::Build(X.getMap(), X.getNumVectors(), true);
704
705 RCP<const Import> importer;
706
707 if (Levels_.size() > 1) {
708 Coarse = Levels_[1];
709 if (Coarse->IsAvailable("Importer"))
710 importer = Coarse->Get< RCP<const Import> >("Importer");
711
712 R = Coarse->Get< RCP<Operator> >("R");
713 P = Coarse->Get< RCP<Operator> >("P");
714
715
716 //if(Coarse->IsAvailable("Pbar"))
717 Pbar = Coarse->Get< RCP<Operator> >("Pbar");
718
719 coarseRhs = MultiVectorFactory::Build(R->getRangeMap(), B.getNumVectors(), true);
720
721 Ac = Coarse->Get< RCP< Operator > >("A");
722
723 ApplyR->start();
724 R->apply(B, *coarseRhs, Teuchos::NO_TRANS, one, zero);
725 //P->apply(B, *coarseRhs, Teuchos::TRANS, one, zero);
726 ApplyR->stop();
727
728 if (doPRrebalance_ || importer.is_null()) {
729 coarseX = MultiVectorFactory::Build(coarseRhs->getMap(), X.getNumVectors(), true);
730
731 } else {
732
733 RCP<TimeMonitor> ITime = rcp(new TimeMonitor(*this, prefix + "Solve : import (total)" , Timings0));
734 RCP<TimeMonitor> ILevelTime = rcp(new TimeMonitor(*this, prefix + "Solve : import" + levelSuffix1, Timings0));
735
736 // Import: range map of R --> domain map of rebalanced Ac (before subcomm replacement)
737 RCP<MultiVector> coarseTmp = MultiVectorFactory::Build(importer->getTargetMap(), coarseRhs->getNumVectors());
738 coarseTmp->doImport(*coarseRhs, *importer, Xpetra::INSERT);
739 coarseRhs.swap(coarseTmp);
740
741 coarseX = MultiVectorFactory::Build(importer->getTargetMap(), X.getNumVectors(), true);
742 }
743
744 if (Coarse->IsAvailable("PreSmoother"))
745 preSmoo_coarse = Coarse->Get< RCP<SmootherBase> >("PreSmoother");
746 if (Coarse->IsAvailable("PostSmoother"))
747 postSmoo_coarse = Coarse->Get< RCP<SmootherBase> >("PostSmoother");
748 }
749
750 // ==========================================================
751
752
753 MagnitudeType prevNorm = STS::magnitude(STS::one()), curNorm = STS::magnitude(STS::one());
754 rate_ = 1.0;
755
756 for (LO i = 1; i <= nIts; i++) {
757#ifdef HAVE_MUELU_DEBUG
758 if (A->getDomainMap()->isCompatible(*(X.getMap())) == false) {
759 std::ostringstream ss;
760 ss << "Level " << startLevel << ": level A's domain map is not compatible with X";
761 throw Exceptions::Incompatible(ss.str());
762 }
763
764 if (A->getRangeMap()->isCompatible(*(B.getMap())) == false) {
765 std::ostringstream ss;
766 ss << "Level " << startLevel << ": level A's range map is not compatible with B";
767 throw Exceptions::Incompatible(ss.str());
768 }
769#endif
770 }
771
772 bool emptyFineSolve = true;
773
774 RCP< MultiVector > fineX;
775 fineX = MultiVectorFactory::Build(X.getMap(), X.getNumVectors(), true);
776
777 //Synchronize_center->start();
778 //communicator->barrier();
779 //Synchronize_center->stop();
780
781 Concurrent->start();
782
783 // NOTE: we need to check using IsAvailable before Get here to avoid building default smoother
784 if (Fine->IsAvailable("PreSmoother")) {
785 RCP<SmootherBase> preSmoo = Fine->Get< RCP<SmootherBase> >("PreSmoother");
786 CompFine->start();
787 preSmoo->Apply(*fineX, B, zeroGuess);
788 CompFine->stop();
789 emptyFineSolve = false;
790 }
791 if (Fine->IsAvailable("PostSmoother")) {
792 RCP<SmootherBase> postSmoo = Fine->Get< RCP<SmootherBase> >("PostSmoother");
793 CompFine->start();
794 postSmoo->Apply(*fineX, B, zeroGuess);
795 CompFine->stop();
796
797 emptyFineSolve = false;
798 }
799 if (emptyFineSolve == true) {
800 GetOStream(Warnings1) << "No fine grid smoother" << std::endl;
801 // Fine grid smoother is identity
802 fineX->update(one, B, zero);
803 }
804
805 if (Levels_.size() > 1) {
806 // NOTE: we need to check using IsAvailable before Get here to avoid building default smoother
807 if (Coarse->IsAvailable("PreSmoother")) {
808 CompCoarse->start();
809 preSmoo_coarse->Apply(*coarseX, *coarseRhs, zeroGuess);
810 CompCoarse->stop();
811 emptyCoarseSolve = false;
812
813 }
814 if (Coarse->IsAvailable("PostSmoother")) {
815 CompCoarse->start();
816 postSmoo_coarse->Apply(*coarseX, *coarseRhs, zeroGuess);
817 CompCoarse->stop();
818 emptyCoarseSolve = false;
819
820 }
821 if (emptyCoarseSolve == true) {
822 GetOStream(Warnings1) << "No coarse grid solver" << std::endl;
823 // Coarse operator is identity
824 coarseX->update(one, *coarseRhs, zero);
825 }
826 Concurrent->stop();
827 //Synchronize_end->start();
828 //communicator->barrier();
829 //Synchronize_end->stop();
830
831 if (!doPRrebalance_ && !importer.is_null()) {
832 RCP<TimeMonitor> ITime = rcp(new TimeMonitor(*this, prefix + "Solve : export (total)" , Timings0));
833 RCP<TimeMonitor> ILevelTime = rcp(new TimeMonitor(*this, prefix + "Solve : export" + levelSuffix1, Timings0));
834
835 // Import: range map of rebalanced Ac (before subcomm replacement) --> domain map of P
836 RCP<MultiVector> coarseTmp = MultiVectorFactory::Build(importer->getSourceMap(), coarseX->getNumVectors());
837 coarseTmp->doExport(*coarseX, *importer, Xpetra::INSERT);
838 coarseX.swap(coarseTmp);
839 }
840
841 ApplyPbar->start();
842 Pbar->apply(*coarseX, *coarseX_prolonged, Teuchos::NO_TRANS, one, zero);
843 ApplyPbar->stop();
844 }
845
846 ApplySum->start();
847 X.update(1.0, *fineX, 1.0, *coarseX_prolonged, 0.0);
848 ApplySum->stop();
849
850 CompTime->stop();
851
852 //communicator->barrier();
853
854 return (tol > 0 ? Unconverged : Undefined);
855}
856#else
857 // ---------------------------------------- Iterate -------------------------------------------------------
858 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
859 ReturnType Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Iterate(const MultiVector& B, MultiVector& X, ConvData conv,
860 bool InitialGuessIsZero, LO startLevel) {
861 LO nIts = conv.maxIts_;
862 MagnitudeType tol = conv.tol_;
863
864 // These timers work as follows. "iterateTime" records total time spent in
865 // iterate. "levelTime" records time on a per level basis. The label is
866 // crafted to mimic the per-level messages used in Monitors. Note that a
867 // given level is timed with a TimeMonitor instead of a Monitor or
868 // SubMonitor. This is mainly because I want to time each level
869 // separately, and Monitors/SubMonitors print "(total) xx yy zz" ,
870 // "(sub,total) xx yy zz", respectively, which is subject to
871 // misinterpretation. The per-level TimeMonitors are stopped/started
872 // manually before/after a recursive call to Iterate. A side artifact to
873 // this approach is that the counts for intermediate level timers are twice
874 // the counts for the finest and coarsest levels.
875
876 RCP<Level> Fine = Levels_[startLevel];
877
878 std::string label = FormattingHelper::getColonLabel(Fine->getObjectLabel());
879 std::string prefix = label + this->ShortClassName() + ": ";
880 std::string levelSuffix = " (level=" + toString(startLevel) + ")";
881 std::string levelSuffix1 = " (level=" + toString(startLevel+1) + ")";
882
883 bool useStackedTimer = !Teuchos::TimeMonitor::getStackedTimer().is_null();
884
885 RCP<Monitor> iterateTime;
886 RCP<TimeMonitor> iterateTime1;
887 if (startLevel == 0)
888 iterateTime = rcp(new Monitor(*this, "Solve", label, (nIts == 1) ? None : Runtime0, Timings0));
889 else if (!useStackedTimer)
890 iterateTime1 = rcp(new TimeMonitor(*this, prefix + "Solve (total, level=" + toString(startLevel) + ")", Timings0));
891
892 std::string iterateLevelTimeLabel = prefix + "Solve" + levelSuffix;
893 RCP<TimeMonitor> iterateLevelTime = rcp(new TimeMonitor(*this, iterateLevelTimeLabel, Timings0));
894
895 bool zeroGuess = InitialGuessIsZero;
896
897 RCP<Operator> A = Fine->Get< RCP<Operator> >("A");
898 using namespace Teuchos;
899 RCP<Time> CompCoarse = Teuchos::TimeMonitor::getNewCounter(prefix + "Coarse: Computational Time");
900
901 if (A.is_null()) {
902 // This processor does not have any data for this process on coarser
903 // levels. This can only happen when there are multiple processors and
904 // we use repartitioning.
905 return Undefined;
906 }
907
908 // If we switched the number of vectors, we'd need to reallocate here.
909 // If the number of vectors is unchanged, this is a noop.
910 // NOTE: We need to check against B because the tests in AllocateLevelMultiVectors
911 // will fail on Stokhos Scalar types (due to the so-called 'hidden dimension')
912 const BlockedMultiVector * Bblocked = dynamic_cast<const BlockedMultiVector*>(&B);
913 if(residual_.size() > startLevel &&
914 ( ( Bblocked && !Bblocked->isSameSize(*residual_[startLevel])) ||
915 (!Bblocked && !residual_[startLevel]->isSameSize(B))))
916 DeleteLevelMultiVectors();
917 AllocateLevelMultiVectors(X.getNumVectors());
918
919 // Print residual information before iterating
920 typedef Teuchos::ScalarTraits<typename STS::magnitudeType> STM;
921 MagnitudeType prevNorm = STM::one(), curNorm = STM::one();
922 rate_ = 1.0;
923 if (startLevel == 0 && !isPreconditioner_ &&
924 (IsPrint(Statistics1) || tol > 0)) {
925 // We calculate the residual only if we want to print it out, or if we
926 // want to stop once we achive the tolerance
927 Teuchos::Array<MagnitudeType> rn;
928 rn = Utilities::ResidualNorm(*A, X, B,*residual_[startLevel]);
929
930 if (tol > 0) {
931 bool passed = true;
932 for (LO k = 0; k < rn.size(); k++)
933 if (rn[k] >= tol)
934 passed = false;
935
936 if (passed)
937 return Converged;
938 }
939
940 if (IsPrint(Statistics1))
941 GetOStream(Statistics1) << "iter: "
942 << std::setiosflags(std::ios::left)
943 << std::setprecision(3) << 0 // iter 0
944 << " residual = "
945 << std::setprecision(10) << rn
946 << std::endl;
947 }
948
949 SC one = STS::one(), zero = STS::zero();
950 for (LO i = 1; i <= nIts; i++) {
951#ifdef HAVE_MUELU_DEBUG
952#if 0 // TODO fix me
953 if (A->getDomainMap()->isCompatible(*(X.getMap())) == false) {
954 std::ostringstream ss;
955 ss << "Level " << startLevel << ": level A's domain map is not compatible with X";
956 throw Exceptions::Incompatible(ss.str());
957 }
958
959 if (A->getRangeMap()->isCompatible(*(B.getMap())) == false) {
960 std::ostringstream ss;
961 ss << "Level " << startLevel << ": level A's range map is not compatible with B";
962 throw Exceptions::Incompatible(ss.str());
963 }
964#endif
965#endif
966
967 if (startLevel == as<LO>(Levels_.size())-1) {
968 // On the coarsest level, we do either smoothing (if defined) or a direct solve.
969 RCP<TimeMonitor> CLevelTime = rcp(new TimeMonitor(*this, prefix + "Solve : coarse" + levelSuffix, Timings0));
970
971 bool emptySolve = true;
972
973 // NOTE: we need to check using IsAvailable before Get here to avoid building default smoother
974 if (Fine->IsAvailable("PreSmoother")) {
975 RCP<SmootherBase> preSmoo = Fine->Get< RCP<SmootherBase> >("PreSmoother");
976 CompCoarse->start();
977 preSmoo->Apply(X, B, zeroGuess);
978 CompCoarse->stop();
979 zeroGuess = false;
980 emptySolve = false;
981 }
982 if (Fine->IsAvailable("PostSmoother")) {
983 RCP<SmootherBase> postSmoo = Fine->Get< RCP<SmootherBase> >("PostSmoother");
984 CompCoarse->start();
985 postSmoo->Apply(X, B, zeroGuess);
986 CompCoarse->stop();
987 emptySolve = false;
988 zeroGuess = false;
989 }
990 if (emptySolve == true) {
991 GetOStream(Warnings1) << "No coarse grid solver" << std::endl;
992 // Coarse operator is identity
993 X.update(one, B, zero);
994 }
995
996 } else {
997 // On intermediate levels, we do cycles
998 RCP<Level> Coarse = Levels_[startLevel+1];
999 {
1000 // ============== PRESMOOTHING ==============
1001 RCP<TimeMonitor> STime;
1002 if (!useStackedTimer)
1003 STime = rcp(new TimeMonitor(*this, prefix + "Solve : smoothing (total)" , Timings0));
1004 RCP<TimeMonitor> SLevelTime = rcp(new TimeMonitor(*this, prefix + "Solve : smoothing" + levelSuffix, Timings0));
1005
1006 if (Fine->IsAvailable("PreSmoother")) {
1007 RCP<SmootherBase> preSmoo = Fine->Get< RCP<SmootherBase> >("PreSmoother");
1008 preSmoo->Apply(X, B, zeroGuess);
1009 zeroGuess = false;
1010 }
1011 }
1012
1013 RCP<MultiVector> residual;
1014 {
1015 RCP<TimeMonitor> ATime;
1016 if (!useStackedTimer)
1017 ATime = rcp(new TimeMonitor(*this, prefix + "Solve : residual calculation (total)" , Timings0));
1018 RCP<TimeMonitor> ALevelTime = rcp(new TimeMonitor(*this, prefix + "Solve : residual calculation" + levelSuffix, Timings0));
1019 if (zeroGuess) {
1020 // If there's a pre-smoother, then zeroGuess is false. If there isn't and people still have zeroGuess set,
1021 // then then X still has who knows what, so we need to zero it before we go to the coarse grid.
1022 X.putScalar(zero);
1023 }
1024
1025 Utilities::Residual(*A, X, B,*residual_[startLevel]);
1026 residual = residual_[startLevel];
1027 }
1028
1029 RCP<Operator> P = Coarse->Get< RCP<Operator> >("P");
1030 if (Coarse->IsAvailable("Pbar"))
1031 P = Coarse->Get< RCP<Operator> >("Pbar");
1032
1033 RCP<MultiVector> coarseRhs, coarseX;
1034 // const bool initializeWithZeros = true;
1035 {
1036 // ============== RESTRICTION ==============
1037 RCP<TimeMonitor> RTime;
1038 if (!useStackedTimer)
1039 RTime = rcp(new TimeMonitor(*this, prefix + "Solve : restriction (total)" , Timings0));
1040 RCP<TimeMonitor> RLevelTime = rcp(new TimeMonitor(*this, prefix + "Solve : restriction" + levelSuffix, Timings0));
1041 coarseRhs = coarseRhs_[startLevel];
1042
1043 if (implicitTranspose_) {
1044 P->apply(*residual, *coarseRhs, Teuchos::TRANS, one, zero);
1045
1046 } else {
1047 RCP<Operator> R = Coarse->Get< RCP<Operator> >("R");
1048 R->apply(*residual, *coarseRhs, Teuchos::NO_TRANS, one, zero);
1049 }
1050 }
1051
1052 RCP<const Import> importer;
1053 if (Coarse->IsAvailable("Importer"))
1054 importer = Coarse->Get< RCP<const Import> >("Importer");
1055
1056 coarseX = coarseX_[startLevel];
1057 if (!doPRrebalance_ && !importer.is_null()) {
1058 RCP<TimeMonitor> ITime;
1059 if (!useStackedTimer)
1060 ITime = rcp(new TimeMonitor(*this, prefix + "Solve : import (total)" , Timings0));
1061 RCP<TimeMonitor> ILevelTime = rcp(new TimeMonitor(*this, prefix + "Solve : import" + levelSuffix1, Timings0));
1062
1063 // Import: range map of R --> domain map of rebalanced Ac (before subcomm replacement)
1064 RCP<MultiVector> coarseTmp = coarseImport_[startLevel];
1065 coarseTmp->doImport(*coarseRhs, *importer, Xpetra::INSERT);
1066 coarseRhs.swap(coarseTmp);
1067 }
1068
1069 RCP<Operator> Ac = Coarse->Get< RCP<Operator> >("A");
1070 if (!Ac.is_null()) {
1071 RCP<const Map> origXMap = coarseX->getMap();
1072 RCP<const Map> origRhsMap = coarseRhs->getMap();
1073
1074 // Replace maps with maps with a subcommunicator
1075 coarseRhs->replaceMap(Ac->getRangeMap());
1076 coarseX ->replaceMap(Ac->getDomainMap());
1077
1078 {
1079 iterateLevelTime = Teuchos::null; // stop timing this level
1080
1081 Iterate(*coarseRhs, *coarseX, 1, true, startLevel+1);
1082 // ^^ zero initial guess
1083 if (Cycle_ == WCYCLE && WCycleStartLevel_ >= startLevel)
1084 Iterate(*coarseRhs, *coarseX, 1, false, startLevel+1);
1085 // ^^ nonzero initial guess
1086
1087 iterateLevelTime = rcp(new TimeMonitor(*this, iterateLevelTimeLabel)); // restart timing this level
1088 }
1089 coarseX->replaceMap(origXMap);
1090 coarseRhs->replaceMap(origRhsMap);
1091 }
1092
1093 if (!doPRrebalance_ && !importer.is_null()) {
1094 RCP<TimeMonitor> ITime;
1095 if (!useStackedTimer)
1096 ITime = rcp(new TimeMonitor(*this, prefix + "Solve : export (total)" , Timings0));
1097 RCP<TimeMonitor> ILevelTime = rcp(new TimeMonitor(*this, prefix + "Solve : export" + levelSuffix1, Timings0));
1098
1099 // Import: range map of rebalanced Ac (before subcomm replacement) --> domain map of P
1100 RCP<MultiVector> coarseTmp = coarseExport_[startLevel];
1101 coarseTmp->doExport(*coarseX, *importer, Xpetra::INSERT);
1102 coarseX.swap(coarseTmp);
1103 }
1104
1105 {
1106 // ============== PROLONGATION ==============
1107 RCP<TimeMonitor> PTime;
1108 if (!useStackedTimer)
1109 PTime = rcp(new TimeMonitor(*this, prefix + "Solve : prolongation (total)" , Timings0));
1110 RCP<TimeMonitor> PLevelTime = rcp(new TimeMonitor(*this, prefix + "Solve : prolongation" + levelSuffix, Timings0));
1111 // Update X += P * coarseX
1112 // Note that due to what may be round-off error accumulation, use of the fused kernel
1113 // P->apply(*coarseX, X, Teuchos::NO_TRANS, one, one);
1114 // can in some cases result in slightly higher iteration counts.
1115 if (fuseProlongationAndUpdate_) {
1116 P->apply(*coarseX, X, Teuchos::NO_TRANS, scalingFactor_, one);
1117 } else {
1118 RCP<MultiVector> correction = correction_[startLevel];
1119 P->apply(*coarseX, *correction, Teuchos::NO_TRANS, one, zero);
1120 X.update(scalingFactor_, *correction, one);
1121 }
1122 }
1123
1124 {
1125 // ============== POSTSMOOTHING ==============
1126 RCP<TimeMonitor> STime;
1127 if (!useStackedTimer)
1128 STime = rcp(new TimeMonitor(*this, prefix + "Solve : smoothing (total)" , Timings0));
1129 RCP<TimeMonitor> SLevelTime = rcp(new TimeMonitor(*this, prefix + "Solve : smoothing" + levelSuffix, Timings0));
1130
1131 if (Fine->IsAvailable("PostSmoother")) {
1132 RCP<SmootherBase> postSmoo = Fine->Get< RCP<SmootherBase> >("PostSmoother");
1133 postSmoo->Apply(X, B, false);
1134 }
1135 }
1136 }
1137 zeroGuess = false;
1138
1139 if (startLevel == 0 && !isPreconditioner_ &&
1140 (IsPrint(Statistics1) || tol > 0)) {
1141 // We calculate the residual only if we want to print it out, or if we
1142 // want to stop once we achive the tolerance
1143 Teuchos::Array<MagnitudeType> rn;
1144 rn = Utilities::ResidualNorm(*A, X, B,*residual_[startLevel]);
1145
1146 prevNorm = curNorm;
1147 curNorm = rn[0];
1148 rate_ = as<MagnitudeType>(curNorm / prevNorm);
1149
1150 if (IsPrint(Statistics1))
1151 GetOStream(Statistics1) << "iter: "
1152 << std::setiosflags(std::ios::left)
1153 << std::setprecision(3) << i
1154 << " residual = "
1155 << std::setprecision(10) << rn
1156 << std::endl;
1157
1158 if (tol > 0) {
1159 bool passed = true;
1160 for (LO k = 0; k < rn.size(); k++)
1161 if (rn[k] >= tol)
1162 passed = false;
1163
1164 if (passed)
1165 return Converged;
1166 }
1167 }
1168 }
1169 return (tol > 0 ? Unconverged : Undefined);
1170 }
1171#endif
1172
1173
1174 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1175 void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write(const LO& start, const LO& end, const std::string &suffix) {
1176 LO startLevel = (start != -1 ? start : 0);
1177 LO endLevel = (end != -1 ? end : Levels_.size()-1);
1178
1179 TEUCHOS_TEST_FOR_EXCEPTION(startLevel > endLevel, Exceptions::RuntimeError,
1180 "MueLu::Hierarchy::Write : startLevel must be <= endLevel");
1181
1182 TEUCHOS_TEST_FOR_EXCEPTION(startLevel < 0 || endLevel >= Levels_.size(), Exceptions::RuntimeError,
1183 "MueLu::Hierarchy::Write bad start or end level");
1184
1185 for (LO i = startLevel; i < endLevel + 1; i++) {
1186 RCP<Matrix> A = rcp_dynamic_cast<Matrix>(Levels_[i]-> template Get< RCP< Operator> >("A")), P, R;
1187 if (i > 0) {
1188 P = rcp_dynamic_cast<Matrix>(Levels_[i]-> template Get< RCP< Operator> >("P"));
1189 if (!implicitTranspose_)
1190 R = rcp_dynamic_cast<Matrix>(Levels_[i]-> template Get< RCP< Operator> >("R"));
1191 }
1192
1193 if (!A.is_null()) Xpetra::IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write("A_" + toString(i) + suffix + ".m", *A);
1194 if (!P.is_null()) {
1195 Xpetra::IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write("P_" + toString(i) + suffix + ".m", *P);
1196 }
1197 if (!R.is_null()) {
1198 Xpetra::IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write("R_" + toString(i) + suffix + ".m", *R);
1199 }
1200 }
1201 }
1202
1203 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1204 void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Keep(const std::string & ename, const FactoryBase* factory) {
1205 for (Array<RCP<Level> >::iterator it = Levels_.begin(); it != Levels_.end(); ++it)
1206 (*it)->Keep(ename, factory);
1207 }
1208
1209 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1210 void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Delete(const std::string& ename, const FactoryBase* factory) {
1211 for (Array<RCP<Level> >::iterator it = Levels_.begin(); it != Levels_.end(); ++it)
1212 (*it)->Delete(ename, factory);
1213 }
1214
1215 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1216 void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::AddKeepFlag(const std::string & ename, const FactoryBase* factory, KeepType keep) {
1217 for (Array<RCP<Level> >::iterator it = Levels_.begin(); it != Levels_.end(); ++it)
1218 (*it)->AddKeepFlag(ename, factory, keep);
1219 }
1220
1221 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1222 void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::RemoveKeepFlag(const std::string & ename, const FactoryBase* factory, KeepType keep) {
1223 for (Array<RCP<Level> >::iterator it = Levels_.begin(); it != Levels_.end(); ++it)
1224 (*it)->RemoveKeepFlag(ename, factory, keep);
1225 }
1226
1227 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1228 std::string Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::description() const {
1229 if ( description_.empty() )
1230 {
1231 std::ostringstream out;
1232 out << BaseClass::description();
1233 out << "{#levels = " << GetGlobalNumLevels() << ", complexity = " << GetOperatorComplexity() << "}";
1234 description_ = out.str();
1235 }
1236 return description_;
1237 }
1238
1239 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1240 void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::describe(Teuchos::FancyOStream& out, const Teuchos::EVerbosityLevel tVerbLevel) const {
1241 describe(out, toMueLuVerbLevel(tVerbLevel));
1242 }
1243
1244 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1245 void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::describe(Teuchos::FancyOStream& out, const VerbLevel verbLevel) const {
1246 RCP<Operator> A0 = Levels_[0]->template Get<RCP<Operator> >("A");
1247 RCP<const Teuchos::Comm<int> > comm = A0->getDomainMap()->getComm();
1248
1249 int numLevels = GetNumLevels();
1250 RCP<Operator> Ac = Levels_[numLevels-1]->template Get<RCP<Operator> >("A");
1251 if (Ac.is_null()) {
1252 // It may happen that we do repartition on the last level, but the matrix
1253 // is small enough to satisfy "max coarse size" requirement. Then, even
1254 // though we have the level, the matrix would be null on all but one processors
1255 numLevels--;
1256 }
1257 int root = comm->getRank();
1258
1259#ifdef HAVE_MPI
1260 int smartData = numLevels*comm->getSize() + comm->getRank(), maxSmartData;
1261 reduceAll(*comm, Teuchos::REDUCE_MAX, smartData, Teuchos::ptr(&maxSmartData));
1262 root = maxSmartData % comm->getSize();
1263#endif
1264
1265 // Compute smoother complexity, if needed
1266 double smoother_comp =-1.0;
1267 if (verbLevel & (Statistics0 | Test))
1268 smoother_comp = GetSmootherComplexity();
1269
1270 std::string outstr;
1271 if (comm->getRank() == root && verbLevel & (Statistics0 | Test)) {
1272 std::vector<Xpetra::global_size_t> nnzPerLevel;
1273 std::vector<Xpetra::global_size_t> rowsPerLevel;
1274 std::vector<int> numProcsPerLevel;
1275 bool aborted = false;
1276 for (int i = 0; i < numLevels; i++) {
1277 TEUCHOS_TEST_FOR_EXCEPTION(!(Levels_[i]->IsAvailable("A")) , Exceptions::RuntimeError,
1278 "Operator A is not available on level " << i);
1279
1280 RCP<Operator> A = Levels_[i]->template Get<RCP<Operator> >("A");
1281 TEUCHOS_TEST_FOR_EXCEPTION(A.is_null(), Exceptions::RuntimeError,
1282 "Operator A on level " << i << " is null.");
1283
1284 RCP<Matrix> Am = rcp_dynamic_cast<Matrix>(A);
1285 if (Am.is_null()) {
1286 GetOStream(Warnings0) << "Some level operators are not matrices, statistics calculation aborted" << std::endl;
1287 aborted = true;
1288 break;
1289 }
1290
1291 Xpetra::global_size_t nnz = Am->getGlobalNumEntries();
1292 nnzPerLevel .push_back(nnz);
1293 rowsPerLevel .push_back(Am->getGlobalNumRows());
1294 numProcsPerLevel.push_back(Am->getRowMap()->getComm()->getSize());
1295 }
1296
1297 if (!aborted) {
1298 std::string label = Levels_[0]->getObjectLabel();
1299 std::ostringstream oss;
1300 oss << std::setfill(' ');
1301 oss << "\n--------------------------------------------------------------------------------\n";
1302 oss << "--- Multigrid Summary " << std::setw(28) << std::left << label << "---\n";
1303 oss << "--------------------------------------------------------------------------------" << std::endl;
1304 if (verbLevel & Parameters1)
1305 oss << "Scalar = " << Teuchos::ScalarTraits<Scalar>::name() << std::endl;
1306 oss << "Number of levels = " << numLevels << std::endl;
1307 oss << "Operator complexity = " << std::setprecision(2) << std::setiosflags(std::ios::fixed)
1308 << GetOperatorComplexity() << std::endl;
1309
1310 if(smoother_comp!=-1.0) {
1311 oss << "Smoother complexity = " << std::setprecision(2) << std::setiosflags(std::ios::fixed)
1312 << smoother_comp << std::endl;
1313 }
1314
1315 switch (Cycle_) {
1316 case VCYCLE:
1317 oss << "Cycle type = V" << std::endl;
1318 break;
1319 case WCYCLE:
1320 oss << "Cycle type = W" << std::endl;
1321 if (WCycleStartLevel_ > 0)
1322 oss << "Cycle start level = " << WCycleStartLevel_ << std::endl;
1323 break;
1324 default:
1325 break;
1326 };
1327 oss << std::endl;
1328
1329 Xpetra::global_size_t tt = rowsPerLevel[0];
1330 int rowspacer = 2; while (tt != 0) { tt /= 10; rowspacer++; }
1331 tt = nnzPerLevel[0];
1332 int nnzspacer = 2; while (tt != 0) { tt /= 10; nnzspacer++; }
1333 tt = numProcsPerLevel[0];
1334 int npspacer = 2; while (tt != 0) { tt /= 10; npspacer++; }
1335 oss << "level " << std::setw(rowspacer) << " rows " << std::setw(nnzspacer) << " nnz " << " nnz/row" << std::setw(npspacer) << " c ratio" << " procs" << std::endl;
1336 for (size_t i = 0; i < nnzPerLevel.size(); ++i) {
1337 oss << " " << i << " ";
1338 oss << std::setw(rowspacer) << rowsPerLevel[i];
1339 oss << std::setw(nnzspacer) << nnzPerLevel[i];
1340 oss << std::setprecision(2) << std::setiosflags(std::ios::fixed);
1341 oss << std::setw(9) << as<double>(nnzPerLevel[i]) / rowsPerLevel[i];
1342 if (i) oss << std::setw(9) << as<double>(rowsPerLevel[i-1])/rowsPerLevel[i];
1343 else oss << std::setw(9) << " ";
1344 oss << " " << std::setw(npspacer) << numProcsPerLevel[i] << std::endl;
1345 }
1346 oss << std::endl;
1347 for (int i = 0; i < GetNumLevels(); ++i) {
1348 RCP<SmootherBase> preSmoo, postSmoo;
1349 if (Levels_[i]->IsAvailable("PreSmoother"))
1350 preSmoo = Levels_[i]->template Get< RCP<SmootherBase> >("PreSmoother");
1351 if (Levels_[i]->IsAvailable("PostSmoother"))
1352 postSmoo = Levels_[i]->template Get< RCP<SmootherBase> >("PostSmoother");
1353
1354 if (preSmoo != null && preSmoo == postSmoo)
1355 oss << "Smoother (level " << i << ") both : " << preSmoo->description() << std::endl;
1356 else {
1357 oss << "Smoother (level " << i << ") pre : "
1358 << (preSmoo != null ? preSmoo->description() : "no smoother") << std::endl;
1359 oss << "Smoother (level " << i << ") post : "
1360 << (postSmoo != null ? postSmoo->description() : "no smoother") << std::endl;
1361 }
1362
1363 oss << std::endl;
1364 }
1365
1366 outstr = oss.str();
1367 }
1368 }
1369
1370#ifdef HAVE_MPI
1371 RCP<const Teuchos::MpiComm<int> > mpiComm = rcp_dynamic_cast<const Teuchos::MpiComm<int> >(comm);
1372 MPI_Comm rawComm = (*mpiComm->getRawMpiComm())();
1373
1374 int strLength = outstr.size();
1375 MPI_Bcast(&strLength, 1, MPI_INT, root, rawComm);
1376 if (comm->getRank() != root)
1377 outstr.resize(strLength);
1378 MPI_Bcast(&outstr[0], strLength, MPI_CHAR, root, rawComm);
1379#endif
1380
1381 out << outstr;
1382 }
1383
1384 // NOTE: at some point this should be replaced by a friend operator <<
1385 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1386 void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::print(std::ostream& out, const VerbLevel verbLevel) const {
1387 Teuchos::OSTab tab2(out);
1388 for (int i = 0; i < GetNumLevels(); ++i)
1389 Levels_[i]->print(out, verbLevel);
1390 }
1391
1392 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1393 void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::IsPreconditioner(const bool flag) {
1394 isPreconditioner_ = flag;
1395 }
1396
1397 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1398 void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::DumpCurrentGraph(int currLevel) const {
1399 if (GetProcRankVerbose() != 0)
1400 return;
1401#if defined(HAVE_MUELU_BOOST) && defined(HAVE_MUELU_BOOST_FOR_REAL) && defined(BOOST_VERSION) && (BOOST_VERSION >= 104400)
1402
1403 BoostGraph graph;
1404
1405 BoostProperties dp;
1406 dp.property("label", boost::get(boost::vertex_name, graph));
1407 dp.property("id", boost::get(boost::vertex_index, graph));
1408 dp.property("label", boost::get(boost::edge_name, graph));
1409 dp.property("color", boost::get(boost::edge_color, graph));
1410
1411 // create local maps
1412 std::map<const FactoryBase*, BoostVertex> vindices;
1413 typedef std::map<std::pair<BoostVertex,BoostVertex>, std::string> emap; emap edges;
1414
1415 static int call_id=0;
1416
1417 RCP<Operator> A = Levels_[0]->template Get<RCP<Operator> >("A");
1418 int rank = A->getDomainMap()->getComm()->getRank();
1419
1420 // printf("[%d] CMS: ----------------------\n",rank);
1421 for (int i = currLevel; i <= currLevel+1 && i < GetNumLevels(); i++) {
1422 edges.clear();
1423 Levels_[i]->UpdateGraph(vindices, edges, dp, graph);
1424
1425 for (emap::const_iterator eit = edges.begin(); eit != edges.end(); eit++) {
1426 std::pair<BoostEdge, bool> boost_edge = boost::add_edge(eit->first.first, eit->first.second, graph);
1427 // printf("[%d] CMS: Hierarchy, adding edge (%d->%d) %d\n",rank,(int)eit->first.first,(int)eit->first.second,(int)boost_edge.second);
1428 // Because xdot.py views 'Graph' as a keyword
1429 if(eit->second==std::string("Graph")) boost::put("label", dp, boost_edge.first, std::string("Graph_"));
1430 else boost::put("label", dp, boost_edge.first, eit->second);
1431 if (i == currLevel)
1432 boost::put("color", dp, boost_edge.first, std::string("red"));
1433 else
1434 boost::put("color", dp, boost_edge.first, std::string("blue"));
1435 }
1436 }
1437
1438 std::ofstream out(dumpFile_.c_str()+std::string("_")+std::to_string(currLevel)+std::string("_")+std::to_string(call_id)+std::string("_")+ std::to_string(rank) + std::string(".dot"));
1439 boost::write_graphviz_dp(out, graph, dp, std::string("id"));
1440 out.close();
1441 call_id++;
1442#else
1443 GetOStream(Errors) << "Dependency graph output requires boost and MueLu_ENABLE_Boost_for_real" << std::endl;
1444#endif
1445 }
1446
1447 // Enforce that coordinate vector's map is consistent with that of A
1448 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1449 void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::ReplaceCoordinateMap(Level& level) {
1450 RCP<Operator> Ao = level.Get<RCP<Operator> >("A");
1451 RCP<Matrix> A = rcp_dynamic_cast<Matrix>(Ao);
1452 if (A.is_null()) {
1453 GetOStream(Warnings1) << "Hierarchy::ReplaceCoordinateMap: operator is not a matrix, skipping..." << std::endl;
1454 return;
1455 }
1456 if(Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A) != Teuchos::null) {
1457 GetOStream(Warnings1) << "Hierarchy::ReplaceCoordinateMap: operator is a BlockedCrsMatrix, skipping..." << std::endl;
1458 return;
1459 }
1460
1461 typedef Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType,LO,GO,NO> xdMV;
1462
1463 RCP<xdMV> coords = level.Get<RCP<xdMV> >("Coordinates");
1464
1465 if (A->getRowMap()->isSameAs(*(coords->getMap()))) {
1466 GetOStream(Warnings1) << "Hierarchy::ReplaceCoordinateMap: matrix and coordinates maps are same, skipping..." << std::endl;
1467 return;
1468 }
1469
1470 if (A->IsView("stridedMaps") && rcp_dynamic_cast<const StridedMap>(A->getRowMap("stridedMaps")) != Teuchos::null) {
1471 RCP<const StridedMap> stridedRowMap = rcp_dynamic_cast<const StridedMap>(A->getRowMap("stridedMaps"));
1472
1473 // It is better to through an exceptions if maps may be inconsistent, than to ignore it and experience unfathomable breakdowns
1474 TEUCHOS_TEST_FOR_EXCEPTION(stridedRowMap->getStridedBlockId() != -1 || stridedRowMap->getOffset() != 0,
1475 Exceptions::RuntimeError, "Hierarchy::ReplaceCoordinateMap: nontrivial maps (block id = " << stridedRowMap->getStridedBlockId()
1476 << ", offset = " << stridedRowMap->getOffset() << ")");
1477 }
1478
1479 GetOStream(Runtime1) << "Replacing coordinate map" << std::endl;
1480
1481 size_t blkSize = A->GetFixedBlockSize();
1482
1483 RCP<const Map> nodeMap = A->getRowMap();
1484 if (blkSize > 1) {
1485 // Create a nodal map, as coordinates have not been expanded to a DOF map yet.
1486 RCP<const Map> dofMap = A->getRowMap();
1487 GO indexBase = dofMap->getIndexBase();
1488 size_t numLocalDOFs = dofMap->getNodeNumElements();
1489 TEUCHOS_TEST_FOR_EXCEPTION(numLocalDOFs % blkSize, Exceptions::RuntimeError,
1490 "Hierarchy::ReplaceCoordinateMap: block size (" << blkSize << ") is incompatible with the number of local dofs in a row map (" << numLocalDOFs);
1491 ArrayView<const GO> GIDs = dofMap->getNodeElementList();
1492
1493 Array<GO> nodeGIDs(numLocalDOFs/blkSize);
1494 for (size_t i = 0; i < numLocalDOFs; i += blkSize)
1495 nodeGIDs[i/blkSize] = (GIDs[i] - indexBase)/blkSize + indexBase;
1496
1497 Xpetra::global_size_t INVALID = Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid();
1498 nodeMap = MapFactory::Build(dofMap->lib(), INVALID, nodeGIDs(), indexBase, dofMap->getComm());
1499 } else {
1500 // blkSize == 1
1501 // Check whether the length of vectors fits to the size of A
1502 // If yes, make sure that the maps are matching
1503 // If no, throw a warning but do not touch the Coordinates
1504 if(coords->getLocalLength() != A->getRowMap()->getNodeNumElements()) {
1505 GetOStream(Warnings) << "Coordinate vector does not match row map of matrix A!" << std::endl;
1506 return;
1507 }
1508 }
1509
1510 Array<ArrayView<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> > coordDataView;
1511 std::vector<ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> > coordData;
1512 for (size_t i = 0; i < coords->getNumVectors(); i++) {
1513 coordData.push_back(coords->getData(i));
1514 coordDataView.push_back(coordData[i]());
1515 }
1516
1517 RCP<xdMV> newCoords = Xpetra::MultiVectorFactory<typename Teuchos::ScalarTraits<Scalar>::coordinateType,LO,GO,NO>::Build(nodeMap, coordDataView(), coords->getNumVectors());
1518 level.Set("Coordinates", newCoords);
1519 }
1520
1521 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1522 void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::AllocateLevelMultiVectors(int numvecs, bool forceMapCheck) {
1523 int N = Levels_.size();
1524 if( ( (sizeOfAllocatedLevelMultiVectors_ == numvecs && residual_.size() == N) || numvecs<=0 ) && !forceMapCheck) return;
1525
1526 // If, somehow, we changed the number of levels, delete everything first
1527 if(residual_.size() != N) {
1528 DeleteLevelMultiVectors();
1529
1530 residual_.resize(N);
1531 coarseRhs_.resize(N);
1532 coarseX_.resize(N);
1533 coarseImport_.resize(N);
1534 coarseExport_.resize(N);
1535 correction_.resize(N);
1536 }
1537
1538 for(int i=0; i<N; i++) {
1539 RCP<Operator> A = Levels_[i]->template Get< RCP<Operator> >("A");
1540 if(!A.is_null()) {
1541 // This dance is because we allow A to have a BlockedMap and X/B to have (compatible) non-blocked map
1542 RCP<const BlockedCrsMatrix> A_as_blocked = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(A);
1543 RCP<const Map> Arm = A->getRangeMap();
1544 RCP<const Map> Adm = A->getDomainMap();
1545 if(!A_as_blocked.is_null()) {
1546 Adm = A_as_blocked->getFullDomainMap();
1547 }
1548
1549 if (residual_[i].is_null() || !residual_[i]->getMap()->isSameAs(*Arm))
1550 // This is zero'd by default since it is filled via an operator apply
1551 residual_[i] = MultiVectorFactory::Build(Arm, numvecs, true);
1552 if (correction_[i].is_null() || !correction_[i]->getMap()->isSameAs(*Adm))
1553 correction_[i] = MultiVectorFactory::Build(Adm, numvecs, false);
1554 }
1555
1556 if(i+1<N) {
1557 // This is zero'd by default since it is filled via an operator apply
1558 if(implicitTranspose_) {
1559 RCP<Operator> P = Levels_[i+1]->template Get< RCP<Operator> >("P");
1560 if(!P.is_null()) {
1561 RCP<const Map> map = P->getDomainMap();
1562 if (coarseRhs_[i].is_null() || !coarseRhs_[i]->getMap()->isSameAs(*map))
1563 coarseRhs_[i] = MultiVectorFactory::Build(map, numvecs, true);
1564 }
1565 } else {
1566 RCP<Operator> R = Levels_[i+1]->template Get< RCP<Operator> >("R");
1567 if (!R.is_null()) {
1568 RCP<const Map> map = R->getRangeMap();
1569 if (coarseRhs_[i].is_null() || !coarseRhs_[i]->getMap()->isSameAs(*map))
1570 coarseRhs_[i] = MultiVectorFactory::Build(map, numvecs, true);
1571 }
1572 }
1573
1574
1575 RCP<const Import> importer;
1576 if(Levels_[i+1]->IsAvailable("Importer"))
1577 importer = Levels_[i+1]->template Get< RCP<const Import> >("Importer");
1578 if (doPRrebalance_ || importer.is_null()) {
1579 RCP<const Map> map = coarseRhs_[i]->getMap();
1580 if (coarseX_[i].is_null() || !coarseX_[i]->getMap()->isSameAs(*map))
1581 coarseX_[i] = MultiVectorFactory::Build(map, numvecs, true);
1582 } else {
1583 RCP<const Map> map;
1584 map = importer->getTargetMap();
1585 if (coarseImport_[i].is_null() || !coarseImport_[i]->getMap()->isSameAs(*map)) {
1586 coarseImport_[i] = MultiVectorFactory::Build(map, numvecs,false);
1587 coarseX_[i] = MultiVectorFactory::Build(map,numvecs,false);
1588 }
1589 map = importer->getSourceMap();
1590 if (coarseExport_[i].is_null() || !coarseExport_[i]->getMap()->isSameAs(*map))
1591 coarseExport_[i] = MultiVectorFactory::Build(map, numvecs,false);
1592 }
1593 }
1594 }
1595 sizeOfAllocatedLevelMultiVectors_ = numvecs;
1596 }
1597
1598
1599template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1600void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::DeleteLevelMultiVectors() {
1601 if(sizeOfAllocatedLevelMultiVectors_==0) return;
1602 residual_.resize(0);
1603 coarseRhs_.resize(0);
1604 coarseX_.resize(0);
1605 coarseImport_.resize(0);
1606 coarseExport_.resize(0);
1607 correction_.resize(0);
1608 sizeOfAllocatedLevelMultiVectors_ = 0;
1609}
1610
1611
1612} //namespace MueLu
1613
1614#endif // MUELU_HIERARCHY_DEF_HPP
MueLu_BoostGraphviz.hpp
MueLu_Hierarchy_decl.hpp
MueLu_Level.hpp
MueLu_Monitor.hpp
MueLu_SmootherBase.hpp
MueLu_SmootherFactoryBase.hpp
MueLu::Describable::description
virtual std::string description() const
Return a simple one-line description of this object.
Definition: MueLu_Describable.cpp:61
MueLu::Exceptions::Incompatible
Exception throws to report incompatible objects (like maps).
Definition: MueLu_Exceptions.hpp:82
MueLu::Exceptions::RuntimeError
Exception throws to report errors in the internal logical of the program.
Definition: MueLu_Exceptions.hpp:70
MueLu::FactoryBase
Base class for factories (e.g., R, P, and A_coarse).
Definition: MueLu_FactoryBase.hpp:60
MueLu::FactoryManagerBase
Class that provides default factories within Needs class.
Definition: MueLu_FactoryManagerBase.hpp:63
MueLu::Hierarchy
Provides methods to build a multigrid hierarchy and apply multigrid cycles.
Definition: MueLu_Hierarchy_decl.hpp:103
MueLu::Hierarchy::AddLevel
void AddLevel(const RCP< Level > &level)
Add a level at the end of the hierarchy.
Definition: MueLu_Hierarchy_def.hpp:125
MueLu::Hierarchy::GetSmootherComplexity
double GetSmootherComplexity() const
Definition: MueLu_Hierarchy_def.hpp:196
MueLu::Hierarchy::Write
void Write(const LO &start=-1, const LO &end=-1, const std::string &suffix="")
Print matrices in the multigrid hierarchy to file.
Definition: MueLu_Hierarchy_def.hpp:1175
MueLu::Hierarchy::GetLevel
RCP< Level > & GetLevel(const int levelID=0)
Retrieve a certain level from hierarchy.
Definition: MueLu_Hierarchy_def.hpp:149
MueLu::Hierarchy::CheckLevel
void CheckLevel(Level &level, int levelID)
Helper function.
Definition: MueLu_Hierarchy_def.hpp:236
MueLu::Hierarchy::description
std::string description() const
Return a simple one-line description of this object.
Definition: MueLu_Hierarchy_def.hpp:1228
MueLu::Hierarchy::IsPreconditioner
void IsPreconditioner(const bool flag)
Definition: MueLu_Hierarchy_def.hpp:1393
MueLu::Hierarchy::Levels_
Array< RCP< Level > > Levels_
Container for Level objects.
Definition: MueLu_Hierarchy_decl.hpp:359
MueLu::Hierarchy::Setup
bool Setup(int coarseLevelID, const RCP< const FactoryManagerBase > fineLevelManager, const RCP< const FactoryManagerBase > coarseLevelManager, const RCP< const FactoryManagerBase > nextLevelManager=Teuchos::null)
Multi-level setup phase: build a new level of the hierarchy.
Definition: MueLu_Hierarchy_def.hpp:290
MueLu::Hierarchy::MagnitudeType
STS::magnitudeType MagnitudeType
Definition: MueLu_Hierarchy_decl.hpp:108
MueLu::Hierarchy::describe
void describe(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the Hierarchy with some verbosity level to a FancyOStream object.
Definition: MueLu_Hierarchy_def.hpp:1245
MueLu::Hierarchy::DumpCurrentGraph
void DumpCurrentGraph(int level) const
Definition: MueLu_Hierarchy_def.hpp:1398
MueLu::Hierarchy::SetMatvecParams
void SetMatvecParams(RCP< ParameterList > matvecParams)
Definition: MueLu_Hierarchy_def.hpp:246
MueLu::Hierarchy::lib_
Xpetra::UnderlyingLib lib_
Epetra/Tpetra mode.
Definition: MueLu_Hierarchy_decl.hpp:396
MueLu::Hierarchy::Clear
void Clear(int startLevel=0)
Clear impermanent data from previous setup.
Definition: MueLu_Hierarchy_def.hpp:638
MueLu::Hierarchy::Iterate
ReturnType Iterate(const MultiVector &B, MultiVector &X, ConvData conv=ConvData(), bool InitialGuessIsZero=false, LO startLevel=0)
Apply the multigrid preconditioner.
Definition: MueLu_Hierarchy_def.hpp:859
MueLu::Hierarchy::GetOperatorComplexity
double GetOperatorComplexity() const
Definition: MueLu_Hierarchy_def.hpp:173
MueLu::Hierarchy::AllocateLevelMultiVectors
void AllocateLevelMultiVectors(int numvecs, bool forceMapCheck=false)
Definition: MueLu_Hierarchy_def.hpp:1522
MueLu::Hierarchy::print
void print(std::ostream &out=std::cout, const VerbLevel verbLevel=(MueLu::Parameters|MueLu::Statistics0)) const
Hierarchy::print is local hierarchy function, thus the statistics can be different from global ones.
Definition: MueLu_Hierarchy_def.hpp:1386
MueLu::Hierarchy::Delete
void Delete(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Call Level::Delete(ename, factory) for each level of the Hierarchy.
Definition: MueLu_Hierarchy_def.hpp:1210
MueLu::Hierarchy::GetGlobalNumLevels
int GetGlobalNumLevels() const
Definition: MueLu_Hierarchy_def.hpp:161
MueLu::Hierarchy::Keep
void Keep(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Call Level::Keep(ename, factory) for each level of the Hierarchy.
Definition: MueLu_Hierarchy_def.hpp:1204
MueLu::Hierarchy::AddKeepFlag
void AddKeepFlag(const std::string &ename, const FactoryBase *factory=NoFactory::get(), KeepType keep=MueLu::Keep)
Call Level::AddKeepFlag for each level of the Hierarchy.
Definition: MueLu_Hierarchy_def.hpp:1216
MueLu::Hierarchy::AddNewLevel
void AddNewLevel()
Add a new level at the end of the hierarchy.
Definition: MueLu_Hierarchy_def.hpp:142
MueLu::Hierarchy::RemoveKeepFlag
void RemoveKeepFlag(const std::string &ename, const FactoryBase *factory, KeepType keep=MueLu::All)
Call Level::RemoveKeepFlag for each level of the Hierarchy.
Definition: MueLu_Hierarchy_def.hpp:1222
MueLu::Hierarchy::ReplaceCoordinateMap
void ReplaceCoordinateMap(Level &level)
Definition: MueLu_Hierarchy_def.hpp:1449
MueLu::Level
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
MueLu::Level::IsAvailable
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need's value has been saved.
Definition: MueLu_Level.hpp:334
MueLu::Level::SetComm
void SetComm(RCP< const Teuchos::Comm< int > > const &comm)
Definition: MueLu_Level.hpp:382
MueLu::Level::setlib
void setlib(Xpetra::UnderlyingLib lib2)
Definition: MueLu_Level.hpp:379
MueLu::Level::GetPreviousLevel
RCP< Level > & GetPreviousLevel()
Previous level.
Definition: MueLu_Level.hpp:132
MueLu::Level::Release
void Release(const FactoryBase &factory)
Decrement the storage counter for all the inputs of a factory.
Definition: MueLu_Level.cpp:143
MueLu::Level::GetComm
RCP< const Teuchos::Comm< int > > GetComm() const
Definition: MueLu_Level.hpp:383
MueLu::Level::GetLevelID
int GetLevelID() const
Return level number.
Definition: MueLu_Level.cpp:76
MueLu::Level::Get
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access)....
Definition: MueLu_Level.hpp:192
MueLu::Level::Set
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
Definition: MueLu_Level.hpp:157
MueLu::Level::Request
void Request(const FactoryBase &factory)
Increment the storage counter for all the inputs of a factory.
Definition: MueLu_Level.cpp:136
MueLu::Level::lib
Xpetra::UnderlyingLib lib()
Definition: MueLu_Level.hpp:380
MueLu::Monitor
Timer to be used in non-factories.
Definition: MueLu_Monitor.hpp:107
MueLu::PerfUtils::PrintMatrixInfo
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
Definition: MueLu_PerfUtils_def.hpp:127
MueLu::SetFactoryManager
An exception safe way to call the method 'Level::SetFactoryManager()'.
Definition: MueLu_HierarchyUtils_decl.hpp:65
MueLu::TimeMonitor
Integrates Teuchos::TimeMonitor with MueLu verbosity system.
Definition: MueLu_TimeMonitor.hpp:68
MueLu::TopRAPFactory::Build
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
Definition: MueLu_TopRAPFactory_def.hpp:55
MueLu::Utilities::ResidualNorm
static Teuchos::Array< Magnitude > ResidualNorm(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &X, const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &RHS)
Definition: MueLu_Utilities_decl.hpp:218
MueLu::Utilities::Residual
static RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Residual(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &X, const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &RHS)
Definition: MueLu_Utilities_decl.hpp:220
MueLu::Utilities::SetRandomSeed
static void SetRandomSeed(const Teuchos::Comm< int > &comm)
Definition: MueLu_Utilities_decl.hpp:233
MueLu
Namespace for MueLu classes and methods.
Definition: MueLu_BrickAggregationFactory_decl.hpp:76
MueLu::Warnings0
@ Warnings0
Important warning messages (one line)
Definition: MueLu_VerbosityLevel.hpp:57
MueLu::Statistics2
@ Statistics2
Print even more statistics.
Definition: MueLu_VerbosityLevel.hpp:72
MueLu::Warnings
@ Warnings
Print all warning messages.
Definition: MueLu_VerbosityLevel.hpp:89
MueLu::Errors
@ Errors
Errors.
Definition: MueLu_VerbosityLevel.hpp:55
MueLu::Statistics1
@ Statistics1
Print more statistics.
Definition: MueLu_VerbosityLevel.hpp:71
MueLu::Timings0
@ Timings0
High level timing information (use Teuchos::TimeMonitor::summarize() to print)
Definition: MueLu_VerbosityLevel.hpp:74
MueLu::Runtime0
@ Runtime0
One-liner description of what is happening.
Definition: MueLu_VerbosityLevel.hpp:62
MueLu::Runtime1
@ Runtime1
Description of what is happening (more verbose)
Definition: MueLu_VerbosityLevel.hpp:63
MueLu::Warnings1
@ Warnings1
Additional warnings.
Definition: MueLu_VerbosityLevel.hpp:59
MueLu::Statistics0
@ Statistics0
Print statistics that do not involve significant additional computation.
Definition: MueLu_VerbosityLevel.hpp:70
MueLu::Parameters1
@ Parameters1
Print class parameters (more parameters, more verbose)
Definition: MueLu_VerbosityLevel.hpp:68
MueLu::KeepType
short KeepType
Definition: MueLu_KeepType.hpp:63
MueLu::toString
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
Definition: MueLu_Utilities_decl.hpp:1167
MueLu::toMueLuVerbLevel
VerbLevel toMueLuVerbLevel(const Teuchos::EVerbosityLevel verbLevel)
Translate Teuchos verbosity level to MueLu verbosity level.
Definition: MueLu_VerbosityLevel.cpp:54
MueLu::FormattingHelper::getColonLabel
static std::string getColonLabel(const std::string &label)
Helper function for object label.
Definition: MueLu_Monitor.hpp:66
MueLu::Hierarchy::ConvData
Data struct for defining stopping criteria of multigrid iteration.
Definition: MueLu_Hierarchy_decl.hpp:111