Skip to content

Commit

Permalink
Merge Pull Request #9627 from cgcgcg/Trilinos/mlDefaults
Browse files Browse the repository at this point in the history 
Automatically Merged using Trilinos Pull Request AutoTester
PR Title: Ifpack2 Chebyshev: Add CG option for computation of spectral radius
PR Author: cgcgcg
  • Loading branch information
@trilinos-autotester
trilinos-autotester authored Sep 1, 2021
2 parents 049b079 + c9148be commit a24500e
Show file tree
Hide file tree
Showing 5 changed files with 228 additions and 16 deletions.
30 changes: 30 additions & 0 deletions packages/ifpack2/src/Ifpack2_Details_Chebyshev_decl.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -482,6 +482,9 @@ class Chebyshev : public Teuchos::Describable {
//! Whether the iteration vectors of the power method should be saved.
bool eigKeepVectors_;

//! Type of eigen-analysis, i.e. power method or cg.
std::string eigenAnalysisType_;

//! Iteration vectors of the power method. Can be saved for the purpose of multiple setups.
Teuchos::RCP<V> eigVector_;
Teuchos::RCP<V> eigVector2_;
Expand Down Expand Up @@ -693,6 +696,33 @@ class Chebyshev : public Teuchos::Describable {
ST
powerMethod (const op_type& A, const V& D_inv, const int numIters);

/// \brief Use the cg method to estimate the maximum eigenvalue
/// of A*D_inv, given an initial guess vector x.
///
/// \param A [in] The Operator to use.
/// \param D_inv [in] Vector to use as implicit right scaling of A.
/// \param numIters [in] Maximum number of iterations of the power
/// method.
/// \param x [in/out] On input: Initial guess Vector for the cg
/// method. Its Map must be the same as that of the domain Map of
/// A. This method may use this Vector as scratch space.
///
/// \return Estimate of the maximum eigenvalue of A*D_inv.
ST
cgMethodWithInitGuess (const op_type& A, const V& D_inv, const int numIters, V& x);

/// \brief Use the cg method to estimate the maximum eigenvalue
/// of A*D_inv.
///
/// \param A [in] The Operator to use.
/// \param D_inv [in] Vector to use as implicit right scaling of A.
/// \param numIters [in] Maximum number of iterations of the power
/// method.
///
/// \return Estimate of the maximum eigenvalue of A*D_inv.
ST
cgMethod (const op_type& A, const V& D_inv, const int numIters);

//! The maximum infinity norm of all the columns of X.
static MT maxNormInf (const MV& X);

Expand Down
161 changes: 154 additions & 7 deletions packages/ifpack2/src/Ifpack2_Details_Chebyshev_def.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -58,6 +58,8 @@
#include "Teuchos_FancyOStream.hpp"
#include "Teuchos_oblackholestream.hpp"
#include "Tpetra_Details_residual.hpp"
#include "Teuchos_LAPACK.hpp"
#include "Ifpack2_Details_LapackSupportsScalar.hpp"
#include <cmath>
#include <iostream>

Expand Down Expand Up @@ -218,6 +220,51 @@ class PositivizeVector {

} // namespace (anonymous)


template<class ScalarType, const bool lapackSupportsScalarType = LapackSupportsScalar<ScalarType>::value>
struct LapackHelper{
static
ScalarType
tri_diag_spectral_radius(Teuchos::ArrayRCP<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> diag,
Teuchos::ArrayRCP<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> offdiag) {
throw std::runtime_error("LAPACK does not support the scalar type.");
}
};


template<class ScalarType>
struct LapackHelper<ScalarType,true> {
static
ScalarType
tri_diag_spectral_radius(Teuchos::ArrayRCP<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> diag,
Teuchos::ArrayRCP<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> offdiag) {
using STS = Teuchos::ScalarTraits<ScalarType>;
using MagnitudeType = typename STS::magnitudeType;
int info = 0;
const int N = diag.size ();
ScalarType scalar_dummy;
std::vector<MagnitudeType> mag_dummy(4*N);
char char_N = 'N';

// lambdaMin = one;
ScalarType lambdaMax = STS::one();
if( N > 2 ) {
Teuchos::LAPACK<int,ScalarType> lapack;
lapack.PTEQR (char_N, N, diag.getRawPtr (), offdiag.getRawPtr (),
&scalar_dummy, 1, &mag_dummy[0], &info);
TEUCHOS_TEST_FOR_EXCEPTION
(info < 0, std::logic_error, "Ifpack2::Details::LapackHelper::tri_diag_spectral_radius:"
"LAPACK's _PTEQR failed with info = "
<< info << " < 0. This suggests there might be a bug in the way Ifpack2 "
"is calling LAPACK. Please report this to the Ifpack2 developers.");
// lambdaMin = Teuchos::as<ScalarType> (diag[N-1]);
lambdaMax = Teuchos::as<ScalarType> (diag[0]);
}
return lambdaMax;
}
};


template<class ScalarType, class MV>
void Chebyshev<ScalarType, MV>::checkInputMatrix () const
{
Expand Down Expand Up @@ -288,6 +335,7 @@ Chebyshev (Teuchos::RCP<const row_matrix_type> A) :
eigMaxIters_ (10),
eigRelTolerance_(Teuchos::ScalarTraits<MT>::zero ()),
eigKeepVectors_(false),
eigenAnalysisType_("power method"),
eigNormalizationFreq_(1),
zeroStartingSolution_ (true),
assumeMatrixUnchanged_ (false),
Expand Down Expand Up @@ -318,6 +366,7 @@ Chebyshev (Teuchos::RCP<const row_matrix_type> A,
eigMaxIters_ (10),
eigRelTolerance_(Teuchos::ScalarTraits<MT>::zero ()),
eigKeepVectors_(false),
eigenAnalysisType_("power method"),
eigNormalizationFreq_(1),
zeroStartingSolution_ (true),
assumeMatrixUnchanged_ (false),
Expand Down Expand Up @@ -649,13 +698,10 @@ setParameters (Teuchos::ParameterList& plist)
eigenAnalysisType = plist.get<std::string> ("eigen-analysis: type");
TEUCHOS_TEST_FOR_EXCEPTION(
eigenAnalysisType != "power-method" &&
eigenAnalysisType != "power method",
eigenAnalysisType != "power method" &&
eigenAnalysisType != "cg",
std::invalid_argument,
"Ifpack2::Chebyshev: This class supports the ML parameter \"eigen-"
"analysis: type\" for backwards compatibility. However, Ifpack2 "
"currently only supports the \"power-method\" option for this "
"parameter. This imitates Ifpack, which only implements the power "
"method for eigenanalysis.");
"Ifpack2::Chebyshev: Ifpack2 only supports \"power method\" and \"cg\" for \"eigen-analysis: type\".");
}

// We've validated all the parameters, so it's safe now to "commit" them.
Expand All @@ -670,6 +716,7 @@ setParameters (Teuchos::ParameterList& plist)
eigRelTolerance_ = eigRelTolerance;
eigKeepVectors_ = eigKeepVectors;
eigNormalizationFreq_ = eigNormalizationFreq;
eigenAnalysisType_ = eigenAnalysisType;
zeroStartingSolution_ = zeroStartingSolution;
assumeMatrixUnchanged_ = assumeMatrixUnchanged;
textbookAlgorithm_ = textbookAlgorithm;
Expand Down Expand Up @@ -852,7 +899,11 @@ Chebyshev<ScalarType, MV>::compute ()
// most important one if using Chebyshev as a smoother.
if (! assumeMatrixUnchanged_ ||
(! computedEigenvalueEstimates && STS::isnaninf (userLambdaMax_))) {
const ST computedLambdaMax = powerMethod (*A_, *D_, eigMaxIters_);
ST computedLambdaMax;
if ((eigenAnalysisType_ == "power method") || (eigenAnalysisType_ == "power-method"))
computedLambdaMax = powerMethod (*A_, *D_, eigMaxIters_);
else
computedLambdaMax = cgMethod (*A_, *D_, eigMaxIters_);
TEUCHOS_TEST_FOR_EXCEPTION(
STS::isnaninf (computedLambdaMax),
std::runtime_error,
Expand Down Expand Up @@ -1548,6 +1599,102 @@ powerMethod (const op_type& A, const V& D_inv, const int numIters)
return lambdaMax;
}


template<class ScalarType, class MV>
typename Chebyshev<ScalarType, MV>::ST
Chebyshev<ScalarType, MV>::
cgMethodWithInitGuess (const op_type& A,
const V& D_inv,
const int numIters,
V& r)
{
using std::endl;
using STS = Teuchos::ScalarTraits<ST>;
using MagnitudeType = typename STS::magnitudeType;
if (debug_) {
*out_ << " cgMethodWithInitGuess:" << endl;
}

const ST one = STS::one();
ST beta, betaOld = one, pAp, pApOld = one, alpha, rHz, rHzOld, rHzOld2 = one, lambdaMax;
// ST lambdaMin;
Teuchos::ArrayRCP<MagnitudeType> diag, offdiag;
Teuchos::RCP<V> p, z, Ap;
diag.resize(numIters);
offdiag.resize(numIters-1);

p = rcp(new V(A.getRangeMap ()));
z = rcp(new V(A.getRangeMap ()));
Ap = rcp(new V(A.getRangeMap ()));

// Tpetra::Details::residual (A, x, *b, *r);
solve (*p, D_inv, r);
rHz = r.dot (*p);

for (int iter = 0; iter < numIters; ++iter) {
if (debug_) {
*out_ << " Iteration " << iter << endl;
}
A.apply (*p, *Ap);
pAp = p->dot (*Ap);
alpha = rHz/pAp;
r.update (-alpha, *Ap, one);
rHzOld = rHz;
solve (*z, D_inv, r);
rHz = r.dot (*z);
beta = rHz / rHzOld;
p->update (one, *z, beta);
if (iter>0) {
diag[iter] = STS::real((betaOld*betaOld * pApOld + pAp) / rHzOld);
offdiag[iter-1] = -STS::real((betaOld * pApOld / (sqrt(rHzOld * rHzOld2))));
if (debug_) {
*out_ << " diag[" << iter << "] = " << diag[iter] << endl;
*out_ << " offdiag["<< iter-1 << "] = " << offdiag[iter-1] << endl;
}
} else {
diag[iter] = STS::real(pAp/rHzOld);
if (debug_) {
*out_ << " diag[" << iter << "] = " << diag[iter] << endl;
}
}
rHzOld2 = rHzOld;
betaOld = beta;
pApOld = pAp;
}

lambdaMax = LapackHelper<ST>::tri_diag_spectral_radius(diag, offdiag);

return lambdaMax;
}


template<class ScalarType, class MV>
typename Chebyshev<ScalarType, MV>::ST
Chebyshev<ScalarType, MV>::
cgMethod (const op_type& A, const V& D_inv, const int numIters)
{
using std::endl;
if (debug_) {
*out_ << "cgMethod:" << endl;
}

Teuchos::RCP<V> r;
if (eigVector_.is_null()) {
r = rcp(new V(A.getDomainMap ()));
if (eigKeepVectors_)
eigVector_ = r;
// For the first pass, just let the pseudorandom number generator
// fill x with whatever values it wants; don't try to make its
// entries nonnegative.
computeInitialGuessForPowerMethod (*r, false);
} else
r = eigVector_;

ST lambdaMax = cgMethodWithInitGuess (A, D_inv, numIters, *r);

return lambdaMax;
}

template<class ScalarType, class MV>
Teuchos::RCP<const typename Chebyshev<ScalarType, MV>::row_matrix_type>
Chebyshev<ScalarType, MV>::getMatrix () const {
Expand Down
24 changes: 24 additions & 0 deletions packages/ifpack2/test/unit_tests/Ifpack2_UnitTestChebyshev.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -151,6 +151,30 @@ TEUCHOS_UNIT_TEST_TEMPLATE_3_DECL(Ifpack2Chebyshev, Test0, Scalar, LocalOrdinal,
Teuchos::ArrayRCP<const Scalar> yview = y.get1dView();
TEST_COMPARE_FLOATING_ARRAYS(yview, halfs(), tol);
}

crsmatrix = tif_utest::create_test_matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>(rowmap, -one);
Scalar n = Teuchos::as<Scalar>(rowmap->getGlobalNumElements());
Scalar expectedLambdaMax = one-std::cos(Teuchos::ScalarTraits<Scalar>::pi()*n/(n+1));

prec.setMatrix(crsmatrix);

params.set("debug", true);
params.remove("chebyshev: max eigenvalue");
params.remove("chebyshev: min eigenvalue");

params.set("eigen-analysis: type", "power method");
params.set("chebyshev: eigenvalue max iterations",30);
prec.setParameters(params);
prec.compute();

TEST_FLOATING_EQUALITY(prec.getLambdaMaxForApply(),expectedLambdaMax,2e-2);

params.set("eigen-analysis: type", "cg");
params.set("chebyshev: eigenvalue max iterations",10);
prec.setParameters(params);
prec.compute();

TEST_FLOATING_EQUALITY(prec.getLambdaMaxForApply(),expectedLambdaMax,2e-2);
}

#define UNIT_TEST_GROUP_SC_LO_GO(Scalar,LocalOrdinal,GlobalOrdinal) \
Expand Down
25 changes: 16 additions & 9 deletions packages/muelu/src/Interface/MueLu_ML2MueLuParameterTranslator.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,9 +46,12 @@

#include "MueLu_ConfigDefs.hpp"
#if defined(HAVE_MUELU_ML)
#include <ml_ValidateParameters.h>
#include <ml_MultiLevelPreconditioner.h> // for default values
#include <ml_RefMaxwell.h>
# include <ml_config.h>
# if defined(HAVE_ML_EPETRA) && defined(HAVE_ML_TEUCHOS)
# include <ml_ValidateParameters.h>
# include <ml_MultiLevelPreconditioner.h> // for default values
# include <ml_RefMaxwell.h>
# endif
#endif

#include <MueLu_ML2MueLuParameterTranslator.hpp>
Expand Down Expand Up @@ -161,6 +164,8 @@ namespace MueLu {
else { mueluss << "<Parameter name=\"chebyshev: degree\" type=\"int\" value=\"2\"/>" << std::endl; }
if ( paramList.isParameter("smoother: Chebyshev alpha") ) { mueluss << "<Parameter name=\"chebyshev: ratio eigenvalue\" type=\"double\" value=\"" << paramList.get<double>("smoother: Chebyshev alpha") << "\"/>" << std::endl; adaptingParamList.remove("smoother: Chebyshev alpha",false); }
else { mueluss << "<Parameter name=\"chebyshev: ratio eigenvalue\" type=\"double\" value=\"20\"/>" << std::endl; adaptingParamList.remove("smoother: Chebyshev alpha",false); }
if ( paramList.isParameter("eigen-analysis: type") ) { mueluss << "<Parameter name=\"eigen-analysis: type\" type=\"string\" value=\"" << paramList.get<std::string>("eigen-analysis: type") << "\"/>" << std::endl; adaptingParamList.remove("eigen-analysis: type",false); }
else { mueluss << "<Parameter name=\"eigen-analysis: type\" type=\"string\" value=\"cg\"/>" << std::endl; }
}

// MLS
Expand All @@ -171,6 +176,8 @@ namespace MueLu {
if ( paramList.isParameter("smoother: MLS alpha") ) { mueluss << "<Parameter name=\"chebyshev: ratio eigenvalue\" type=\"double\" value=\"" << paramList.get<double>("smoother: MLS alpha") << "\"/>" << std::endl; adaptingParamList.remove("smoother: MLS alpha",false); }
else if ( paramList.isParameter("smoother: Chebyshev alpha") ) { mueluss << "<Parameter name=\"chebyshev: ratio eigenvalue\" type=\"double\" value=\"" << paramList.get<double>("smoother: Chebyshev alpha") << "\"/>" << std::endl; adaptingParamList.remove("smoother: Chebyshev alpha",false); }
else { mueluss << "<Parameter name=\"chebyshev: ratio eigenvalue\" type=\"double\" value=\"20\"/>" << std::endl; }
if ( paramList.isParameter("eigen-analysis: type") ) { mueluss << "<Parameter name=\"eigen-analysis: type\" type=\"string\" value=\"" << paramList.get<std::string>("eigen-analysis: type") << "\"/>" << std::endl; adaptingParamList.remove("eigen-analysis: type",false); }
else { mueluss << "<Parameter name=\"eigen-analysis: type\" type=\"string\" value=\"cg\"/>" << std::endl; }
}

// parameters for ILU based preconditioners
Expand Down Expand Up @@ -198,7 +205,7 @@ namespace MueLu {

RCP<Teuchos::FancyOStream> out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout)); // TODO: use internal out (GetOStream())

#if defined(HAVE_MUELU_ML)
#if defined(HAVE_MUELU_ML) && defined(HAVE_ML_EPETRA) && defined(HAVE_ML_TEUCHOS)

// TODO alternative with standard parameterlist from ML user guide?

Expand All @@ -218,9 +225,9 @@ namespace MueLu {
#else
if (defaultVals != "") {
// If no validator available: issue a warning and set parameter value to false in the output list
*out << "Warning: MueLu_ENABLE_ML=OFF and/or MueLu_ENABLE_Epetra=OFF. No ML default values available." << std::endl;
*out << "Warning: MueLu_ENABLE_ML=OFF, ML_ENABLE_Epetra=OFF or ML_ENABLE_TEUCHOS=OFF. No ML default values available." << std::endl;
}
#endif // HAVE_MUELU_ML
#endif // HAVE_MUELU_ML && HAVE_ML_EPETRA && HAVE_ML_TEUCHOS

//
// Move smoothers/aggregation/coarse parameters to sublists
Expand All @@ -240,17 +247,17 @@ namespace MueLu {
bool validate = paramList.get("ML validate parameter list", true); /* true = default in ML */
if (validate && defaultVals!="refmaxwell") {

#if defined(HAVE_MUELU_ML)
#if defined(HAVE_MUELU_ML) && defined(HAVE_ML_EPETRA) && defined(HAVE_ML_TEUCHOS)
// Validate parameter list using ML validator
int depth = paramList.get("ML validate depth", 5); /* 5 = default in ML */
TEUCHOS_TEST_FOR_EXCEPTION(! ML_Epetra::ValidateMLPParameters(paramList, depth), Exceptions::RuntimeError,
"ERROR: ML's Teuchos::ParameterList contains incorrect parameter!");
#else
// If no validator available: issue a warning and set parameter value to false in the output list
*out << "Warning: MueLu_ENABLE_ML=OFF and/or MueLu_ENABLE_Epetra=OFF. The parameter list cannot be validated." << std::endl;
*out << "Warning: MueLu_ENABLE_ML=OFF, ML_ENABLE_Epetra=OFF or ML_ENABLE_TEUCHOS=OFF. The parameter list cannot be validated." << std::endl;
paramList.set("ML validate parameter list", false);

#endif // HAVE_MUELU_ML
#endif // HAVE_MUELU_ML && HAVE_ML_EPETRA && HAVE_ML_TEUCHOS
} // if(validate)
} // scope

Expand Down
Loading

0 comments on commit a24500e

Please sign in to comment.