52 #include "Teuchos_CommandLineProcessor.hpp" 53 #include "Teuchos_ParameterList.hpp" 54 #include "Teuchos_StandardCatchMacros.hpp" 55 #include "Teuchos_StandardCatchMacros.hpp" 62 #ifdef HAVE_BELOS_TRIUTILS 63 #include "Trilinos_Util_iohb.h" 72 int main(
int argc,
char *argv[]) {
74 typedef std::complex<double> ST;
75 typedef ScalarTraits<ST> SCT;
76 typedef SCT::magnitudeType MT;
82 ST zero = SCT::zero();
84 Teuchos::GlobalMPISession session(&argc, &argv, NULL);
85 int MyPID = session.getRank();
94 bool norm_failure =
false;
95 bool proc_verbose =
false;
96 bool userandomrhs =
true;
102 int maxsubspace = 50;
103 int maxrestarts = 15;
104 std::string outersolver(
"Block Gmres");
105 std::string filename(
"mhd1280b.cua");
106 std::string polyType(
"Arnoldi");
110 Teuchos::CommandLineProcessor cmdp(
false,
true);
111 cmdp.setOption(
"verbose",
"quiet",&verbose,
"Print messages and results.");
112 cmdp.setOption(
"use-random-rhs",
"use-rhs",&userandomrhs,
"Use linear system RHS or random RHS to generate polynomial.");
113 cmdp.setOption(
"frequency",&frequency,
"Solvers frequency for printing residuals (#iters).");
114 cmdp.setOption(
"filename",&filename,
"Filename for test matrix. Acceptable file extensions: *.hb,*.mtx,*.triU,*.triS");
115 cmdp.setOption(
"outersolver",&outersolver,
"Name of outer solver to be used with GMRES poly");
116 cmdp.setOption(
"tol",&tol,
"Relative residual tolerance used by GMRES solver.");
117 cmdp.setOption(
"poly-tol",&polytol,
"Relative residual tolerance used to construct the GMRES polynomial.");
118 cmdp.setOption(
"poly-type",&polyType,
"Polynomial type (Roots, Arnoldi, or Gmres).");
119 cmdp.setOption(
"num-rhs",&numrhs,
"Number of right-hand sides to be solved for.");
120 cmdp.setOption(
"block-size",&blocksize,
"Block size used by GMRES.");
121 cmdp.setOption(
"max-iters",&maxiters,
"Maximum number of iterations per linear system (-1 = adapted to problem/block size).");
122 cmdp.setOption(
"max-degree",&maxdegree,
"Maximum degree of the GMRES polynomial.");
123 cmdp.setOption(
"max-subspace",&maxsubspace,
"Maximum number of blocks the solver can use for the subspace.");
124 cmdp.setOption(
"max-restarts",&maxrestarts,
"Maximum number of restarts allowed for GMRES solver.");
125 if (cmdp.parse(argc,argv) != Teuchos::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
172 RCP<MyMultiVec<ST> > soln = rcp(
new MyMultiVec<ST>(dim,numrhs) );
174 MVT::MvRandom( *soln );
175 OPT::Apply( *A, *soln, *rhs );
176 MVT::MvInit( *soln, zero );
180 RCP<Belos::LinearProblem<ST,MV,OP> > problem =
182 problem->setInitResVec( rhs );
183 bool set = problem->setProblem();
186 std::cout << std::endl <<
"ERROR: Belos::LinearProblem failed to set up correctly!" << std::endl;
194 maxiters = dim/blocksize - 1;
196 ParameterList belosList;
197 belosList.set(
"Num Blocks", maxsubspace);
198 belosList.set(
"Block Size", blocksize );
199 belosList.set(
"Maximum Iterations", maxiters );
200 belosList.set(
"Maximum Restarts", maxrestarts );
201 belosList.set(
"Convergence Tolerance", tol );
206 belosList.set(
"Output Frequency", frequency );
208 belosList.set(
"Verbosity", verbosity );
210 ParameterList polyList;
211 polyList.set(
"Maximum Degree", maxdegree );
212 polyList.set(
"Polynomial Tolerance", polytol );
213 polyList.set(
"Polynomial Type", polyType );
214 polyList.set(
"Verbosity", verbosity );
215 polyList.set(
"Random RHS", userandomrhs );
216 if ( outersolver !=
"" ) {
217 polyList.set(
"Outer Solver", outersolver );
218 polyList.set(
"Outer Solver Params", belosList );
239 std::cout << std::endl << std::endl;
240 std::cout <<
"Dimension of matrix: " << dim << std::endl;
241 std::cout <<
"Number of right-hand sides: " << numrhs << std::endl;
242 std::cout <<
"Block size used by solver: " << blocksize << std::endl;
243 std::cout <<
"Max number of Gmres iterations: " << maxiters << std::endl;
244 std::cout <<
"Relative residual tolerance: " << tol << std::endl;
245 std::cout << std::endl;
254 RCP<MyMultiVec<ST> > temp = rcp(
new MyMultiVec<ST>(dim,numrhs) );
255 OPT::Apply( *A, *soln, *temp );
256 MVT::MvAddMv( one, *rhs, -one, *temp, *temp );
257 std::vector<MT> norm_num(numrhs), norm_denom(numrhs);
258 MVT::MvNorm( *temp, norm_num );
259 MVT::MvNorm( *rhs, norm_denom );
260 for (
int i=0; i<numrhs; ++i) {
262 std::cout <<
"Relative residual "<<i<<
" : " << norm_num[i] / norm_denom[i] << std::endl;
263 if ( norm_num[i] / norm_denom[i] > tol ) {
269 const std::vector<MT> residualLog = debugTest.
getLogResNorm();
270 if (numrhs==1 && proc_verbose && residualLog.size())
272 std::cout <<
"Absolute residual 2-norm [ " << residualLog.size() <<
" ] : ";
273 for (
unsigned int i=0; i<residualLog.size(); i++)
274 std::cout << residualLog[i] <<
" ";
275 std::cout << std::endl;
276 std::cout <<
"Final abs 2-norm / rhs 2-norm : " << residualLog[residualLog.size()-1] / norm_denom[0] << std::endl;
288 std::cout <<
"End Result: TEST PASSED" << std::endl;
291 std::cout <<
"End Result: TEST FAILED" << std::endl;
294 TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success);
296 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...
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.
int main(int argc, char *argv[])
Alternative run-time polymorphic interface for operators.
Declaration and definition of Belos::GmresPolySolMgr (hybrid block GMRES linear solver).
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
The GMRES polynomial can be created in conjunction with any standard preconditioner.
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.
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...
Simple example of a user's defined Belos::Operator class.