Teko: Enable preconditioned sub-solves in block Jacobi and Gauss-Seidel

packages/teko/cmake/Dependencies.cmake

@@ -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)

packages/teko/src/Teko_BlockInvDiagonalStrategy.cpp

@@ -45,6 +45,7 @@
 */
 
 #include "Teko_BlockInvDiagonalStrategy.hpp"
+#include "Teko_InverseFactory.hpp"
 
 namespace Teko {
 
@@ -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
  */
@@ -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;
 }

packages/teko/src/Teko_BlockInvDiagonalStrategy.hpp

@@ -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
@@ -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();

packages/teko/src/Teko_GaussSeidelPreconditionerFactory.cpp

@@ -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;
@@ -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

packages/teko/src/Teko_JacobiPreconditionerFactory.cpp

@@ -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);
@@ -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

packages/teko/tests/src/tBlockJacobiPreconditionerFactory_tpetra.cpp

@@ -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 {
 
@@ -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",
@@ -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

packages/teko/tests/src/tBlockJacobiPreconditionerFactory_tpetra.hpp

@@ -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_;