51#include "Teuchos_CommandLineProcessor.hpp"
52#include "Teuchos_ParameterList.hpp"
53#include "Teuchos_StandardCatchMacros.hpp"
60#ifdef HAVE_BELOS_TRIUTILS
61#include "Trilinos_Util_iohb.h"
68using namespace Teuchos;
70int main(
int argc,
char *argv[]) {
72 typedef std::complex<double> ST;
73 typedef ScalarTraits<ST> SCT;
74 typedef SCT::magnitudeType MT;
80 ST zero = SCT::zero();
83 bool norm_failure =
false;
85 Teuchos::GlobalMPISession session(&argc, &argv, NULL);
87 int MyPID = session.getRank();
95 bool proc_verbose =
false;
99 std::string filename(
"mhd1280b.cua");
102 CommandLineProcessor cmdp(
false,
true);
103 cmdp.setOption(
"verbose",
"quiet",&verbose,
"Print messages and results.");
104 cmdp.setOption(
"frequency",&frequency,
"Solvers frequency for printing residuals (#iters).");
105 cmdp.setOption(
"filename",&filename,
"Filename for Harwell-Boeing test matrix.");
106 cmdp.setOption(
"tol",&tol,
"Relative residual tolerance used by pseudo-block TFQMR solver.");
107 cmdp.setOption(
"num-rhs",&numrhs,
"Number of right-hand sides to be solved for.");
108 if (cmdp.parse(argc,argv) != CommandLineProcessor::PARSE_SUCCESSFUL) {
112 proc_verbose = verbose && (MyPID==0);
120#ifndef HAVE_BELOS_TRIUTILS
121 std::cout <<
"This test requires Triutils. Please configure with --enable-triutils." << std::endl;
123 std::cout <<
"End Result: TEST FAILED" << std::endl;
134 info = readHB_newmat_double(filename.c_str(),&dim,&dim2,&nnz,
135 &colptr,&rowind,&dvals);
136 if (info == 0 || nnz < 0) {
138 std::cout <<
"Error reading '" << filename <<
"'" << std::endl;
139 std::cout <<
"End Result: TEST FAILED" << std::endl;
145 for (
int ii=0; ii<nnz; ii++) {
146 cvals[ii] = ST(dvals[ii*2],dvals[ii*2+1]);
149 RCP< MyBetterOperator<ST> > A
155 int maxits = dim/blocksize;
157 ParameterList belosList;
158 belosList.set(
"Maximum Iterations", maxits );
159 belosList.set(
"Convergence Tolerance", tol );
164 belosList.set(
"Output Frequency", frequency );
173 RCP<MyMultiVec<ST> > soln = rcp(
new MyMultiVec<ST>(dim,numrhs) );
175 MVT::MvRandom( *soln );
176 OPT::Apply( *A, *soln, *rhs );
177 MVT::MvInit( *soln, zero );
181 RCP<Belos::LinearProblem<ST,MV,OP> > problem =
183 bool set = problem->setProblem();
186 std::cout << std::endl <<
"ERROR: Belos::LinearProblem failed to set up correctly!" << std::endl;
200 std::cout << std::endl << std::endl;
201 std::cout <<
"Dimension of matrix: " << dim << std::endl;
202 std::cout <<
"Number of right-hand sides: " << numrhs << std::endl;
203 std::cout <<
"Block size used by solver: " << blocksize << std::endl;
204 std::cout <<
"Max number of pseudo-block TFQMR iterations: " << maxits << std::endl;
205 std::cout <<
"Relative residual tolerance: " << tol << std::endl;
206 std::cout << std::endl;
215 RCP<MyMultiVec<ST> > temp = rcp(
new MyMultiVec<ST>(dim,numrhs) );
216 OPT::Apply( *A, *soln, *temp );
217 MVT::MvAddMv( one, *rhs, -one, *temp, *temp );
218 std::vector<MT> norm_num(numrhs), norm_denom(numrhs);
219 MVT::MvNorm( *temp, norm_num );
220 MVT::MvNorm( *rhs, norm_denom );
221 for (
int i=0; i<numrhs; ++i) {
223 std::cout <<
"Relative residual "<<i<<
" : " << norm_num[i] / norm_denom[i] << std::endl;
224 if ( norm_num[i] / norm_denom[i] > tol ) {
238 std::cout <<
"End Result: TEST PASSED" << std::endl;
241 std::cout <<
"End Result: TEST FAILED" << std::endl;
244 TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success);
246 return ( success ? EXIT_SUCCESS : EXIT_FAILURE );
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::PseudoBlockTFQMRSolMgr provides a solver manager for the pseudo-block TFQMR linear solver.
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.
The Belos::PseudoBlockTFQMRSolMgr provides a powerful and fully-featured solver manager over the pseu...
ReturnType solve() override
This method performs possibly repeated calls to the underlying linear solver's iterate() routine unti...
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[])