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"
71using namespace Teuchos;
73int 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();
96 bool norm_failure =
false;
97 bool proc_verbose =
false;
102 int maxrestarts = 15;
104 std::string filename(
"mhd1280b.cua");
105 std::string ortho(
"ICGS");
108 CommandLineProcessor cmdp(
false,
true);
109 cmdp.setOption(
"verbose",
"quiet",&verbose,
"Print messages and results.");
110 cmdp.setOption(
"debug",
"no-debug",&debug,
"Print debug messages.");
111 cmdp.setOption(
"pseudo",
"regular",&pseudo,
"Use pseudo-block GMRES to solve the linear systems.");
112 cmdp.setOption(
"frequency",&frequency,
"Solvers frequency for printing residuals (#iters).");
113 cmdp.setOption(
"filename",&filename,
"Filename for Harwell-Boeing test matrix.");
114 cmdp.setOption(
"tol",&tol,
"Relative residual tolerance used by GMRES solver.");
115 cmdp.setOption(
"num-rhs",&numrhs,
"Number of right-hand sides to be solved for.");
116 cmdp.setOption(
"num-restarts",&maxrestarts,
"Maximum number of restarts allowed for the GMRES solver.");
117 cmdp.setOption(
"blocksize",&blocksize,
"Block size used by GMRES.");
118 cmdp.setOption(
"subspace-length",&length,
"Maximum dimension of block-subspace used by GMRES solver.");
119 cmdp.setOption(
"ortho-type",&ortho,
"Orthogonalization type, either DGKS, ICGS or IMGS");
120 if (cmdp.parse(argc,argv) != CommandLineProcessor::PARSE_SUCCESSFUL) {
124 proc_verbose = verbose && (MyPID==0);
131#ifndef HAVE_BELOS_TRIUTILS
132 std::cout <<
"This test requires Triutils. Please configure with --enable-triutils." << std::endl;
134 std::cout <<
"End Result: TEST FAILED" << std::endl;
145 info = readHB_newmat_double(filename.c_str(),&dim,&dim2,&nnz,
146 &colptr,&rowind,&dvals);
147 if (info == 0 || nnz < 0) {
149 std::cout <<
"Error reading '" << filename <<
"'" << std::endl;
150 std::cout <<
"End Result: TEST FAILED" << std::endl;
156 for (
int ii=0; ii<nnz; ii++) {
157 cvals[ii] = ST(dvals[ii*2],dvals[ii*2+1]);
160 RCP< MyBetterOperator<ST> > A
166 int maxits = dim/blocksize;
168 ParameterList belosList;
169 belosList.set(
"Num Blocks", length );
170 belosList.set(
"Block Size", blocksize );
171 belosList.set(
"Maximum Iterations", maxits );
172 belosList.set(
"Maximum Restarts", maxrestarts );
173 belosList.set(
"Convergence Tolerance", tol );
174 belosList.set(
"Orthogonalization", ortho );
180 belosList.set(
"Verbosity", verbosity );
182 belosList.set(
"Output Frequency", frequency );
191 RCP<MyMultiVec<ST> > soln = rcp(
new MyMultiVec<ST>(dim,numrhs) );
193 MVT::MvRandom( *soln );
194 OPT::Apply( *A, *soln, *rhs );
195 MVT::MvInit( *soln, zero );
199 RCP<Belos::LinearProblem<ST,MV,OP> > problem =
201 bool set = problem->setProblem();
204 std::cout << std::endl <<
"ERROR: Belos::LinearProblem failed to set up correctly!" << std::endl;
217 Teuchos::RCP< Belos::SolverManager<ST,MV,OP> > solver;
223 solver->setDebugStatusTest( Teuchos::rcp(&debugTest,
false) );
229 std::cout << std::endl << std::endl;
230 std::cout <<
"Dimension of matrix: " << dim << std::endl;
231 std::cout <<
"Number of right-hand sides: " << numrhs << std::endl;
232 std::cout <<
"Block size used by solver: " << blocksize << std::endl;
233 std::cout <<
"Max number of Gmres iterations: " << maxits << std::endl;
234 std::cout <<
"Relative residual tolerance: " << tol << std::endl;
235 std::cout << std::endl;
244 RCP<MyMultiVec<ST> > temp = rcp(
new MyMultiVec<ST>(dim,numrhs) );
245 OPT::Apply( *A, *soln, *temp );
246 MVT::MvAddMv( one, *rhs, -one, *temp, *temp );
247 std::vector<MT> norm_num(numrhs), norm_denom(numrhs);
248 MVT::MvNorm( *temp, norm_num );
249 MVT::MvNorm( *rhs, norm_denom );
250 for (
int i=0; i<numrhs; ++i) {
252 std::cout <<
"Relative residual "<<i<<
" : " << norm_num[i] / norm_denom[i] << std::endl;
253 if ( norm_num[i] / norm_denom[i] > tol ) {
259 const std::vector<MT> residualLog = debugTest.
getLogResNorm();
260 if (numrhs==1 && proc_verbose && residualLog.size())
262 std::cout <<
"Absolute residual 2-norm [ " << residualLog.size() <<
" ] : ";
263 for (
unsigned int i=0; i<residualLog.size(); i++)
264 std::cout << residualLog[i] <<
" ";
265 std::cout << std::endl;
266 std::cout <<
"Final abs 2-norm / rhs 2-norm : " << residualLog[residualLog.size()-1] / norm_denom[0] << std::endl;
278 std::cout <<
"End Result: TEST PASSED" << std::endl;
281 std::cout <<
"End Result: TEST FAILED" << std::endl;
284 TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success);
286 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::Operator class.
Simple example of a user's defined Belos::MultiVec class.
ReturnType
Whether the Belos solve converged for all linear systems.
std::string Belos_Version()
int main(int argc, char *argv[])