[Trilinos-Users] Example program for solving a sparse linear problem Ax=y with Epetra, Ifpack, Belos - asking for code review
Joachim Wuttke
j.wuttke at fz-juelich.de
Wed Jun 21 16:25:38 EDT 2017
In the extant documentation, I found no simple, self-contained, complete code example
that shows how to solve a sparse linear problem Ax=y using Trilinos.
With considerable effort, and with substantial help from this mailing list, I finally
succeeded to write an example program of my own. Perhaps it could be useful for
other newcomers. I would be delighted if you see a place for it in your tutorial.
First, however, I would like to ask you to review my code for correctness, and for
making best use of Trilinos.
Also, my numeric example seems to be badly chosen, because it is ill conditioned.
What could be a better, but equally simple example?
- Joachim
//==================================================================================================
//
// belos_sparse_solver.cpp
//
// Example program that shows how to solve a sparse linear problem A x = y
// using the libraries Epetra, Ifpack, and Belos from Trilinos.
//
// (C) Joachim Wuttke, Forschungszentrum Jülich GmbH, 2017
// using code fragments from Epetra, Ifpack, and Belos tutorials,
// written with considerable help from the trilinos-users mailing list.
// Licence: CC-BY-SA
//
// Compilation:
// mpic++ -I/usr/include/trilinos -ltrilinos_belos -ltrilinos_teuchoscomm
// -ltrilinos_teuchoscore -ltrilinos_teuchosnumerics -ltrilinos_teuchosparameterlist
// -ltrilinos_epetra -ltrilinos_belosepetra -ltrilinos_ifpack -o belos_sparse_solver.cpp
// Execution:
// mpirun -n 4 a.out
// Alternatively without MPI:
// Compile with g++ instead of mpic++, and simply run a.out without mpirun wrapper
//
//==================================================================================================
#include "Epetra_CrsMatrix.h"
#include "Epetra_Vector.h"
#include "Epetra_MultiVector.h"
#include "Epetra_LinearProblem.h"
#include "Epetra_Export.h"
#include "Teuchos_ParameterList.hpp"
#include "Teuchos_RCP.hpp"
#include "Ifpack.h"
#include "Ifpack_AdditiveSchwarz.h"
#include "BelosLinearProblem.hpp"
#include "BelosBlockGmresSolMgr.hpp"
#include "BelosEpetraAdapter.hpp"
#ifdef HAVE_MPI
#include "Epetra_MpiComm.h"
#else
#include "Epetra_SerialComm.h"
#endif
#include <vector>
//==================================================================================================
// Parameterization
//==================================================================================================
Teuchos::ParameterList my_preconditioner_parameters()
{
Teuchos::ParameterList ret;
ret.set("fact: drop tolerance", 1e-9);
ret.set("fact: level-of-fill", 1);
ret.set("schwarz: combine mode", "Add"); //explained in Epetra_CombineMode.h
return ret;
}
Teuchos::ParameterList my_solver_parameters()
{
Teuchos::ParameterList ret;
ret.set("Block Size", 1); // Blocksize to be used by iterative solver
ret.set("Maximum Iterations", 1550); // Maximum number of iterations allowed
ret.set("Convergence Tolerance", 1e-12);// Relative convergence tolerance requested
using namespace Belos;
ret.set ("Verbosity", Errors+Warnings+TimingDetails+FinalSummary );
return ret;
}
//==================================================================================================
// Set A and y of exemplary linear problem
//==================================================================================================
// We choose a linear problem for which we know the solution x = (0, 1, ..., N-1).
//! Sets row number jj of the sparse NxN matrix A.
void my_row(const int jj, const int N, std::vector<int>& indices, std::vector<double>& values)
{
// Specifically, in this example, A is a tridiagonal band matrix (-1, 2, 1).
if (jj>0) {
values.push_back(-1.);
indices.push_back(jj-1);
}
values.push_back(2.);
indices.push_back(jj);
if (jj<N-1) {
values.push_back(-1.);
indices.push_back(jj+1);
}
}
//! Entry number jj of the rhs vector y = (-1, 0, ..., 0, N).
double my_rhs(const int jj, const int N)
{
if (jj==0)
return -1;
else if (jj==N-1)
return N;
else
return 0;
}
//! Entry number jj of expected solution x = (0, 1, ..., N-1).
double expected_solution(const int jj)
{
return jj;
}
//==================================================================================================
// Auxiliary
//==================================================================================================
void print_solution(const Epetra_Vector& Solution, const int N)
{
std::cout << "Solution:";
for (int i=0; i<N; ++i) {
if (i>10 && i<N-10)
continue;
std::cout << " " << Solution[i];
}
std::cout << "\n";
std::cout << "Error:";
for (int i=0; i<N; ++i) {
if (i>10 && i<N-10)
continue;
std::cout << " " << Solution[i]-expected_solution(i);
}
std::cout << "\n";
double max_error = 0;
for (int i=0; i<N; ++i)
max_error = std::max(std::abs(Solution[i]-expected_solution(i)), max_error);
std::cout << "Max error: " << max_error << "\n";
}
//==================================================================================================
// Main program
//==================================================================================================
//! Solve the linear problem A x = y, with A, y as given above. Print the solution x.
int main(int argc, char *argv[])
{
using Teuchos::RCP;
typedef Epetra_Vector SV;
typedef Epetra_MultiVector MV;
typedef Epetra_Operator OP;
#ifdef EPETRA_NO_32BIT_GLOBAL_INDICES
std::cerr << "Unexpectedly, global indices are not 32 bit\n";
exit(1);
#endif
#ifdef HAVE_MPI
MPI_Init(&argc,&argv);
Epetra_MpiComm Comm( MPI_COMM_WORLD );
#else
Epetra_SerialComm Comm;
#endif
// The problem is defined on a 2D grid, global size is N * N.
int N = 120000;
Epetra_Map map(N, 0 , Comm);
int myNumRows = map.NumMyElements(); // approximately N divided by number of processes
int* myRows = nullptr;
myRows = map.MyGlobalElements();
// Set the sparse matrix A.
RCP<Epetra_CrsMatrix> A { new Epetra_CrsMatrix(Copy, map, 0) };
for (int j=0; j<myNumRows; ++j) {
std::vector<int> indices;
std::vector<double> values;
my_row(myRows[j], N, indices, values);
int err = 0;
err |= A->IndicesAreLocal();
err |= A->IndicesAreContiguous();
err |= (A->InsertGlobalValues(myRows[j], indices.size(), &values[0], &indices[0]) < 0);
TEUCHOS_TEST_FOR_EXCEPTION(err, std::runtime_error, "filling A failed");
}
A->FillComplete();
// Set the right-hand side vector y of the linear problem A x = y.
RCP<SV> RHS { new SV(A->OperatorDomainMap(), 1) };
for (int j=0; j<myNumRows; ++j)
(*RHS)[j] = my_rhs(myRows[j], N);
// Construct the preconditioner.
int err = 0;
RCP<Ifpack_Preconditioner> Prec { Ifpack().Create("ILU", &*A, 1) };
err |= Prec == Teuchos::null;
Teuchos::ParameterList preconditionerParameters { my_preconditioner_parameters() };
err |= Prec->SetParameters(preconditionerParameters);
err |= Prec->Initialize();
err |= Prec->Compute();
// Setup the linear problem, consisting of the numeric input and preconditioner.
RCP<SV> LHS { new SV(A->OperatorDomainMap(), 1) };
RCP<Belos::LinearProblem<double,MV,OP>> problem {
new Belos::LinearProblem<double,MV,OP>(A, LHS, RHS) };
RCP<Belos::EpetraPrecOp> PrecOp { new Belos::EpetraPrecOp(Prec) };
problem->setRightPrec(PrecOp);
err |= !problem->setProblem();
TEUCHOS_TEST_FOR_EXCEPTION(err, std::runtime_error, "Some error in Belos construction");
// Iterative solve.
RCP<Teuchos::ParameterList> solverParams{ new Teuchos::ParameterList(my_solver_parameters()) };
Belos::BlockGmresSolMgr<double,MV,OP> belosSolver(problem, solverParams);
Belos::ReturnType ret = belosSolver.solve(); // is enum { Belos::Converged, Belos::Unconverged }
// Collect solution from all processors into processor 0.
Epetra_Map TargetMap(-1, Comm.MyPID()==0 ? N : 0, 0, Comm);
Epetra_Export Exporter(A->OperatorDomainMap(), TargetMap);
const Teuchos::RCP<MV> pLHS { problem->getLHS() };
Epetra_Vector Solution(TargetMap);
Solution.Export(*pLHS, Exporter, Add);
// Print result.
if (Comm.MyPID() == 0) {
if (ret == Belos::Converged)
std::cout << "Belos converged.\n";
else
std::cout << "Belos did not converge.\n";
std::cout << "\nPreconditioner:" << *Prec << "/Preconditioner\n";
print_solution(Solution, N);
}
#ifdef HAVE_MPI
MPI_Finalize() ;
#endif
return 0;
}
-------------- next part --------------
A non-text attachment was scrubbed...
Name: smime.p7s
Type: application/pkcs7-signature
Size: 5110 bytes
Desc: S/MIME Cryptographic Signature
URL: <https://trilinos.org/pipermail/trilinos-users/attachments/20170621/1e7f5cc1/attachment-0001.p7s>
More information about the Trilinos-Users
mailing list