73int main(
int argc,
char *argv[]) {
75 typedef std::complex<double> ST;
76 typedef ScalarTraits<ST> SCT;
77 typedef SCT::magnitudeType MT;
83 ST zero = SCT::zero();
85 Teuchos::GlobalMPISession session(&argc, &argv, NULL);
86 int MyPID = session.getRank();
96 bool norm_failure =
false;
97 bool proc_verbose =
false;
102 int maxrestarts = 15;
104 std::string filename(
"mhd1280b.cua");
105 std::string ortho(
"ICGS");
108 CommandLineProcessor cmdp(
false,
true);
109 cmdp.setOption(
"verbose",
"quiet",&verbose,
"Print messages and results.");
110 cmdp.setOption(
"debug",
"no-debug",&debug,
"Print debug messages.");
111 cmdp.setOption(
"pseudo",
"regular",&pseudo,
"Use pseudo-block GMRES to solve the linear systems.");
112 cmdp.setOption(
"frequency",&frequency,
"Solvers frequency for printing residuals (#iters).");
113 cmdp.setOption(
"filename",&filename,
"Filename for Harwell-Boeing test matrix.");
114 cmdp.setOption(
"tol",&tol,
"Relative residual tolerance used by GMRES solver.");
115 cmdp.setOption(
"num-rhs",&numrhs,
"Number of right-hand sides to be solved for.");
116 cmdp.setOption(
"num-restarts",&maxrestarts,
"Maximum number of restarts allowed for the GMRES solver.");
117 cmdp.setOption(
"blocksize",&blocksize,
"Block size used by GMRES.");
118 cmdp.setOption(
"subspace-length",&length,
"Maximum dimension of block-subspace used by GMRES solver.");
119 cmdp.setOption(
"ortho-type",&ortho,
"Orthogonalization type, either DGKS, ICGS or IMGS");
120 if (cmdp.parse(argc,argv) != CommandLineProcessor::PARSE_SUCCESSFUL) {
124 proc_verbose = verbose && (MyPID==0);
131#ifndef HAVE_BELOS_TRIUTILS
132 std::cout <<
"This test requires Triutils. Please configure with --enable-triutils." << std::endl;
134 std::cout <<
"End Result: TEST FAILED" << std::endl;
145 info = readHB_newmat_double(filename.c_str(),&dim,&dim2,&nnz,
146 &colptr,&rowind,&dvals);
147 if (info == 0 || nnz < 0) {
149 std::cout <<
"Error reading '" << filename <<
"'" << std::endl;
150 std::cout <<
"End Result: TEST FAILED" << std::endl;
156 for (
int ii=0; ii<nnz; ii++) {
157 cvals[ii] = ST(dvals[ii*2],dvals[ii*2+1]);
160 RCP< MyBetterOperator<ST> > A
166 int maxits = dim/blocksize;
168 ParameterList belosList;
169 belosList.set(
"Num Blocks", length );
170 belosList.set(
"Block Size", blocksize );
171 belosList.set(
"Maximum Iterations", maxits );
172 belosList.set(
"Maximum Restarts", maxrestarts );
173 belosList.set(
"Convergence Tolerance", tol );
174 belosList.set(
"Orthogonalization", ortho );
180 belosList.set(
"Verbosity", verbosity );
182 belosList.set(
"Output Frequency", frequency );
191 RCP<MyMultiVec<ST> > soln = rcp(
new MyMultiVec<ST>(dim,numrhs) );
193 MVT::MvRandom( *soln );
194 OPT::Apply( *A, *soln, *rhs );
195 MVT::MvInit( *soln, zero );
199 RCP<Belos::LinearProblem<ST,MV,OP> > problem =
201 bool set = problem->setProblem();
204 std::cout << std::endl <<
"ERROR: Belos::LinearProblem failed to set up correctly!" << std::endl;
217 Teuchos::RCP< Belos::SolverManager<ST,MV,OP> > solver;
223 solver->setDebugStatusTest( Teuchos::rcp(&debugTest,
false) );
229 std::cout << std::endl << std::endl;
230 std::cout <<
"Dimension of matrix: " << dim << std::endl;
231 std::cout <<
"Number of right-hand sides: " << numrhs << std::endl;
232 std::cout <<
"Block size used by solver: " << blocksize << std::endl;
233 std::cout <<
"Max number of Gmres iterations: " << maxits << std::endl;
234 std::cout <<
"Relative residual tolerance: " << tol << std::endl;
235 std::cout << std::endl;
244 RCP<MyMultiVec<ST> > temp = rcp(
new MyMultiVec<ST>(dim,numrhs) );
245 OPT::Apply( *A, *soln, *temp );
246 MVT::MvAddMv( one, *rhs, -one, *temp, *temp );
247 std::vector<MT> norm_num(numrhs), norm_denom(numrhs);
248 MVT::MvNorm( *temp, norm_num );
249 MVT::MvNorm( *rhs, norm_denom );
250 for (
int i=0; i<numrhs; ++i) {
252 std::cout <<
"Relative residual "<<i<<
" : " << norm_num[i] / norm_denom[i] << std::endl;
253 if ( norm_num[i] / norm_denom[i] > tol ) {
259 const std::vector<MT> residualLog = debugTest.
getLogResNorm();
260 if (numrhs==1 && proc_verbose && residualLog.size())
262 std::cout <<
"Absolute residual 2-norm [ " << residualLog.size() <<
" ] : ";
263 for (
unsigned int i=0; i<residualLog.size(); i++)
264 std::cout << residualLog[i] <<
" ";
265 std::cout << std::endl;
266 std::cout <<
"Final abs 2-norm / rhs 2-norm : " << residualLog[residualLog.size()-1] / norm_denom[0] << std::endl;
278 std::cout <<
"End Result: TEST PASSED" << std::endl;
281 std::cout <<
"End Result: TEST FAILED" << std::endl;
284 TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success);
286 return ( success ? EXIT_SUCCESS : EXIT_FAILURE );