[Trilinos-Users] [EXTERNAL] GMRES inconsistent when using multiple threads (Belos+Tpetra, since 12.12.1)
Christopher Thiele
ct37 at rice.edu
Thu Nov 9 15:34:37 EST 2017
Judging by the output of Teuchos::TimeMonitor::summarize() at the end,
it seems to be using PseudoBlockGmresSolMgr. I just used
Belos::SolverFactory with the "GMRES" argument, so this should be the
default. I attached the code below.
Christopher
// solver_driver.cc
// Import linear systems from matrix market files and solve them with
Trilinos.
#include <string>
#include <iostream>
// Trilinos headers
#include <Teuchos_TimeMonitor.hpp>
#include <Teuchos_DefaultComm.hpp>
#include <Tpetra_Vector.hpp>
#include <Tpetra_CrsMatrix.hpp>
#include <BelosTpetraAdapter.hpp>
#include <BelosSolverManager.hpp>
#include <BelosSolverFactory.hpp>
#include <Ifpack2_Preconditioner.hpp>
#include <Ifpack2_Factory.hpp>
#include <KokkosCompat_ClassicNodeAPI_Wrapper.hpp>
#include <MatrixMarket_Tpetra.hpp>
// Type definitions for the Trilinos classes that we will use.
typedef Kokkos::Compat::KokkosOpenMPWrapperNode Node;
typedef Tpetra::Operator<double, int, int, Node> Operator;
typedef Tpetra::Map<int, int, Node> Map;
typedef Tpetra::CrsMatrix<double, int, int, Node> Matrix;
typedef Tpetra::Vector<double, int, int, Node> Vector;
typedef Tpetra::MultiVector<double, int, int, Node> MultiVector;
int main(int argc, char **argv)
{
// Initialize MPI.
Teuchos::GlobalMPISession mpi_session(&argc, &argv, NULL);
// We only initialized MPI for good measure and to reproduce the
environment
// of an actual Trilinos-based solver.
// However, running the driver with more than one process is rather
pointless,
// because the proper distribution of matrices and vectors over MPI ranks
// cannot be restored from matrix market files.
if (mpi_session.getNProc() > 1) {
if (mpi_session.getRank() == 0) {
std::cerr << "Error: Do not run this driver with more than one
MPI rank!"
<< std::endl;
}
return 1;
}
// Define command line arguments and default values.
int threads = 1;
std::string solver_name = "GMRES";
std::string preconditioner_name = "Gauss-Seidel";
bool use_preconditioner = true;
int max_iter = 9999;
int gmres_restart_length = 30;
double tol = 1.0e-6;
std::string matrix_file = "";
std::string rhs_file = "";
std::string initial_guess_file = "";
Teuchos::CommandLineProcessor parser(false);
parser.setDocString(
"Known solver names that work are: CG, GMRES, BiCGStab\n"
"Known preconditioner names that work are: "
"Jacobi, Gauss-Seidel, Symmetric Gauss-Seidel\n\n"
"The multi-threaded Gauss-Seidel preconditioners \"MT Gauss-Seidel\"\n"
"and \"MT Symmetric Gauss-Seidel\" will only work if they have been\n"
"enabled in the Trilinos installation.\n\n"
"You will have to provide at least the number of threads, a matrix,
and a\n"
"right-hand side vector. All other arguments are optional.\n\n"
"Example: mpirun -n 1 ./solver_driver --threads=4 --matrix=A.mm
--rhs=b.mm"
);
parser.setOption("threads", &threads, "Number of threads", true);
parser.setOption("solver", &solver_name, "Name of the linear solver",
false);
parser.setOption("preconditioner",
&preconditioner_name,
"Name of the preconditioner",
false);
parser.setOption("enable-preconditioner",
"disable-preconditioner",
&use_preconditioner,
"Enable or disable the preconditioner");
parser.setOption("max_iter",
&max_iter,
"Maximum linear solver iterations",
false);
parser.setOption("gmres_restart_length",
&gmres_restart_length,
"GMRES restart length",
false);
parser.setOption("tol", &tol, "Relative tolerance", false);
parser.setOption("matrix",
&matrix_file,
"The matrix market file for the system matrix",
true);
parser.setOption("rhs",
&rhs_file,
"The matrix market file for the right-hand side vector",
true);
parser.setOption("initial_guess",
&initial_guess_file,
"The matrix market file for the initial guess "
"(default: zero vector)",
false);
// Parse command line arguments.
if (parser.parse(argc, argv)
!= Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
std::cerr << "Error: Could not parse command line arguments!"
<< std::endl
<< "Run with --help for instructions."
<< std::endl;
return 1;
}
// Initialize Kokkos (for OpenMP).
Kokkos::OpenMP::initialize(threads, 1);
// Import linear system.
Teuchos::RCP<const Teuchos::Comm<int> > comm
= Teuchos::DefaultComm<int>::getComm();
std::cout << "Loading matrix: " << matrix_file << std::endl;
Teuchos::RCP<const Matrix> matrix
= Tpetra::MatrixMarket::Reader<Matrix>::readSparseFile(
matrix_file,
comm
);
Teuchos::RCP<const Map> domain_map = matrix->getDomainMap();
Teuchos::RCP<const Map> range_map = matrix->getRangeMap();
std::cout << "Loading right-hand side vector: " << rhs_file << std::endl;
Teuchos::RCP<const Vector> rhs
= Tpetra::MatrixMarket::Reader<Matrix>::readVectorFile(
rhs_file,
comm,
range_map
);
Teuchos::RCP<Vector> initial_guess;
// Load initial guess if provided, otherwise use zero vector.
if (initial_guess_file != "") {
std::cout << "Loading initial guess: " << initial_guess_file <<
std::endl;
initial_guess = Tpetra::MatrixMarket::Reader<Matrix>::readVectorFile(
initial_guess_file,
comm,
domain_map
);
}
else {
std::cout << "Using zero vector as initial guess." << std::endl;
initial_guess = Teuchos::rcp(new Vector(domain_map));
}
// Create the linear solver and preconditioner.
std::cout << "Creating linear solver and preconditioner." << std::endl;
Teuchos::ParameterList solver_params;
Teuchos::ParameterList preconditioner_params;
if (solver_name == "GMRES") {
solver_params.set("Num Blocks", gmres_restart_length);
}
// Make sure we limit the solver by iterations rather than restarts.
solver_params.set("Verbosity", Belos::MsgType::Errors
| Belos::MsgType::IterationDetails
| Belos::MsgType::FinalSummary
| Belos::MsgType::StatusTestDetails);
solver_params.set("Maximum Restarts", 99999);
solver_params.set("Maximum Iterations", max_iter);
solver_params.set("Convergence Tolerance", tol);
solver_params.set("Explicit Residual Scaling", "Norm of Initial
Residual");
preconditioner_params.set("relaxation: type", preconditioner_name);
preconditioner_params.set("relaxation: sweeps", 1);
Belos::SolverFactory<double, MultiVector, Operator> solver_factory;
Ifpack2::Factory preconditioner_factory;
Teuchos::RCP<Belos::SolverManager<double, MultiVector, Operator> >
solver;
Teuchos::RCP<Ifpack2::Preconditioner<double, int, int> > preconditioner;
solver = solver_factory.create(solver_name,
Teuchos::rcpFromRef(solver_params));
if (use_preconditioner) {
preconditioner = preconditioner_factory.create("RELAXATION", matrix);
preconditioner->setParameters(preconditioner_params);
preconditioner->compute();
}
// Solve the linear system.
std::cout << "Solving linear system." << std::endl;
Belos::LinearProblem<double, MultiVector, Operator> problem(matrix,
initial_guess,
rhs);
if (use_preconditioner) problem.setRightPrec(preconditioner);
problem.setProblem();
solver->setProblem(Teuchos::rcpFromRef(problem));
const Belos::ReturnType status = solver->solve();
if (status == Belos::ReturnType::Converged) {
std::cout << "Linear solver converged." << std::endl;
}
else {
std::cerr << "Error: Linear solver did not converge!" << std::endl;
}
// Print time measurements.
Teuchos::TimeMonitor::summarize();
return 0;
}
On 11/09/2017 01:52 PM, Thornquist, Heidi K wrote:
> Hi Christopher,
>
> That is odd behavior. Does your driver use the BelosBlockGmresSolMgr or BelosPseudoBlockGmresSolMgr? If you don’t mind sharing your solver driver, that would be helpful too!
>
> Many thanks!
>
> Heidi
>
>
> On 11/9/17, 11:28 AM, "Trilinos-Users on behalf of Christopher Thiele" <trilinos-users-bounces at trilinos.org on behalf of ct37 at rice.edu> wrote:
>
> Hello,
>
> I noticed some strange behavior in Belos' GMRES solver when using Tpetra
> and multiple threads. I wrote a little tool (solver_driver) that just
> reads in a matrix market file and solves a linear system to reproduce
> it. Here are some results with the rma10 matrix from the University of
> Florida collection (the choice is rather arbitrary):
>
> ./solver_driver_12.10.1 --threads=2 --matrix=rma10.mtx --rhs=rma10_b.mtx
> --tol=1.0e-1 --solver=GMRES --disable-preconditioner
> ---> Converges in just 3 iterations, independent of the number of
> threads etc. This seems reasonable given the tolerance of 1e-1.
>
> ./solver_driver_12.12.1 --threads=2 --matrix=rma10.mtx --rhs=rma10_b.mtx
> --tol=1.0e-1 --solver=GMRES --disable-preconditioner
> ---> Converges in either 61 or 62 iterations, but the final residual
> changes in the order of 1e-4 with each run.
>
> OMP_PROC_BIND=true ./solver_driver_12.12.1 --threads=2
> --matrix=rma10.mtx --rhs=rma10_b.mtx --tol=1.0e-1 --solver=GMRES
> --disable-preconditioner
> ---> Converges in 121 iterations, and the final residual changes in the
> order of 1e-3.
>
> mpirun -n 1 ./solver_driver_12.12.1 --threads=2 --matrix=rma10.mtx
> --rhs=rma10_b.mtx --tol=1.0e-1 --solver=GMRES --disable-preconditioner
> ---> Seems to solve the problem, consistently reproduces the result
> obtained with 12.10.1
>
> I also tested the current master branch from GitHub, and it shows the
> same behavior as the 12.12.1 release. The problem does not occur with
> other solvers like BiCGStab or CG. I added some information about my
> system and the Trilinos configuration below. I can also share the code
> (~200 loc), but it is really just a driver for the solvers.
>
> Best regards,
> Christopher
>
>
>
> System:
>
> OS: RHEL 7 (x84_64)
> Compiler: GCC 4.8.5
> MPI: Open MPI 1.10.6
> CPU: Xeon E3-1220 V2 (4 cores)
> Memory: 8 GB
>
>
>
> Trilinos configuration script:
>
> #!/bin/bash
>
> set -e
>
> SOURCE_BASE=/path/to/source
> INSTALL_DIR=/path/to/installation
>
> TRILINOS_PKG_CONF=" \
> -DTrilinos_ENABLE_Teuchos=ON \
> -DTrilinos_ENABLE_TeuchosParser=ON \
> -DTrilinos_ENABLE_Belos=ON \
> -DTrilinos_ENABLE_Ifpack2=ON \
> -DTrilinos_ENABLE_MueLu=OFF \
> -DTrilinos_ENABLE_Tpetra=ON \
> -DTrilinos_ENABLE_Teko=OFF \
> -DTrilinos_ENABLE_Zoltan2=OFF \
> -DTrilinos_ASSERT_MISSING_PACKAGES=OFF \
> "
>
> TP_PKG_CONF=" \
> -DTPL_ENABLE_MPI=ON \
> -DTPL_ENABLE_CUDA=OFF \
> -DTPL_ENABLE_MKL=ON \
> -DTPL_MKL_LIBRARIES=/path/to/mkl/lib/intel64/libmkl_rt.so \
> -DTPL_MKL_INCLUDE_DIRS=/path/to/mkl/include \
> "
>
> cmake \
> -D MPI_C_COMPILER=mpicc \
> -D MPI_CXX_COMPILER=mpicxx \
> -D MPI_Fortran_COMPILER=mpifort \
> -D CMAKE_C_COMPILER=mpicc \
> -D CMAKE_CXX_COMPILER=mpicxx \
> -D CMAKE_Fortran_COMPILER=mpifort \
> -D CMAKE_CXX_FLAGS="-g -O3 -march=native" \
> -D CMAKE_BUILD_TYPE=RELEASE \
> -D Trilinos_ENABLE_TESTS=OFF \
> -D Trilinos_ENABLE_EXPLICIT_INSTANTIATION=ON \
> -D Trilinos_CXX11_FLAGS="-std=c++11" \
> -D Kokkos_ENABLE_Pthread=OFF \
> -D Kokkos_ENABLE_OpenMP=ON \
> -D Kokkos_ENABLE_Cuda=OFF \
> -D Kokkos_ENABLE_Cuda_UVM=OFF \
> -D Trilinos_ENABLE_OpenMP=ON \
> -D Tpetra_INST_SERIAL=ON \
> -D Tpetra_INST_CUDA=OFF \
> -D Tpetra_INST_OPENMP=ON \
> -D Tpetra_INST_INT_INT=ON \
> -D Tpetra_INST_INT_LONG=OFF \
> -D Tpetra_INST_COMPLEX_DOUBLE=OFF \
> -D Teuchos_ENABLE_LONG_LONG_INT=OFF \
> -D CMAKE_INSTALL_PREFIX=${INSTALL_DIR} \
> ${TRILINOS_PKG_CONF} \
> ${TP_PKG_CONF} \
> ${SOURCE_BASE}
> _______________________________________________
> Trilinos-Users mailing list
> Trilinos-Users at trilinos.org
> https://trilinos.org/mailman/listinfo/trilinos-users
>
>
More information about the Trilinos-Users
mailing list