53#include "Teuchos_CommandLineProcessor.hpp"
54#include "Teuchos_ParameterList.hpp"
55#include "Teuchos_StandardCatchMacros.hpp"
64using namespace Teuchos;
66int main(
int argc,
char *argv[]) {
68 typedef std::complex<double> ST;
69 typedef ScalarTraits<ST> SCT;
70 typedef SCT::magnitudeType MT;
76 ST zero = SCT::zero();
78 Teuchos::GlobalMPISession session(&argc, &argv, NULL);
79 int MyPID = session.getRank();
87 bool norm_failure =
false;
88 bool proc_verbose =
false;
95 std::string ortho(
"DGKS");
98 CommandLineProcessor cmdp(
false,
true);
99 cmdp.setOption(
"verbose",
"quiet",&verbose,
"Print messages and results.");
100 cmdp.setOption(
"pseudo",
"regular",&pseudo,
"Use pseudo-block GMRES to solve the linear systems.");
101 cmdp.setOption(
"frequency",&frequency,
"Solvers frequency for printing residuals (#iters).");
102 cmdp.setOption(
"tol",&tol,
"Relative residual tolerance used by GMRES solver.");
103 cmdp.setOption(
"num-rhs",&numrhs,
"Number of right-hand sides to be solved for.");
104 cmdp.setOption(
"num-restarts",&maxrestarts,
"Maximum number of restarts allowed for the GMRES solver.");
105 cmdp.setOption(
"blocksize",&blocksize,
"Block size used by GMRES.");
106 cmdp.setOption(
"subspace-length",&length,
"Maximum dimension of block-subspace used by GMRES solver.");
107 cmdp.setOption(
"ortho-type",&ortho,
"Orthogonalization type, either DGKS, ICGS or IMGS");
108 if (cmdp.parse(argc,argv) != CommandLineProcessor::PARSE_SUCCESSFUL) {
112 proc_verbose = verbose && (MyPID==0);
119#ifndef HAVE_BELOS_TRIUTILS
120 std::cout <<
"This test requires Triutils. Please configure with --enable-triutils." << std::endl;
122 std::cout <<
"End Result: TEST FAILED" << std::endl;
131 std::vector<ST> diag( dim, (ST)4.0 );
132 RCP< MyOperator<ST> > A
138 int maxits = dim/blocksize;
140 ParameterList belosList;
141 belosList.set(
"Num Blocks", length );
142 belosList.set(
"Block Size", blocksize );
143 belosList.set(
"Maximum Iterations", maxits );
144 belosList.set(
"Maximum Restarts", maxrestarts );
145 belosList.set(
"Convergence Tolerance", tol );
146 belosList.set(
"Orthogonalization", ortho );
151 belosList.set(
"Output Frequency", frequency );
160 RCP<MyMultiVec<ST> > soln = rcp(
new MyMultiVec<ST>(dim,numrhs) );
162 MVT::MvInit( *rhs, 1.0 );
163 MVT::MvInit( *soln, zero );
167 RCP<Belos::LinearProblem<ST,MV,OP> > problem =
169 bool set = problem->setProblem();
172 std::cout << std::endl <<
"ERROR: Belos::LinearProblem failed to set up correctly!" << std::endl;
185 Teuchos::RCP< Belos::SolverManager<ST,MV,OP> > solver;
191 solver->setDebugStatusTest( Teuchos::rcp(&debugTest,
false) );
197 std::cout << std::endl << std::endl;
198 std::cout <<
"Dimension of matrix: " << dim << std::endl;
199 std::cout <<
"Number of right-hand sides: " << numrhs << std::endl;
200 std::cout <<
"Block size used by solver: " << blocksize << std::endl;
201 std::cout <<
"Max number of Gmres iterations: " << maxits << std::endl;
202 std::cout <<
"Relative residual tolerance: " << tol << std::endl;
203 std::cout << std::endl;
212 RCP<MyMultiVec<ST> > temp = rcp(
new MyMultiVec<ST>(dim,numrhs) );
213 OPT::Apply( *A, *soln, *temp );
214 MVT::MvAddMv( one, *rhs, -one, *temp, *temp );
215 std::vector<MT> norm_num(numrhs), norm_denom(numrhs);
216 MVT::MvNorm( *temp, norm_num );
217 MVT::MvNorm( *rhs, norm_denom );
218 for (
int i=0; i<numrhs; ++i) {
220 std::cout <<
"Relative residual "<<i<<
" : " << norm_num[i] / norm_denom[i] << std::endl;
221 if ( norm_num[i] / norm_denom[i] > tol ) {
227 const std::vector<MT> residualLog = debugTest.
getLogResNorm();
228 if (numrhs==1 && proc_verbose && residualLog.size())
230 std::cout <<
"Absolute residual 2-norm [ " << residualLog.size() <<
" ] : ";
231 for (
unsigned int i=0; i<residualLog.size(); i++)
232 std::cout << residualLog[i] <<
" ";
233 std::cout << std::endl;
234 std::cout <<
"Final abs 2-norm / rhs 2-norm : " << residualLog[residualLog.size()-1] / norm_denom[0] << std::endl;
240 std::cout <<
"End Result: TEST PASSED" << std::endl;
243 std::cout <<
"End Result: TEST FAILED" << std::endl;
246 TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success);
248 return ( success ? EXIT_SUCCESS : EXIT_FAILURE );
The Belos::BlockGmresSolMgr provides a solver manager for the BlockGmres linear solver.
Belos header file which uses auto-configuration information to include necessary C++ headers.
Class which describes the linear problem to be solved by the iterative solver.
The Belos::PseudoBlockGmresSolMgr provides a solver manager for the BlockGmres linear solver.
Belos::StatusTest debugging class for storing the absolute residual norms generated during a solve.
Interface to Block GMRES and Flexible GMRES.
A linear system to solve, and its associated information.
Traits class which defines basic operations on multivectors.
Interface for multivectors used by Belos' linear solvers.
Class which defines basic traits for the operator type.
Alternative run-time polymorphic interface for operators.
Interface to standard and "pseudoblock" GMRES.
A Belos::StatusTest debugging class for storing the absolute residual norms generated during a solve.
const std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > & getLogResNorm() const
Returns the log of the absolute residual norm from the iteration.
Simple example of a user's defined Belos::MultiVec class.
Simple example of a user's defined Belos::Operator class.
ReturnType
Whether the Belos solve converged for all linear systems.
std::string Belos_Version()
int main(int argc, char *argv[])