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 );