MueLu Version of the Day
Loading...
Searching...
No Matches
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
61
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"
72#include "MueLu_SmootherFactory.hpp"
74
75#include "Teuchos_TimeMonitor.hpp"
76
77
78
79namespace MueLu {
80
81 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
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>
94 : Hierarchy()
95 {
96 setObjectLabel(label);
97 Levels_[0]->setObjectLabel(label);
98 }
99
100 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
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>
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>
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>
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>
157 return Levels_.size();
158 }
159
160 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
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>
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>
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>
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>
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>
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
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
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;
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 if (coarseFact.is_null())
407 coarseFact = rcp(new TopSmootherFactory(coarseLevelManager, "CoarseSolver"));
408 level.Request(*coarseFact);
409 }
410
411 GetOStream(Runtime0) << std::endl;
412 PrintMonitor m0(*this, "Level " + Teuchos::toString(coarseLevelID), static_cast<MsgType>(Runtime0 | Test));
413
414 // Build coarse level hierarchy
415 RCP<Operator> Ac = Teuchos::null;
416 TopRAPFactory coarseRAPFactory(fineLevelManager, coarseLevelManager);
417
418 if (level.IsAvailable("A")) {
419 Ac = level.Get<RCP<Operator> >("A");
420 } else if (!isFinestLevel) {
421 // We only build here, the release is done later
422 coarseRAPFactory.Build(*level.GetPreviousLevel(), level);
423 }
424
425 bool setLastLevelviaMaxCoarseSize = false;
426 if (level.IsAvailable("A"))
427 Ac = level.Get<RCP<Operator> >("A");
428 RCP<Matrix> Acm = rcp_dynamic_cast<Matrix>(Ac);
429
430 // Record the communicator on the level
431 if (!Ac.is_null())
432 level.SetComm(Ac->getDomainMap()->getComm());
433
434 // Test if we reach the end of the hierarchy
435 bool isOrigLastLevel = isLastLevel;
436 if (isLastLevel) {
437 // Last level as we have achieved the max limit
438 isLastLevel = true;
439
440 } else if (Ac.is_null()) {
441 // Last level for this processor, as it does not belong to the next
442 // subcommunicator. Other processors may continue working on the
443 // hierarchy
444 isLastLevel = true;
445
446 } else {
447 if (!Acm.is_null() && Acm->getGlobalNumRows() <= maxCoarseSize_) {
448 // Last level as the size of the coarse matrix became too small
449 GetOStream(Runtime0) << "Max coarse size (<= " << maxCoarseSize_ << ") achieved" << std::endl;
450 isLastLevel = true;
451 if (Acm->getGlobalNumRows() != 0) setLastLevelviaMaxCoarseSize = true;
452 }
453 }
454
455 if (!Ac.is_null() && !isFinestLevel) {
456 RCP<Operator> A = Levels_[coarseLevelID-1]->template Get< RCP<Operator> >("A");
457 RCP<Matrix> Am = rcp_dynamic_cast<Matrix>(A);
458
459 const double maxCoarse2FineRatio = 0.8;
460 if (!Acm.is_null() && !Am.is_null() && Acm->getGlobalNumRows() > maxCoarse2FineRatio * Am->getGlobalNumRows()) {
461 // We could abort here, but for now we simply notify user.
462 // Couple of additional points:
463 // - if repartitioning is delayed until level K, but the aggregation
464 // procedure stagnates between levels K-1 and K. In this case,
465 // repartitioning could enable faster coarsening once again, but the
466 // hierarchy construction will abort due to the stagnation check.
467 // - if the matrix is small enough, we could move it to one processor.
468 GetOStream(Warnings0) << "Aggregation stagnated. Please check your matrix and/or adjust your configuration file."
469 << "Possible fixes:\n"
470 << " - reduce the maximum number of levels\n"
471 << " - enable repartitioning\n"
472 << " - increase the minimum coarse size." << std::endl;
473
474 }
475 }
476
477 if (isLastLevel) {
478 if (!isOrigLastLevel) {
479 // We did not expect to finish this early so we did request a smoother.
480 // We need a coarse solver instead. Do the magic.
481 level.Release(*smootherFact);
482 if (coarseFact.is_null())
483 coarseFact = rcp(new TopSmootherFactory(coarseLevelManager, "CoarseSolver"));
484 level.Request(*coarseFact);
485 }
486
487 // Do the actual build, if we have any data.
488 // NOTE: this is not a great check, we may want to call Build() regardless.
489 if (!Ac.is_null())
490 coarseFact->Build(level);
491
492 // Once the dirty deed is done, release stuff. The smoother has already
493 // been released.
494 level.Release(*coarseFact);
495
496 } else {
497 // isLastLevel = false => isOrigLastLevel = false, meaning that we have
498 // requested the smoother. Now we need to build it and to release it.
499 // We don't need to worry about the coarse solver, as we didn't request it.
500 if (!Ac.is_null())
501 smootherFact->Build(level);
502
503 level.Release(*smootherFact);
504 }
505
506 if (isLastLevel == true) {
507 int actualNumLevels = nextLevelID;
508 if (isOrigLastLevel == false) {
509 // Earlier in the function, we constructed the next coarse level, and requested data for the that level,
510 // assuming that we are not at the coarsest level. Now, we changed our mind, so we have to release those.
511 Levels_[nextLevelID]->Release(TopRAPFactory(coarseLevelManager, nextLevelManager));
512
513 // We truncate/resize the hierarchy and possibly remove the last created level if there is
514 // something wrong with it as indicated by its P not being valid. This might happen
515 // if the global number of aggregates turns out to be zero
516
517
518 if (!setLastLevelviaMaxCoarseSize) {
519 if (Levels_[nextLevelID-1]->IsAvailable("P")) {
520 if (Levels_[nextLevelID-1]->template Get<RCP<Matrix> >("P") == Teuchos::null) actualNumLevels = nextLevelID-1;
521 }
522 else actualNumLevels = nextLevelID-1;
523 }
524 }
525 if (actualNumLevels == nextLevelID-1) {
526 // Didn't expect to finish early so we requested smoother but need coarse solver instead.
527 Levels_[nextLevelID-2]->Release(*smootherFact);
528
529 if (Levels_[nextLevelID-2]->IsAvailable("PreSmoother") ) Levels_[nextLevelID-2]->RemoveKeepFlag("PreSmoother" ,NoFactory::get());
530 if (Levels_[nextLevelID-2]->IsAvailable("PostSmoother")) Levels_[nextLevelID-2]->RemoveKeepFlag("PostSmoother",NoFactory::get());
531 if (coarseFact.is_null())
532 coarseFact = rcp(new TopSmootherFactory(coarseLevelManager, "CoarseSolver"));
533 Levels_[nextLevelID-2]->Request(*coarseFact);
534 if ( !(Levels_[nextLevelID-2]->template Get<RCP<Matrix> >("A").is_null() ))
535 coarseFact->Build( *(Levels_[nextLevelID-2]));
536 Levels_[nextLevelID-2]->Release(*coarseFact);
537 }
538 Levels_.resize(actualNumLevels);
539 }
540
541 // I think this is the proper place for graph so that it shows every dependence
542 if (isDumpingEnabled_ && ( (dumpLevel_ > 0 && coarseLevelID == dumpLevel_) || dumpLevel_ == -1 ) )
543 DumpCurrentGraph(coarseLevelID);
544
545 if (!isFinestLevel) {
546 // Release the hierarchy data
547 // We release so late to help blocked solvers, as the smoothers for them need A blocks
548 // which we construct in RAPFactory
549 level.Release(coarseRAPFactory);
550 }
551
552 if (oldRank != -1)
553 SetProcRankVerbose(oldRank);
554
555 return isLastLevel;
556 }
557
558 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
560 int numLevels = Levels_.size();
561 TEUCHOS_TEST_FOR_EXCEPTION(levelManagers_.size() != numLevels, Exceptions::RuntimeError,
562 "Hierarchy::SetupRe: " << Levels_.size() << " levels, but " << levelManagers_.size() << " level factory managers");
563
564 const int startLevel = 0;
565 Clear(startLevel);
566
567#ifdef HAVE_MUELU_DEBUG
568 // Reset factories' data used for debugging
569 for (int i = 0; i < numLevels; i++)
570 levelManagers_[i]->ResetDebugData();
571
572#endif
573
574 int levelID;
575 for (levelID = startLevel; levelID < numLevels;) {
576 bool r = Setup(levelID,
577 (levelID != 0 ? levelManagers_[levelID-1] : Teuchos::null),
578 levelManagers_[levelID],
579 (levelID+1 != numLevels ? levelManagers_[levelID+1] : Teuchos::null));
580 levelID++;
581 if (r) break;
582 }
583 // We may construct fewer levels for some reason, make sure we continue
584 // doing that in the future
585 Levels_ .resize(levelID);
586 levelManagers_.resize(levelID);
587
588 int sizeOfVecs = sizeOfAllocatedLevelMultiVectors_;
589
590 AllocateLevelMultiVectors(sizeOfVecs, true);
591
592 // since the # of levels, etc. may have changed, force re-determination of description during next call to description()
593 ResetDescription();
594
595 describe(GetOStream(Statistics0), GetVerbLevel());
596 }
597
598 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
599 void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Setup(const FactoryManagerBase& manager, int startLevel, int numDesiredLevels) {
600 // Use MueLu::BaseClass::description() to avoid printing "{numLevels = 1}" (numLevels is increasing...)
601 PrintMonitor m0(*this, "Setup (" + this->MueLu::BaseClass::description() + ")", Runtime0);
602
603 Clear(startLevel);
604
605 // Check Levels_[startLevel] exists.
606 TEUCHOS_TEST_FOR_EXCEPTION(Levels_.size() <= startLevel, Exceptions::RuntimeError,
607 "MueLu::Hierarchy::Setup(): fine level (" << startLevel << ") does not exist");
608
609 TEUCHOS_TEST_FOR_EXCEPTION(numDesiredLevels <= 0, Exceptions::RuntimeError,
610 "Constructing non-positive (" << numDesiredLevels << ") number of levels does not make sense.");
611
612 // Check for fine level matrix A
613 TEUCHOS_TEST_FOR_EXCEPTION(!Levels_[startLevel]->IsAvailable("A"), Exceptions::RuntimeError,
614 "MueLu::Hierarchy::Setup(): fine level (" << startLevel << ") has no matrix A! "
615 "Set fine level matrix A using Level.Set()");
616
617 RCP<Operator> A = Levels_[startLevel]->template Get<RCP<Operator> >("A");
618 lib_ = A->getDomainMap()->lib();
619
620 if (IsPrint(Statistics2)) {
621 RCP<Matrix> Amat = rcp_dynamic_cast<Matrix>(A);
622
623 if (!Amat.is_null()) {
624 RCP<ParameterList> params = rcp(new ParameterList());
625 params->set("printLoadBalancingInfo", true);
626 params->set("printCommInfo", true);
627
628 GetOStream(Statistics2) << PerfUtils::PrintMatrixInfo(*Amat, "A0", params);
629 } else {
630 GetOStream(Warnings1) << "Fine level operator is not a matrix, statistics are not available" << std::endl;
631 }
632 }
633
634 RCP<const FactoryManagerBase> rcpmanager = rcpFromRef(manager);
635
636 const int lastLevel = startLevel + numDesiredLevels - 1;
637 GetOStream(Runtime0) << "Setup loop: startLevel = " << startLevel << ", lastLevel = " << lastLevel
638 << " (stop if numLevels = " << numDesiredLevels << " or Ac.size() < " << maxCoarseSize_ << ")" << std::endl;
639
640 // Setup multigrid levels
641 int iLevel = 0;
642 if (numDesiredLevels == 1) {
643 iLevel = 0;
644 Setup(startLevel, Teuchos::null, rcpmanager, Teuchos::null); // setup finest==coarsest level (first and last managers are Teuchos::null)
645
646 } else {
647 bool bIsLastLevel = Setup(startLevel, Teuchos::null, rcpmanager, rcpmanager); // setup finest level (level 0) (first manager is Teuchos::null)
648 if (bIsLastLevel == false) {
649 for (iLevel = startLevel + 1; iLevel < lastLevel; iLevel++) {
650 bIsLastLevel = Setup(iLevel, rcpmanager, rcpmanager, rcpmanager); // setup intermediate levels
651 if (bIsLastLevel == true)
652 break;
653 }
654 if (bIsLastLevel == false)
655 Setup(lastLevel, rcpmanager, rcpmanager, Teuchos::null); // setup coarsest level (last manager is Teuchos::null)
656 }
657 }
658
659 // TODO: some check like this should be done at the beginning of the routine
660 TEUCHOS_TEST_FOR_EXCEPTION(iLevel != Levels_.size() - 1, Exceptions::RuntimeError,
661 "MueLu::Hierarchy::Setup(): number of level");
662
663 // TODO: this is not exception safe: manager will still hold default
664 // factories if you exit this function with an exception
665 manager.Clean();
666
667 describe(GetOStream(Statistics0), GetVerbLevel());
668 }
669
670 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
672 if (startLevel < GetNumLevels())
673 GetOStream(Runtime0) << "Clearing old data (if any)" << std::endl;
674
675 for (int iLevel = startLevel; iLevel < GetNumLevels(); iLevel++)
676 Levels_[iLevel]->Clear();
677 }
678
679 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
681 GetOStream(Runtime0) << "Clearing old data (expert)" << std::endl;
682 for (int iLevel = 0; iLevel < GetNumLevels(); iLevel++)
683 Levels_[iLevel]->ExpertClear();
684 }
685
686#if defined(HAVE_MUELU_EXPERIMENTAL) && defined(HAVE_MUELU_ADDITIVE_VARIANT)
687 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
688 ConvergenceStatus Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Iterate(const MultiVector& B, MultiVector& X, ConvData conv,
689 bool InitialGuessIsZero, LO startLevel) {
690 LO nIts = conv.maxIts_;
691 MagnitudeType tol = conv.tol_;
692
693 std::string prefix = this->ShortClassName() + ": ";
694 std::string levelSuffix = " (level=" + toString(startLevel) + ")";
695 std::string levelSuffix1 = " (level=" + toString(startLevel+1) + ")";
696
697 using namespace Teuchos;
698 RCP<Time> CompTime = Teuchos::TimeMonitor::getNewCounter(prefix + "Computational Time (total)");
699 RCP<Time> Concurrent = Teuchos::TimeMonitor::getNewCounter(prefix + "Concurrent portion");
700 RCP<Time> ApplyR = Teuchos::TimeMonitor::getNewCounter(prefix + "R: Computational Time");
701 RCP<Time> ApplyPbar = Teuchos::TimeMonitor::getNewCounter(prefix + "Pbar: Computational Time");
702 RCP<Time> CompFine = Teuchos::TimeMonitor::getNewCounter(prefix + "Fine: Computational Time");
703 RCP<Time> CompCoarse = Teuchos::TimeMonitor::getNewCounter(prefix + "Coarse: Computational Time");
704 RCP<Time> ApplySum = Teuchos::TimeMonitor::getNewCounter(prefix + "Sum: Computational Time");
705 RCP<Time> Synchronize_beginning = Teuchos::TimeMonitor::getNewCounter(prefix + "Synchronize_beginning");
706 RCP<Time> Synchronize_center = Teuchos::TimeMonitor::getNewCounter(prefix + "Synchronize_center");
707 RCP<Time> Synchronize_end = Teuchos::TimeMonitor::getNewCounter(prefix + "Synchronize_end");
708
709 RCP<Level> Fine = Levels_[0];
710 RCP<Level> Coarse;
711
712 RCP<Operator> A = Fine->Get< RCP<Operator> >("A");
713 Teuchos::RCP< const Teuchos::Comm< int > > communicator = A->getDomainMap()->getComm();
714
715
716 //Synchronize_beginning->start();
717 //communicator->barrier();
718 //Synchronize_beginning->stop();
719
720 CompTime->start();
721
722 SC one = STS::one(), zero = STS::zero();
723
724 bool zeroGuess = InitialGuessIsZero;
725
726 // ======= UPFRONT DEFINITION OF COARSE VARIABLES ===========
727
728 //RCP<const Map> origMap;
729 RCP< Operator > P;
730 RCP< Operator > Pbar;
731 RCP< Operator > R;
732 RCP< MultiVector > coarseRhs, coarseX;
733 RCP< Operator > Ac;
734 RCP<SmootherBase> preSmoo_coarse, postSmoo_coarse;
735 bool emptyCoarseSolve = true;
736 RCP<MultiVector> coarseX_prolonged = MultiVectorFactory::Build(X.getMap(), X.getNumVectors(), true);
737
738 RCP<const Import> importer;
739
740 if (Levels_.size() > 1) {
741 Coarse = Levels_[1];
742 if (Coarse->IsAvailable("Importer"))
743 importer = Coarse->Get< RCP<const Import> >("Importer");
744
745 R = Coarse->Get< RCP<Operator> >("R");
746 P = Coarse->Get< RCP<Operator> >("P");
747
748
749 //if(Coarse->IsAvailable("Pbar"))
750 Pbar = Coarse->Get< RCP<Operator> >("Pbar");
751
752 coarseRhs = MultiVectorFactory::Build(R->getRangeMap(), B.getNumVectors(), true);
753
754 Ac = Coarse->Get< RCP< Operator > >("A");
755
756 ApplyR->start();
757 R->apply(B, *coarseRhs, Teuchos::NO_TRANS, one, zero);
758 //P->apply(B, *coarseRhs, Teuchos::TRANS, one, zero);
759 ApplyR->stop();
760
761 if (doPRrebalance_ || importer.is_null()) {
762 coarseX = MultiVectorFactory::Build(coarseRhs->getMap(), X.getNumVectors(), true);
763
764 } else {
765
766 RCP<TimeMonitor> ITime = rcp(new TimeMonitor(*this, prefix + "Solve : import (total)" , Timings0));
767 RCP<TimeMonitor> ILevelTime = rcp(new TimeMonitor(*this, prefix + "Solve : import" + levelSuffix1, Timings0));
768
769 // Import: range map of R --> domain map of rebalanced Ac (before subcomm replacement)
770 RCP<MultiVector> coarseTmp = MultiVectorFactory::Build(importer->getTargetMap(), coarseRhs->getNumVectors());
771 coarseTmp->doImport(*coarseRhs, *importer, Xpetra::INSERT);
772 coarseRhs.swap(coarseTmp);
773
774 coarseX = MultiVectorFactory::Build(importer->getTargetMap(), X.getNumVectors(), true);
775 }
776
777 if (Coarse->IsAvailable("PreSmoother"))
778 preSmoo_coarse = Coarse->Get< RCP<SmootherBase> >("PreSmoother");
779 if (Coarse->IsAvailable("PostSmoother"))
780 postSmoo_coarse = Coarse->Get< RCP<SmootherBase> >("PostSmoother");
781 }
782
783 // ==========================================================
784
785
786 MagnitudeType prevNorm = STS::magnitude(STS::one()), curNorm = STS::magnitude(STS::one());
787 rate_ = 1.0;
788
789 for (LO i = 1; i <= nIts; i++) {
790#ifdef HAVE_MUELU_DEBUG
791 if (A->getDomainMap()->isCompatible(*(X.getMap())) == false) {
792 std::ostringstream ss;
793 ss << "Level " << startLevel << ": level A's domain map is not compatible with X";
794 throw Exceptions::Incompatible(ss.str());
795 }
796
797 if (A->getRangeMap()->isCompatible(*(B.getMap())) == false) {
798 std::ostringstream ss;
799 ss << "Level " << startLevel << ": level A's range map is not compatible with B";
800 throw Exceptions::Incompatible(ss.str());
801 }
802#endif
803 }
804
805 bool emptyFineSolve = true;
806
807 RCP< MultiVector > fineX;
808 fineX = MultiVectorFactory::Build(X.getMap(), X.getNumVectors(), true);
809
810 //Synchronize_center->start();
811 //communicator->barrier();
812 //Synchronize_center->stop();
813
814 Concurrent->start();
815
816 // NOTE: we need to check using IsAvailable before Get here to avoid building default smoother
817 if (Fine->IsAvailable("PreSmoother")) {
818 RCP<SmootherBase> preSmoo = Fine->Get< RCP<SmootherBase> >("PreSmoother");
819 CompFine->start();
820 preSmoo->Apply(*fineX, B, zeroGuess);
821 CompFine->stop();
822 emptyFineSolve = false;
823 }
824 if (Fine->IsAvailable("PostSmoother")) {
825 RCP<SmootherBase> postSmoo = Fine->Get< RCP<SmootherBase> >("PostSmoother");
826 CompFine->start();
827 postSmoo->Apply(*fineX, B, zeroGuess);
828 CompFine->stop();
829
830 emptyFineSolve = false;
831 }
832 if (emptyFineSolve == true) {
833 GetOStream(Warnings1) << "No fine grid smoother" << std::endl;
834 // Fine grid smoother is identity
835 fineX->update(one, B, zero);
836 }
837
838 if (Levels_.size() > 1) {
839 // NOTE: we need to check using IsAvailable before Get here to avoid building default smoother
840 if (Coarse->IsAvailable("PreSmoother")) {
841 CompCoarse->start();
842 preSmoo_coarse->Apply(*coarseX, *coarseRhs, zeroGuess);
843 CompCoarse->stop();
844 emptyCoarseSolve = false;
845
846 }
847 if (Coarse->IsAvailable("PostSmoother")) {
848 CompCoarse->start();
849 postSmoo_coarse->Apply(*coarseX, *coarseRhs, zeroGuess);
850 CompCoarse->stop();
851 emptyCoarseSolve = false;
852
853 }
854 if (emptyCoarseSolve == true) {
855 GetOStream(Warnings1) << "No coarse grid solver" << std::endl;
856 // Coarse operator is identity
857 coarseX->update(one, *coarseRhs, zero);
858 }
859 Concurrent->stop();
860 //Synchronize_end->start();
861 //communicator->barrier();
862 //Synchronize_end->stop();
863
864 if (!doPRrebalance_ && !importer.is_null()) {
865 RCP<TimeMonitor> ITime = rcp(new TimeMonitor(*this, prefix + "Solve : export (total)" , Timings0));
866 RCP<TimeMonitor> ILevelTime = rcp(new TimeMonitor(*this, prefix + "Solve : export" + levelSuffix1, Timings0));
867
868 // Import: range map of rebalanced Ac (before subcomm replacement) --> domain map of P
869 RCP<MultiVector> coarseTmp = MultiVectorFactory::Build(importer->getSourceMap(), coarseX->getNumVectors());
870 coarseTmp->doExport(*coarseX, *importer, Xpetra::INSERT);
871 coarseX.swap(coarseTmp);
872 }
873
874 ApplyPbar->start();
875 Pbar->apply(*coarseX, *coarseX_prolonged, Teuchos::NO_TRANS, one, zero);
876 ApplyPbar->stop();
877 }
878
879 ApplySum->start();
880 X.update(1.0, *fineX, 1.0, *coarseX_prolonged, 0.0);
881 ApplySum->stop();
882
883 CompTime->stop();
884
885 //communicator->barrier();
886
888}
889#else
890 // ---------------------------------------- Iterate -------------------------------------------------------
891 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
893 bool InitialGuessIsZero, LO startLevel) {
894 LO nIts = conv.maxIts_;
895 MagnitudeType tol = conv.tol_;
896
897 // These timers work as follows. "iterateTime" records total time spent in
898 // iterate. "levelTime" records time on a per level basis. The label is
899 // crafted to mimic the per-level messages used in Monitors. Note that a
900 // given level is timed with a TimeMonitor instead of a Monitor or
901 // SubMonitor. This is mainly because I want to time each level
902 // separately, and Monitors/SubMonitors print "(total) xx yy zz" ,
903 // "(sub,total) xx yy zz", respectively, which is subject to
904 // misinterpretation. The per-level TimeMonitors are stopped/started
905 // manually before/after a recursive call to Iterate. A side artifact to
906 // this approach is that the counts for intermediate level timers are twice
907 // the counts for the finest and coarsest levels.
908
909 RCP<Level> Fine = Levels_[startLevel];
910
911 std::string label = FormattingHelper::getColonLabel(Fine->getObjectLabel());
912 std::string prefix = label + this->ShortClassName() + ": ";
913 std::string levelSuffix = " (level=" + toString(startLevel) + ")";
914 std::string levelSuffix1 = " (level=" + toString(startLevel+1) + ")";
915
916 bool useStackedTimer = !Teuchos::TimeMonitor::getStackedTimer().is_null();
917
918 RCP<Monitor> iterateTime;
919 RCP<TimeMonitor> iterateTime1;
920 if (startLevel == 0)
921 iterateTime = rcp(new Monitor(*this, "Solve", label, (nIts == 1) ? None : Runtime0, Timings0));
922 else if (!useStackedTimer)
923 iterateTime1 = rcp(new TimeMonitor(*this, prefix + "Solve (total, level=" + toString(startLevel) + ")", Timings0));
924
925 std::string iterateLevelTimeLabel = prefix + "Solve" + levelSuffix;
926 RCP<TimeMonitor> iterateLevelTime = rcp(new TimeMonitor(*this, iterateLevelTimeLabel, Timings0));
927
928 bool zeroGuess = InitialGuessIsZero;
929
930 RCP<Operator> A = Fine->Get< RCP<Operator> >("A");
931 using namespace Teuchos;
932 RCP<Time> CompCoarse = Teuchos::TimeMonitor::getNewCounter(prefix + "Coarse: Computational Time");
933
934 if (A.is_null()) {
935 // This processor does not have any data for this process on coarser
936 // levels. This can only happen when there are multiple processors and
937 // we use repartitioning.
939 }
940
941 // If we switched the number of vectors, we'd need to reallocate here.
942 // If the number of vectors is unchanged, this is a noop.
943 // NOTE: We need to check against B because the tests in AllocateLevelMultiVectors
944 // will fail on Stokhos Scalar types (due to the so-called 'hidden dimension')
945 const BlockedMultiVector * Bblocked = dynamic_cast<const BlockedMultiVector*>(&B);
946 if(residual_.size() > startLevel &&
947 ( ( Bblocked && !Bblocked->isSameSize(*residual_[startLevel])) ||
948 (!Bblocked && !residual_[startLevel]->isSameSize(B))))
949 DeleteLevelMultiVectors();
950 AllocateLevelMultiVectors(X.getNumVectors());
951
952 // Print residual information before iterating
953 typedef Teuchos::ScalarTraits<typename STS::magnitudeType> STM;
954 MagnitudeType prevNorm = STM::one();
955 rate_ = 1.0;
956 if (IsCalculationOfResidualRequired(startLevel, conv)) {
957 ConvergenceStatus convergenceStatus = ComputeResidualAndPrintHistory(*A, X, B, Teuchos::ScalarTraits<LO>::zero(), startLevel, conv, prevNorm);
958 if (convergenceStatus == MueLu::ConvergenceStatus::Converged)
959 return convergenceStatus;
960 }
961
962 SC one = STS::one(), zero = STS::zero();
963 for (LO iteration = 1; iteration <= nIts; iteration++) {
964#ifdef HAVE_MUELU_DEBUG
965#if 0 // TODO fix me
966 if (A->getDomainMap()->isCompatible(*(X.getMap())) == false) {
967 std::ostringstream ss;
968 ss << "Level " << startLevel << ": level A's domain map is not compatible with X";
969 throw Exceptions::Incompatible(ss.str());
970 }
971
972 if (A->getRangeMap()->isCompatible(*(B.getMap())) == false) {
973 std::ostringstream ss;
974 ss << "Level " << startLevel << ": level A's range map is not compatible with B";
975 throw Exceptions::Incompatible(ss.str());
976 }
977#endif
978#endif
979
980 if (startLevel == as<LO>(Levels_.size())-1) {
981 // On the coarsest level, we do either smoothing (if defined) or a direct solve.
982 RCP<TimeMonitor> CLevelTime = rcp(new TimeMonitor(*this, prefix + "Solve : coarse" + levelSuffix, Timings0));
983
984 bool emptySolve = true;
985
986 // NOTE: we need to check using IsAvailable before Get here to avoid building default smoother
987 if (Fine->IsAvailable("PreSmoother")) {
988 RCP<SmootherBase> preSmoo = Fine->Get< RCP<SmootherBase> >("PreSmoother");
989 CompCoarse->start();
990 preSmoo->Apply(X, B, zeroGuess);
991 CompCoarse->stop();
992 zeroGuess = false;
993 emptySolve = false;
994 }
995 if (Fine->IsAvailable("PostSmoother")) {
996 RCP<SmootherBase> postSmoo = Fine->Get< RCP<SmootherBase> >("PostSmoother");
997 CompCoarse->start();
998 postSmoo->Apply(X, B, zeroGuess);
999 CompCoarse->stop();
1000 emptySolve = false;
1001 zeroGuess = false;
1002 }
1003 if (emptySolve == true) {
1004 GetOStream(Warnings1) << "No coarse grid solver" << std::endl;
1005 // Coarse operator is identity
1006 X.update(one, B, zero);
1007 }
1008
1009 } else {
1010 // On intermediate levels, we do cycles
1011 RCP<Level> Coarse = Levels_[startLevel+1];
1012 {
1013 // ============== PRESMOOTHING ==============
1014 RCP<TimeMonitor> STime;
1015 if (!useStackedTimer)
1016 STime = rcp(new TimeMonitor(*this, prefix + "Solve : smoothing (total)" , Timings0));
1017 RCP<TimeMonitor> SLevelTime = rcp(new TimeMonitor(*this, prefix + "Solve : smoothing" + levelSuffix, Timings0));
1018
1019 if (Fine->IsAvailable("PreSmoother")) {
1020 RCP<SmootherBase> preSmoo = Fine->Get< RCP<SmootherBase> >("PreSmoother");
1021 preSmoo->Apply(X, B, zeroGuess);
1022 zeroGuess = false;
1023 }
1024 }
1025
1026 RCP<MultiVector> residual;
1027 {
1028 RCP<TimeMonitor> ATime;
1029 if (!useStackedTimer)
1030 ATime = rcp(new TimeMonitor(*this, prefix + "Solve : residual calculation (total)" , Timings0));
1031 RCP<TimeMonitor> ALevelTime = rcp(new TimeMonitor(*this, prefix + "Solve : residual calculation" + levelSuffix, Timings0));
1032 if (zeroGuess) {
1033 // If there's a pre-smoother, then zeroGuess is false. If there isn't and people still have zeroGuess set,
1034 // then then X still has who knows what, so we need to zero it before we go to the coarse grid.
1035 X.putScalar(zero);
1036 }
1037
1038 Utilities::Residual(*A, X, B,*residual_[startLevel]);
1039 residual = residual_[startLevel];
1040 }
1041
1042 RCP<Operator> P = Coarse->Get< RCP<Operator> >("P");
1043 if (Coarse->IsAvailable("Pbar"))
1044 P = Coarse->Get< RCP<Operator> >("Pbar");
1045
1046 RCP<MultiVector> coarseRhs, coarseX;
1047 // const bool initializeWithZeros = true;
1048 {
1049 // ============== RESTRICTION ==============
1050 RCP<TimeMonitor> RTime;
1051 if (!useStackedTimer)
1052 RTime = rcp(new TimeMonitor(*this, prefix + "Solve : restriction (total)" , Timings0));
1053 RCP<TimeMonitor> RLevelTime = rcp(new TimeMonitor(*this, prefix + "Solve : restriction" + levelSuffix, Timings0));
1054 coarseRhs = coarseRhs_[startLevel];
1055
1056 if (implicitTranspose_) {
1057 P->apply(*residual, *coarseRhs, Teuchos::TRANS, one, zero);
1058
1059 } else {
1060 RCP<Operator> R = Coarse->Get< RCP<Operator> >("R");
1061 R->apply(*residual, *coarseRhs, Teuchos::NO_TRANS, one, zero);
1062 }
1063 }
1064
1065 RCP<const Import> importer;
1066 if (Coarse->IsAvailable("Importer"))
1067 importer = Coarse->Get< RCP<const Import> >("Importer");
1068
1069 coarseX = coarseX_[startLevel];
1070 if (!doPRrebalance_ && !importer.is_null()) {
1071 RCP<TimeMonitor> ITime;
1072 if (!useStackedTimer)
1073 ITime = rcp(new TimeMonitor(*this, prefix + "Solve : import (total)" , Timings0));
1074 RCP<TimeMonitor> ILevelTime = rcp(new TimeMonitor(*this, prefix + "Solve : import" + levelSuffix1, Timings0));
1075
1076 // Import: range map of R --> domain map of rebalanced Ac (before subcomm replacement)
1077 RCP<MultiVector> coarseTmp = coarseImport_[startLevel];
1078 coarseTmp->doImport(*coarseRhs, *importer, Xpetra::INSERT);
1079 coarseRhs.swap(coarseTmp);
1080 }
1081
1082 RCP<Operator> Ac = Coarse->Get< RCP<Operator> >("A");
1083 if (!Ac.is_null()) {
1084 RCP<const Map> origXMap = coarseX->getMap();
1085 RCP<const Map> origRhsMap = coarseRhs->getMap();
1086
1087 // Replace maps with maps with a subcommunicator
1088 coarseRhs->replaceMap(Ac->getRangeMap());
1089 coarseX ->replaceMap(Ac->getDomainMap());
1090
1091 {
1092 iterateLevelTime = Teuchos::null; // stop timing this level
1093
1094 Iterate(*coarseRhs, *coarseX, 1, true, startLevel+1);
1095 // ^^ zero initial guess
1096 if (Cycle_ == WCYCLE && WCycleStartLevel_ >= startLevel)
1097 Iterate(*coarseRhs, *coarseX, 1, false, startLevel+1);
1098 // ^^ nonzero initial guess
1099
1100 iterateLevelTime = rcp(new TimeMonitor(*this, iterateLevelTimeLabel)); // restart timing this level
1101 }
1102 coarseX->replaceMap(origXMap);
1103 coarseRhs->replaceMap(origRhsMap);
1104 }
1105
1106 if (!doPRrebalance_ && !importer.is_null()) {
1107 RCP<TimeMonitor> ITime;
1108 if (!useStackedTimer)
1109 ITime = rcp(new TimeMonitor(*this, prefix + "Solve : export (total)" , Timings0));
1110 RCP<TimeMonitor> ILevelTime = rcp(new TimeMonitor(*this, prefix + "Solve : export" + levelSuffix1, Timings0));
1111
1112 // Import: range map of rebalanced Ac (before subcomm replacement) --> domain map of P
1113 RCP<MultiVector> coarseTmp = coarseExport_[startLevel];
1114 coarseTmp->doExport(*coarseX, *importer, Xpetra::INSERT);
1115 coarseX.swap(coarseTmp);
1116 }
1117
1118 {
1119 // ============== PROLONGATION ==============
1120 RCP<TimeMonitor> PTime;
1121 if (!useStackedTimer)
1122 PTime = rcp(new TimeMonitor(*this, prefix + "Solve : prolongation (total)" , Timings0));
1123 RCP<TimeMonitor> PLevelTime = rcp(new TimeMonitor(*this, prefix + "Solve : prolongation" + levelSuffix, Timings0));
1124 // Update X += P * coarseX
1125 // Note that due to what may be round-off error accumulation, use of the fused kernel
1126 // P->apply(*coarseX, X, Teuchos::NO_TRANS, one, one);
1127 // can in some cases result in slightly higher iteration counts.
1128 if (fuseProlongationAndUpdate_) {
1129 P->apply(*coarseX, X, Teuchos::NO_TRANS, scalingFactor_, one);
1130 } else {
1131 RCP<MultiVector> correction = correction_[startLevel];
1132 P->apply(*coarseX, *correction, Teuchos::NO_TRANS, one, zero);
1133 X.update(scalingFactor_, *correction, one);
1134 }
1135 }
1136
1137 {
1138 // ============== POSTSMOOTHING ==============
1139 RCP<TimeMonitor> STime;
1140 if (!useStackedTimer)
1141 STime = rcp(new TimeMonitor(*this, prefix + "Solve : smoothing (total)" , Timings0));
1142 RCP<TimeMonitor> SLevelTime = rcp(new TimeMonitor(*this, prefix + "Solve : smoothing" + levelSuffix, Timings0));
1143
1144 if (Fine->IsAvailable("PostSmoother")) {
1145 RCP<SmootherBase> postSmoo = Fine->Get< RCP<SmootherBase> >("PostSmoother");
1146 postSmoo->Apply(X, B, false);
1147 }
1148 }
1149 }
1150 zeroGuess = false;
1151
1152
1153 if (IsCalculationOfResidualRequired(startLevel, conv)) {
1154 ConvergenceStatus convergenceStatus = ComputeResidualAndPrintHistory(*A, X, B, iteration, startLevel, conv, prevNorm);
1155 if (convergenceStatus == MueLu::ConvergenceStatus::Converged)
1156 return convergenceStatus;
1157 }
1158 }
1160 }
1161#endif
1162
1163
1164 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1165 void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write(const LO& start, const LO& end, const std::string &suffix) {
1166 LO startLevel = (start != -1 ? start : 0);
1167 LO endLevel = (end != -1 ? end : Levels_.size()-1);
1168
1169 TEUCHOS_TEST_FOR_EXCEPTION(startLevel > endLevel, Exceptions::RuntimeError,
1170 "MueLu::Hierarchy::Write : startLevel must be <= endLevel");
1171
1172 TEUCHOS_TEST_FOR_EXCEPTION(startLevel < 0 || endLevel >= Levels_.size(), Exceptions::RuntimeError,
1173 "MueLu::Hierarchy::Write bad start or end level");
1174
1175 for (LO i = startLevel; i < endLevel + 1; i++) {
1176 RCP<Matrix> A = rcp_dynamic_cast<Matrix>(Levels_[i]-> template Get< RCP< Operator> >("A")), P, R;
1177 if (i > 0) {
1178 P = rcp_dynamic_cast<Matrix>(Levels_[i]-> template Get< RCP< Operator> >("P"));
1179 if (!implicitTranspose_)
1180 R = rcp_dynamic_cast<Matrix>(Levels_[i]-> template Get< RCP< Operator> >("R"));
1181 }
1182
1183 if (!A.is_null()) Xpetra::IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write("A_" + toString(i) + suffix + ".m", *A);
1184 if (!P.is_null()) {
1185 Xpetra::IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write("P_" + toString(i) + suffix + ".m", *P);
1186 }
1187 if (!R.is_null()) {
1188 Xpetra::IO<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write("R_" + toString(i) + suffix + ".m", *R);
1189 }
1190 }
1191 }
1192
1193 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1194 void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Keep(const std::string & ename, const FactoryBase* factory) {
1195 for (Array<RCP<Level> >::iterator it = Levels_.begin(); it != Levels_.end(); ++it)
1196 (*it)->Keep(ename, factory);
1197 }
1198
1199 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1200 void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Delete(const std::string& ename, const FactoryBase* factory) {
1201 for (Array<RCP<Level> >::iterator it = Levels_.begin(); it != Levels_.end(); ++it)
1202 (*it)->Delete(ename, factory);
1203 }
1204
1205 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1207 for (Array<RCP<Level> >::iterator it = Levels_.begin(); it != Levels_.end(); ++it)
1208 (*it)->AddKeepFlag(ename, factory, keep);
1209 }
1210
1211 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1213 for (Array<RCP<Level> >::iterator it = Levels_.begin(); it != Levels_.end(); ++it)
1214 (*it)->RemoveKeepFlag(ename, factory, keep);
1215 }
1216
1217 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1219 if ( description_.empty() )
1220 {
1221 std::ostringstream out;
1222 out << BaseClass::description();
1223 out << "{#levels = " << GetGlobalNumLevels() << ", complexity = " << GetOperatorComplexity() << "}";
1224 description_ = out.str();
1225 }
1226 return description_;
1227 }
1228
1229 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1230 void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::describe(Teuchos::FancyOStream& out, const Teuchos::EVerbosityLevel tVerbLevel) const {
1231 describe(out, toMueLuVerbLevel(tVerbLevel));
1232 }
1233
1234 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1235 void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::describe(Teuchos::FancyOStream& out, const VerbLevel verbLevel) const {
1236 RCP<Operator> A0 = Levels_[0]->template Get<RCP<Operator> >("A");
1237 RCP<const Teuchos::Comm<int> > comm = A0->getDomainMap()->getComm();
1238
1239 int numLevels = GetNumLevels();
1240 RCP<Operator> Ac = Levels_[numLevels-1]->template Get<RCP<Operator> >("A");
1241 if (Ac.is_null()) {
1242 // It may happen that we do repartition on the last level, but the matrix
1243 // is small enough to satisfy "max coarse size" requirement. Then, even
1244 // though we have the level, the matrix would be null on all but one processors
1245 numLevels--;
1246 }
1247 int root = comm->getRank();
1248
1249#ifdef HAVE_MPI
1250 int smartData = numLevels*comm->getSize() + comm->getRank(), maxSmartData;
1251 reduceAll(*comm, Teuchos::REDUCE_MAX, smartData, Teuchos::ptr(&maxSmartData));
1252 root = maxSmartData % comm->getSize();
1253#endif
1254
1255 // Compute smoother complexity, if needed
1256 double smoother_comp =-1.0;
1257 if (verbLevel & (Statistics0 | Test))
1258 smoother_comp = GetSmootherComplexity();
1259
1260 std::string outstr;
1261 if (comm->getRank() == root && verbLevel & (Statistics0 | Test)) {
1262 std::vector<Xpetra::global_size_t> nnzPerLevel;
1263 std::vector<Xpetra::global_size_t> rowsPerLevel;
1264 std::vector<int> numProcsPerLevel;
1265 bool someOpsNotMatrices = false;
1266 const Xpetra::global_size_t INVALID = Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid();
1267 for (int i = 0; i < numLevels; i++) {
1268 TEUCHOS_TEST_FOR_EXCEPTION(!(Levels_[i]->IsAvailable("A")) , Exceptions::RuntimeError,
1269 "Operator A is not available on level " << i);
1270
1271 RCP<Operator> A = Levels_[i]->template Get<RCP<Operator> >("A");
1272 TEUCHOS_TEST_FOR_EXCEPTION(A.is_null(), Exceptions::RuntimeError,
1273 "Operator A on level " << i << " is null.");
1274
1275 RCP<Matrix> Am = rcp_dynamic_cast<Matrix>(A);
1276 if (Am.is_null()) {
1277 someOpsNotMatrices = true;
1278 nnzPerLevel .push_back(INVALID);
1279 rowsPerLevel .push_back(A->getDomainMap()->getGlobalNumElements());
1280 numProcsPerLevel.push_back(A->getDomainMap()->getComm()->getSize());
1281 } else {
1282 LO storageblocksize=Am->GetStorageBlockSize();
1283 Xpetra::global_size_t nnz = Am->getGlobalNumEntries()*storageblocksize*storageblocksize;
1284 nnzPerLevel .push_back(nnz);
1285 rowsPerLevel .push_back(Am->getGlobalNumRows()*storageblocksize);
1286 numProcsPerLevel.push_back(Am->getRowMap()->getComm()->getSize());
1287 }
1288 }
1289 if (someOpsNotMatrices)
1290 GetOStream(Warnings0) << "Some level operators are not matrices, statistics calculation are incomplete" << std::endl;
1291
1292 {
1293 std::string label = Levels_[0]->getObjectLabel();
1294 std::ostringstream oss;
1295 oss << std::setfill(' ');
1296 oss << "\n--------------------------------------------------------------------------------\n";
1297 oss << "--- Multigrid Summary " << std::setw(28) << std::left << label << "---\n";
1298 oss << "--------------------------------------------------------------------------------" << std::endl;
1299 if (verbLevel & Parameters1)
1300 oss << "Scalar = " << Teuchos::ScalarTraits<Scalar>::name() << std::endl;
1301 oss << "Number of levels = " << numLevels << std::endl;
1302 oss << "Operator complexity = " << std::setprecision(2) << std::setiosflags(std::ios::fixed);
1303 if (!someOpsNotMatrices)
1304 oss << GetOperatorComplexity() << std::endl;
1305 else
1306 oss << "not available (Some operators in hierarchy are not matrices.)" << std::endl;
1307
1308 if(smoother_comp!=-1.0) {
1309 oss << "Smoother complexity = " << std::setprecision(2) << std::setiosflags(std::ios::fixed)
1310 << smoother_comp << std::endl;
1311 }
1312
1313 switch (Cycle_) {
1314 case VCYCLE:
1315 oss << "Cycle type = V" << std::endl;
1316 break;
1317 case WCYCLE:
1318 oss << "Cycle type = W" << std::endl;
1319 if (WCycleStartLevel_ > 0)
1320 oss << "Cycle start level = " << WCycleStartLevel_ << std::endl;
1321 break;
1322 default:
1323 break;
1324 };
1325 oss << std::endl;
1326
1327 Xpetra::global_size_t tt = rowsPerLevel[0];
1328 int rowspacer = 2; while (tt != 0) { tt /= 10; rowspacer++; }
1329 for (size_t i = 0; i < nnzPerLevel.size(); ++i) {
1330 tt = nnzPerLevel[i];
1331 if (tt != INVALID)
1332 break;
1333 tt = 100; // This will get used if all levels are operators.
1334 }
1335 int nnzspacer = 2; while (tt != 0) { tt /= 10; nnzspacer++; }
1336 tt = numProcsPerLevel[0];
1337 int npspacer = 2; while (tt != 0) { tt /= 10; npspacer++; }
1338 oss << "level " << std::setw(rowspacer) << " rows " << std::setw(nnzspacer) << " nnz " << " nnz/row" << std::setw(npspacer) << " c ratio" << " procs" << std::endl;
1339 for (size_t i = 0; i < nnzPerLevel.size(); ++i) {
1340 oss << " " << i << " ";
1341 oss << std::setw(rowspacer) << rowsPerLevel[i];
1342 if (nnzPerLevel[i] != INVALID) {
1343 oss << std::setw(nnzspacer) << nnzPerLevel[i];
1344 oss << std::setprecision(2) << std::setiosflags(std::ios::fixed);
1345 oss << std::setw(9) << as<double>(nnzPerLevel[i]) / rowsPerLevel[i];
1346 } else {
1347 oss << std::setw(nnzspacer) << "Operator";
1348 oss << std::setprecision(2) << std::setiosflags(std::ios::fixed);
1349 oss << std::setw(9) << " ";
1350 }
1351 if (i) oss << std::setw(9) << as<double>(rowsPerLevel[i-1])/rowsPerLevel[i];
1352 else oss << std::setw(9) << " ";
1353 oss << " " << std::setw(npspacer) << numProcsPerLevel[i] << std::endl;
1354 }
1355 oss << std::endl;
1356 for (int i = 0; i < GetNumLevels(); ++i) {
1357 RCP<SmootherBase> preSmoo, postSmoo;
1358 if (Levels_[i]->IsAvailable("PreSmoother"))
1359 preSmoo = Levels_[i]->template Get< RCP<SmootherBase> >("PreSmoother");
1360 if (Levels_[i]->IsAvailable("PostSmoother"))
1361 postSmoo = Levels_[i]->template Get< RCP<SmootherBase> >("PostSmoother");
1362
1363 if (preSmoo != null && preSmoo == postSmoo)
1364 oss << "Smoother (level " << i << ") both : " << preSmoo->description() << std::endl;
1365 else {
1366 oss << "Smoother (level " << i << ") pre : "
1367 << (preSmoo != null ? preSmoo->description() : "no smoother") << std::endl;
1368 oss << "Smoother (level " << i << ") post : "
1369 << (postSmoo != null ? postSmoo->description() : "no smoother") << std::endl;
1370 }
1371
1372 oss << std::endl;
1373 }
1374
1375 outstr = oss.str();
1376 }
1377 }
1378
1379#ifdef HAVE_MPI
1380 RCP<const Teuchos::MpiComm<int> > mpiComm = rcp_dynamic_cast<const Teuchos::MpiComm<int> >(comm);
1381 MPI_Comm rawComm = (*mpiComm->getRawMpiComm())();
1382
1383 int strLength = outstr.size();
1384 MPI_Bcast(&strLength, 1, MPI_INT, root, rawComm);
1385 if (comm->getRank() != root)
1386 outstr.resize(strLength);
1387 MPI_Bcast(&outstr[0], strLength, MPI_CHAR, root, rawComm);
1388#endif
1389
1390 out << outstr;
1391 }
1392
1393 // NOTE: at some point this should be replaced by a friend operator <<
1394 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1395 void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::print(std::ostream& out, const VerbLevel verbLevel) const {
1396 Teuchos::OSTab tab2(out);
1397 for (int i = 0; i < GetNumLevels(); ++i)
1398 Levels_[i]->print(out, verbLevel);
1399 }
1400
1401 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1403 isPreconditioner_ = flag;
1404 }
1405
1406 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1408 if (GetProcRankVerbose() != 0)
1409 return;
1410#if defined(HAVE_MUELU_BOOST) && defined(HAVE_MUELU_BOOST_FOR_REAL) && defined(BOOST_VERSION) && (BOOST_VERSION >= 104400)
1411
1412 BoostGraph graph;
1413
1414 BoostProperties dp;
1415 dp.property("label", boost::get(boost::vertex_name, graph));
1416 dp.property("id", boost::get(boost::vertex_index, graph));
1417 dp.property("label", boost::get(boost::edge_name, graph));
1418 dp.property("color", boost::get(boost::edge_color, graph));
1419
1420 // create local maps
1421 std::map<const FactoryBase*, BoostVertex> vindices;
1422 typedef std::map<std::pair<BoostVertex,BoostVertex>, std::string> emap; emap edges;
1423
1424 static int call_id=0;
1425
1426 RCP<Operator> A = Levels_[0]->template Get<RCP<Operator> >("A");
1427 int rank = A->getDomainMap()->getComm()->getRank();
1428
1429 // printf("[%d] CMS: ----------------------\n",rank);
1430 for (int i = currLevel; i <= currLevel+1 && i < GetNumLevels(); i++) {
1431 edges.clear();
1432 Levels_[i]->UpdateGraph(vindices, edges, dp, graph);
1433
1434 for (emap::const_iterator eit = edges.begin(); eit != edges.end(); eit++) {
1435 std::pair<BoostEdge, bool> boost_edge = boost::add_edge(eit->first.first, eit->first.second, graph);
1436 // printf("[%d] CMS: Hierarchy, adding edge (%d->%d) %d\n",rank,(int)eit->first.first,(int)eit->first.second,(int)boost_edge.second);
1437 // Because xdot.py views 'Graph' as a keyword
1438 if(eit->second==std::string("Graph")) boost::put("label", dp, boost_edge.first, std::string("Graph_"));
1439 else boost::put("label", dp, boost_edge.first, eit->second);
1440 if (i == currLevel)
1441 boost::put("color", dp, boost_edge.first, std::string("red"));
1442 else
1443 boost::put("color", dp, boost_edge.first, std::string("blue"));
1444 }
1445 }
1446
1447 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"));
1448 boost::write_graphviz_dp(out, graph, dp, std::string("id"));
1449 out.close();
1450 call_id++;
1451#else
1452 GetOStream(Errors) << "Dependency graph output requires boost and MueLu_ENABLE_Boost_for_real" << std::endl;
1453#endif
1454 }
1455
1456 // Enforce that coordinate vector's map is consistent with that of A
1457 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1459 RCP<Operator> Ao = level.Get<RCP<Operator> >("A");
1460 RCP<Matrix> A = rcp_dynamic_cast<Matrix>(Ao);
1461 if (A.is_null()) {
1462 GetOStream(Runtime1) << "Hierarchy::ReplaceCoordinateMap: operator is not a matrix, skipping..." << std::endl;
1463 return;
1464 }
1465 if(Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A) != Teuchos::null) {
1466 GetOStream(Runtime1) << "Hierarchy::ReplaceCoordinateMap: operator is a BlockedCrsMatrix, skipping..." << std::endl;
1467 return;
1468 }
1469
1470 typedef Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType,LO,GO,NO> xdMV;
1471
1472 RCP<xdMV> coords = level.Get<RCP<xdMV> >("Coordinates");
1473
1474 if (A->getRowMap()->isSameAs(*(coords->getMap()))) {
1475 GetOStream(Runtime1) << "Hierarchy::ReplaceCoordinateMap: matrix and coordinates maps are same, skipping..." << std::endl;
1476 return;
1477 }
1478
1479 if (A->IsView("stridedMaps") && rcp_dynamic_cast<const StridedMap>(A->getRowMap("stridedMaps")) != Teuchos::null) {
1480 RCP<const StridedMap> stridedRowMap = rcp_dynamic_cast<const StridedMap>(A->getRowMap("stridedMaps"));
1481
1482 // It is better to through an exceptions if maps may be inconsistent, than to ignore it and experience unfathomable breakdowns
1483 TEUCHOS_TEST_FOR_EXCEPTION(stridedRowMap->getStridedBlockId() != -1 || stridedRowMap->getOffset() != 0,
1484 Exceptions::RuntimeError, "Hierarchy::ReplaceCoordinateMap: nontrivial maps (block id = " << stridedRowMap->getStridedBlockId()
1485 << ", offset = " << stridedRowMap->getOffset() << ")");
1486 }
1487
1488 GetOStream(Runtime1) << "Replacing coordinate map" << std::endl;
1489 TEUCHOS_TEST_FOR_EXCEPTION(A->GetFixedBlockSize() % A->GetStorageBlockSize() != 0, Exceptions::RuntimeError, "Hierarchy::ReplaceCoordinateMap: Storage block size does not evenly divide fixed block size");
1490
1491 size_t blkSize = A->GetFixedBlockSize() / A->GetStorageBlockSize();
1492
1493 RCP<const Map> nodeMap = A->getRowMap();
1494 if (blkSize > 1) {
1495 // Create a nodal map, as coordinates have not been expanded to a DOF map yet.
1496 RCP<const Map> dofMap = A->getRowMap();
1497 GO indexBase = dofMap->getIndexBase();
1498 size_t numLocalDOFs = dofMap->getLocalNumElements();
1499 TEUCHOS_TEST_FOR_EXCEPTION(numLocalDOFs % blkSize, Exceptions::RuntimeError,
1500 "Hierarchy::ReplaceCoordinateMap: block size (" << blkSize << ") is incompatible with the number of local dofs in a row map (" << numLocalDOFs);
1501 ArrayView<const GO> GIDs = dofMap->getLocalElementList();
1502
1503 Array<GO> nodeGIDs(numLocalDOFs/blkSize);
1504 for (size_t i = 0; i < numLocalDOFs; i += blkSize)
1505 nodeGIDs[i/blkSize] = (GIDs[i] - indexBase)/blkSize + indexBase;
1506
1507 Xpetra::global_size_t INVALID = Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid();
1508 nodeMap = MapFactory::Build(dofMap->lib(), INVALID, nodeGIDs(), indexBase, dofMap->getComm());
1509 } else {
1510 // blkSize == 1
1511 // Check whether the length of vectors fits to the size of A
1512 // If yes, make sure that the maps are matching
1513 // If no, throw a warning but do not touch the Coordinates
1514 if(coords->getLocalLength() != A->getRowMap()->getLocalNumElements()) {
1515 GetOStream(Warnings) << "Coordinate vector does not match row map of matrix A!" << std::endl;
1516 return;
1517 }
1518 }
1519
1520 Array<ArrayView<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> > coordDataView;
1521 std::vector<ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> > coordData;
1522 for (size_t i = 0; i < coords->getNumVectors(); i++) {
1523 coordData.push_back(coords->getData(i));
1524 coordDataView.push_back(coordData[i]());
1525 }
1526
1527 RCP<xdMV> newCoords = Xpetra::MultiVectorFactory<typename Teuchos::ScalarTraits<Scalar>::coordinateType,LO,GO,NO>::Build(nodeMap, coordDataView(), coords->getNumVectors());
1528 level.Set("Coordinates", newCoords);
1529 }
1530
1531 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1533 int N = Levels_.size();
1534 if( ( (sizeOfAllocatedLevelMultiVectors_ == numvecs && residual_.size() == N) || numvecs<=0 ) && !forceMapCheck) return;
1535
1536 // If, somehow, we changed the number of levels, delete everything first
1537 if(residual_.size() != N) {
1538 DeleteLevelMultiVectors();
1539
1540 residual_.resize(N);
1541 coarseRhs_.resize(N);
1542 coarseX_.resize(N);
1543 coarseImport_.resize(N);
1544 coarseExport_.resize(N);
1545 correction_.resize(N);
1546 }
1547
1548 for(int i=0; i<N; i++) {
1549 RCP<Operator> A = Levels_[i]->template Get< RCP<Operator> >("A");
1550 if(!A.is_null()) {
1551 // This dance is because we allow A to have a BlockedMap and X/B to have (compatible) non-blocked map
1552 RCP<const BlockedCrsMatrix> A_as_blocked = Teuchos::rcp_dynamic_cast<const BlockedCrsMatrix>(A);
1553 RCP<const Map> Arm = A->getRangeMap();
1554 RCP<const Map> Adm = A->getDomainMap();
1555 if(!A_as_blocked.is_null()) {
1556 Adm = A_as_blocked->getFullDomainMap();
1557 }
1558
1559 if (residual_[i].is_null() || !residual_[i]->getMap()->isSameAs(*Arm))
1560 // This is zero'd by default since it is filled via an operator apply
1561 residual_[i] = MultiVectorFactory::Build(Arm, numvecs, true);
1562 if (correction_[i].is_null() || !correction_[i]->getMap()->isSameAs(*Adm))
1563 correction_[i] = MultiVectorFactory::Build(Adm, numvecs, false);
1564 }
1565
1566 if(i+1<N) {
1567 // This is zero'd by default since it is filled via an operator apply
1568 if(implicitTranspose_) {
1569 RCP<Operator> P = Levels_[i+1]->template Get< RCP<Operator> >("P");
1570 if(!P.is_null()) {
1571 RCP<const Map> map = P->getDomainMap();
1572 if (coarseRhs_[i].is_null() || !coarseRhs_[i]->getMap()->isSameAs(*map))
1573 coarseRhs_[i] = MultiVectorFactory::Build(map, numvecs, true);
1574 }
1575 } else {
1576 RCP<Operator> R = Levels_[i+1]->template Get< RCP<Operator> >("R");
1577 if (!R.is_null()) {
1578 RCP<const Map> map = R->getRangeMap();
1579 if (coarseRhs_[i].is_null() || !coarseRhs_[i]->getMap()->isSameAs(*map))
1580 coarseRhs_[i] = MultiVectorFactory::Build(map, numvecs, true);
1581 }
1582 }
1583
1584
1585 RCP<const Import> importer;
1586 if(Levels_[i+1]->IsAvailable("Importer"))
1587 importer = Levels_[i+1]->template Get< RCP<const Import> >("Importer");
1588 if (doPRrebalance_ || importer.is_null()) {
1589 RCP<const Map> map = coarseRhs_[i]->getMap();
1590 if (coarseX_[i].is_null() || !coarseX_[i]->getMap()->isSameAs(*map))
1591 coarseX_[i] = MultiVectorFactory::Build(map, numvecs, true);
1592 } else {
1593 RCP<const Map> map;
1594 map = importer->getTargetMap();
1595 if (coarseImport_[i].is_null() || !coarseImport_[i]->getMap()->isSameAs(*map)) {
1596 coarseImport_[i] = MultiVectorFactory::Build(map, numvecs,false);
1597 coarseX_[i] = MultiVectorFactory::Build(map,numvecs,false);
1598 }
1599 map = importer->getSourceMap();
1600 if (coarseExport_[i].is_null() || !coarseExport_[i]->getMap()->isSameAs(*map))
1601 coarseExport_[i] = MultiVectorFactory::Build(map, numvecs,false);
1602 }
1603 }
1604 }
1605 sizeOfAllocatedLevelMultiVectors_ = numvecs;
1606 }
1607
1608
1609template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1611 if(sizeOfAllocatedLevelMultiVectors_==0) return;
1612 residual_.resize(0);
1613 coarseRhs_.resize(0);
1614 coarseX_.resize(0);
1615 coarseImport_.resize(0);
1616 coarseExport_.resize(0);
1617 correction_.resize(0);
1618 sizeOfAllocatedLevelMultiVectors_ = 0;
1619}
1620
1621
1622template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1624 const LO startLevel, const ConvData& conv) const
1625{
1626 return (startLevel == 0 && !isPreconditioner_ && (IsPrint(Statistics1) || conv.tol_ > 0));
1627}
1628
1629
1630template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1632 const Teuchos::Array<MagnitudeType>& residualNorm, const MagnitudeType convergenceTolerance) const
1633{
1635
1636 if (convergenceTolerance > Teuchos::ScalarTraits<MagnitudeType>::zero())
1637 {
1638 bool passed = true;
1639 for (LO k = 0; k < residualNorm.size(); k++)
1640 if (residualNorm[k] >= convergenceTolerance)
1641 passed = false;
1642
1643 if (passed)
1644 convergenceStatus = ConvergenceStatus::Converged;
1645 else
1646 convergenceStatus = ConvergenceStatus::Unconverged;
1647 }
1648
1649 return convergenceStatus;
1650}
1651
1652
1653template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1655 const LO iteration, const Teuchos::Array<MagnitudeType>& residualNorm) const
1656{
1657 GetOStream(Statistics1) << "iter: "
1658 << std::setiosflags(std::ios::left)
1659 << std::setprecision(3) << std::setw(4) << iteration
1660 << " residual = "
1661 << std::setprecision(10) << residualNorm
1662 << std::endl;
1663}
1664
1665template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1667 const Operator& A, const MultiVector& X, const MultiVector& B, const LO iteration,
1668 const LO startLevel, const ConvData& conv, MagnitudeType& previousResidualNorm)
1669{
1670 Teuchos::Array<MagnitudeType> residualNorm;
1671 residualNorm = Utilities::ResidualNorm(A, X, B, *residual_[startLevel]);
1672
1673 const MagnitudeType currentResidualNorm = residualNorm[0];
1674 rate_ = currentResidualNorm / previousResidualNorm;
1675 previousResidualNorm = currentResidualNorm;
1676
1677 if (IsPrint(Statistics1))
1678 PrintResidualHistory(iteration, residualNorm);
1679
1680 return IsConverged(residualNorm, conv.tol_);
1681}
1682
1683
1684} //namespace MueLu
1685
1686#endif // MUELU_HIERARCHY_DEF_HPP
virtual std::string description() const
Return a simple one-line description of this object.
Exception throws to report incompatible objects (like maps).
Exception throws to report errors in the internal logical of the program.
Base class for factories (e.g., R, P, and A_coarse).
Class that provides default factories within Needs class.
Provides methods to build a multigrid hierarchy and apply multigrid cycles.
void AddLevel(const RCP< Level > &level)
Add a level at the end of the hierarchy.
double GetSmootherComplexity() const
void Write(const LO &start=-1, const LO &end=-1, const std::string &suffix="")
Print matrices in the multigrid hierarchy to file.
RCP< Level > & GetLevel(const int levelID=0)
Retrieve a certain level from hierarchy.
void CheckLevel(Level &level, int levelID)
Helper function.
std::string description() const
Return a simple one-line description of this object.
void IsPreconditioner(const bool flag)
Array< RCP< Level > > Levels_
Container for Level objects.
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.
STS::magnitudeType MagnitudeType
void describe(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the Hierarchy with some verbosity level to a FancyOStream object.
ConvergenceStatus IsConverged(const Teuchos::Array< MagnitudeType > &residualNorm, const MagnitudeType convergenceTolerance) const
Decide if the multigrid iteration is converged.
ConvergenceStatus Iterate(const MultiVector &B, MultiVector &X, ConvData conv=ConvData(), bool InitialGuessIsZero=false, LO startLevel=0)
Apply the multigrid preconditioner.
void DumpCurrentGraph(int level) const
void SetMatvecParams(RCP< ParameterList > matvecParams)
Xpetra::UnderlyingLib lib_
Epetra/Tpetra mode.
void Clear(int startLevel=0)
Clear impermanent data from previous setup.
bool IsCalculationOfResidualRequired(const LO startLevel, const ConvData &conv) const
Decide if the residual needs to be computed.
ConvergenceStatus ComputeResidualAndPrintHistory(const Operator &A, const MultiVector &X, const MultiVector &B, const LO iteration, const LO startLevel, const ConvData &conv, MagnitudeType &previousResidualNorm)
Compute the residual norm and print it depending on the verbosity level.
double GetOperatorComplexity() const
void PrintResidualHistory(const LO iteration, const Teuchos::Array< MagnitudeType > &residualNorm) const
Print residualNorm for this iteration to the screen.
void AllocateLevelMultiVectors(int numvecs, bool forceMapCheck=false)
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.
void Delete(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Call Level::Delete(ename, factory) for each level of the Hierarchy.
int GetGlobalNumLevels() const
void Keep(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Call Level::Keep(ename, factory) for each level of the Hierarchy.
void AddKeepFlag(const std::string &ename, const FactoryBase *factory=NoFactory::get(), KeepType keep=MueLu::Keep)
Call Level::AddKeepFlag for each level of the Hierarchy.
void AddNewLevel()
Add a new level at the end of the hierarchy.
void RemoveKeepFlag(const std::string &ename, const FactoryBase *factory, KeepType keep=MueLu::All)
Call Level::RemoveKeepFlag for each level of the Hierarchy.
void ReplaceCoordinateMap(Level &level)
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need's value has been saved.
void SetComm(RCP< const Teuchos::Comm< int > > const &comm)
void setlib(Xpetra::UnderlyingLib lib2)
RCP< Level > & GetPreviousLevel()
Previous level.
void Release(const FactoryBase &factory)
Decrement the storage counter for all the inputs of a factory.
RCP< const Teuchos::Comm< int > > GetComm() const
int GetLevelID() const
Return level number.
Definition: MueLu_Level.cpp:76
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())
void Request(const FactoryBase &factory)
Increment the storage counter for all the inputs of a factory.
Xpetra::UnderlyingLib lib()
Timer to be used in non-factories.
static const NoFactory * get()
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
An exception safe way to call the method 'Level::SetFactoryManager()'.
Integrates Teuchos::TimeMonitor with MueLu verbosity system.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
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)
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)
static void SetRandomSeed(const Teuchos::Comm< int > &comm)
Namespace for MueLu classes and methods.
@ Warnings0
Important warning messages (one line)
@ Statistics2
Print even more statistics.
@ Warnings
Print all warning messages.
@ Errors
Errors.
@ Statistics1
Print more statistics.
@ Timings0
High level timing information (use Teuchos::TimeMonitor::summarize() to print)
@ Runtime0
One-liner description of what is happening.
@ Runtime1
Description of what is happening (more verbose)
@ Warnings1
Additional warnings.
@ Statistics0
Print statistics that do not involve significant additional computation.
@ Parameters1
Print class parameters (more parameters, more verbose)
short KeepType
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
VerbLevel toMueLuVerbLevel(const Teuchos::EVerbosityLevel verbLevel)
Translate Teuchos verbosity level to MueLu verbosity level.
static std::string getColonLabel(const std::string &label)
Helper function for object label.
Data struct for defining stopping criteria of multigrid iteration.