62 #ifdef HAVE_BELOS_TRIUTILS
63 #include "Trilinos_Util_iohb.h"
70 using namespace Teuchos;
72 int main(
int argc,
char *argv[]) {
74 typedef std::complex<double> ST;
76 typedef SCT::magnitudeType MT;
82 ST zero = SCT::zero();
94 bool norm_failure =
false;
95 bool proc_verbose =
false;
96 bool userandomrhs =
true;
102 int maxsubspace = 50;
103 int maxrestarts = 15;
104 std::string outersolver(
"Block Gmres");
105 std::string
filename(
"mhd1280b.cua");
106 std::string polyType(
"Arnoldi");
111 cmdp.
setOption(
"verbose",
"quiet",&verbose,
"Print messages and results.");
112 cmdp.
setOption(
"use-random-rhs",
"use-rhs",&userandomrhs,
"Use linear system RHS or random RHS to generate polynomial.");
113 cmdp.
setOption(
"frequency",&frequency,
"Solvers frequency for printing residuals (#iters).");
114 cmdp.
setOption(
"filename",&filename,
"Filename for test matrix. Acceptable file extensions: *.hb,*.mtx,*.triU,*.triS");
115 cmdp.
setOption(
"outersolver",&outersolver,
"Name of outer solver to be used with GMRES poly");
116 cmdp.
setOption(
"tol",&tol,
"Relative residual tolerance used by GMRES solver.");
117 cmdp.
setOption(
"poly-tol",&polytol,
"Relative residual tolerance used to construct the GMRES polynomial.");
118 cmdp.
setOption(
"poly-type",&polyType,
"Polynomial type (Roots, Arnoldi, or Gmres).");
119 cmdp.
setOption(
"num-rhs",&numrhs,
"Number of right-hand sides to be solved for.");
120 cmdp.
setOption(
"block-size",&blocksize,
"Block size used by GMRES.");
121 cmdp.
setOption(
"max-iters",&maxiters,
"Maximum number of iterations per linear system (-1 = adapted to problem/block size).");
122 cmdp.
setOption(
"max-degree",&maxdegree,
"Maximum degree of the GMRES polynomial.");
123 cmdp.
setOption(
"max-subspace",&maxsubspace,
"Maximum number of blocks the solver can use for the subspace.");
124 cmdp.
setOption(
"max-restarts",&maxrestarts,
"Maximum number of restarts allowed for GMRES solver.");
129 proc_verbose = verbose && (MyPID==0);
136 #ifndef HAVE_BELOS_TRIUTILS
137 std::cout <<
"This test requires Triutils. Please configure with --enable-triutils." << std::endl;
139 std::cout <<
"End Result: TEST FAILED" << std::endl;
150 info = readHB_newmat_double(filename.c_str(),&dim,&dim2,&nnz,
151 &colptr,&rowind,&dvals);
152 if (info == 0 || nnz < 0) {
154 std::cout <<
"Error reading '" << filename <<
"'" << std::endl;
155 std::cout <<
"End Result: TEST FAILED" << std::endl;
161 for (
int ii=0; ii<nnz; ii++) {
162 cvals[ii] = ST(dvals[ii*2],dvals[ii*2+1]);
174 MVT::MvRandom( *soln );
175 OPT::Apply( *A, *soln, *rhs );
176 MVT::MvInit( *soln, zero );
182 problem->setInitResVec( rhs );
183 bool set = problem->setProblem();
186 std::cout << std::endl <<
"ERROR: Belos::LinearProblem failed to set up correctly!" << std::endl;
194 maxiters = dim/blocksize - 1;
197 belosList.
set(
"Num Blocks", maxsubspace);
198 belosList.
set(
"Block Size", blocksize );
199 belosList.
set(
"Maximum Iterations", maxiters );
200 belosList.
set(
"Maximum Restarts", maxrestarts );
201 belosList.
set(
"Convergence Tolerance", tol );
206 belosList.
set(
"Output Frequency", frequency );
208 belosList.
set(
"Verbosity", verbosity );
211 polyList.
set(
"Maximum Degree", maxdegree );
212 polyList.
set(
"Polynomial Tolerance", polytol );
213 polyList.
set(
"Polynomial Type", polyType );
214 polyList.
set(
"Verbosity", verbosity );
215 polyList.
set(
"Random RHS", userandomrhs );
216 if ( outersolver !=
"" ) {
217 polyList.
set(
"Outer Solver", outersolver );
218 polyList.
set(
"Outer Solver Params", belosList );
239 std::cout << std::endl << std::endl;
240 std::cout <<
"Dimension of matrix: " << dim << std::endl;
241 std::cout <<
"Number of right-hand sides: " << numrhs << std::endl;
242 std::cout <<
"Block size used by solver: " << blocksize << std::endl;
243 std::cout <<
"Max number of Gmres iterations: " << maxiters << std::endl;
244 std::cout <<
"Relative residual tolerance: " << tol << std::endl;
245 std::cout << std::endl;
252 std::cout <<
"Belos::GmresPolySolMgr returned achievedTol: " << solver->achievedTol() << std::endl << std::endl;
256 OPT::Apply( *A, *soln, *temp );
257 MVT::MvAddMv( one, *rhs, -one, *temp, *temp );
258 std::vector<MT> norm_num(numrhs), norm_denom(numrhs);
259 MVT::MvNorm( *temp, norm_num );
260 MVT::MvNorm( *rhs, norm_denom );
261 for (
int i=0; i<numrhs; ++i) {
263 std::cout <<
"Relative residual "<<i<<
" : " << norm_num[i] / norm_denom[i] << std::endl;
264 if ( norm_num[i] / norm_denom[i] > tol ) {
270 const std::vector<MT> residualLog = debugTest.
getLogResNorm();
271 if (numrhs==1 && proc_verbose && residualLog.size())
273 std::cout <<
"Absolute residual 2-norm [ " << residualLog.size() <<
" ] : ";
274 for (
unsigned int i=0; i<residualLog.size(); i++)
275 std::cout << residualLog[i] <<
" ";
276 std::cout << std::endl;
277 std::cout <<
"Final abs 2-norm / rhs 2-norm : " << residualLog[residualLog.size()-1] / norm_denom[0] << std::endl;
289 std::cout <<
"End Result: TEST PASSED" << std::endl;
292 std::cout <<
"End Result: TEST FAILED" << std::endl;
297 return ( success ? EXIT_SUCCESS : EXIT_FAILURE );
std::string Belos_Version()
int main(int argc, char *argv[])
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
const std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > & getLogResNorm() const
Returns the log of the absolute residual norm from the iteration.
A Belos::StatusTest debugging class for storing the absolute residual norms generated during a solve...
Belos::StatusTest debugging class for storing the absolute residual norms generated during a solve...
Traits class which defines basic operations on multivectors.
Simple example of a user's defined Belos::MultiVec class.
Alternative run-time polymorphic interface for operators.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void setOption(const char option_true[], const char option_false[], bool *option_val, const char documentation[]=NULL)
Declaration and definition of Belos::GmresPolySolMgr (hybrid block GMRES linear solver).
#define TEUCHOS_STANDARD_CATCH_STATEMENTS(VERBOSE, ERR_STREAM, SUCCESS_FLAG)
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
EParseCommandLineReturn parse(int argc, char *argv[], std::ostream *errout=&std::cerr) const
The GMRES polynomial can be created in conjunction with any standard preconditioner.
ReturnType
Whether the Belos solve converged for all linear systems.
Interface for multivectors used by Belos' linear solvers.
Class which defines basic traits for the operator type.
Belos header file which uses auto-configuration information to include necessary C++ headers...
Simple example of a user's defined Belos::Operator class.