Skip to content

Commit

Permalink
Merge Pull Request #9793 from mayrmt/Trilinos/muelu-refactor-residual…
Browse files Browse the repository at this point in the history 
…-printout

Automatically Merged using Trilinos Pull Request AutoTester
PR Title: MueLu: refactor residual printout in Hierarchy::Iterate()
PR Author: mayrmt
  • Loading branch information
@trilinos-autotester
trilinos-autotester authored Oct 12, 2021
2 parents 73b5b34 + 222aeb7 commit 7644811
Show file tree
Hide file tree
Showing 4 changed files with 113 additions and 78 deletions.
26 changes: 24 additions & 2 deletions packages/muelu/src/MueCentral/MueLu_Hierarchy_decl.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -80,7 +80,7 @@

namespace MueLu {

enum ReturnType {
enum class ConvergenceStatus {
Converged,
Unconverged,
Undefined
Expand Down Expand Up @@ -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);

/*!
Expand Down Expand Up @@ -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<MagnitudeType>& residualNorm,
const MagnitudeType convergenceTolerance) const;

//! Print \c residualNorm for this \c iteration to the screen
void PrintResidualHistory(const LO iteration,
const Teuchos::Array<MagnitudeType>& 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<RCP<Level> > Levels_;

Expand Down
151 changes: 82 additions & 69 deletions packages/muelu/src/MueCentral/MueLu_Hierarchy_def.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -652,7 +652,7 @@ namespace MueLu {

#if defined(HAVE_MUELU_EXPERIMENTAL) && defined(HAVE_MUELU_ADDITIVE_VARIANT)
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
ReturnType Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Iterate(const MultiVector& B, MultiVector& X, ConvData conv,
ConvergenceStatus Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Iterate(const MultiVector& B, MultiVector& X, ConvData conv,
bool InitialGuessIsZero, LO startLevel) {
LO nIts = conv.maxIts_;
MagnitudeType tol = conv.tol_;
Expand Down Expand Up @@ -851,12 +851,12 @@ namespace MueLu {

//communicator->barrier();

return (tol > 0 ? Unconverged : Undefined);
return (tol > 0 ? ConvergenceStatus::Unconverged : ConvergenceStatus::Undefined);
}
#else
// ---------------------------------------- Iterate -------------------------------------------------------
template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
ReturnType Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Iterate(const MultiVector& B, MultiVector& X, ConvData conv,
ConvergenceStatus Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Iterate(const MultiVector& B, MultiVector& X, ConvData conv,
bool InitialGuessIsZero, LO startLevel) {
LO nIts = conv.maxIts_;
MagnitudeType tol = conv.tol_;
Expand Down Expand Up @@ -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.
Expand All @@ -918,36 +918,13 @@ namespace MueLu {

// Print residual information before iterating
typedef Teuchos::ScalarTraits<typename STS::magnitudeType> 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<MagnitudeType> 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<LO>::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) {
Expand Down Expand Up @@ -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<MagnitudeType> rn;
rn = Utilities::ResidualNorm(*A, X, B,*residual_[startLevel]);

prevNorm = curNorm;
curNorm = rn[0];
rate_ = as<MagnitudeType>(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

Expand Down Expand Up @@ -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<const FactoryBase*, BoostVertex> vindices;
typedef std::map<std::pair<BoostVertex,BoostVertex>, std::string> emap; emap edges;

static int call_id=0;

RCP<Operator> A = Levels_[0]->template Get<RCP<Operator> >("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<BoostEdge, bool> 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);
Expand All @@ -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();
Expand Down Expand Up @@ -1609,6 +1560,68 @@ void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::DeleteLevelMultiVecto
}


template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
bool Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::IsCalculationOfResidualRequired(
const LO startLevel, const ConvData& conv) const
{
return (startLevel == 0 && !isPreconditioner_ && (IsPrint(Statistics1) || conv.tol_ > 0));
}


template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
ConvergenceStatus Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::IsConverged(
const Teuchos::Array<MagnitudeType>& residualNorm, const MagnitudeType convergenceTolerance) const
{
ConvergenceStatus convergenceStatus = ConvergenceStatus::Undefined;

if (convergenceTolerance > Teuchos::ScalarTraits<MagnitudeType>::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 <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
void Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::PrintResidualHistory(
const LO iteration, const Teuchos::Array<MagnitudeType>& 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 <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
ConvergenceStatus Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::ComputeResidualAndPrintHistory(
const Operator& A, const MultiVector& X, const MultiVector& B, const LO iteration,
const LO startLevel, const ConvData& conv, MagnitudeType& previousResidualNorm)
{
Teuchos::Array<MagnitudeType> 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
10 changes: 5 additions & 5 deletions packages/muelu/test/convergence/Convergence.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<LO,MagnitudeType>(maxIts, tol));
MueLu::ConvergenceStatus ret = H->Iterate(*B, *X, std::pair<LO,MagnitudeType>(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
Expand Down
4 changes: 2 additions & 2 deletions packages/muelu/test/scaling/Driver.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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");
Expand All @@ -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
Expand Down

0 comments on commit 7644811

Please sign in to comment.