Belos Package Browser (Single Doxygen Collection)  Development
test_bl_cg_real_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  typedef double ST;
74 
75  typedef ScalarTraits<ST> SCT;
76  typedef SCT::magnitudeType MT;
77  typedef Belos::MultiVec<ST> MV;
78  typedef Belos::Operator<ST> OP;
79  typedef Belos::MultiVecTraits<ST,MV> MVT;
81  ST one = SCT::one();
82  ST zero = SCT::zero();
83 
84  Teuchos::GlobalMPISession session(&argc, &argv, NULL);
85  //
86  using Teuchos::RCP;
87  using Teuchos::rcp;
88 
89  bool success = false;
90  bool verbose = false;
91  try {
92  int info = 0;
93  int MyPID = 0;
94  bool pseudo = false; // use pseudo block CG to solve this linear system.
95  bool norm_failure = false;
96  bool proc_verbose = false;
97  bool use_single_red = false;
98  bool combineConvInner = false;
99  int frequency = -1; // how often residuals are printed by solver
100  int blocksize = 1;
101  int numrhs = 1;
102  std::string filename("A.hb");
103  MT tol = 1.0e-5; // relative residual tolerance
104 
105  CommandLineProcessor cmdp(false,true);
106  cmdp.setOption("verbose","quiet",&verbose,"Print messages and results.");
107  cmdp.setOption("pseudo","regular",&pseudo,"Use pseudo-block CG to solve the linear systems.");
108  cmdp.setOption("frequency",&frequency,"Solvers frequency for printing residuals (#iters).");
109  cmdp.setOption("filename",&filename,"Filename for Harwell-Boeing test matrix.");
110  cmdp.setOption("tol",&tol,"Relative residual tolerance used by CG solver.");
111  cmdp.setOption("num-rhs",&numrhs,"Number of right-hand sides to be solved for.");
112  cmdp.setOption("blocksize",&blocksize,"Block size used by CG .");
113  cmdp.setOption("use-single-red","use-standard-red",&use_single_red,"Use single-reduction variant of CG iteration.");
114  cmdp.setOption("combine-conv-inner","separate-conv-inner",&combineConvInner,"Combine convergence detection and inner products of single-reduce CG iteration.");
115  if (cmdp.parse(argc,argv) != CommandLineProcessor::PARSE_SUCCESSFUL) {
116  return -1;
117  }
118 
119  proc_verbose = verbose && (MyPID==0); /* Only print on the zero processor */
120  if (proc_verbose) {
121  std::cout << Belos::Belos_Version() << std::endl << std::endl;
122  }
123  if (!verbose)
124  frequency = -1; // reset frequency if test is not verbose
125 
126 
127 #ifndef HAVE_BELOS_TRIUTILS
128  std::cout << "This test requires Triutils. Please configure with --enable-triutils." << std::endl;
129  if (MyPID==0) {
130  std::cout << "End Result: TEST FAILED" << std::endl;
131  }
132  return -1;
133 #endif
134 
135  // Get the data from the HB file
136  int dim,dim2,nnz;
137  int *colptr,*rowind;
138  ST *cvals;
139  nnz = -1;
140  info = readHB_newmat_double(filename.c_str(),&dim,&dim2,&nnz,
141  &colptr,&rowind,&cvals);
142  if (info == 0 || nnz < 0) {
143  if (MyPID==0) {
144  std::cout << "Error reading '" << filename << "'" << std::endl;
145  std::cout << "End Result: TEST FAILED" << std::endl;
146  }
147  return -1;
148  }
149  // Build the problem matrix
150  RCP< MyBetterOperator<ST> > A
151  = rcp( new MyBetterOperator<ST>(dim,colptr,nnz,rowind,cvals) );
152  // for (int j=0; j<nnz; j++)
153  // std::cout << cvals[j] << std::endl;
154  // A->Print(std::cout);
155  //
156  // ********Other information used by block solver***********
157  // *****************(can be user specified)******************
158  //
159  int maxits = dim/blocksize; // maximum number of iterations to run
160  //
161  ParameterList belosList;
162  belosList.set( "Block Size", blocksize ); // Blocksize to be used by iterative solver
163  belosList.set( "Maximum Iterations", maxits ); // Maximum number of iterations allowed
164  belosList.set( "Convergence Tolerance", tol ); // Relative convergence tolerance requested
165  if (blocksize == 1) {
166  if (use_single_red)
167  belosList.set( "Use Single Reduction", use_single_red ); // Use single reduction CG iteration
168  if (combineConvInner)
169  belosList.set( "Fold Convergence Detection Into Allreduce", combineConvInner );
170  }
171  if (verbose) {
172  belosList.set( "Verbosity", Belos::Errors + Belos::Warnings +
174  if (frequency > 0)
175  belosList.set( "Output Frequency", frequency );
176  }
177  else
178  belosList.set( "Verbosity", Belos::Errors + Belos::Warnings );
179  //
180  // Construct the right-hand side and solution multivectors.
181  // NOTE: The right-hand side will be constructed such that the solution is
182  // a vectors of one.
183  //
184  RCP<MyMultiVec<ST> > soln = rcp( new MyMultiVec<ST>(dim,numrhs) );
185  RCP<MyMultiVec<ST> > rhs = rcp( new MyMultiVec<ST>(dim,numrhs) );
186  MVT::MvRandom( *soln );
187  OPT::Apply( *A, *soln, *rhs );
188  MVT::MvInit( *soln, zero );
189  //
190  // Construct an unpreconditioned linear problem instance.
191  //
192  RCP<Belos::LinearProblem<ST,MV,OP> > problem =
193  rcp( new Belos::LinearProblem<ST,MV,OP>( A, soln, rhs ) );
194  bool set = problem->setProblem();
195  if (set == false) {
196  if (proc_verbose)
197  std::cout << std::endl << "ERROR: Belos::LinearProblem failed to set up correctly!" << std::endl;
198  return -1;
199  }
200 
201  //
202  // *******************************************************************
203  // *************Start the block CG iteration***********************
204  // *******************************************************************
205  //
206  Teuchos::RCP< Belos::SolverManager<ST,MV,OP> > solver;
207  if (pseudo)
208  solver = Teuchos::rcp( new Belos::PseudoBlockCGSolMgr<ST,MV,OP>( problem, Teuchos::rcp(&belosList,false) ) );
209  else
210  solver = Teuchos::rcp( new Belos::BlockCGSolMgr<ST,MV,OP>( problem, Teuchos::rcp(&belosList,false) ) );
211 
212  //
213  // **********Print out information about problem*******************
214  //
215  if (proc_verbose) {
216  std::cout << std::endl << std::endl;
217  std::cout << "Dimension of matrix: " << dim << std::endl;
218  std::cout << "Number of right-hand sides: " << numrhs << std::endl;
219  std::cout << "Block size used by solver: " << blocksize << std::endl;
220  std::cout << "Max number of CG iterations: " << maxits << std::endl;
221  std::cout << "Relative residual tolerance: " << tol << std::endl;
222  std::cout << std::endl;
223  }
224  //
225  // Perform solve
226  //
227  Belos::ReturnType ret = solver->solve();
228  //
229  // Compute actual residuals.
230  //
231  RCP<MyMultiVec<ST> > temp = rcp( new MyMultiVec<ST>(dim,numrhs) );
232  OPT::Apply( *A, *soln, *temp );
233  MVT::MvAddMv( one, *rhs, -one, *temp, *temp );
234  std::vector<MT> norm_num(numrhs), norm_denom(numrhs);
235  MVT::MvNorm( *temp, norm_num );
236  MVT::MvNorm( *rhs, norm_denom );
237  for (int i=0; i<numrhs; ++i) {
238  if (proc_verbose)
239  std::cout << "Relative residual "<<i<<" : " << norm_num[i] / norm_denom[i] << std::endl;
240  if ( norm_num[i] / norm_denom[i] > tol ) {
241  norm_failure = true;
242  }
243  }
244 
245  // Test achievedTol output
246  MT ach_tol = solver->achievedTol();
247  if (proc_verbose)
248  std::cout << "Achieved tol : "<<ach_tol<<std::endl;
249 
250  // Clean up.
251  delete [] colptr;
252  delete [] rowind;
253  delete [] cvals;
254 
255  success = ret==Belos::Converged && !norm_failure;
256 
257  if (success) {
258  if (proc_verbose)
259  std::cout << "End Result: TEST PASSED" << std::endl;
260  } else {
261  if (proc_verbose)
262  std::cout << "End Result: TEST FAILED" << std::endl;
263  }
264  }
265  TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success);
266 
267  return ( success ? EXIT_SUCCESS : EXIT_FAILURE );
268 } // end test_bl_cg_complex_hb.cpp
int main(int argc, char *argv[])
std::string Belos_Version()
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:65
Alternative run-time polymorphic interface for operators.
A linear system to solve, and its associated information.
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:155
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.