‍Data files. The drivers read Matrix-Market files from examples/data/, which ships nslhs_test.mm and nsrhs_test.mm. All three example programs (both BuildPreconditioner drivers and StridedSolve) additionally read ../data/nsjac.mm, which is not currently in examples/data/ — supply your own Navier–Stokes Jacobian there before running them.

All examples are built by the package's CMake when Teko's tests/examples are enabled; the binaries land next to their source directories in the build tree.

BuildPreconditioner — the canonical driver

examples/BuildPreconditioner/ contains two drivers that read a strided Navier–Stokes matrix, build a block preconditioner, and solve with Belos.

**example-driver.cpp** — the native Tpetra path. It splits the monolithic matrix with a StridedTpetraOperator ({2, 1} = 2 velocity + 1 pressure per node), builds an LSC preconditioner via InvLSCStrategy, wraps it in a TpetraBlockPreconditioner, and hands it to Belos::BlockGmresSolMgr as a right preconditioner.
**example-driver-belos.cpp** — the Thyra/Stratimikos path. It assembles a 2×2 block operator with Thyra::block2x2, builds an InverseLibrary from a hand-built parameter list (a "Gauss-Seidel" entry over Ifpack2 block solves), calls Teko::buildInverse, and solves with Belos. This is the driver dissected in Getting Started.

Run (from the example's build directory):

$ mpirun -np 2 ./Teko_example-driver-belos.exe

Full source (Thyra/Stratimikos path):

// @HEADER
// *****************************************************************************
//      Teko: A package for block and physics based preconditioning
//
// Copyright 2010 NTESS and the Teko contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER
 
// Teuchos includes /*@ \label{lned:being-includes} @*/
#include "Teuchos_ConfigDefs.hpp"
#include "Teuchos_GlobalMPISession.hpp"
#include "Teuchos_RCP.hpp"
#include "Teuchos_XMLParameterListHelpers.hpp"
#include "Teuchos_DefaultComm.hpp"
#include "Teuchos_AbstractFactoryStd.hpp"
 
#include "mpi.h"
 
#include "Tpetra_CrsMatrix.hpp"
#include "Tpetra_Core.hpp"
#include "MatrixMarket_Tpetra.hpp"
 
// Teko-Package includes
#include "Teko_Utilities.hpp"
#include "Teko_InverseFactory.hpp"
#include "Teko_InverseLibrary.hpp"
#include "Teko_LSCPreconditionerFactory.hpp"
#include "Teko_InvLSCStrategy.hpp"
#include "Teko_SIMPLEPreconditionerFactory.hpp"
 
#include "Thyra_TpetraLinearOp.hpp"
#include "Thyra_TpetraVector.hpp"
#include "Thyra_TpetraVectorSpace.hpp"
 
#include "Stratimikos_DefaultLinearSolverBuilder.hpp"
 
// Belos includes
#include "BelosConfigDefs.hpp"
#include "BelosLinearProblem.hpp"
#include "BelosBlockGmresSolMgr.hpp"
#include "BelosThyraAdapter.hpp"  // Requires Stratimikos...
 
#include <iostream> /*@ \label{lned:end-includes} @*/
 
// for simplicity
using Teuchos::RCP;
using Teuchos::rcp;
using Teuchos::rcpFromRef;
 
RCP<Teuchos::ParameterList> buildLibPL();
 
int main(int argc, char* argv[]) {
  typedef double ST;
  typedef Thyra::MultiVectorBase<ST> MV;
  typedef Thyra::LinearOpBase<ST> OP;
 
  typedef Tpetra::Vector<ST> TP_Vec;
  typedef Tpetra::CrsMatrix<ST> TP_Crs;
  typedef Tpetra::Operator<ST> TP_Op;
 
  typedef TP_Vec::local_ordinal_type LO;
  typedef TP_Vec::global_ordinal_type GO;
  typedef TP_Vec::node_type NT;
 
  // calls MPI_Init and MPI_Finalize
  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
 
  RCP<Stratimikos::DefaultLinearSolverBuilder> linearSolverBuilder =
      Teuchos::rcp(new Stratimikos::DefaultLinearSolverBuilder);
 
  // Build the Tpetra matrices and vectors
 
  // read in the CRS matrix
  RCP<TP_Crs> crsMat = Tpetra::MatrixMarket::Reader<TP_Crs>::readSparseFile(
      "../data/nsjac.mm", Tpetra::getDefaultComm());
 
  RCP<TP_Crs> zeroCrsMat = rcp(new TP_Crs(*crsMat, Teuchos::Copy));
  zeroCrsMat->setAllToScalar(0.0);
 
  RCP<TP_Op> Mat     = crsMat;
  RCP<TP_Op> zeroMat = zeroCrsMat;
 
  // Allocate some right handside vectors
  RCP<TP_Vec> x0_tp = rcp(new TP_Vec(Mat->getDomainMap()));
  RCP<TP_Vec> x1_tp = rcp(new TP_Vec(Mat->getDomainMap()));
  RCP<TP_Vec> b0_tp = rcp(new TP_Vec(Mat->getRangeMap()));
  RCP<TP_Vec> b1_tp = rcp(new TP_Vec(Mat->getRangeMap()));
  b0_tp->randomize();
  b1_tp->randomize();
 
  RCP<const Thyra::TpetraVectorSpace<ST, LO, GO, NT> > domain =
      Thyra::tpetraVectorSpace<ST>(Mat->getDomainMap());
  RCP<const Thyra::TpetraVectorSpace<ST, LO, GO, NT> > range =
      Thyra::tpetraVectorSpace<ST>(Mat->getRangeMap());
 
  // Build Teko compatible matrices and vectors
 
  // convert them to teko compatible sub vectors
  Teko::MultiVector x0_th = Thyra::tpetraVector(domain, x0_tp);
  Teko::MultiVector x1_th = Thyra::tpetraVector(domain, x1_tp);
  Teko::MultiVector b0_th = Thyra::tpetraVector(range, b0_tp);
  Teko::MultiVector b1_th = Thyra::tpetraVector(range, b1_tp);
  std::vector<Teko::MultiVector> x_vec;
  x_vec.push_back(x0_th);
  x_vec.push_back(x1_th);
  std::vector<Teko::MultiVector> b_vec;
  b_vec.push_back(b0_th);
  b_vec.push_back(b1_th);
 
  Teko::MultiVector x =
      Teko::buildBlockedMultiVector(x_vec);  // these will be used in the Teko solve
  Teko::MultiVector b = Teko::buildBlockedMultiVector(b_vec);
 
  // Build the Teko compatible linear system
  Teko::LinearOp thMat  = Thyra::tpetraLinearOp<double>(range, domain, Mat);
  Teko::LinearOp thZero = Thyra::tpetraLinearOp<double>(range, domain, zeroMat);
  Teko::LinearOp A =
      Thyra::block2x2(thMat, thZero, thZero, thMat);  // build an upper triangular 2x2
 
  // Build the preconditioner
 
  // build an InverseLibrary
  RCP<Teko::InverseLibrary> invLib =
      Teko::InverseLibrary::buildFromParameterList(*buildLibPL(), linearSolverBuilder);
 
  // build the inverse factory needed by the example preconditioner
  RCP<Teko::InverseFactory> inverse = invLib->getInverseFactory("Gauss-Seidel");
 
  // build the preconditioner from the jacobian
  Teko::LinearOp prec = Teko::buildInverse(*inverse, A);
 
  // Setup the Belos solver
 
  Teuchos::ParameterList belosList;
  belosList.set("Num Blocks", 200);              // Maximum number of blocks in Krylov factorization
  belosList.set("Block Size", 1);                // Blocksize to be used by iterative solver
  belosList.set("Maximum Iterations", 200);      // Maximum number of iterations allowed
  belosList.set("Maximum Restarts", 1);          // Maximum number of restarts allowed
  belosList.set("Convergence Tolerance", 1e-5);  // Relative convergence tolerance requested
  belosList.set(
      "Verbosity",
      33);  // Belos::Errors + Belos::Warnings + Belos::TimingDetails + Belos::StatusTestDetails );
  belosList.set("Output Frequency", 1);
  belosList.set("Output Style", 1);
 
  RCP<Belos::LinearProblem<double, MV, OP> > problem =
      rcp(new Belos::LinearProblem<double, MV, OP>(A, x, b));
  problem->setLeftPrec(prec);
  problem->setProblem();  // should check the return type!!!
 
  RCP<Belos::SolverManager<double, MV, OP> > solver =
      rcp(new Belos::BlockGmresSolMgr<double, MV, OP>(problem, rcpFromRef(belosList)));
 
  //
  // Perform solve
  //
  Belos::ReturnType ret = solver->solve();
 
  if (ret != Belos::Converged) {
    std::cout << std::endl << "ERROR:  Belos did not converge!" << std::endl;
    return -1;
  }
 
  return 0;
}
 
RCP<Teuchos::ParameterList> buildLibPL() {
  RCP<Teuchos::ParameterList> pl = rcp(new Teuchos::ParameterList());
 
  {
    Teuchos::ParameterList& sub_jac = pl->sublist("Jacobi");
    sub_jac.set("Type", "Block Jacobi");
    sub_jac.set("Inverse Type", "Ifpack2");
 
    Teuchos::ParameterList& sub_gs = pl->sublist("Gauss-Seidel");
    sub_gs.set("Type", "Block Gauss-Seidel");
    sub_gs.set("Use Upper Triangle", true);
    sub_gs.set("Inverse Type", "Ifpack2");
  }
  return pl;
}

Full source (native Tpetra path):

// @HEADER
// *****************************************************************************
//      Teko: A package for block and physics based preconditioning
//
// Copyright 2010 NTESS and the Teko contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER
 
// Teuchos includes
#include "Teuchos_ConfigDefs.hpp" /*@ \label{lned:being-includes} @*/
#include "Teuchos_RCP.hpp"
 
// Tpetra includes
#include "Tpetra_Core.hpp"
#include "Tpetra_CrsMatrix.hpp"
#include "Tpetra_Vector.hpp"
#include "Tpetra_MultiVector.hpp"
#include "MatrixMarket_Tpetra.hpp"
 
// Teko-Package includes
#include "Teko_Utilities.hpp"
#include "Teko_InverseFactory.hpp"
#include "Teko_InverseLibrary.hpp"
#include "Teko_StridedTpetraOperator.hpp"
#include "Teko_TpetraBlockPreconditioner.hpp"
#include "Teko_LSCPreconditionerFactory.hpp"
#include "Teko_InvLSCStrategy.hpp"
#include "Teko_SIMPLEPreconditionerFactory.hpp"
#include "Teko_ConfigDefs.hpp"
 
// Belos includes
#include "BelosConfigDefs.hpp"
#include "BelosLinearProblem.hpp"
#include "BelosBlockGmresSolMgr.hpp"
#include "BelosTpetraAdapter.hpp"
 
#include <iostream>
#include <vector> /*@ \label{lned:end-includes} @*/
 
// for simplicity
using Teuchos::RCP;
using Teuchos::rcp;
 
void run_driver() {
  using ST = double;
  using LO = Teko::LO;
  using GO = Teko::GO;
  using NT = Teko::NT;
 
  using map_t = Tpetra::Map<LO, GO, NT>;
  using crs_t = Tpetra::CrsMatrix<ST, LO, GO, NT>;
  using vec_t = Tpetra::Vector<ST, LO, GO, NT>;
  using mv_t  = Tpetra::MultiVector<ST, LO, GO, NT>;
  using op_t  = Tpetra::Operator<ST, LO, GO, NT>;
 
  auto comm = Tpetra::getDefaultComm();
 
  // Read in the matrix
  RCP<crs_t> A = Tpetra::MatrixMarket::Reader<crs_t>::readSparseFile("../data/nsjac.mm", comm);
 
  // Read in the RHS vector
  RCP<const map_t> rangeMap  = A->getRangeMap();
  RCP<const map_t> domainMap = A->getDomainMap();
 
  RCP<vec_t> b = Tpetra::MatrixMarket::Reader<crs_t>::readVectorFile("../data/nsrhs_test.mm", comm,
                                                                     rangeMap, false, false);
 
  // Allocate solution vector
  RCP<vec_t> x = rcp(new vec_t(domainMap));
 
  // Break apart the strided linear system
 
  std::vector<int> vec(2);
  vec[0] = 2;
  vec[1] = 1; /*@ \label{lned:define-strided} @*/
 
  Teuchos::RCP<Teko::TpetraHelpers::StridedTpetraOperator> sA =
      Teuchos::rcp(new Teko::TpetraHelpers::StridedTpetraOperator(vec, A));
 
  // Build the preconditioner /*@ \label{lned:construct-prec} @*/
 
  RCP<Teko::InverseLibrary> invLib =
      Teko::InverseLibrary::buildFromStratimikos(); /*@ \label{lned:define-inv-params} @*/
 
  RCP<Teko::InverseFactory> inverse /*@ \label{lned:define-inv-fact} @*/
      = invLib->getInverseFactory("Amesos2");
 
  RCP<Teko::NS::LSCStrategy> strategy =
      rcp(new Teko::NS::InvLSCStrategy(inverse, true)); /*@ \label{lned:const-prec-strategy} @*/
 
  RCP<Teko::BlockPreconditionerFactory> precFact =
      rcp(new Teko::NS::LSCPreconditionerFactory(strategy)); /*@ \label{lned:const-prec-fact} @*/
 
  Teko::TpetraHelpers::TpetraBlockPreconditioner prec(
      precFact); /*@ \label{lned:const-tpetra-prec} @*/
  prec.buildPreconditioner(sA);
 
  // Build and solve the linear system
 
  RCP<mv_t> X = x;
  RCP<mv_t> B = b;
 
  using problem_t        = Belos::LinearProblem<ST, mv_t, op_t>;
  RCP<problem_t> problem = rcp(new problem_t(A, X, B)); /*@ \label{lned:belos-solve} @*/
 
  RCP<op_t> precOp = Teuchos::rcp(&prec, false);
  problem->setRightPrec(precOp);
 
  const bool set = problem->setProblem();
  TEUCHOS_TEST_FOR_EXCEPTION(!set, std::runtime_error,
                             "Belos::LinearProblem::setProblem() failed.");
 
  using solver_t = Belos::BlockGmresSolMgr<ST, mv_t, op_t>;
 
  Teuchos::ParameterList belosList;
  belosList.set("Maximum Iterations", 1000);
  belosList.set("Convergence Tolerance", 1e-5);
  belosList.set("Num Blocks", 1000);
  belosList.set("Verbosity",
                Belos::Errors + Belos::Warnings + Belos::IterationDetails + Belos::FinalSummary);
  belosList.set("Output Frequency", 10);
 
  solver_t solver(problem, rcpFromRef(belosList));
  Belos::ReturnType result = solver.solve();
 
  TEUCHOS_TEST_FOR_EXCEPTION(result != Belos::Converged, std::runtime_error,
                             "Belos solver failed to converge.");
}
 
int main(int argc, char* argv[]) {
  Tpetra::ScopeGuard tpetraScope(&argc, &argv);
  { run_driver(); }
  return 0;
}

StridedSolve — SIMPLE on a strided operator, solver chosen at runtime

examples/StridedSolve/strided_solve.cpp reads ../data/nsjac.mm, splits it with a strided operator, and builds a SIMPLE preconditioner (Teko::NS::SIMPLEPreconditionerFactory(inverse, alpha)). The inner solver name is read from solverparams.xml (an ordinary Stratimikos solver list), so you can swap Amesos2/Belos/… without recompiling. A good template for "drive the inner solver from XML while building the block structure in code."

Full source:

// @HEADER
// *****************************************************************************
//      Teko: A package for block and physics based preconditioning
//
// Copyright 2010 NTESS and the Teko contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER
 
#include <sys/types.h>
#include <unistd.h>
 
#include "Teuchos_ConfigDefs.hpp"
#include "Teuchos_FancyOStream.hpp"
#include "Teuchos_RCP.hpp"
#include "Teuchos_XMLParameterListHelpers.hpp"
#include "Teuchos_CommHelpers.hpp"
 
// Thyra includes
#include "Thyra_LinearOpBase.hpp"
#include "Thyra_TpetraLinearOp.hpp"
#include "Thyra_TpetraThyraWrappers.hpp"
 
// Tpetra includes
#include "Tpetra_Core.hpp"
#include "Tpetra_Map.hpp"
#include "Tpetra_CrsMatrix.hpp"
#include "Tpetra_Vector.hpp"
#include "Tpetra_MultiVector.hpp"
#include "MatrixMarket_Tpetra.hpp"
 
// Teko includes
#include "Teko_Utilities.hpp"
#include "Teko_InverseFactory.hpp"
#include "Teko_SIMPLEPreconditionerFactory.hpp"
#include "Teko_LSCPreconditionerFactory.hpp"
#include "Teko_StridedTpetraOperator.hpp"
#include "Teko_TpetraBlockPreconditioner.hpp"
#include "Teko_ConfigDefs.hpp"
 
// Belos includes
#include "BelosConfigDefs.hpp"
#include "BelosLinearProblem.hpp"
#include "BelosBlockGmresSolMgr.hpp"
#include "BelosTpetraAdapter.hpp"
 
#include <iostream>
#include <fstream>
#include <cmath>
#include <string>
#include <vector>
 
using Teuchos::FancyOStream;
using Teuchos::ParameterList;
using Teuchos::RCP;
using Teuchos::rcp;
 
void run_driver(std::string solveName) {
  RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
 
  using ST = double;
  using LO = Teko::LO;
  using GO = Teko::GO;
  using NT = Teko::NT;
 
  using map_t = Tpetra::Map<LO, GO, NT>;
  using crs_t = Tpetra::CrsMatrix<ST, LO, GO, NT>;
  using vec_t = Tpetra::Vector<ST, LO, GO, NT>;
  using mv_t  = Tpetra::MultiVector<ST, LO, GO, NT>;
  using op_t  = Tpetra::Operator<ST, LO, GO, NT>;
 
  auto comm = Tpetra::getDefaultComm();
 
  RCP<FancyOStream> fos = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
  fos->setOutputToRootOnly(0);
 
  *fos << "Using \"" << solveName << "\" for approximate solve" << std::endl;
 
  const int numProc = comm->getSize();
  const int myPID   = comm->getRank();
 
  std::cout << "MPI_PID = " << myPID << ", UNIX_PID = " << getpid() << std::endl;
 
  *out << "Approaching Barrier: proc = " << numProc << ", pid = " << myPID << std::endl;
  Teuchos::barrier(*comm);
 
  RCP<Teuchos::ParameterList> paramList = Teuchos::getParametersFromXmlFile("solverparams.xml");
 
  *fos << "Reading matrix market file" << std::endl;
 
  RCP<crs_t> A = Tpetra::MatrixMarket::Reader<crs_t>::readSparseFile("../data/nsjac.mm", comm);
 
  RCP<const map_t> rangeMap  = A->getRangeMap();
  RCP<const map_t> domainMap = A->getDomainMap();
 
  RCP<vec_t> b = Tpetra::MatrixMarket::Reader<crs_t>::readVectorFile("../data/nsrhs_test.mm", comm,
                                                                     rangeMap, false, false);
 
  RCP<vec_t> x = Tpetra::MatrixMarket::Reader<crs_t>::readVectorFile("../data/nslhs_test.mm", comm,
                                                                     domainMap, false, false);
 
  *fos << "Building strided operator" << std::endl;
 
  std::vector<int> vars(2);
  vars[0] = 2;
  vars[1] = 1;
 
  Teuchos::RCP<Teko::TpetraHelpers::StridedTpetraOperator> sA =
      Teuchos::rcp(new Teko::TpetraHelpers::StridedTpetraOperator(vars, A));
 
  double alpha                      = 0.9;
  RCP<Teko::InverseFactory> inverse = Teko::invFactoryFromParamList(*paramList, solveName);
 
  RCP<Teko::BlockPreconditionerFactory> precFact =
      rcp(new Teko::NS::SIMPLEPreconditionerFactory(inverse, alpha));
 
  *fos << "Preconditioner factory built" << std::endl;
 
  Teko::TpetraHelpers::TpetraBlockPreconditioner prec(precFact);
  prec.buildPreconditioner(sA);
 
  *fos << "Preconditioner built" << std::endl;
 
  RCP<mv_t> X = x;
  RCP<mv_t> B = b;
 
  using problem_t        = Belos::LinearProblem<ST, mv_t, op_t>;
  RCP<problem_t> problem = rcp(new problem_t(A, X, B));
 
  RCP<op_t> precOp = Teuchos::rcp(&prec, false);
  problem->setRightPrec(precOp);
 
  const bool set = problem->setProblem();
  TEUCHOS_TEST_FOR_EXCEPTION(!set, std::runtime_error,
                             "Belos::LinearProblem::setProblem() failed.");
 
  *fos << "Setting solver parameters" << std::endl;
 
  Teuchos::ParameterList belosList;
  belosList.set("Maximum Iterations", 1000);
  belosList.set("Convergence Tolerance", 1e-5);
  belosList.set("Num Blocks", 50);
  belosList.set("Verbosity",
                Belos::Errors + Belos::Warnings + Belos::IterationDetails + Belos::FinalSummary);
  belosList.set("Output Frequency", 10);
 
  using solver_t = Belos::BlockGmresSolMgr<ST, mv_t, op_t>;
  solver_t solver(problem, rcpFromRef(belosList));
 
  *fos << "Solving" << std::endl;
 
  Belos::ReturnType result = solver.solve();
 
  TEUCHOS_TEST_FOR_EXCEPTION(result != Belos::Converged, std::runtime_error,
                             "Belos solver failed to converge.");
 
  *fos << "Solve converged" << std::endl;
}
 
int main(int argc, char* argv[]) {
  Tpetra::ScopeGuard tpetraScope(&argc, &argv);
  {
    std::string solveName = "Amesos2";
    if (argc > 1) solveName = argv[1];
 
    run_driver(solveName);
  }
  return 0;
}

Arbitrary Tpetra blocking — field blocks from monolithic GIDs

The strided examples assume a regular nodal ordering such as [u v p] at every mesh node. For application orderings that are not regular strides, use Teko::TpetraHelpers::BlockedTpetraOperator. It takes an explicit list of monolithic global IDs for each Teko block:

#include "Teko_BlockedTpetraOperator.hpp"
 
std::vector<std::vector<Teko::GO>> blockGIDs(3);
const auto baseMap = A->getDomainMap();
 
for (size_t lid = 0; lid < baseMap->getLocalNumElements(); ++lid) {
  const Teko::GO gid = baseMap->getGlobalElement(static_cast<Teko::LO>(lid));
  if (gidBelongsToVelocity(gid)) {
    blockGIDs[0].push_back(gid);
  } else if (gidBelongsToPressure(gid)) {
    blockGIDs[1].push_back(gid);
  } else {
    blockGIDs[2].push_back(gid);
  }
}
 
auto blockedA = Teuchos::rcp(
    new Teko::TpetraHelpers::BlockedTpetraOperator(blockGIDs, A, "A"));

The outer vector entry is the Teko block number (blockGIDs[0] becomes block 0, and so on). Each inner vector is this MPI rank's owned monolithic GID list for that block. The wrapped operator is assumed to be square, with matching domain and range maps. The unit test tests/src/Tpetra/tBlockedTpetraOperator.cpp is the current in-tree reference for constructing these lists, checking BlockedTpetraOperator::testAgainstFullOperator, rebuilding after matrix value changes, and applying reorder managers. A small runnable example built from that test would be a useful addition to examples/.

AddMultiplyPrecs — composing preconditioners

examples/AddMultiplyPrecs/Driver.cpp builds a block Gauss–Seidel factory and a block Jacobi factory over the same diagonal strategy, then composes them:

RCP<Teko::InverseFactory> inverse =
    Teko::invFactoryFromParamList(*paramList, "Amesos2");
RCP<Teko::BlockInvDiagonalStrategy> strategy =
    rcp(new Teko::InvFactoryDiagStrategy(inverse));
 
RCP<Teko::BlockPreconditionerFactory> GSFactory =
    rcp(new Teko::GaussSeidelPreconditionerFactory(Teko::GS_UseLowerTriangle, strategy));
RCP<Teko::BlockPreconditionerFactory> JacobiFactory =
    rcp(new Teko::JacobiPreconditionerFactory(strategy));
 
// Product of the two (MultPreconditionerFactory); AddPreconditionerFactory forms the sum
RCP<Teko::BlockPreconditionerFactory> MasterFactory =
    rcp(new Teko::MultPreconditionerFactory(GSFactory, JacobiFactory));

The equivalent parameter-driven form uses the "Block Add" / "Block Multiply" block types with "Preconditioner A" / "Preconditioner B".

Full source:

// @HEADER
// *****************************************************************************
//      Teko: A package for block and physics based preconditioning
//
// Copyright 2010 NTESS and the Teko contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER
 
#include "Teuchos_ConfigDefs.hpp"
#include "Teuchos_FancyOStream.hpp"
#include "Teuchos_RCP.hpp"
#include "Teuchos_XMLParameterListHelpers.hpp"
#include "Teuchos_DefaultComm.hpp"
#include "Teuchos_ParameterList.hpp"
#include "Teuchos_CommHelpers.hpp"
 
// Thyra includes
#include "Thyra_LinearOpBase.hpp"
#include "Thyra_TpetraLinearOp.hpp"
#include "Thyra_TpetraThyraWrappers.hpp"
 
// Tpetra includes
#include "Tpetra_Core.hpp"
#include "Tpetra_Map.hpp"
#include "Tpetra_CrsMatrix.hpp"
#include "Tpetra_Vector.hpp"
#include "Tpetra_MultiVector.hpp"
#include "MatrixMarket_Tpetra.hpp"
 
// Teko includes
#include "Teko_Utilities.hpp"
#include "Teko_InverseFactory.hpp"
#include "Teko_JacobiPreconditionerFactory.hpp"
#include "Teko_GaussSeidelPreconditionerFactory.hpp"
#include "Teko_BlockInvDiagonalStrategy.hpp"
#include "Teko_StridedTpetraOperator.hpp"
#include "Teko_TpetraBlockPreconditioner.hpp"
#include "Teko_AddPreconditionerFactory.hpp"
#include "Teko_MultPreconditionerFactory.hpp"
#include "Teko_ConfigDefs.hpp"
 
// Belos includes
#include "BelosConfigDefs.hpp"
#include "BelosLinearProblem.hpp"
#include "BelosBlockGmresSolMgr.hpp"
#include "BelosTpetraAdapter.hpp"
 
#include <iostream>
#include <vector>
 
using Teuchos::FancyOStream;
using Teuchos::ParameterList;
using Teuchos::RCP;
using Teuchos::rcp;
 
void run_driver() {
  RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
 
  using ST = double;
  using LO = Teko::LO;
  using GO = Teko::GO;
  using NT = Teko::NT;
 
  using map_t = Tpetra::Map<LO, GO, NT>;
  using crs_t = Tpetra::CrsMatrix<ST, LO, GO, NT>;
  using vec_t = Tpetra::Vector<ST, LO, GO, NT>;
  using mv_t  = Tpetra::MultiVector<ST, LO, GO, NT>;
  using op_t  = Tpetra::Operator<ST, LO, GO, NT>;
 
  auto comm = Tpetra::getDefaultComm();
 
  RCP<FancyOStream> fos = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
  fos->setOutputToRootOnly(0);
 
  const int numProc = comm->getSize();
  const int myPID   = comm->getRank();
 
  *out << "Approaching Barrier: proc = " << numProc << ", pid = " << myPID << std::endl;
  Teuchos::barrier(*comm);
 
  RCP<Teuchos::ParameterList> paramList = Teuchos::getParametersFromXmlFile("solverparams.xml");
 
  *fos << "Reading matrix market files" << std::endl;
 
  RCP<crs_t> A = Tpetra::MatrixMarket::Reader<crs_t>::readSparseFile("./modified.mm", comm);
 
  RCP<const map_t> rangeMap  = A->getRangeMap();
  RCP<const map_t> domainMap = A->getDomainMap();
 
  RCP<vec_t> b = Tpetra::MatrixMarket::Reader<crs_t>::readVectorFile("./rhs_test.mm", comm,
                                                                     rangeMap, false, false);
 
  RCP<vec_t> x = Tpetra::MatrixMarket::Reader<crs_t>::readVectorFile("./lhs_test.mm", comm,
                                                                     domainMap, false, false);
 
  *fos << "Building strided operator" << std::endl;
 
  std::vector<int> vars(2);
  vars[0] = 2;
  vars[1] = 1;
 
  Teuchos::RCP<Teko::TpetraHelpers::StridedTpetraOperator> sA =
      Teuchos::rcp(new Teko::TpetraHelpers::StridedTpetraOperator(vars, A));
 
  RCP<Teko::InverseFactory> inverse = Teko::invFactoryFromParamList(*paramList, "Amesos2");
 
  RCP<Teko::BlockInvDiagonalStrategy> strategy = rcp(new Teko::InvFactoryDiagStrategy(inverse));
 
  RCP<Teko::BlockPreconditionerFactory> GSFactory =
      rcp(new Teko::GaussSeidelPreconditionerFactory(Teko::GS_UseLowerTriangle, strategy));
 
  RCP<Teko::BlockPreconditionerFactory> JacobiFactory =
      rcp(new Teko::JacobiPreconditionerFactory(strategy));
 
#ifdef ADD_PREC
  RCP<Teko::BlockPreconditionerFactory> MasterFactory =
      rcp(new Teko::AddPreconditionerFactory(GSFactory, JacobiFactory));
#else
  RCP<Teko::BlockPreconditionerFactory> MasterFactory =
      rcp(new Teko::MultPreconditionerFactory(GSFactory, JacobiFactory));
#endif
 
  Teko::TpetraHelpers::TpetraBlockPreconditioner MyPreconditioner(MasterFactory);
  MyPreconditioner.buildPreconditioner(sA);
 
  RCP<mv_t> X = x;
  RCP<mv_t> B = rcp(new mv_t(b->getMap(), 1));
 
  using problem_t        = Belos::LinearProblem<ST, mv_t, op_t>;
  RCP<problem_t> problem = rcp(new problem_t(A, X, B));
 
  RCP<op_t> precOp = Teuchos::rcp(&MyPreconditioner, false);
  problem->setRightPrec(precOp);
 
  const bool set = problem->setProblem();
  TEUCHOS_TEST_FOR_EXCEPTION(!set, std::runtime_error,
                             "Belos::LinearProblem::setProblem() failed.");
 
  Teuchos::ParameterList belosList;
  belosList.set("Maximum Iterations", 100);
  belosList.set("Convergence Tolerance", 1e-5);
  belosList.set("Num Blocks", 50);
  belosList.set("Verbosity",
                Belos::Errors + Belos::Warnings + Belos::IterationDetails + Belos::FinalSummary);
  belosList.set("Output Frequency", 10);
 
  using solver_t = Belos::BlockGmresSolMgr<ST, mv_t, op_t>;
  solver_t solver(problem, rcpFromRef(belosList));
 
  *fos << "Starting Belos solve" << std::endl;
 
  Belos::ReturnType result = solver.solve();
 
  TEUCHOS_TEST_FOR_EXCEPTION(result != Belos::Converged, std::runtime_error,
                             "Belos solver failed to converge.");
 
  *fos << "Solve converged" << std::endl;
}
 
int main(int argc, char* argv[]) {
  Tpetra::ScopeGuard tpetraScope(&argc, &argv);
  { run_driver(); }
  return 0;
}

step1 — writing your own BlockPreconditionerFactory

examples/BuildPreconditioner/step1/ is the tutorial for extending Teko. To create a custom block preconditioner you subclass Teko::BlockPreconditionerFactory and override one method, buildPreconditionerOperator, using the block-operator utilities in Teko_Utilities.hpp:

Teko::LinearOp ExamplePreconditionerFactory::buildPreconditionerOperator(
    Teko::BlockedLinearOp& blockOp, Teko::BlockPreconditionerState& state) const {
  // 1. Pull out the sub-blocks
  const Teko::LinearOp A_00 = Teko::getBlock(0, 0, blockOp);
  const Teko::LinearOp A_01 = Teko::getBlock(0, 1, blockOp);
  const Teko::LinearOp A_10 = Teko::getBlock(1, 0, blockOp);
  const Teko::LinearOp A_11 = Teko::getBlock(1, 1, blockOp);
 
  // 2. Cheap inverse of the (1,1) block, and an explicit (0,0) approximation
  const Teko::LinearOp invH = Teko::getInvDiagonalOp(A_11);
  const Teko::LinearOp P    = Teko::explicitAdd(A_00, Teko::scale(alpha_, A_01));
  const Teko::LinearOp invP = Teko::buildInverse(*inverse_, P);  // uses the injected inverse
 
  // 3. Assemble a block lower-triangular inverse operator
  Teko::BlockedLinearOp L = Teko::zeroBlockedOp(blockOp);
  Teko::setBlock(1, 0, L, A_10);
  Teko::endBlockFill(L);
 
  std::vector<Teko::LinearOp> invDiag = {invP, invH};
  return Teko::createBlockLowerTriInverseOp(L, invDiag);
}

Key utilities on display (all in Teko_Utilities.hpp): blockRowCount/blockColCount, getBlock, getInvDiagonalOp, explicitAdd, scale, buildInverse, zeroBlockedOp / setBlock / endBlockFill, and createBlockLowerTriInverseOp. The companion example-test.cpp shows how to assemble four Tpetra sub-blocks, wrap them with Thyra::tpetraLinearOp, combine with Teko::block2x2, and apply the resulting preconditioner. See Advanced Topics for how to register such a factory so it can be selected by a "Type" string, and for the state/rebuild pattern that extends this example to cache inverse and explicit operators across nonlinear or time-step updates.

Full source (the custom factory):

// @HEADER
// *****************************************************************************
//      Teko: A package for block and physics based preconditioning
//
// Copyright 2010 NTESS and the Teko contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER
 
#include "Teko_BlockPreconditionerFactory.hpp"
#include "Teko_Utilities.hpp"
#include "Teko_InverseFactory.hpp"
#include "Teko_BlockLowerTriInverseOp.hpp"
#include "Teko_BlockUpperTriInverseOp.hpp"
 
using Teuchos::RCP;
 
// Declaration of the preconditioner factory
class ExamplePreconditionerFactory /*@ \label{lne1:begin-decl} @*/
    : public Teko::BlockPreconditionerFactory {
 public:
  // Constructor
  ExamplePreconditionerFactory(const RCP<const Teko::InverseFactory>& inverse, double alpha);
 
  // Function inherited from Teko::BlockPreconditionerFactory
  Teko::LinearOp buildPreconditionerOperator(Teko::BlockedLinearOp& blo,
                                             Teko::BlockPreconditionerState& state) const;
 
 protected:
  // class members
  RCP<const Teko::InverseFactory> inverse_;
  double alpha_;
 
}; /*@ \label{lne1:end-decl} @*/
 
// Constructor definition
ExamplePreconditionerFactory /*@ \label{lne1:begin-constructor} @*/
    ::ExamplePreconditionerFactory(const RCP<const Teko::InverseFactory>& inverse, double alpha)
    : inverse_(inverse), alpha_(alpha) {} /*@ \label{lne1:end-constructor} @*/
 
// Use the factory to build the preconditioner (this is where the work goes)
Teko::LinearOp ExamplePreconditionerFactory /*@ \label{lne1:begin-bpo} @*/
    ::buildPreconditionerOperator(Teko::BlockedLinearOp& blockOp,
                                  Teko::BlockPreconditionerState& state) const {
  int rows = Teko::blockRowCount(blockOp); /*@ \label{lne1:begin-extraction} @*/
  int cols = Teko::blockColCount(blockOp);
 
  TEUCHOS_ASSERT(rows == 2);  // sanity checks
  TEUCHOS_ASSERT(cols == 2);
 
  // extract subblocks
  const Teko::LinearOp A_00 = Teko::getBlock(0, 0, blockOp);
  const Teko::LinearOp A_01 = Teko::getBlock(0, 1, blockOp);
  const Teko::LinearOp A_10 = Teko::getBlock(1, 0, blockOp);
  const Teko::LinearOp A_11 = Teko::getBlock(1, 1, blockOp); /*@ \label{lne1:end-extraction} @*/
 
  // get inverse of diag(A11)
  const Teko::LinearOp invH = Teko::getInvDiagonalOp(A_11); /*@ \label{lne1:invH} @*/
 
  // build 0,0 block in the preconditioner
  const Teko::LinearOp P =
      Teko::explicitAdd(A_00, Teko::scale(alpha_, A_01)); /*@ \label{lne1:P} @*/
  const Teko::LinearOp invP =
      Teko::buildInverse(*inverse_, P);  // build inverse P /*@ \label{lne1:invP} @*/
 
  // build lower triangular inverse matrix
  Teko::BlockedLinearOp L = Teko::zeroBlockedOp(blockOp); /*@ \label{lne1:begin-trisolve} @*/
  Teko::setBlock(1, 0, L, A_10);
  Teko::endBlockFill(L);
 
  std::vector<Teko::LinearOp> invDiag(
      2);  // vector storing inverses /*@ \label{lne1:begin-invdiags} @*/
  invDiag[0] = invP;
  invDiag[1] = invH; /*@ \label{lne1:end-invdiags} @*/
 
  Teko::LinearOp invTildeA =
      Teko::createBlockLowerTriInverseOp(L, invDiag); /*@ \label{lne1:invLower} @*/
 
  // return fully constructed preconditioner
  return invTildeA; /*@ \label{lne1:end-trisolve} @*/
} /*@ \label{lne1:end-bpo} @*/

Full source (the driver that exercises it):

// @HEADER
// *****************************************************************************
//      Teko: A package for block and physics based preconditioning
//
// Copyright 2010 NTESS and the Teko contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER
 
// Teuchos includes /*@ \label{lnet:being-includes} @*/
#include "Teuchos_ConfigDefs.hpp"
#include "Teuchos_RCP.hpp"
 
// Thyra includes
#include "Thyra_TpetraLinearOp.hpp"
#include "Thyra_TpetraThyraWrappers.hpp"
 
// Tpetra includes
#include "Tpetra_Core.hpp"
#include "Tpetra_Map.hpp"
#include "Tpetra_CrsMatrix.hpp"
#include "Tpetra_Vector.hpp"
#include "Tpetra_MultiVector.hpp"
 
// Teko-Package includes
#include "Teko_Utilities.hpp"
#include "Teko_InverseFactory.hpp"
#include "Teko_InverseLibrary.hpp"
#include "Teko_TpetraOperatorWrapper.hpp"
#include "Teko_TpetraBlockPreconditioner.hpp"
#include "Teko_ConfigDefs.hpp"
 
#include "ExamplePreconditionerFactory.cpp"
 
#include <iostream>
#include <vector> /*@ \label{lnet:end-includes} @*/
 
// for simplicity
using Teuchos::RCP;
using Teuchos::rcp;
 
// utility function to construct Tpetra operators
RCP<Tpetra::CrsMatrix<Teko::ST, Teko::LO, Teko::GO, Teko::NT> > build2x2(
    double a, double b, double c, double d, const RCP<const Teuchos::Comm<int> >& comm) {
  using ST = Teko::ST;
  using LO = Teko::LO;
  using GO = Teko::GO;
  using NT = Teko::NT;
 
  using map_t = Tpetra::Map<LO, GO, NT>;
  using crs_t = Tpetra::CrsMatrix<ST, LO, GO, NT>;
 
  RCP<const map_t> map = rcp(new map_t(2, 0, comm));
  RCP<crs_t> matrix    = rcp(new crs_t(map, 2));
 
  Teuchos::Array<GO> indices(2);
  Teuchos::Array<ST> values(2);
 
  indices[0] = 0;
  indices[1] = 1;
 
  // build first row
  if (map->isNodeGlobalElement(0)) {
    values[0] = a;
    values[1] = b;
    matrix->insertGlobalValues(0, indices(), values());
  }
 
  // build second row
  if (map->isNodeGlobalElement(1)) {
    values[0] = c;
    values[1] = d;
    matrix->insertGlobalValues(1, indices(), values());
  }
 
  matrix->fillComplete();
 
  return matrix;
}
 
void run_test() {
  using ST = Teko::ST;
  using LO = Teko::LO;
  using GO = Teko::GO;
  using NT = Teko::NT;
 
  using vec_t = Tpetra::Vector<ST, LO, GO, NT>;
  using mv_t  = Tpetra::MultiVector<ST, LO, GO, NT>;
 
  auto comm = Tpetra::getDefaultComm();
 
  // build the sub blocks
  auto mat_00         = build2x2(1, 2, 2, 1, comm);
  Teko::LinearOp A_00 = Thyra::tpetraLinearOp<ST, LO, GO, NT>(
      Thyra::tpetraVectorSpace<ST, LO, GO, NT>(mat_00->getRangeMap()),
      Thyra::tpetraVectorSpace<ST, LO, GO, NT>(mat_00->getDomainMap()), mat_00);
 
  auto mat_01         = build2x2(0, -1, 3, 4, comm);
  Teko::LinearOp A_01 = Thyra::tpetraLinearOp<ST, LO, GO, NT>(
      Thyra::tpetraVectorSpace<ST, LO, GO, NT>(mat_01->getRangeMap()),
      Thyra::tpetraVectorSpace<ST, LO, GO, NT>(mat_01->getDomainMap()), mat_01);
 
  auto mat_10         = build2x2(1, 6, -3, 2, comm);
  Teko::LinearOp A_10 = Thyra::tpetraLinearOp<ST, LO, GO, NT>(
      Thyra::tpetraVectorSpace<ST, LO, GO, NT>(mat_10->getRangeMap()),
      Thyra::tpetraVectorSpace<ST, LO, GO, NT>(mat_10->getDomainMap()), mat_10);
 
  auto mat_11         = build2x2(2, 1, 1, 2, comm);
  Teko::LinearOp A_11 = Thyra::tpetraLinearOp<ST, LO, GO, NT>(
      Thyra::tpetraVectorSpace<ST, LO, GO, NT>(mat_11->getRangeMap()),
      Thyra::tpetraVectorSpace<ST, LO, GO, NT>(mat_11->getDomainMap()), mat_11);
 
  // build the Tpetra operator wrapper
  Teuchos::RCP<Teko::TpetraHelpers::TpetraOperatorWrapper> A = Teuchos::rcp(
      new Teko::TpetraHelpers::TpetraOperatorWrapper(Teko::block2x2(A_00, A_01, A_10, A_11)));
 
  // build the Tpetra vectors
  RCP<vec_t> b = rcp(new vec_t(A->getRangeMap()));
  RCP<vec_t> x = rcp(new vec_t(A->getDomainMap()));
 
  // build the RHS vector
  b->replaceGlobalValue(0, 1.0);
  b->replaceGlobalValue(1, 2.0);
  b->replaceGlobalValue(2, 3.0);
  b->replaceGlobalValue(3, 4.0);
 
  // Build the preconditioner
 
  RCP<Teko::InverseLibrary> invLib = Teko::InverseLibrary::buildFromStratimikos();
 
  RCP<const Teko::InverseFactory> inverse = invLib->getInverseFactory("Amesos2");
 
  RCP<Teko::BlockPreconditionerFactory> precFact =
      rcp(new ExamplePreconditionerFactory(inverse, 0.9));
 
  Teko::TpetraHelpers::TpetraBlockPreconditioner prec(precFact);
  prec.buildPreconditioner(A);
 
  // apply the preconditioner
  RCP<mv_t> B = b;
  RCP<mv_t> X = x;
 
  prec.apply(*B, *X);
 
  x->describe(*(Teuchos::VerboseObjectBase::getDefaultOStream()),
              Teuchos::EVerbosityLevel::VERB_EXTREME);
}
 
int main(int argc, char* argv[]) {
  Tpetra::ScopeGuard tpetraScope(&argc, &argv);
  { run_test(); }
  return 0;
}

Table of Contents

BuildPreconditioner — the canonical driver

StridedSolve — SIMPLE on a strided operator, solver chosen at runtime

Arbitrary Tpetra blocking — field blocks from monolithic GIDs

AddMultiplyPrecs — composing preconditioners

step1 — writing your own BlockPreconditionerFactory