51#include "Teuchos_CommandLineProcessor.hpp"
52#include "Teuchos_ParameterList.hpp"
59#ifdef HAVE_BELOS_TRIUTILS
60#include "Trilinos_Util_iohb.h"
67using namespace Teuchos;
69int main(
int argc,
char *argv[]) {
71 typedef std::complex<double> ST;
72 typedef ScalarTraits<ST> SCT;
73 typedef SCT::magnitudeType MT;
79 ST zero = SCT::zero();
82 bool norm_failure =
false;
84 Teuchos::GlobalMPISession session(&argc, &argv, NULL);
85 int MyPID = session.getRank();
93 bool proc_verbose =
false;
97 std::string filename(
"mhd1280b.cua");
100 CommandLineProcessor cmdp(
false,
true);
101 cmdp.setOption(
"verbose",
"quiet",&verbose,
"Print messages and results.");
102 cmdp.setOption(
"frequency",&frequency,
"Solvers frequency for printing residuals (#iters).");
103 cmdp.setOption(
"filename",&filename,
"Filename for Harwell-Boeing test matrix.");
104 cmdp.setOption(
"tol",&tol,
"Relative residual tolerance used by GCRODR solver.");
105 cmdp.setOption(
"num-rhs",&numrhs,
"Number of right-hand sides to be solved for.");
106 if (cmdp.parse(argc,argv) != CommandLineProcessor::PARSE_SUCCESSFUL) {
110 proc_verbose = verbose && (MyPID==0);
117#ifndef HAVE_BELOS_TRIUTILS
118 std::cout <<
"This test requires Triutils. Please configure with --enable-triutils." << std::endl;
120 std::cout <<
"End Result: TEST FAILED" << std::endl;
130 info = readHB_newmat_double(filename.c_str(),&dim,&dim2,&nnz,
131 &colptr,&rowind,&dvals);
132 if (info == 0 || nnz < 0) {
134 std::cout <<
"Error reading '" << filename <<
"'" << std::endl;
135 std::cout <<
"End Result: TEST FAILED" << std::endl;
141 for (
int ii=0; ii<nnz; ii++) {
142 cvals[ii] = ST(dvals[ii*2],dvals[ii*2+1]);
145 RCP< MyBetterOperator<ST> > A
149 int maxits = dim/blocksize;
151 int numRecycledBlocks = 20;
152 int numIters1, numIters2, numIters3;
153 ParameterList belosList;
154 belosList.set(
"Maximum Iterations", maxits );
155 belosList.set(
"Convergence Tolerance", tol );
157 belosList.set(
"Num Blocks", numBlocks );
158 belosList.set(
"Num Recycled Blocks", numRecycledBlocks );
162 RCP<MyMultiVec<ST> > soln = rcp(
new MyMultiVec<ST>(dim,numrhs) );
164 MVT::MvRandom( *soln );
165 OPT::Apply( *A, *soln, *rhs );
166 MVT::MvInit( *soln, zero );
168 RCP<Belos::LinearProblem<ST,MV,OP> > problem =
170 bool set = problem->setProblem();
173 std::cout << std::endl <<
"ERROR: Belos::LinearProblem failed to set up correctly!" << std::endl;
182 std::cout << std::endl << std::endl;
183 std::cout <<
"Dimension of matrix: " << dim << std::endl;
184 std::cout <<
"Number of right-hand sides: " << numrhs << std::endl;
185 std::cout <<
"Block size used by solver: " << blocksize << std::endl;
186 std::cout <<
"Max number of GCRODR iterations: " << maxits << std::endl;
187 std::cout <<
"Relative residual tolerance: " << tol << std::endl;
188 std::cout << std::endl;
193 numIters1=solver.getNumIters();
195 RCP<MyMultiVec<ST> > temp = rcp(
new MyMultiVec<ST>(dim,numrhs) );
196 std::vector<MT> norm_num(numrhs), norm_denom(numrhs);
197 OPT::Apply( *A, *soln, *temp );
198 MVT::MvAddMv( one, *rhs, -one, *temp, *temp );
199 MVT::MvNorm( *temp, norm_num );
200 MVT::MvNorm( *rhs, norm_denom );
201 for (
int i=0; i<numrhs; ++i) {
203 std::cout <<
"Relative residual "<<i<<
" : " << norm_num[i] / norm_denom[i] << std::endl;
204 if ( norm_num[i] / norm_denom[i] > tol ) {
209 MVT::MvInit( *soln, zero );
211 ret = solver.solve();
212 numIters2=solver.getNumIters();
215 MVT::MvInit( *soln, zero );
217 ret = solver.solve();
218 numIters3=solver.getNumIters();
225 if ( ret!=
Belos::Converged || norm_failure || numIters1 < numIters2 || numIters2 < numIters3 ) {
228 std::cout <<
"End Result: TEST FAILED" << std::endl;
232 std::cout <<
"End Result: TEST PASSED" << std::endl;
235 TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success);
237 return ( success ? EXIT_SUCCESS : EXIT_FAILURE );
Belos header file which uses auto-configuration information to include necessary C++ headers.
Declaration and definition of Belos::GCRODRSolMgr, which implements the GCRODR (recycling GMRES) solv...
Class which describes the linear problem to be solved by the iterative solver.
Implementation of the GCRODR (Recycling GMRES) iterative linear solver.
A linear system to solve, and its associated information.
Traits class which defines basic operations on multivectors.
Interface for multivectors used by Belos' linear solvers.
Class which defines basic traits for the operator type.
Alternative run-time polymorphic interface for operators.
Simple example of a user's defined Belos::Operator class.
Simple example of a user's defined Belos::MultiVec class.
ReturnType
Whether the Belos solve converged for all linear systems.
std::string Belos_Version()
int main(int argc, char *argv[])