The provided example uses PCPGSolMgr with an ML preconditioner.
The provided example uses PCPGSolMgr with an ML preconditioner.
#include "ml_include.h"
#include "Epetra_Map.h"
#ifdef EPETRA_MPI
#include "Epetra_MpiComm.h"
#else
#include "Epetra_SerialComm.h"
#endif
#include "Epetra_Vector.h"
#include "Epetra_CrsMatrix.h"
#include "EpetraExt_MatrixMatrix.h"
#include "ml_MultiLevelPreconditioner.h"
#include "BelosEpetraAdapter.hpp"
#include "Teuchos_CommandLineProcessor.hpp"
#include "Teuchos_ParameterList.hpp"
#include "Teuchos_StandardCatchMacros.hpp"
int main(int argc, char *argv[]) {
int MyPID = 0;
int numProc = 1;
#ifdef EPETRA_MPI
MPI_Init(&argc,&argv);
Epetra_MpiComm Comm(MPI_COMM_WORLD);
MyPID = Comm.MyPID();
numProc = Comm.NumProc();
#else
Epetra_SerialComm Comm;
#endif
typedef double ST;
typedef Teuchos::ScalarTraits<ST> SCT;
typedef SCT::magnitudeType MT;
typedef Epetra_MultiVector MV;
typedef Epetra_Operator OP;
using Teuchos::ParameterList;
using Teuchos::RCP;
using Teuchos::rcp;
bool verbose = false;
bool success = true;
try {
bool proc_verbose = false;
int frequency = -1;
int blocksize = 1;
int numrhs = 1;
int maxiters = 30;
int maxDeflate = 8;
int maxSave = 16;
std::string ortho("ICGS");
MT tol = 1.0e-8;
Teuchos::CommandLineProcessor cmdp(false,true);
cmdp.setOption("verbose","quiet",&verbose,"Print messages and results");
cmdp.setOption("frequency",&frequency,"Solvers frequency for printing residuals (#iters)");
cmdp.setOption("tol",&tol,"Relative residual tolerance used by PCPG solver");
cmdp.setOption("num-rhs",&numrhs,"Number of right-hand sides to be solved for");
cmdp.setOption("max-iters",&maxiters,"Maximum number of iterations per linear system (-1 = adapted to problem/block size)");
cmdp.setOption("num-deflate",&maxDeflate,"Number of vectors deflated from the linear system");
cmdp.setOption("num-save",&maxSave,"Number of vectors saved from old Krylov subspaces");
cmdp.setOption("ortho-type",&ortho,"Orthogonalization type, either DGKS, ICGS or IMGS");
if (cmdp.parse(argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
return -1;
}
if (!verbose)
frequency = -1;
int numElePerDirection = 14*numProc;
int num_time_step = 4;
int numNodes = (numElePerDirection - 1)*(numElePerDirection - 1);
RCP<Epetra_Map> Map = rcp(new Epetra_Map(numNodes, 0, Comm) );
RCP<Epetra_CrsMatrix> Stiff = rcp(new Epetra_CrsMatrix(Epetra_DataAccess::Copy, *Map, 0) );
RCP<Epetra_CrsMatrix> Mass = rcp(new Epetra_CrsMatrix(Epetra_DataAccess::Copy, *Map, 0) );
RCP<Epetra_Vector> vecLHS = rcp( new Epetra_Vector(*Map) );
RCP<Epetra_Vector> vecRHS = rcp( new Epetra_Vector(*Map) );
RCP<Epetra_MultiVector> LHS, RHS;
double ko = 8.0/3.0, k1 = -1.0/3.0;
double h = 1.0/(double) numElePerDirection;
double mo = h*h*4.0/9.0, m1 = h*h/9.0, m2 = h*h/36.0;
double pi = 4.0*atan(1.0), valueLHS;
int lid, node, pos, iX, iY;
for(lid = Map->MinLID(); lid <= Map->MaxLID(); lid++){
node = Map->GID(lid);
iX = node % (numElePerDirection-1);
iY = ( node - iX )/(numElePerDirection-1);
pos = node;
Stiff->InsertGlobalValues(node, 1, &ko, &pos);
Mass->InsertGlobalValues(node, 1, &mo, &pos);
valueLHS = sin( pi*h*((double) iX+1) )*cos( 2.0 * pi*h*((double) iY+1) );
vecLHS->ReplaceGlobalValues( 1, &valueLHS, &node);
if (iY > 0) {
pos = iX + (iY-1)*(numElePerDirection-1);
Stiff->InsertGlobalValues(node, 1, &k1, &pos);
Mass->InsertGlobalValues(node, 1, &m1, &pos);
}
if (iY < numElePerDirection-2) {
pos = iX + (iY+1)*(numElePerDirection-1);
Stiff->InsertGlobalValues(node, 1, &k1, &pos);
Mass->InsertGlobalValues(node, 1, &m1, &pos);
}
if (iX > 0) {
pos = iX-1 + iY*(numElePerDirection-1);
Stiff->InsertGlobalValues(node, 1, &k1, &pos);
Mass->InsertGlobalValues(node, 1, &m1, &pos);
if (iY > 0) {
pos = iX-1 + (iY-1)*(numElePerDirection-1);
Stiff->InsertGlobalValues(node, 1, &k1, &pos);
Mass->InsertGlobalValues(node, 1, &m2, &pos);
}
if (iY < numElePerDirection-2) {
pos = iX-1 + (iY+1)*(numElePerDirection-1);
Stiff->InsertGlobalValues(node, 1, &k1, &pos);
Mass->InsertGlobalValues(node, 1, &m2, &pos);
}
}
if (iX < numElePerDirection - 2) {
pos = iX+1 + iY*(numElePerDirection-1);
Stiff->InsertGlobalValues(node, 1, &k1, &pos);
Mass->InsertGlobalValues(node, 1, &m1, &pos);
if (iY > 0) {
pos = iX+1 + (iY-1)*(numElePerDirection-1);
Stiff->InsertGlobalValues(node, 1, &k1, &pos);
Mass->InsertGlobalValues(node, 1, &m2, &pos);
}
if (iY < numElePerDirection-2) {
pos = iX+1 + (iY+1)*(numElePerDirection-1);
Stiff->InsertGlobalValues(node, 1, &k1, &pos);
Mass->InsertGlobalValues(node, 1, &m2, &pos);
}
}
}
Stiff->FillComplete();
Stiff->OptimizeStorage();
Mass->FillComplete();
Mass->OptimizeStorage();
double one = 1.0, hdt = .00005;
RCP<Epetra_CrsMatrix> A = rcp(new Epetra_CrsMatrix(*Stiff) );
int err = EpetraExt::MatrixMatrix::Add(*Mass, false, one,*A,hdt);
if (err != 0) {
std::cout << "err "<<err<<" from MatrixMatrix::Add "<<std::endl;
return(err);
}
A->FillComplete();
A->OptimizeStorage();
hdt = -hdt;
RCP<Epetra_CrsMatrix> B = rcp(new Epetra_CrsMatrix(*Stiff) );
EpetraExt::MatrixMatrix::Add(*Mass, false, one,*B,hdt);
if (err != 0) {
std::cout << "err "<<err<<" from MatrixMatrix::Add "<<std::endl;
return(err);
}
B->FillComplete();
B->OptimizeStorage();
B->Multiply(false, *vecLHS, *vecRHS);
proc_verbose = verbose && (MyPID==0);
LHS = Teuchos::rcp_implicit_cast<Epetra_MultiVector>(vecLHS);
RHS = Teuchos::rcp_implicit_cast<Epetra_MultiVector>(vecRHS);
ParameterList MLList;
ML_Epetra::SetDefaults("SA",MLList);
MLList.set("smoother: type","Chebyshev");
MLList.set("smoother: sweeps",3);
MLList.set("smoother: pre or post", "both");
#ifdef HAVE_ML_AMESOS
MLList.set("coarse: type","Amesos-KLU");
#else
MLList.set("coarse: type","Jacobi");
puts("Warning: Iterative coarse grid solve");
#endif
RCP<Epetra_Operator> Prec = rcp( new ML_Epetra::MultiLevelPreconditioner(*A, MLList) );
assert(Prec != Teuchos::null);
RCP<Belos::EpetraPrecOp> belosPrec = rcp( new Belos::EpetraPrecOp( Prec ) );
const int NumGlobalElements = RHS->GlobalLength();
if (maxiters == -1)
maxiters = NumGlobalElements/blocksize - 1;
ParameterList belosList;
belosList.set( "Block Size", blocksize );
belosList.set( "Maximum Iterations", maxiters );
belosList.set( "Convergence Tolerance", tol );
belosList.set( "Num Deflated Blocks", maxDeflate );
belosList.set( "Num Saved Blocks", maxSave );
belosList.set( "Orthogonalization", ortho );
if (numrhs > 1) {
belosList.set( "Show Maximum Residual Norm Only", true );
}
if (verbose) {
if (frequency > 0)
belosList.set( "Output Frequency", frequency );
}
else
RCP<Belos::LinearProblem<double,MV,OP> > problem
problem->setLeftPrec( belosPrec );
bool set = problem->setProblem();
if (set == false) {
if (proc_verbose)
std::cout << std::endl << "ERROR: Belos::LinearProblem failed to set up correctly!" << std::endl;
return -1;
}
RCP< Belos::SolverManager<double,MV,OP> > solver
if (proc_verbose) {
std::cout << std::endl << std::endl;
std::cout << "Dimension of matrix: " << NumGlobalElements << std::endl;
std::cout << "Number of right-hand sides: " << numrhs << std::endl;
std::cout << "Block size used by solver: " << blocksize << std::endl;
std::cout << "Maximum number of iterations allowed: " << maxiters << std::endl;
std::cout << "Relative residual tolerance: " << tol << std::endl;
std::cout << std::endl;
}
bool badRes;
for( int time_step = 0; time_step < num_time_step; time_step++){
if( time_step ){
B->Multiply(false, *LHS, *RHS);
set = problem->setProblem(LHS,RHS);
if (set == false) {
if (proc_verbose)
std::cout << std::endl << "ERROR: Belos::LinearProblem failed to set up correctly!" << std::endl;
return -1;
}
}
std::vector<double> rhs_norm( numrhs );
MVT::MvNorm( *RHS, rhs_norm );
std::cout << " RHS norm is ... " << rhs_norm[0] << std::endl;
badRes = false;
std::vector<double> actual_resids( numrhs );
Epetra_MultiVector resid(*Map, numrhs);
OPT::Apply( *A, *LHS, resid );
MVT::MvAddMv( -1.0, resid, 1.0, *RHS, resid );
MVT::MvNorm( resid, actual_resids );
MVT::MvNorm( *RHS, rhs_norm );
std::cout << " RHS norm is ... " << rhs_norm[0] << std::endl;
if (proc_verbose) {
std::cout<< "---------- Actual Residuals (normalized) ----------"<<std::endl<<std::endl;
for ( int i=0; i<numrhs; i++) {
double actRes = actual_resids[i]/rhs_norm[i];
std::cout<<"Problem "<<i<<" : \t"<< actRes <<std::endl;
if (actRes > tol) badRes = true;
}
}
success = false;
break;
}
}
if (proc_verbose) {
if (success)
std::cout << std::endl << "SUCCESS: Belos converged!" << std::endl;
else
std::cout << std::endl << "ERROR: Belos did not converge!" << std::endl;
}
}
TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success);
#ifdef EPETRA_MPI
MPI_Finalize();
#endif
return success ? EXIT_SUCCESS : EXIT_FAILURE;
}
Belos header file which uses auto-configuration information to include necessary C++ headers.
Class which describes the linear problem to be solved by the iterative solver.
Declaration and definition of Belos::PCPGSolMgr (PCPG iterative linear solver).
A linear system to solve, and its associated information.
Traits class which defines basic operations on multivectors.
Class which defines basic traits for the operator type.
PCPG iterative linear solver.
ReturnType
Whether the Belos solve converged for all linear systems.