53 #include "Teuchos_CommandLineProcessor.hpp" 54 #include "Teuchos_ParameterList.hpp" 55 #include "Teuchos_StandardCatchMacros.hpp" 56 #include "Teuchos_StandardCatchMacros.hpp" 63 #ifdef HAVE_BELOS_TRIUTILS 64 #include "Trilinos_Util_iohb.h" 73 int main(
int argc,
char *argv[]) {
76 typedef std::complex<double> ST;
78 typedef std::complex<double> ST;
80 std::cout <<
"Not compiled with std::complex support." << std::endl;
81 std::cout <<
"End Result: TEST FAILED" << std::endl;
85 typedef ScalarTraits<ST> SCT;
86 typedef SCT::magnitudeType MT;
92 ST zero = SCT::zero();
94 Teuchos::GlobalMPISession session(&argc, &argv, NULL);
95 int MyPID = session.getRank();
100 bool success =
false;
101 bool verbose =
false;
104 bool norm_failure =
false;
105 bool proc_verbose =
false;
110 int maxrestarts = 15;
112 std::string filename(
"mhd1280b.cua");
115 CommandLineProcessor cmdp(
false,
true);
116 cmdp.setOption(
"verbose",
"quiet",&verbose,
"Print messages and results.");
117 cmdp.setOption(
"pseudo",
"regular",&pseudo,
"Use pseudo-block GMRES to solve the linear systems.");
118 cmdp.setOption(
"frequency",&frequency,
"Solvers frequency for printing residuals (#iters).");
119 cmdp.setOption(
"filename",&filename,
"Filename for Harwell-Boeing test matrix.");
120 cmdp.setOption(
"tol",&tol,
"Relative residual tolerance used by GMRES solver.");
121 cmdp.setOption(
"num-rhs",&numrhs,
"Number of right-hand sides to be solved for.");
122 cmdp.setOption(
"num-restarts",&maxrestarts,
"Maximum number of restarts allowed for the GMRES solver.");
123 cmdp.setOption(
"blocksize",&blocksize,
"Block size used by GMRES.");
124 cmdp.setOption(
"subspace-length",&length,
"Maximum dimension of block-subspace used by GMRES solver.");
125 if (cmdp.parse(argc,argv) != 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
171 int maxits = dim/blocksize;
173 ParameterList belosList;
174 belosList.set(
"Num Blocks", length );
175 belosList.set(
"Block Size", blocksize );
176 belosList.set(
"Maximum Iterations", maxits );
177 belosList.set(
"Maximum Restarts", maxrestarts );
178 belosList.set(
"Convergence Tolerance", tol );
183 belosList.set(
"Output Frequency", frequency );
192 RCP<MyMultiVec<ST> > soln = rcp(
new MyMultiVec<ST>(dim,numrhs) );
194 MVT::MvRandom( *soln );
195 OPT::Apply( *A, *soln, *rhs );
196 MVT::MvInit( *soln, zero );
200 RCP<Belos::LinearProblem<ST,MV,OP> > problem =
202 bool set = problem->setProblem();
205 std::cout << std::endl <<
"ERROR: Belos::LinearProblem failed to set up correctly!" << std::endl;
218 Teuchos::RCP< Belos::SolverManager<ST,MV,OP> > solver;
224 solver->setDebugStatusTest( Teuchos::rcp(&debugTest,
false) );
230 std::cout << std::endl << std::endl;
231 std::cout <<
"Dimension of matrix: " << dim << std::endl;
232 std::cout <<
"Number of right-hand sides: " << numrhs << std::endl;
233 std::cout <<
"Block size used by solver: " << blocksize << std::endl;
234 std::cout <<
"Max number of Gmres iterations: " << maxits << std::endl;
235 std::cout <<
"Relative residual tolerance: " << tol << std::endl;
236 std::cout << std::endl;
245 RCP<MyMultiVec<ST> > temp = rcp(
new MyMultiVec<ST>(dim,numrhs) );
246 OPT::Apply( *A, *soln, *temp );
247 MVT::MvAddMv( one, *rhs, -one, *temp, *temp );
248 std::vector<MT> norm_num(numrhs), norm_denom(numrhs);
249 MVT::MvNorm( *temp, norm_num );
250 MVT::MvNorm( *rhs, norm_denom );
251 for (
int i=0; i<numrhs; ++i) {
253 std::cout <<
"Relative residual "<<i<<
" : " << norm_num[i] / norm_denom[i] << std::endl;
254 if ( norm_num[i] / norm_denom[i] > tol ) {
260 const std::vector<MT> residualLog = debugTest.
getLogResNorm();
261 if (numrhs==1 && proc_verbose && residualLog.size())
263 std::cout <<
"Absolute residual 2-norm [ " << residualLog.size() <<
" ] : ";
264 for (
unsigned int i=0; i<residualLog.size(); i++)
265 std::cout << residualLog[i] <<
" ";
266 std::cout << std::endl;
267 std::cout <<
"Final abs 2-norm / rhs 2-norm : " << residualLog[residualLog.size()-1] / norm_denom[0] << std::endl;
279 std::cout <<
"End Result: TEST PASSED" << std::endl;
282 std::cout <<
"End Result: TEST FAILED" << std::endl;
285 TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success);
287 return ( success ? EXIT_SUCCESS : EXIT_FAILURE );
std::string Belos_Version()
A Belos::StatusTest debugging class for storing the absolute residual norms generated during a solve...
Interface to Block GMRES and Flexible GMRES.
Belos::StatusTest debugging class for storing the absolute residual norms generated during a solve...
Traits class which defines basic operations on multivectors.
Simple example of a user's defined Belos::MultiVec class.
Alternative run-time polymorphic interface for operators.
The Belos::BlockGmresSolMgr provides a solver manager for the BlockGmres linear solver.
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
Interface to standard and "pseudoblock" GMRES.
ReturnType
Whether the Belos solve converged for all linear systems.
Interface for multivectors used by Belos' linear solvers.
Class which defines basic traits for the operator type.
The Belos::PseudoBlockGmresSolMgr provides a solver manager for the BlockGmres linear solver...
const std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > & getLogResNorm() const
Returns the log of the absolute residual norm from the iteration.
Belos header file which uses auto-configuration information to include necessary C++ headers...
int main(int argc, char *argv[])
Simple example of a user's defined Belos::Operator class.