52#include "Teuchos_CommandLineProcessor.hpp"
53#include "Teuchos_ParameterList.hpp"
54#include "Teuchos_StandardCatchMacros.hpp"
61#ifdef HAVE_BELOS_TRIUTILS
62#include "Trilinos_Util_iohb.h"
69using namespace Teuchos;
71int main(
int argc,
char *argv[]) {
73 typedef std::complex<double> ST;
74 typedef ScalarTraits<ST> SCT;
75 typedef SCT::magnitudeType MT;
81 ST zero = SCT::zero();
83 Teuchos::GlobalMPISession session(&argc, &argv, NULL);
94 bool norm_failure =
false;
95 bool proc_verbose =
false;
96 bool use_single_red =
false;
100 std::string filename(
"mhd1280b.cua");
103 CommandLineProcessor cmdp(
false,
true);
104 cmdp.setOption(
"verbose",
"quiet",&verbose,
"Print messages and results.");
105 cmdp.setOption(
"pseudo",
"regular",&pseudo,
"Use pseudo-block CG to solve the linear systems.");
106 cmdp.setOption(
"frequency",&frequency,
"Solvers frequency for printing residuals (#iters).");
107 cmdp.setOption(
"filename",&filename,
"Filename for Harwell-Boeing test matrix.");
108 cmdp.setOption(
"tol",&tol,
"Relative residual tolerance used by CG solver.");
109 cmdp.setOption(
"num-rhs",&numrhs,
"Number of right-hand sides to be solved for.");
110 cmdp.setOption(
"blocksize",&blocksize,
"Block size used by CG .");
111 cmdp.setOption(
"use-single-red",
"use-standard-red",&use_single_red,
"Use single-reduction variant of CG iteration.");
112 if (cmdp.parse(argc,argv) != CommandLineProcessor::PARSE_SUCCESSFUL) {
116 proc_verbose = verbose && (MyPID==0);
124#ifndef HAVE_BELOS_TRIUTILS
125 std::cout <<
"This test requires Triutils. Please configure with --enable-triutils." << std::endl;
127 std::cout <<
"End Result: TEST FAILED" << std::endl;
138 info = readHB_newmat_double(filename.c_str(),&dim,&dim2,&nnz,
139 &colptr,&rowind,&dvals);
140 if (info == 0 || nnz < 0) {
142 std::cout <<
"Error reading '" << filename <<
"'" << std::endl;
143 std::cout <<
"End Result: TEST FAILED" << std::endl;
149 for (
int ii=0; ii<nnz; ii++) {
150 cvals[ii] = ST(dvals[ii*2],dvals[ii*2+1]);
153 RCP< MyBetterOperator<ST> > A
159 int maxits = dim/blocksize;
161 ParameterList belosList;
162 belosList.set(
"Block Size", blocksize );
163 belosList.set(
"Maximum Iterations", maxits );
164 belosList.set(
"Convergence Tolerance", tol );
165 if ((blocksize == 1) && use_single_red)
166 belosList.set(
"Use Single Reduction", use_single_red );
171 belosList.set(
"Output Frequency", frequency );
180 RCP<MyMultiVec<ST> > soln = rcp(
new MyMultiVec<ST>(dim,numrhs) );
182 MVT::MvRandom( *soln );
183 OPT::Apply( *A, *soln, *rhs );
184 MVT::MvInit( *soln, zero );
188 RCP<Belos::LinearProblem<ST,MV,OP> > problem =
190 bool set = problem->setProblem();
193 std::cout << std::endl <<
"ERROR: Belos::LinearProblem failed to set up correctly!" << std::endl;
202 Teuchos::RCP< Belos::SolverManager<ST,MV,OP> > solver;
212 std::cout << std::endl << std::endl;
213 std::cout <<
"Dimension of matrix: " << dim << std::endl;
214 std::cout <<
"Number of right-hand sides: " << numrhs << std::endl;
215 std::cout <<
"Block size used by solver: " << blocksize << std::endl;
216 std::cout <<
"Max number of CG iterations: " << maxits << std::endl;
217 std::cout <<
"Relative residual tolerance: " << tol << std::endl;
218 std::cout << std::endl;
227 RCP<MyMultiVec<ST> > temp = rcp(
new MyMultiVec<ST>(dim,numrhs) );
228 OPT::Apply( *A, *soln, *temp );
229 MVT::MvAddMv( one, *rhs, -one, *temp, *temp );
230 std::vector<MT> norm_num(numrhs), norm_denom(numrhs);
231 MVT::MvNorm( *temp, norm_num );
232 MVT::MvNorm( *rhs, norm_denom );
233 for (
int i=0; i<numrhs; ++i) {
235 std::cout <<
"Relative residual "<<i<<
" : " << norm_num[i] / norm_denom[i] << std::endl;
236 if ( norm_num[i] / norm_denom[i] > tol ) {
242 MT ach_tol = solver->achievedTol();
244 std::cout <<
"Achieved tol : "<<ach_tol<<std::endl;
256 std::cout <<
"End Result: TEST PASSED" << std::endl;
259 std::cout <<
"End Result: TEST FAILED" << std::endl;
262 TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success);
264 return ( success ? EXIT_SUCCESS : EXIT_FAILURE );
The Belos::BlockCGSolMgr provides a solver manager for the BlockCG 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::PseudoBlockCGSolMgr provides a solver manager for the BlockCG linear solver.
The Belos::BlockCGSolMgr provides a powerful and fully-featured solver manager over the CG and BlockC...
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::PseudoBlockCGSolMgr provides a powerful and fully-featured solver manager over the pseudo-...
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[])