72int main(
int argc,
char *argv[]) {
74 typedef std::complex<double> ST;
75 typedef ScalarTraits<ST> SCT;
76 typedef SCT::magnitudeType MT;
82 ST zero = SCT::zero();
84 Teuchos::GlobalMPISession session(&argc, &argv, NULL);
85 int MyPID = session.getRank();
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");
110 Teuchos::CommandLineProcessor cmdp(
false,
true);
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.");
125 if (cmdp.parse(argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
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]);
165 RCP< MyBetterOperator<ST> > A
172 RCP<MyMultiVec<ST> > soln = rcp(
new MyMultiVec<ST>(dim,numrhs) );
174 MVT::MvRandom( *soln );
175 OPT::Apply( *A, *soln, *rhs );
176 MVT::MvInit( *soln, zero );
180 RCP<Belos::LinearProblem<ST,MV,OP> > problem =
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;
196 ParameterList belosList;
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 );
210 ParameterList polyList;
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;
255 RCP<MyMultiVec<ST> > temp = rcp(
new MyMultiVec<ST>(dim,numrhs) );
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;
295 TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success);
297 return ( success ? EXIT_SUCCESS : EXIT_FAILURE );