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[]) {
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();
95 bool norm_failure =
false;
96 bool proc_verbose =
false;
101 int maxrestarts = 15;
103 std::string filename(
"mhd1280b.cua");
106 CommandLineProcessor cmdp(
false,
true);
107 cmdp.setOption(
"verbose",
"quiet",&verbose,
"Print messages and results.");
108 cmdp.setOption(
"pseudo",
"regular",&pseudo,
"Use pseudo-block GMRES to solve the linear systems.");
109 cmdp.setOption(
"frequency",&frequency,
"Solvers frequency for printing residuals (#iters).");
110 cmdp.setOption(
"filename",&filename,
"Filename for Harwell-Boeing test matrix.");
111 cmdp.setOption(
"tol",&tol,
"Relative residual tolerance used by GMRES solver.");
112 cmdp.setOption(
"num-rhs",&numrhs,
"Number of right-hand sides to be solved for.");
113 cmdp.setOption(
"num-restarts",&maxrestarts,
"Maximum number of restarts allowed for the GMRES solver.");
114 cmdp.setOption(
"blocksize",&blocksize,
"Block size used by GMRES.");
115 cmdp.setOption(
"subspace-length",&length,
"Maximum dimension of block-subspace used by GMRES solver.");
116 if (cmdp.parse(argc,argv) != CommandLineProcessor::PARSE_SUCCESSFUL) {
120 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;
141 info = readHB_newmat_double(filename.c_str(),&dim,&dim2,&nnz,
142 &colptr,&rowind,&dvals);
143 if (info == 0 || nnz < 0) {
145 std::cout <<
"Error reading '" << filename <<
"'" << std::endl;
146 std::cout <<
"End Result: TEST FAILED" << std::endl;
152 for (
int ii=0; ii<nnz; ii++) {
153 cvals[ii] = ST(dvals[ii*2],dvals[ii*2+1]);
156 RCP< MyBetterOperator<ST> > A
162 int maxits = dim/blocksize;
164 ParameterList belosList;
165 belosList.set(
"Num Blocks", length );
166 belosList.set(
"Block Size", blocksize );
167 belosList.set(
"Maximum Iterations", maxits );
168 belosList.set(
"Maximum Restarts", maxrestarts );
169 belosList.set(
"Convergence Tolerance", tol );
174 belosList.set(
"Output Frequency", frequency );
183 RCP<MyMultiVec<ST> > soln = rcp(
new MyMultiVec<ST>(dim,numrhs) );
185 MVT::MvRandom( *soln );
186 OPT::Apply( *A, *soln, *rhs );
187 MVT::MvInit( *soln, zero );
191 RCP<Belos::LinearProblem<ST,MV,OP> > problem =
193 bool set = problem->setProblem();
196 std::cout << std::endl <<
"ERROR: Belos::LinearProblem failed to set up correctly!" << std::endl;
209 Teuchos::RCP< Belos::SolverManager<ST,MV,OP> > solver;
215 solver->setDebugStatusTest( Teuchos::rcp(&debugTest,
false) );
221 std::cout << std::endl << std::endl;
222 std::cout <<
"Dimension of matrix: " << dim << std::endl;
223 std::cout <<
"Number of right-hand sides: " << numrhs << std::endl;
224 std::cout <<
"Block size used by solver: " << blocksize << std::endl;
225 std::cout <<
"Max number of Gmres iterations: " << maxits << std::endl;
226 std::cout <<
"Relative residual tolerance: " << tol << std::endl;
227 std::cout << std::endl;
236 RCP<MyMultiVec<ST> > temp = rcp(
new MyMultiVec<ST>(dim,numrhs) );
237 OPT::Apply( *A, *soln, *temp );
238 MVT::MvAddMv( one, *rhs, -one, *temp, *temp );
239 std::vector<MT> norm_num(numrhs), norm_denom(numrhs);
240 MVT::MvNorm( *temp, norm_num );
241 MVT::MvNorm( *rhs, norm_denom );
242 for (
int i=0; i<numrhs; ++i) {
244 std::cout <<
"Relative residual "<<i<<
" : " << norm_num[i] / norm_denom[i] << std::endl;
245 if ( norm_num[i] / norm_denom[i] > tol ) {
251 const std::vector<MT> residualLog = debugTest.
getLogResNorm();
252 if (numrhs==1 && proc_verbose && residualLog.size())
254 std::cout <<
"Absolute residual 2-norm [ " << residualLog.size() <<
" ] : ";
255 for (
unsigned int i=0; i<residualLog.size(); i++)
256 std::cout << residualLog[i] <<
" ";
257 std::cout << std::endl;
258 std::cout <<
"Final abs 2-norm / rhs 2-norm : " << residualLog[residualLog.size()-1] / norm_denom[0] << std::endl;
270 std::cout <<
"End Result: TEST PASSED" << std::endl;
273 std::cout <<
"End Result: TEST FAILED" << std::endl;
276 TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success);
278 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.