71int main(
int argc,
char *argv[]) {
75 typedef ScalarTraits<ST> SCT;
76 typedef SCT::magnitudeType MT;
82 ST zero = SCT::zero();
84 Teuchos::GlobalMPISession session(&argc, &argv, NULL);
95 bool norm_failure =
false;
96 bool proc_verbose =
false;
97 bool use_single_red =
false;
98 bool combineConvInner =
false;
102 std::string filename(
"A.hb");
105 CommandLineProcessor cmdp(
false,
true);
106 cmdp.setOption(
"verbose",
"quiet",&verbose,
"Print messages and results.");
107 cmdp.setOption(
"pseudo",
"regular",&pseudo,
"Use pseudo-block CG to solve the linear systems.");
108 cmdp.setOption(
"frequency",&frequency,
"Solvers frequency for printing residuals (#iters).");
109 cmdp.setOption(
"filename",&filename,
"Filename for Harwell-Boeing test matrix.");
110 cmdp.setOption(
"tol",&tol,
"Relative residual tolerance used by CG solver.");
111 cmdp.setOption(
"num-rhs",&numrhs,
"Number of right-hand sides to be solved for.");
112 cmdp.setOption(
"blocksize",&blocksize,
"Block size used by CG .");
113 cmdp.setOption(
"use-single-red",
"use-standard-red",&use_single_red,
"Use single-reduction variant of CG iteration.");
114 cmdp.setOption(
"combine-conv-inner",
"separate-conv-inner",&combineConvInner,
"Combine convergence detection and inner products of single-reduce CG iteration.");
115 if (cmdp.parse(argc,argv) != CommandLineProcessor::PARSE_SUCCESSFUL) {
119 proc_verbose = verbose && (MyPID==0);
127#ifndef HAVE_BELOS_TRIUTILS
128 std::cout <<
"This test requires Triutils. Please configure with --enable-triutils." << std::endl;
130 std::cout <<
"End Result: TEST FAILED" << std::endl;
140 info = readHB_newmat_double(filename.c_str(),&dim,&dim2,&nnz,
141 &colptr,&rowind,&cvals);
142 if (info == 0 || nnz < 0) {
144 std::cout <<
"Error reading '" << filename <<
"'" << std::endl;
145 std::cout <<
"End Result: TEST FAILED" << std::endl;
150 RCP< MyBetterOperator<ST> > A
159 int maxits = dim/blocksize;
161 ParameterList belosList;
162 belosList.set(
"Block Size", blocksize );
163 belosList.set(
"Maximum Iterations", maxits );
164 belosList.set(
"Convergence Tolerance", tol );
165 if (blocksize == 1) {
167 belosList.set(
"Use Single Reduction", use_single_red );
168 if (combineConvInner)
169 belosList.set(
"Fold Convergence Detection Into Allreduce", combineConvInner );
175 belosList.set(
"Output Frequency", frequency );
184 RCP<MyMultiVec<ST> > soln = rcp(
new MyMultiVec<ST>(dim,numrhs) );
186 MVT::MvRandom( *soln );
187 OPT::Apply( *A, *soln, *rhs );
188 MVT::MvInit( *soln, zero );
192 RCP<Belos::LinearProblem<ST,MV,OP> > problem =
194 bool set = problem->setProblem();
197 std::cout << std::endl <<
"ERROR: Belos::LinearProblem failed to set up correctly!" << std::endl;
206 Teuchos::RCP< Belos::SolverManager<ST,MV,OP> > solver;
216 std::cout << std::endl << std::endl;
217 std::cout <<
"Dimension of matrix: " << dim << std::endl;
218 std::cout <<
"Number of right-hand sides: " << numrhs << std::endl;
219 std::cout <<
"Block size used by solver: " << blocksize << std::endl;
220 std::cout <<
"Max number of CG iterations: " << maxits << std::endl;
221 std::cout <<
"Relative residual tolerance: " << tol << std::endl;
222 std::cout << std::endl;
231 RCP<MyMultiVec<ST> > temp = rcp(
new MyMultiVec<ST>(dim,numrhs) );
232 OPT::Apply( *A, *soln, *temp );
233 MVT::MvAddMv( one, *rhs, -one, *temp, *temp );
234 std::vector<MT> norm_num(numrhs), norm_denom(numrhs);
235 MVT::MvNorm( *temp, norm_num );
236 MVT::MvNorm( *rhs, norm_denom );
237 for (
int i=0; i<numrhs; ++i) {
239 std::cout <<
"Relative residual "<<i<<
" : " << norm_num[i] / norm_denom[i] << std::endl;
240 if ( norm_num[i] / norm_denom[i] > tol ) {
246 MT ach_tol = solver->achievedTol();
248 std::cout <<
"Achieved tol : "<<ach_tol<<std::endl;
259 std::cout <<
"End Result: TEST PASSED" << std::endl;
262 std::cout <<
"End Result: TEST FAILED" << std::endl;
265 TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success);
267 return ( success ? EXIT_SUCCESS : EXIT_FAILURE );