Skip to content

Commit

Permalink
Merge pull request #12675 from cgcgcg/tekoBlockIterative
Browse files Browse the repository at this point in the history 
Teko: Enable preconditioned sub-solves in block Jacobi and Gauss-Seidel
  • Loading branch information
@cgcgcg
cgcgcg authored Feb 13, 2024
2 parents c2e8cef + 4b8ba7a commit ffdcd00
Show file tree
Hide file tree
Showing 7 changed files with 188 additions and 27 deletions.
2 changes: 1 addition & 1 deletion packages/teko/cmake/Dependencies.cmake
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
SET(LIB_REQUIRED_DEP_PACKAGES Teuchos Thyra Stratimikos Anasazi Tpetra ThyraTpetraAdapters)
SET(LIB_OPTIONAL_DEP_PACKAGES Epetra Isorropia Ifpack2 Ifpack AztecOO Amesos Amesos2 Belos ThyraEpetraAdapters ThyraEpetraExtAdapters ML)
SET(TEST_REQUIRED_DEP_PACKAGES Ifpack2)
SET(TEST_REQUIRED_DEP_PACKAGES Belos Ifpack2)
SET(TEST_OPTIONAL_DEP_PACKAGES)
SET(LIB_REQUIRED_DEP_TPLS)
SET(LIB_OPTIONAL_DEP_TPLS)
Expand Down
66 changes: 49 additions & 17 deletions packages/teko/src/Teko_BlockInvDiagonalStrategy.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,6 +45,7 @@
*/

#include "Teko_BlockInvDiagonalStrategy.hpp"
#include "Teko_InverseFactory.hpp"

namespace Teko {

Expand All @@ -65,6 +66,21 @@ InvFactoryDiagStrategy::InvFactoryDiagStrategy(
defaultInvFact_ = defaultFact;
}

InvFactoryDiagStrategy::InvFactoryDiagStrategy(
const std::vector<Teuchos::RCP<InverseFactory> >& inverseFactories,
const std::vector<Teuchos::RCP<InverseFactory> >& preconditionerFactories,
const Teuchos::RCP<InverseFactory>& defaultInverseFact,
const Teuchos::RCP<InverseFactory>& defaultPreconditionerFact) {
invDiagFact_ = inverseFactories;
precDiagFact_ = preconditionerFactories;

if (defaultInverseFact == Teuchos::null)
defaultInvFact_ = invDiagFact_[0];
else
defaultInvFact_ = defaultInverseFact;
defaultPrecFact_ = defaultPreconditionerFact;
}

/** returns an (approximate) inverse of the diagonal blocks of A
* where A is closely related to the original source for invD0 and invD1
*/
Expand All @@ -73,35 +89,51 @@ void InvFactoryDiagStrategy::getInvD(const BlockedLinearOp& A, BlockPrecondition
Teko_DEBUG_SCOPE("InvFactoryDiagStrategy::getInvD", 10);

// loop over diagonals, build an inverse operator for each
int diagCnt = A->productRange()->numBlocks();
int invCnt = invDiagFact_.size();
size_t diagCnt = A->productRange()->numBlocks();

Teko_DEBUG_MSG("# diags = " << diagCnt << ", # inverses = " << invCnt, 6);
Teko_DEBUG_MSG("# diags = " << diagCnt << ", # inverses = " << invDiagFact_.size(), 6);

const std::string opPrefix = "JacobiDiagOp";
if (diagCnt <= invCnt) {
for (int i = 0; i < diagCnt; i++)
invDiag.push_back(buildInverse(*invDiagFact_[i], getBlock(i, i, A), state, opPrefix, i));
} else {
for (int i = 0; i < invCnt; i++)
invDiag.push_back(buildInverse(*invDiagFact_[i], getBlock(i, i, A), state, opPrefix, i));

for (int i = invCnt; i < diagCnt; i++)
invDiag.push_back(buildInverse(*defaultInvFact_, getBlock(i, i, A), state, opPrefix, i));
for (size_t i = 0; i < diagCnt; i++) {
auto precFact = ((i < precDiagFact_.size()) && (!precDiagFact_[i].is_null()))
? precDiagFact_[i]
: defaultPrecFact_;
auto invFact = (i < invDiagFact_.size()) ? invDiagFact_[i] : defaultInvFact_;
invDiag.push_back(buildInverse(*invFact, precFact, getBlock(i, i, A), state, opPrefix, i));
}
}

LinearOp InvFactoryDiagStrategy::buildInverse(const InverseFactory& invFact, const LinearOp& matrix,
LinearOp InvFactoryDiagStrategy::buildInverse(const InverseFactory& invFact,
Teuchos::RCP<InverseFactory>& precFact,
const LinearOp& matrix,
BlockPreconditionerState& state,
const std::string& opPrefix, int i) const {
std::stringstream ss;
ss << opPrefix << "_" << i;

ModifiableLinearOp& invOp = state.getModifiableOp(ss.str());
ModifiableLinearOp& invOp = state.getModifiableOp(ss.str());
ModifiableLinearOp& precOp = state.getModifiableOp("prec_" + ss.str());

if (precFact != Teuchos::null) {
if (precOp == Teuchos::null) {
precOp = precFact->buildInverse(matrix);
state.addModifiableOp("prec_" + ss.str(), precOp);
} else {
Teko::rebuildInverse(*precFact, matrix, precOp);
}
}

if (invOp == Teuchos::null)
invOp = Teko::buildInverse(invFact, matrix);
else
rebuildInverse(invFact, matrix, invOp);
if (precOp.is_null())
invOp = Teko::buildInverse(invFact, matrix);
else
invOp = Teko::buildInverse(invFact, matrix, precOp);
else {
if (precOp.is_null())
Teko::rebuildInverse(invFact, matrix, invOp);
else
Teko::rebuildInverse(invFact, matrix, precOp, invOp);
}

return invOp;
}
Expand Down
13 changes: 11 additions & 2 deletions packages/teko/src/Teko_BlockInvDiagonalStrategy.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -139,6 +139,12 @@ class InvFactoryDiagStrategy : public BlockInvDiagonalStrategy {
InvFactoryDiagStrategy(const std::vector<Teuchos::RCP<InverseFactory> > &factories,
const Teuchos::RCP<InverseFactory> &defaultFact = Teuchos::null);

InvFactoryDiagStrategy(
const std::vector<Teuchos::RCP<InverseFactory> > &inverseFactories,
const std::vector<Teuchos::RCP<InverseFactory> > &preconditionerFactories,
const Teuchos::RCP<InverseFactory> &defaultInverseFact = Teuchos::null,
const Teuchos::RCP<InverseFactory> &defaultPreconditionerFact = Teuchos::null);

virtual ~InvFactoryDiagStrategy() {}

/** returns an (approximate) inverse of the diagonal blocks of A
Expand All @@ -157,11 +163,14 @@ class InvFactoryDiagStrategy : public BlockInvDiagonalStrategy {
protected:
// stored inverse operators
std::vector<Teuchos::RCP<InverseFactory> > invDiagFact_;
std::vector<Teuchos::RCP<InverseFactory> > precDiagFact_;
Teuchos::RCP<InverseFactory> defaultInvFact_;
Teuchos::RCP<InverseFactory> defaultPrecFact_;

//! Conveinence function for building inverse operators
LinearOp buildInverse(const InverseFactory &invFact, const LinearOp &matrix,
BlockPreconditionerState &state, const std::string &opPrefix, int i) const;
LinearOp buildInverse(const InverseFactory &invFact, Teuchos::RCP<InverseFactory> &precFact,
const LinearOp &matrix, BlockPreconditionerState &state,
const std::string &opPrefix, int i) const;

private:
InvFactoryDiagStrategy();
Expand Down
36 changes: 33 additions & 3 deletions packages/teko/src/Teko_GaussSeidelPreconditionerFactory.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -102,19 +102,25 @@ void GaussSeidelPreconditionerFactory::initializeFromParameterList(
pl.print(DEBUG_STREAM);
Teko_DEBUG_MSG_END();

const std::string inverse_type = "Inverse Type";
const std::string inverse_type = "Inverse Type";
const std::string preconditioner_type = "Preconditioner Type";
std::vector<RCP<InverseFactory> > inverses;
std::vector<RCP<InverseFactory> > preconditioners;

RCP<const InverseLibrary> invLib = getInverseLibrary();

// get string specifying default inverse
std::string invStr = "Amesos";
std::string invStr = "Amesos";
std::string precStr = "None";
if (pl.isParameter(inverse_type)) invStr = pl.get<std::string>(inverse_type);
if (pl.isParameter(preconditioner_type)) precStr = pl.get<std::string>(preconditioner_type);
if (pl.isParameter("Use Upper Triangle"))
solveType_ = pl.get<bool>("Use Upper Triangle") ? GS_UseUpperTriangle : GS_UseLowerTriangle;

Teko_DEBUG_MSG("GSPrecFact: Building default inverse \"" << invStr << "\"", 5);
RCP<InverseFactory> defaultInverse = invLib->getInverseFactory(invStr);
RCP<InverseFactory> defaultPrec;
if (precStr != "None") defaultPrec = invLib->getInverseFactory(precStr);

// now check individual solvers
Teuchos::ParameterList::ConstIterator itr;
Expand Down Expand Up @@ -144,14 +150,38 @@ void GaussSeidelPreconditionerFactory::initializeFromParameterList(
inverses[position - 1] = invLib->getInverseFactory(invStr2);
} else
inverses[position - 1] = invLib->getInverseFactory(invStr2);
} else if (fieldName.compare(0, preconditioner_type.length(), preconditioner_type) == 0 &&
fieldName != preconditioner_type) {
int position = -1;
std::string preconditioner, type;

// figure out position
std::stringstream ss(fieldName);
ss >> preconditioner >> type >> position;

if (position <= 0) {
Teko_DEBUG_MSG("Gauss-Seidel \"Preconditioner Type\" must be a (strictly) positive integer",
1);
}

// inserting preconditioner factory into vector
std::string precStr2 = pl.get<std::string>(fieldName);
Teko_DEBUG_MSG(
"GSPrecFact: Building preconditioner " << position << " \"" << precStr2 << "\"", 5);
if (position > (int)preconditioners.size()) {
preconditioners.resize(position, defaultPrec);
preconditioners[position - 1] = invLib->getInverseFactory(precStr2);
} else
preconditioners[position - 1] = invLib->getInverseFactory(precStr2);
}
}

// use default inverse
if (inverses.size() == 0) inverses.push_back(defaultInverse);

// based on parameter type build a strategy
invOpsStrategy_ = rcp(new InvFactoryDiagStrategy(inverses, defaultInverse));
invOpsStrategy_ =
rcp(new InvFactoryDiagStrategy(inverses, preconditioners, defaultInverse, defaultPrec));
}

} // namespace Teko
34 changes: 31 additions & 3 deletions packages/teko/src/Teko_JacobiPreconditionerFactory.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -119,14 +119,19 @@ void JacobiPreconditionerFactory::initializeFromParameterList(const Teuchos::Par
pl.print(DEBUG_STREAM);
Teko_DEBUG_MSG_END();

const std::string inverse_type = "Inverse Type";
const std::string inverse_type = "Inverse Type";
const std::string preconditioner_type = "Preconditioner Type";
std::vector<RCP<InverseFactory> > inverses;
std::vector<RCP<InverseFactory> > preconditioners;

RCP<const InverseLibrary> invLib = getInverseLibrary();

// get string specifying default inverse
std::string invStr = "Amesos";
std::string invStr = "Amesos";
std::string precStr = "None";
if (pl.isParameter(inverse_type)) invStr = pl.get<std::string>(inverse_type);
RCP<InverseFactory> defaultPrec;
if (precStr != "None") defaultPrec = invLib->getInverseFactory(precStr);

Teko_DEBUG_MSG("JacobiPrecFact: Building default inverse \"" << invStr << "\"", 5);
RCP<InverseFactory> defaultInverse = invLib->getInverseFactory(invStr);
Expand Down Expand Up @@ -160,14 +165,37 @@ void JacobiPreconditionerFactory::initializeFromParameterList(const Teuchos::Par
inverses[position - 1] = invLib->getInverseFactory(invStr2);
} else
inverses[position - 1] = invLib->getInverseFactory(invStr2);
} else if (fieldName.compare(0, preconditioner_type.length(), preconditioner_type) == 0 &&
fieldName != preconditioner_type) {
int position = -1;
std::string preconditioner, type;

// figure out position
std::stringstream ss(fieldName);
ss >> preconditioner >> type >> position;

if (position <= 0) {
Teko_DEBUG_MSG("Jacobi \"Preconditioner Type\" must be a (strictly) positive integer", 1);
}

// inserting preconditioner factory into vector
std::string precStr2 = pl.get<std::string>(fieldName);
Teko_DEBUG_MSG(
"JacobiPrecFact: Building preconditioner " << position << " \"" << precStr2 << "\"", 5);
if (position > (int)preconditioners.size()) {
preconditioners.resize(position, defaultPrec);
preconditioners[position - 1] = invLib->getInverseFactory(precStr2);
} else
preconditioners[position - 1] = invLib->getInverseFactory(precStr2);
}
}

// use default inverse
if (inverses.size() == 0) inverses.push_back(defaultInverse);

// based on parameter type build a strategy
invOpsStrategy_ = rcp(new InvFactoryDiagStrategy(inverses, defaultInverse));
invOpsStrategy_ =
rcp(new InvFactoryDiagStrategy(inverses, preconditioners, defaultInverse, defaultPrec));
}

} // namespace Teko
63 changes: 62 additions & 1 deletion packages/teko/tests/src/tBlockJacobiPreconditionerFactory_tpetra.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -74,12 +74,16 @@
#include "Trilinos_Util_CrsMatrixGallery.h"

#include <vector>

#include "Teko_StratimikosFactory.hpp"
#include "Teko_Utilities.hpp"
#include "Teko_TpetraHelpers.hpp"
#include "Thyra_TpetraLinearOp.hpp"
#include "Tpetra_CrsMatrix.hpp"

#include "Teuchos_AbstractFactoryStd.hpp"

#include "Stratimikos_DefaultLinearSolverBuilder.hpp"

namespace Teko {
namespace Test {

Expand Down Expand Up @@ -185,6 +189,13 @@ int tBlockJacobiPreconditionerFactory_tpetra::runTest(int verbosity, std::ostrea
failcount += status ? 0 : 1;
totalrun++;

status = test_iterativeSolves(verbosity, failstrm);
Teko_TEST_MSG(stdstrm, 1, " \"iterativeSolves\" ... PASSED",
" \"iterativeSolves\" ... FAILED");
allTests &= status;
failcount += status ? 0 : 1;
totalrun++;

status = allTests;
if (verbosity >= 10) {
Teko_TEST_MSG(failstrm, 0, "tBlockJacobiPreconditionedFactory...PASSED",
Expand Down Expand Up @@ -413,5 +424,55 @@ bool tBlockJacobiPreconditionerFactory_tpetra::test_isCompatible(int verbosity,
return allPassed;
}

bool tBlockJacobiPreconditionerFactory_tpetra::test_iterativeSolves(int verbosity,
std::ostream& os) {
using Thyra::zero;

bool status = false;
bool allPassed = true;

// allocate new linear operator
const RCP<Thyra::PhysicallyBlockedLinearOpBase<ST> > blkOp = Thyra::defaultBlockedLinearOp<ST>();
blkOp->beginBlockFill(3, 3);
blkOp->setBlock(0, 0, F_);
blkOp->setBlock(0, 1, Bt_);
blkOp->setBlock(1, 0, B_);
blkOp->setBlock(1, 1, C_);
blkOp->setBlock(1, 2, B_);
blkOp->setBlock(2, 1, Bt_);
blkOp->setBlock(2, 2, F_);
blkOp->endBlockFill();

// build stratimikos factory, adding Teko's version
Stratimikos::DefaultLinearSolverBuilder stratFactory;
stratFactory.setPreconditioningStrategyFactory(
Teuchos::abstractFactoryStd<Thyra::PreconditionerFactoryBase<double>,
Teko::StratimikosFactory>(),
"Teko");
RCP<ParameterList> params = Teuchos::rcp(new ParameterList(*stratFactory.getValidParameters()));
ParameterList& tekoList = params->sublist("Preconditioner Types").sublist("Teko");
tekoList.set("Write Block Operator", false);
tekoList.set("Test Block Operator", false);
tekoList.set("Strided Blocking", "1 1");
tekoList.set("Inverse Type", "BGS");
ParameterList& ifl = tekoList.sublist("Inverse Factory Library");
ifl.sublist("BGS").set("Type", "Block Gauss-Seidel");
ifl.sublist("BGS").set("Inverse Type", "Belos");
ifl.sublist("BGS").set("Preconditioner Type", "Ifpack2");

RCP<Teko::InverseLibrary> invLib = Teko::InverseLibrary::buildFromParameterList(ifl);
RCP<const Teko::InverseFactory> invFact = invLib->getInverseFactory("BGS");

Teko::ModifiableLinearOp inv = Teko::buildInverse(*invFact, blkOp);
TEST_ASSERT(!inv.is_null(), "Constructed preconditioner is null");

if (!inv.is_null()) {
Teko::rebuildInverse(*invFact, blkOp, inv);
TEST_ASSERT(!inv.is_null(), "Constructed preconditioner is null");
}

return allPassed;
}

} // end namespace Test
} // end namespace Teko
1 change: 1 addition & 0 deletions packages/teko/tests/src/tBlockJacobiPreconditionerFactory_tpetra.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -73,6 +73,7 @@ class tBlockJacobiPreconditionerFactory_tpetra : public UnitTest {
bool test_initializePrec(int verbosity, std::ostream& os);
bool test_uninitializePrec(int verbosity, std::ostream& os);
bool test_isCompatible(int verbosity, std::ostream& os);
bool test_iterativeSolves(int verbosity, std::ostream& os);

protected:
double tolerance_;
Expand Down

0 comments on commit ffdcd00

Please sign in to comment.