Belos Package Browser (Single Doxygen Collection)  Development
test_bl_cg_complex_hb.cpp
Go to the documentation of this file.
1 //@HEADER
2 // ************************************************************************
3 //
4 // Belos: Block Linear Solvers Package
5 // Copyright 2004 Sandia Corporation
6 //
7 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ************************************************************************
40 //@HEADER
41 //
42 // This driver reads a problem from a Harwell-Boeing (HB) file.
43 // The right-hand-side from the HB file is used instead of random vectors.
44 // The initial guesses are all set to zero.
45 //
46 // NOTE: No preconditioner is used in this case.
47 //
48 #include "BelosConfigDefs.hpp"
49 #include "BelosLinearProblem.hpp"
50 #include "BelosBlockCGSolMgr.hpp"
52 #include "Teuchos_CommandLineProcessor.hpp"
53 #include "Teuchos_ParameterList.hpp"
54 #include "Teuchos_StandardCatchMacros.hpp"
55 
56 #ifdef HAVE_MPI
57 #include <mpi.h>
58 #endif
59 
60 // I/O for Harwell-Boeing files
61 #ifdef HAVE_BELOS_TRIUTILS
62 #include "Trilinos_Util_iohb.h"
63 #endif
64 
65 #include "MyMultiVec.hpp"
66 #include "MyBetterOperator.hpp"
67 #include "MyOperator.hpp"
68 
69 using namespace Teuchos;
70 
71 int main(int argc, char *argv[]) {
72  //
73 #ifdef HAVE_COMPLEX
74  typedef std::complex<double> ST;
75 #elif HAVE_COMPLEX_H
76  typedef std::complex<double> ST;
77 #else
78  std::cout << "Not compiled with std::complex support." << std::endl;
79  std::cout << "End Result: TEST FAILED" << std::endl;
80  return EXIT_FAILURE;
81 #endif
82 
83  typedef ScalarTraits<ST> SCT;
84  typedef SCT::magnitudeType MT;
85  typedef Belos::MultiVec<ST> MV;
86  typedef Belos::Operator<ST> OP;
87  typedef Belos::MultiVecTraits<ST,MV> MVT;
89  ST one = SCT::one();
90  ST zero = SCT::zero();
91 
92  Teuchos::GlobalMPISession session(&argc, &argv, NULL);
93  //
94  using Teuchos::RCP;
95  using Teuchos::rcp;
96 
97  bool success = false;
98  bool verbose = false;
99  try {
100  int info = 0;
101  int MyPID = 0;
102  bool pseudo = false; // use pseudo block CG to solve this linear system.
103  bool norm_failure = false;
104  bool proc_verbose = false;
105  bool use_single_red = false;
106  int frequency = -1; // how often residuals are printed by solver
107  int blocksize = 1;
108  int numrhs = 1;
109  std::string filename("mhd1280b.cua");
110  MT tol = 1.0e-5; // relative residual tolerance
111 
112  CommandLineProcessor cmdp(false,true);
113  cmdp.setOption("verbose","quiet",&verbose,"Print messages and results.");
114  cmdp.setOption("pseudo","regular",&pseudo,"Use pseudo-block CG to solve the linear systems.");
115  cmdp.setOption("frequency",&frequency,"Solvers frequency for printing residuals (#iters).");
116  cmdp.setOption("filename",&filename,"Filename for Harwell-Boeing test matrix.");
117  cmdp.setOption("tol",&tol,"Relative residual tolerance used by CG solver.");
118  cmdp.setOption("num-rhs",&numrhs,"Number of right-hand sides to be solved for.");
119  cmdp.setOption("blocksize",&blocksize,"Block size used by CG .");
120  cmdp.setOption("use-single-red","use-standard-red",&use_single_red,"Use single-reduction variant of CG iteration.");
121  if (cmdp.parse(argc,argv) != CommandLineProcessor::PARSE_SUCCESSFUL) {
122  return -1;
123  }
124 
125  proc_verbose = verbose && (MyPID==0); /* Only print on the zero processor */
126  if (proc_verbose) {
127  std::cout << Belos::Belos_Version() << std::endl << std::endl;
128  }
129  if (!verbose)
130  frequency = -1; // reset frequency if test is not verbose
131 
132 
133 #ifndef HAVE_BELOS_TRIUTILS
134  std::cout << "This test requires Triutils. Please configure with --enable-triutils." << std::endl;
135  if (MyPID==0) {
136  std::cout << "End Result: TEST FAILED" << std::endl;
137  }
138  return -1;
139 #endif
140 
141  // Get the data from the HB file
142  int dim,dim2,nnz;
143  MT *dvals;
144  int *colptr,*rowind;
145  ST *cvals;
146  nnz = -1;
147  info = readHB_newmat_double(filename.c_str(),&dim,&dim2,&nnz,
148  &colptr,&rowind,&dvals);
149  if (info == 0 || nnz < 0) {
150  if (MyPID==0) {
151  std::cout << "Error reading '" << filename << "'" << std::endl;
152  std::cout << "End Result: TEST FAILED" << std::endl;
153  }
154  return -1;
155  }
156  // Convert interleaved doubles to std::complex values
157  cvals = new ST[nnz];
158  for (int ii=0; ii<nnz; ii++) {
159  cvals[ii] = ST(dvals[ii*2],dvals[ii*2+1]);
160  }
161  // Build the problem matrix
162  RCP< MyBetterOperator<ST> > A
163  = rcp( new MyBetterOperator<ST>(dim,colptr,nnz,rowind,cvals) );
164  //
165  // ********Other information used by block solver***********
166  // *****************(can be user specified)******************
167  //
168  int maxits = dim/blocksize; // maximum number of iterations to run
169  //
170  ParameterList belosList;
171  belosList.set( "Block Size", blocksize ); // Blocksize to be used by iterative solver
172  belosList.set( "Maximum Iterations", maxits ); // Maximum number of iterations allowed
173  belosList.set( "Convergence Tolerance", tol ); // Relative convergence tolerance requested
174  if ((blocksize == 1) && use_single_red)
175  belosList.set( "Use Single Reduction", use_single_red ); // Use single reduction CG iteration
176  if (verbose) {
177  belosList.set( "Verbosity", Belos::Errors + Belos::Warnings +
179  if (frequency > 0)
180  belosList.set( "Output Frequency", frequency );
181  }
182  else
183  belosList.set( "Verbosity", Belos::Errors + Belos::Warnings );
184  //
185  // Construct the right-hand side and solution multivectors.
186  // NOTE: The right-hand side will be constructed such that the solution is
187  // a vectors of one.
188  //
189  RCP<MyMultiVec<ST> > soln = rcp( new MyMultiVec<ST>(dim,numrhs) );
190  RCP<MyMultiVec<ST> > rhs = rcp( new MyMultiVec<ST>(dim,numrhs) );
191  MVT::MvRandom( *soln );
192  OPT::Apply( *A, *soln, *rhs );
193  MVT::MvInit( *soln, zero );
194  //
195  // Construct an unpreconditioned linear problem instance.
196  //
197  RCP<Belos::LinearProblem<ST,MV,OP> > problem =
198  rcp( new Belos::LinearProblem<ST,MV,OP>( A, soln, rhs ) );
199  bool set = problem->setProblem();
200  if (set == false) {
201  if (proc_verbose)
202  std::cout << std::endl << "ERROR: Belos::LinearProblem failed to set up correctly!" << std::endl;
203  return -1;
204  }
205 
206  //
207  // *******************************************************************
208  // *************Start the block CG iteration***********************
209  // *******************************************************************
210  //
211  Teuchos::RCP< Belos::SolverManager<ST,MV,OP> > solver;
212  if (pseudo)
213  solver = Teuchos::rcp( new Belos::PseudoBlockCGSolMgr<ST,MV,OP>( problem, Teuchos::rcp(&belosList,false) ) );
214  else
215  solver = Teuchos::rcp( new Belos::BlockCGSolMgr<ST,MV,OP>( problem, Teuchos::rcp(&belosList,false) ) );
216 
217  //
218  // **********Print out information about problem*******************
219  //
220  if (proc_verbose) {
221  std::cout << std::endl << std::endl;
222  std::cout << "Dimension of matrix: " << dim << std::endl;
223  std::cout << "Number of right-hand sides: " << numrhs << std::endl;
224  std::cout << "Block size used by solver: " << blocksize << std::endl;
225  std::cout << "Max number of CG iterations: " << maxits << std::endl;
226  std::cout << "Relative residual tolerance: " << tol << std::endl;
227  std::cout << std::endl;
228  }
229  //
230  // Perform solve
231  //
232  Belos::ReturnType ret = solver->solve();
233  //
234  // Compute actual residuals.
235  //
236  RCP<MyMultiVec<ST> > temp = rcp( new MyMultiVec<ST>(dim,numrhs) );
237  OPT::Apply( *A, *soln, *temp );
238  MVT::MvAddMv( one, *rhs, -one, *temp, *temp );
239  std::vector<MT> norm_num(numrhs), norm_denom(numrhs);
240  MVT::MvNorm( *temp, norm_num );
241  MVT::MvNorm( *rhs, norm_denom );
242  for (int i=0; i<numrhs; ++i) {
243  if (proc_verbose)
244  std::cout << "Relative residual "<<i<<" : " << norm_num[i] / norm_denom[i] << std::endl;
245  if ( norm_num[i] / norm_denom[i] > tol ) {
246  norm_failure = true;
247  }
248  }
249 
250  // Test achievedTol output
251  MT ach_tol = solver->achievedTol();
252  if (proc_verbose)
253  std::cout << "Achieved tol : "<<ach_tol<<std::endl;
254 
255  // Clean up.
256  delete [] dvals;
257  delete [] colptr;
258  delete [] rowind;
259  delete [] cvals;
260 
261  success = ret==Belos::Converged && !norm_failure;
262 
263  if (success) {
264  if (proc_verbose)
265  std::cout << "End Result: TEST PASSED" << std::endl;
266  } else {
267  if (proc_verbose)
268  std::cout << "End Result: TEST FAILED" << std::endl;
269  }
270  }
271  TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success);
272 
273  return ( success ? EXIT_SUCCESS : EXIT_FAILURE );
274 } // end test_bl_cg_complex_hb.cpp
std::string Belos_Version()
int main(int argc, char *argv[])
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...
Traits class which defines basic operations on multivectors.
The Belos::PseudoBlockCGSolMgr provides a powerful and fully-featured solver manager over the pseudo-...
Simple example of a user&#39;s defined Belos::MultiVec class.
Definition: MyMultiVec.hpp:62
Alternative run-time polymorphic interface for operators.
A linear system to solve, and its associated information.
const double tol
Class which describes the linear problem to be solved by the iterative solver.
ReturnType
Whether the Belos solve converged for all linear systems.
Definition: BelosTypes.hpp:154
The Belos::BlockCGSolMgr provides a solver manager for the BlockCG linear solver. ...
Interface for multivectors used by Belos&#39; linear solvers.
Class which defines basic traits for the operator type.
Belos header file which uses auto-configuration information to include necessary C++ headers...
Simple example of a user&#39;s defined Belos::Operator class.