Skip to content

Commit

Permalink
Merge Pull Request #9294 from cgcgcg/Trilinos/halfPrecisionMueLu
Browse files Browse the repository at this point in the history 
Automatically Merged using Trilinos Pull Request AutoTester
PR Title: MueLu: half precision preconditioner
PR Author: cgcgcg
  • Loading branch information
@trilinos-autotester
trilinos-autotester authored Jun 17, 2021
2 parents f9c7dfd + 445624b commit 1a08fde
Show file tree
Hide file tree
Showing 8 changed files with 336 additions and 643 deletions.
396 changes: 9 additions & 387 deletions packages/muelu/adapters/stratimikos/Thyra_MueLuPreconditionerFactory_decl.hpp
View file

Large diffs are not rendered by default.

390 changes: 307 additions & 83 deletions packages/muelu/adapters/stratimikos/Thyra_MueLuPreconditionerFactory_def.hpp
View file

Large diffs are not rendered by default.

6 changes: 2 additions & 4 deletions packages/muelu/adapters/stratimikos/Thyra_MueLuRefMaxwellPreconditionerFactory_decl.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -56,6 +56,7 @@
#include "Thyra_BlockedLinearOpBase.hpp"
#include "Thyra_DiagonalLinearOpBase.hpp"
#include "Thyra_XpetraLinearOp.hpp"
#include <Thyra_MueLuPreconditionerFactory.hpp>
#ifdef HAVE_MUELU_TPETRA
#include "Thyra_TpetraLinearOp.hpp"
#include "Thyra_TpetraThyraWrappers.hpp"
Expand Down Expand Up @@ -95,9 +96,6 @@

namespace Thyra {

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
bool replaceWithXpetra(ParameterList& paramList, std::string parameterName);

/** @brief Concrete preconditioner factory subclass for Thyra based on MueLu.
@ingroup MueLuAdapters
Add support for MueLu preconditioners in Thyra. This class provides an interface both
Expand Down Expand Up @@ -283,7 +281,7 @@ namespace Thyra {
// rebuild preconditioner if startingOver == true
// reuse preconditioner if startingOver == false
const bool startingOver = (thyra_precOp.is_null() || !paramList.isParameter("refmaxwell: enable reuse") || !paramList.get<bool>("refmaxwell: enable reuse"));
const bool useHalfPrecision = paramList.get<bool>("refmaxwell: half precision", false) && bIsTpetra;
const bool useHalfPrecision = paramList.get<bool>("half precision", false) && bIsTpetra;

RCP<XpOp> xpPrecOp;
if (startingOver == true) {
Expand Down
169 changes: 1 addition & 168 deletions packages/muelu/adapters/stratimikos/Thyra_MueLuRefMaxwellPreconditionerFactory_def.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -59,173 +59,6 @@ namespace Thyra {
using Teuchos::rcp_dynamic_cast;
using Teuchos::rcp_const_cast;

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
bool replaceWithXpetra(ParameterList& paramList, std::string parameterName) {
typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
typedef Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpOp;
typedef Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpThyUtils;
typedef Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMatWrap;
typedef Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMat;
typedef Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMat;
typedef Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMultVec;
typedef Xpetra::MultiVector<Magnitude,LocalOrdinal,GlobalOrdinal,Node> XpMagMultVec;
typedef Xpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpVec;

typedef Thyra::LinearOpBase<Scalar> ThyLinOpBase;
typedef Thyra::DiagonalLinearOpBase<Scalar> ThyDiagLinOpBase;
typedef Thyra::XpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node> ThyXpOp;
typedef Thyra::SpmdVectorSpaceBase<Scalar> ThyVSBase;

#ifdef HAVE_MUELU_TPETRA
typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> TpCrsMat;
typedef Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tV;
typedef Thyra::TpetraVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> thyTpV;
typedef Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tMV;
typedef Tpetra::MultiVector<Magnitude, LocalOrdinal, GlobalOrdinal, Node> tMagMV;
# if defined(HAVE_TPETRA_INST_DOUBLE) && defined(HAVE_TPETRA_INST_FLOAT)
typedef typename Teuchos::ScalarTraits<Magnitude>::halfPrecision HalfMagnitude;
typedef Tpetra::MultiVector<HalfMagnitude, LocalOrdinal, GlobalOrdinal, Node> tHalfMagMV;
# endif
#endif
#if defined(HAVE_MUELU_EPETRA)
typedef Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> XpEpCrsMat;
#endif

if (paramList.isParameter(parameterName)) {
if (paramList.isType<RCP<XpMat> >(parameterName))
return true;
else if (paramList.isType<RCP<const XpMat> >(parameterName)) {
RCP<const XpMat> constM = paramList.get<RCP<const XpMat> >(parameterName);
paramList.remove(parameterName);
RCP<XpMat> M = rcp_const_cast<XpMat>(constM);
paramList.set<RCP<XpMat> >(parameterName, M);
return true;
}
else if (paramList.isType<RCP<XpMultVec> >(parameterName))
return true;
else if (paramList.isType<RCP<const XpMultVec> >(parameterName)) {
RCP<const XpMultVec> constX = paramList.get<RCP<const XpMultVec> >(parameterName);
paramList.remove(parameterName);
RCP<XpMultVec> X = rcp_const_cast<XpMultVec>(constX);
paramList.set<RCP<XpMultVec> >(parameterName, X);
return true;
}
else if (paramList.isType<RCP<XpMagMultVec> >(parameterName))
return true;
else if (paramList.isType<RCP<const XpMagMultVec> >(parameterName)) {
RCP<const XpMagMultVec> constX = paramList.get<RCP<const XpMagMultVec> >(parameterName);
paramList.remove(parameterName);
RCP<XpMagMultVec> X = rcp_const_cast<XpMagMultVec>(constX);
paramList.set<RCP<XpMagMultVec> >(parameterName, X);
return true;
}
#ifdef HAVE_MUELU_TPETRA
else if (paramList.isType<RCP<TpCrsMat> >(parameterName)) {
RCP<TpCrsMat> tM = paramList.get<RCP<TpCrsMat> >(parameterName);
paramList.remove(parameterName);
RCP<XpCrsMat> xM = rcp_dynamic_cast<XpCrsMat>(tM, true);
paramList.set<RCP<XpCrsMat> >(parameterName, xM);
return true;
} else if (paramList.isType<RCP<tMV> >(parameterName)) {
RCP<tMV> tpetra_X = paramList.get<RCP<tMV> >(parameterName);
paramList.remove(parameterName);
RCP<XpMultVec> X = MueLu::TpetraMultiVector_To_XpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(tpetra_X);
paramList.set<RCP<XpMultVec> >(parameterName, X);
TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(X));
return true;
} else if (paramList.isType<RCP<tMagMV> >(parameterName)) {
RCP<tMagMV> tpetra_X = paramList.get<RCP<tMagMV> >(parameterName);
paramList.remove(parameterName);
RCP<XpMagMultVec> X = MueLu::TpetraMultiVector_To_XpetraMultiVector<Magnitude,LocalOrdinal,GlobalOrdinal,Node>(tpetra_X);
paramList.set<RCP<XpMagMultVec> >(parameterName, X);
TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(X));
return true;
}
# if defined(HAVE_TPETRA_INST_DOUBLE) && defined(HAVE_TPETRA_INST_FLOAT)
else if (paramList.isType<RCP<tHalfMagMV> >(parameterName)) {
RCP<tHalfMagMV> tpetra_hX = paramList.get<RCP<tHalfMagMV> >(parameterName);
paramList.remove(parameterName);
RCP<tMagMV> tpetra_X = rcp(new tMagMV(tpetra_hX->getMap(),tpetra_hX->getNumVectors()));
Tpetra::deep_copy(*tpetra_X,*tpetra_hX);
RCP<XpMagMultVec> X = MueLu::TpetraMultiVector_To_XpetraMultiVector<Magnitude,LocalOrdinal,GlobalOrdinal,Node>(tpetra_X);
paramList.set<RCP<XpMagMultVec> >(parameterName, X);
TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(X));
return true;
}
# endif
#endif
#ifdef HAVE_MUELU_EPETRA
else if (paramList.isType<RCP<Epetra_CrsMatrix> >(parameterName)) {
RCP<Epetra_CrsMatrix> eM = paramList.get<RCP<Epetra_CrsMatrix> >(parameterName);
paramList.remove(parameterName);
RCP<XpEpCrsMat> xeM = rcp(new XpEpCrsMat(eM));
RCP<XpCrsMat> xCrsM = rcp_dynamic_cast<XpCrsMat>(xeM, true);
RCP<XpCrsMatWrap> xwM = rcp(new XpCrsMatWrap(xCrsM));
RCP<XpMat> xM = rcp_dynamic_cast<XpMat>(xwM);
paramList.set<RCP<XpMat> >(parameterName, xM);
return true;
} else if (paramList.isType<RCP<Epetra_MultiVector> >(parameterName)) {
RCP<Epetra_MultiVector> epetra_X = Teuchos::null;
epetra_X = paramList.get<RCP<Epetra_MultiVector> >(parameterName);
paramList.remove(parameterName);
RCP<Xpetra::EpetraMultiVectorT<int,Node> > xpEpX = rcp(new Xpetra::EpetraMultiVectorT<int,Node>(epetra_X));
RCP<Xpetra::MultiVector<Scalar,int,int,Node> > xpEpXMult = rcp_dynamic_cast<Xpetra::MultiVector<Scalar,int,int,Node> >(xpEpX, true);
RCP<XpMultVec> X = rcp_dynamic_cast<XpMultVec>(xpEpXMult, true);
paramList.set<RCP<XpMultVec> >(parameterName, X);
return true;
}
#endif
else if (paramList.isType<RCP<const ThyLinOpBase> >(parameterName)) {
RCP<const ThyLinOpBase> thyM = paramList.get<RCP<const ThyLinOpBase> >(parameterName);
paramList.remove(parameterName);
RCP<const XpCrsMat> crsM = XpThyUtils::toXpetra(thyM);
TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(crsM));
// MueLu needs a non-const object as input
RCP<XpCrsMat> crsMNonConst = rcp_const_cast<XpCrsMat>(crsM);
TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(crsMNonConst));
// wrap as an Xpetra::Matrix that MueLu can work with
RCP<XpMat> M = rcp(new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(crsMNonConst));
paramList.set<RCP<XpMat> >(parameterName, M);
return true;
} else if (paramList.isType<RCP<const ThyDiagLinOpBase> >(parameterName)) {
RCP<const ThyDiagLinOpBase> thyM = paramList.get<RCP<const ThyDiagLinOpBase> >(parameterName);
paramList.remove(parameterName);
RCP<const Thyra::VectorBase<Scalar> > diag = thyM->getDiag();

RCP<const XpVec> xpDiag;
#ifdef HAVE_MUELU_TPETRA
if (!rcp_dynamic_cast<const thyTpV>(diag).is_null()) {
RCP<const tV> tDiag = Thyra::TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getConstTpetraVector(diag);
if (!tDiag.is_null())
xpDiag = Xpetra::toXpetra(tDiag);
}
#endif
#ifdef HAVE_MUELU_EPETRA
if (xpDiag.is_null()) {
RCP<const Epetra_Comm> comm = Thyra::get_Epetra_Comm(*rcp_dynamic_cast<const ThyVSBase>(thyM->range())->getComm());
RCP<const Epetra_Map> map = Thyra::get_Epetra_Map(*(thyM->range()), comm);
if (!map.is_null()) {
RCP<const Epetra_Vector> eDiag = Thyra::get_Epetra_Vector(*map, diag);
RCP<Epetra_Vector> nceDiag = rcp_const_cast<Epetra_Vector>(eDiag);
RCP<Xpetra::EpetraVectorT<int,Node> > xpEpDiag = rcp(new Xpetra::EpetraVectorT<int,Node>(nceDiag));
xpDiag = rcp_dynamic_cast<XpVec>(xpEpDiag, true);
}
}
#endif
TEUCHOS_ASSERT(!xpDiag.is_null());
RCP<XpMat> M = Xpetra::MatrixFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(xpDiag);
paramList.set<RCP<XpMat> >(parameterName, M);
return true;
}
else {
TEUCHOS_TEST_FOR_EXCEPTION(true, MueLu::Exceptions::RuntimeError, "Parameter " << parameterName << " has wrong type.");
return false;
}
} else
return false;
}


// Constructors/initializers/accessors

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
Expand Down Expand Up @@ -317,7 +150,7 @@ namespace Thyra {
// rebuild preconditioner if startingOver == true
// reuse preconditioner if startingOver == false
const bool startingOver = (thyra_precOp.is_null() || !paramList.isParameter("refmaxwell: enable reuse") || !paramList.get<bool>("refmaxwell: enable reuse"));
const bool useHalfPrecision = paramList.get<bool>("refmaxwell: half precision", false) && bIsTpetra;
const bool useHalfPrecision = paramList.get<bool>("half precision", false) && bIsTpetra;

RCP<XpOp> xpPrecOp;
if (startingOver == true) {
Expand Down
11 changes: 10 additions & 1 deletion packages/muelu/doc/UsersGuide/masterList.xml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -146,6 +146,15 @@
<description>Pass parameters to the underlying linear algebra</description>
<comment-ML>not supported by ML</comment-ML>
<visible>false</visible>
</parameter>

<parameter>
<name>half precision</name>
<type>bool</type>
<description>Build half precision preconditioner</description>
<comment-ML>not supported by ML</comment-ML>
<visible>false</visible>
<default>false</default>
</parameter>

</general>
Expand Down Expand Up @@ -942,7 +951,7 @@
both (``both'').</description>
<comment-ML>not supported by ML</comment-ML>
</parameter>



</aggregate_qualities>
Expand Down
2 changes: 2 additions & 0 deletions packages/muelu/doc/UsersGuide/paramlist_hidden.tex
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,8 @@

\cba{matvec params}{\parameterlist}{Pass parameters to the underlying linear algebra}

\cbb{half precision}{bool}{false}{Build half precision preconditioner}

\cbb{smoother: pre or post}{string}{"both"}{Pre- and post-smoother combination. Possible values: "pre" (only pre-smoother), "post" (only post-smoother), "both" (both pre-and post-smoothers), "none" (no smoothing).}

\cbb{smoother: type}{string}{"RELAXATION"}{Smoother type. Possible values: see Table~\ref{tab:smoothers}.}
Expand Down
2 changes: 2 additions & 0 deletions packages/muelu/src/MueCentral/MueLu_Hierarchy_def.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1301,6 +1301,8 @@ namespace MueLu {
oss << "\n--------------------------------------------------------------------------------\n";
oss << "--- Multigrid Summary " << std::setw(28) << std::left << label << "---\n";
oss << "--------------------------------------------------------------------------------" << std::endl;
if (verbLevel & Parameters1)
oss << "Scalar = " << Teuchos::ScalarTraits<Scalar>::name() << std::endl;
oss << "Number of levels = " << numLevels << std::endl;
oss << "Operator complexity = " << std::setprecision(2) << std::setiosflags(std::ios::fixed)
<< GetOperatorComplexity() << std::endl;
Expand Down
3 changes: 3 additions & 0 deletions packages/muelu/src/MueCentral/MueLu_MasterList.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -184,6 +184,7 @@ namespace MueLu {
"<Parameter name=\"parameterlist: syntax\" type=\"string\" value=\"muelu\"/>"
"<Parameter name=\"hierarchy label\" type=\"string\" value=\"\"/>"
"<ParameterList name=\"matvec params\"/>"
"<Parameter name=\"half precision\" type=\"bool\" value=\"false\"/>"
"<Parameter name=\"smoother: pre or post\" type=\"string\" value=\"both\"/>"
"<Parameter name=\"smoother: type\" type=\"string\" value=\"RELAXATION\"/>"
"<Parameter name=\"smoother: pre type\" type=\"string\" value=\"RELAXATION\"/>"
Expand Down Expand Up @@ -570,6 +571,8 @@ namespace MueLu {

("matvec params","matvec params")

("half precision","half precision")

("smoother: pre or post","smoother: pre or post")

("smoother: type","smoother: type")
Expand Down

0 comments on commit 1a08fde

Please sign in to comment.