diff --git a/packages/muelu/src/MueCentral/MueLu_Hierarchy_decl.hpp b/packages/muelu/src/MueCentral/MueLu_Hierarchy_decl.hpp index 44e86eede79d..04831cf38d8e 100644 --- a/packages/muelu/src/MueCentral/MueLu_Hierarchy_decl.hpp +++ b/packages/muelu/src/MueCentral/MueLu_Hierarchy_decl.hpp @@ -80,7 +80,7 @@ namespace MueLu { - enum ReturnType { + enum class ConvergenceStatus { Converged, Unconverged, Undefined @@ -271,7 +271,7 @@ namespace MueLu { @param InitialGuessIsZero Indicates whether the initial guess is zero @param startLevel index of starting level to build multigrid hierarchy (default = 0) */ - ReturnType Iterate(const MultiVector& B, MultiVector& X, ConvData conv = ConvData(), + ConvergenceStatus Iterate(const MultiVector& B, MultiVector& X, ConvData conv = ConvData(), bool InitialGuessIsZero = false, LO startLevel = 0); /*! @@ -355,6 +355,28 @@ namespace MueLu { //! Copy constructor is not implemented. Hierarchy(const Hierarchy &h); + //! Decide if the residual needs to be computed + bool IsCalculationOfResidualRequired(const LO startLevel, const ConvData& conv) const; + + /*! + \brief Decide if the multigrid iteration is converged + + We judge convergence by comparing the current \c residualNorm + to the user given \c convergenceTolerance and then return the + appropriate \c ConvergenceStatus + */ + ConvergenceStatus IsConverged(const Teuchos::Array& residualNorm, + const MagnitudeType convergenceTolerance) const; + + //! Print \c residualNorm for this \c iteration to the screen + void PrintResidualHistory(const LO iteration, + const Teuchos::Array& residualNorm) const; + + //! Compute the residual norm and print it depending on the verbosity level + ConvergenceStatus ComputeResidualAndPrintHistory(const Operator& A, const MultiVector& X, + const MultiVector& B, const LO iteration, + const LO startLevel, const ConvData& conv, MagnitudeType& previousResidualNorm); + //! Container for Level objects Array > Levels_; diff --git a/packages/muelu/src/MueCentral/MueLu_Hierarchy_def.hpp b/packages/muelu/src/MueCentral/MueLu_Hierarchy_def.hpp index 95d43e03dd03..16061e10b7d2 100644 --- a/packages/muelu/src/MueCentral/MueLu_Hierarchy_def.hpp +++ b/packages/muelu/src/MueCentral/MueLu_Hierarchy_def.hpp @@ -652,7 +652,7 @@ namespace MueLu { #if defined(HAVE_MUELU_EXPERIMENTAL) && defined(HAVE_MUELU_ADDITIVE_VARIANT) template - ReturnType Hierarchy::Iterate(const MultiVector& B, MultiVector& X, ConvData conv, + ConvergenceStatus Hierarchy::Iterate(const MultiVector& B, MultiVector& X, ConvData conv, bool InitialGuessIsZero, LO startLevel) { LO nIts = conv.maxIts_; MagnitudeType tol = conv.tol_; @@ -851,12 +851,12 @@ namespace MueLu { //communicator->barrier(); - return (tol > 0 ? Unconverged : Undefined); + return (tol > 0 ? ConvergenceStatus::Unconverged : ConvergenceStatus::Undefined); } #else // ---------------------------------------- Iterate ------------------------------------------------------- template - ReturnType Hierarchy::Iterate(const MultiVector& B, MultiVector& X, ConvData conv, + ConvergenceStatus Hierarchy::Iterate(const MultiVector& B, MultiVector& X, ConvData conv, bool InitialGuessIsZero, LO startLevel) { LO nIts = conv.maxIts_; MagnitudeType tol = conv.tol_; @@ -902,7 +902,7 @@ namespace MueLu { // This processor does not have any data for this process on coarser // levels. This can only happen when there are multiple processors and // we use repartitioning. - return Undefined; + return ConvergenceStatus::Undefined; } // If we switched the number of vectors, we'd need to reallocate here. @@ -918,36 +918,13 @@ namespace MueLu { // Print residual information before iterating typedef Teuchos::ScalarTraits STM; - MagnitudeType prevNorm = STM::one(), curNorm = STM::one(); + MagnitudeType prevNorm = STM::one(); rate_ = 1.0; - if (startLevel == 0 && !isPreconditioner_ && - (IsPrint(Statistics1) || tol > 0)) { - // We calculate the residual only if we want to print it out, or if we - // want to stop once we achive the tolerance - Teuchos::Array rn; - rn = Utilities::ResidualNorm(*A, X, B,*residual_[startLevel]); - - if (tol > 0) { - bool passed = true; - for (LO k = 0; k < rn.size(); k++) - if (rn[k] >= tol) - passed = false; - - if (passed) - return Converged; - } - - if (IsPrint(Statistics1)) - GetOStream(Statistics1) << "iter: " - << std::setiosflags(std::ios::left) - << std::setprecision(3) << 0 // iter 0 - << " residual = " - << std::setprecision(10) << rn - << std::endl; - } + if (IsCalculationOfResidualRequired(startLevel, conv)) + ComputeResidualAndPrintHistory(*A, X, B, Teuchos::ScalarTraits::zero(), startLevel, conv, prevNorm); SC one = STS::one(), zero = STS::zero(); - for (LO i = 1; i <= nIts; i++) { + for (LO iteration = 1; iteration <= nIts; iteration++) { #ifdef HAVE_MUELU_DEBUG #if 0 // TODO fix me if (A->getDomainMap()->isCompatible(*(X.getMap())) == false) { @@ -1136,37 +1113,11 @@ namespace MueLu { } zeroGuess = false; - if (startLevel == 0 && !isPreconditioner_ && - (IsPrint(Statistics1) || tol > 0)) { - // We calculate the residual only if we want to print it out, or if we - // want to stop once we achive the tolerance - Teuchos::Array rn; - rn = Utilities::ResidualNorm(*A, X, B,*residual_[startLevel]); - - prevNorm = curNorm; - curNorm = rn[0]; - rate_ = as(curNorm / prevNorm); - - if (IsPrint(Statistics1)) - GetOStream(Statistics1) << "iter: " - << std::setiosflags(std::ios::left) - << std::setprecision(3) << i - << " residual = " - << std::setprecision(10) << rn - << std::endl; - - if (tol > 0) { - bool passed = true; - for (LO k = 0; k < rn.size(); k++) - if (rn[k] >= tol) - passed = false; - - if (passed) - return Converged; - } - } + + if (IsCalculationOfResidualRequired(startLevel, conv)) + ComputeResidualAndPrintHistory(*A, X, B, iteration, startLevel, conv, prevNorm); } - return (tol > 0 ? Unconverged : Undefined); + return (tol > 0 ? ConvergenceStatus::Unconverged : ConvergenceStatus::Undefined); } #endif @@ -1399,29 +1350,29 @@ namespace MueLu { if (GetProcRankVerbose() != 0) return; #if defined(HAVE_MUELU_BOOST) && defined(HAVE_MUELU_BOOST_FOR_REAL) && defined(BOOST_VERSION) && (BOOST_VERSION >= 104400) - + BoostGraph graph; - + BoostProperties dp; dp.property("label", boost::get(boost::vertex_name, graph)); dp.property("id", boost::get(boost::vertex_index, graph)); dp.property("label", boost::get(boost::edge_name, graph)); dp.property("color", boost::get(boost::edge_color, graph)); - + // create local maps std::map vindices; typedef std::map, std::string> emap; emap edges; - + static int call_id=0; - + RCP A = Levels_[0]->template Get >("A"); int rank = A->getDomainMap()->getComm()->getRank(); - + // printf("[%d] CMS: ----------------------\n",rank); for (int i = currLevel; i <= currLevel+1 && i < GetNumLevels(); i++) { edges.clear(); Levels_[i]->UpdateGraph(vindices, edges, dp, graph); - + for (emap::const_iterator eit = edges.begin(); eit != edges.end(); eit++) { std::pair boost_edge = boost::add_edge(eit->first.first, eit->first.second, graph); // printf("[%d] CMS: Hierarchy, adding edge (%d->%d) %d\n",rank,(int)eit->first.first,(int)eit->first.second,(int)boost_edge.second); @@ -1434,7 +1385,7 @@ namespace MueLu { boost::put("color", dp, boost_edge.first, std::string("blue")); } } - + 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")); boost::write_graphviz_dp(out, graph, dp, std::string("id")); out.close(); @@ -1609,6 +1560,68 @@ void Hierarchy::DeleteLevelMultiVecto } +template +bool Hierarchy::IsCalculationOfResidualRequired( + const LO startLevel, const ConvData& conv) const +{ + return (startLevel == 0 && !isPreconditioner_ && (IsPrint(Statistics1) || conv.tol_ > 0)); +} + + +template +ConvergenceStatus Hierarchy::IsConverged( + const Teuchos::Array& residualNorm, const MagnitudeType convergenceTolerance) const +{ + ConvergenceStatus convergenceStatus = ConvergenceStatus::Undefined; + + if (convergenceTolerance > Teuchos::ScalarTraits::zero()) + { + bool passed = true; + for (LO k = 0; k < residualNorm.size(); k++) + if (residualNorm[k] >= convergenceTolerance) + passed = false; + + if (passed) + convergenceStatus = ConvergenceStatus::Converged; + else + convergenceStatus = ConvergenceStatus::Unconverged; + } + + return convergenceStatus; +} + + +template +void Hierarchy::PrintResidualHistory( + const LO iteration, const Teuchos::Array& residualNorm) const +{ + GetOStream(Statistics1) << "iter: " + << std::setiosflags(std::ios::left) + << std::setprecision(3) << std::setw(4) << iteration + << " residual = " + << std::setprecision(10) << residualNorm + << std::endl; +} + +template +ConvergenceStatus Hierarchy::ComputeResidualAndPrintHistory( + const Operator& A, const MultiVector& X, const MultiVector& B, const LO iteration, + const LO startLevel, const ConvData& conv, MagnitudeType& previousResidualNorm) +{ + Teuchos::Array residualNorm; + residualNorm = Utilities::ResidualNorm(A, X, B, *residual_[startLevel]); + + const MagnitudeType currentResidualNorm = residualNorm[0]; + rate_ = currentResidualNorm / previousResidualNorm; + previousResidualNorm = currentResidualNorm; + + if (IsPrint(Statistics1)) + PrintResidualHistory(iteration, residualNorm); + + return IsConverged(residualNorm, conv.tol_); +} + + } //namespace MueLu #endif // MUELU_HIERARCHY_DEF_HPP diff --git a/packages/muelu/test/convergence/Convergence.cpp b/packages/muelu/test/convergence/Convergence.cpp index d38a83656e26..954035d908a4 100644 --- a/packages/muelu/test/convergence/Convergence.cpp +++ b/packages/muelu/test/convergence/Convergence.cpp @@ -341,21 +341,21 @@ int main_(Teuchos::CommandLineProcessor &clp, Xpetra::UnderlyingLib lib, int arg H->IsPreconditioner(isPrec); if (isPrec == false) { - MueLu::ReturnType ret = H->Iterate(*B, *X, std::pair(maxIts, tol)); + MueLu::ConvergenceStatus ret = H->Iterate(*B, *X, std::pair(maxIts, tol)); double rate = H->GetRate(); if (std::abs(rate-goldRate) < 0.02) { out << xmlFile << ": passed (" << - (ret == MueLu::Converged ? "converged, " : "unconverged, ") << + (ret == MueLu::ConvergenceStatus::Converged ? "converged, " : "unconverged, ") << "expected rate = " << goldRate << ", real rate = " << rate << - (ret == MueLu::Converged ? "" : " (after " + Teuchos::toString(maxIts) + " iterations)") + (ret == MueLu::ConvergenceStatus::Converged ? "" : " (after " + Teuchos::toString(maxIts) + " iterations)") << ")" << std::endl; } else { out << xmlFile << ": failed (" << - (ret == MueLu::Converged ? "converged, " : "unconverged, ") << + (ret == MueLu::ConvergenceStatus::Converged ? "converged, " : "unconverged, ") << "expected rate = " << goldRate << ", real rate = " << rate << - (ret == MueLu::Converged ? "" : " (after " + Teuchos::toString(maxIts) + " iterations)") + (ret == MueLu::ConvergenceStatus::Converged ? "" : " (after " + Teuchos::toString(maxIts) + " iterations)") << ")" << std::endl; // ap: we need to understand what's going on with the convergence rate // At the moment, disable failure state, so that the test passes as long diff --git a/packages/muelu/test/scaling/Driver.cpp b/packages/muelu/test/scaling/Driver.cpp index e84784bd11fb..678c414148c6 100644 --- a/packages/muelu/test/scaling/Driver.cpp +++ b/packages/muelu/test/scaling/Driver.cpp @@ -219,7 +219,7 @@ int main_(Teuchos::CommandLineProcessor &clp, Xpetra::UnderlyingLib& lib, int ar bool printTimings = true; clp.setOption("timings", "notimings", &printTimings, "print timings to screen"); std::string timingsFormat = "table-fixed"; clp.setOption("time-format", &timingsFormat, "timings format (table-fixed | table-scientific | yaml)"); int writeMatricesOPT = -2; clp.setOption("write", &writeMatricesOPT, "write matrices to file (-1 means all; i>=0 means level i)"); - std::string dsolveType = "belos", solveType; clp.setOption("solver", &dsolveType, "solve type: (none | cg | gmres | standalone | matvec)"); + std::string dsolveType = "belos", solveType; clp.setOption("solver", &dsolveType, "solve type: (none | belos | standalone | matvec)"); std::string belosType = "cg"; clp.setOption("belosType", &belosType, "belos solver type: (Pseudoblock CG | Block CG | Pseudoblock GMRES | Block GMRES | ...) see BelosSolverFactory.hpp for exhaustive list of solvers"); double dtol = 1e-12, tol; clp.setOption("tol", &dtol, "solver convergence tolerance"); bool binaryFormat = false; clp.setOption("binary", "ascii", &binaryFormat, "print timings to screen"); @@ -239,7 +239,7 @@ int main_(Teuchos::CommandLineProcessor &clp, Xpetra::UnderlyingLib& lib, int ar int maxIts = 200; clp.setOption("its", &maxIts, "maximum number of solver iterations"); int numVectors = 1; clp.setOption("multivector", &numVectors, "number of rhs to solve simultaneously"); bool scaleResidualHist = true; clp.setOption("scale", "noscale", &scaleResidualHist, "scaled Krylov residual history"); - bool solvePreconditioned = true; clp.setOption("solve-preconditioned","no-solve-preconditioned", &solvePreconditioned, "use MueLu preconditioner in solve"); + bool solvePreconditioned = true; clp.setOption("solve-preconditioned","no-solve-preconditioned", &solvePreconditioned, "use MueLu preconditioner in solve"); bool useStackedTimer = false; clp.setOption("stacked-timer","no-stacked-timer", &useStackedTimer, "use stacked timer"); #ifdef HAVE_MUELU_TPETRA