65int main(
int argc,
char *argv[]) {
67 typedef std::complex<double> ST;
68 typedef ScalarTraits<ST> SCT;
69 typedef SCT::magnitudeType MT;
75 ST zero = SCT::zero();
77 Teuchos::GlobalMPISession session(&argc, &argv, NULL);
78 int MyPID = session.getRank();
86 bool norm_failure =
false;
87 bool proc_verbose =
false;
96 CommandLineProcessor cmdp(
false,
true);
97 cmdp.setOption(
"verbose",
"quiet",&verbose,
"Print messages and results.");
98 cmdp.setOption(
"pseudo",
"regular",&pseudo,
"Use pseudo-block TFQMR to solve the linear systems.");
99 cmdp.setOption(
"frequency",&frequency,
"Solvers frequency for printing residuals (#iters).");
100 cmdp.setOption(
"tol",&tol,
"Relative residual tolerance used by TFQMR solver.");
101 cmdp.setOption(
"num-rhs",&numrhs,
"Number of right-hand sides to be solved for.");
102 cmdp.setOption(
"num-restarts",&maxrestarts,
"Maximum number of restarts allowed for the TFQMR solver.");
103 cmdp.setOption(
"blocksize",&blocksize,
"Block size used by TFQMR.");
104 cmdp.setOption(
"subspace-length",&length,
"Maximum dimension of block-subspace used by TFQMR solver.");
105 if (cmdp.parse(argc,argv) != CommandLineProcessor::PARSE_SUCCESSFUL) {
109 proc_verbose = verbose && (MyPID==0);
116#ifndef HAVE_BELOS_TRIUTILS
117 std::cout <<
"This test requires Triutils. Please configure with --enable-triutils." << std::endl;
119 std::cout <<
"End Result: TEST FAILED" << std::endl;
128 std::vector<ST> diag( dim, (ST)4.0 );
129 RCP< MyOperator<ST> > A
135 int maxits = dim/blocksize;
137 ParameterList belosList;
138 belosList.set(
"Num Blocks", length );
139 belosList.set(
"Block Size", blocksize );
140 belosList.set(
"Maximum Iterations", maxits );
141 belosList.set(
"Maximum Restarts", maxrestarts );
142 belosList.set(
"Convergence Tolerance", tol );
147 belosList.set(
"Output Frequency", frequency );
156 RCP<MyMultiVec<ST> > soln = rcp(
new MyMultiVec<ST>(dim,numrhs) );
158 MVT::MvInit( *rhs, 1.0 );
159 MVT::MvInit( *soln, zero );
163 RCP<Belos::LinearProblem<ST,MV,OP> > problem =
165 bool set = problem->setProblem();
168 std::cout << std::endl <<
"ERROR: Belos::LinearProblem failed to set up correctly!" << std::endl;
177 Teuchos::RCP< Belos::SolverManager<ST,MV,OP> > solver;
187 std::cout << std::endl << std::endl;
188 std::cout <<
"Dimension of matrix: " << dim << std::endl;
189 std::cout <<
"Number of right-hand sides: " << numrhs << std::endl;
190 std::cout <<
"Block size used by solver: " << blocksize << std::endl;
191 std::cout <<
"Max number of TFQMR iterations: " << maxits << std::endl;
192 std::cout <<
"Relative residual tolerance: " << tol << std::endl;
193 std::cout << std::endl;
202 RCP<MyMultiVec<ST> > temp = rcp(
new MyMultiVec<ST>(dim,numrhs) );
203 OPT::Apply( *A, *soln, *temp );
204 MVT::MvAddMv( one, *rhs, -one, *temp, *temp );
205 std::vector<MT> norm_num(numrhs), norm_denom(numrhs);
206 MVT::MvNorm( *temp, norm_num );
207 MVT::MvNorm( *rhs, norm_denom );
208 for (
int i=0; i<numrhs; ++i) {
210 std::cout <<
"Relative residual "<<i<<
" : " << norm_num[i] / norm_denom[i] << std::endl;
211 if ( norm_num[i] / norm_denom[i] > tol ) {
219 std::cout <<
"End Result: TEST PASSED" << std::endl;
222 std::cout <<
"End Result: TEST FAILED" << std::endl;
225 TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success);
227 return ( success ? EXIT_SUCCESS : EXIT_FAILURE );