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"
51 #include "Teuchos_CommandLineProcessor.hpp"
52 #include "Teuchos_ParameterList.hpp"
53 #include "Teuchos_StandardCatchMacros.hpp"
54 
55 #ifdef HAVE_MPI
56 #include <mpi.h>
57 #endif
58 
59 // I/O for Harwell-Boeing files
60 #ifdef HAVE_BELOS_TRIUTILS
61 #include "Trilinos_Util_iohb.h"
62 #endif
63 
64 #include "MyMultiVec.hpp"
65 #include "MyBetterOperator.hpp"
66 #include "MyOperator.hpp"
67 
68 using namespace Teuchos;
69 
70 int main(int argc, char *argv[]) {
71  //
72 #ifdef HAVE_COMPLEX
73  typedef std::complex<double> ST;
74 #elif HAVE_COMPLEX_H
75  typedef std::complex<double> ST;
76 #else
77  std::cout << "Not compiled with std::complex support." << std::endl;
78  std::cout << "End Result: TEST FAILED" << std::endl;
79  return EXIT_FAILURE;
80 #endif
81 
82  typedef ScalarTraits<ST> SCT;
83  typedef SCT::magnitudeType MT;
84  typedef Belos::MultiVec<ST> MV;
85  typedef Belos::Operator<ST> OP;
86  typedef Belos::MultiVecTraits<ST,MV> MVT;
88  ST one = SCT::one();
89  ST zero = SCT::zero();
90 
91  Teuchos::GlobalMPISession session(&argc, &argv, NULL);
92  //
93  using Teuchos::RCP;
94  using Teuchos::rcp;
95 
96  bool success = false;
97  bool verbose = false;
98  try {
99  int info = 0;
100  int MyPID = 0;
101  bool norm_failure = false;
102  bool proc_verbose = false;
103  int frequency = -1; // how often residuals are printed by solver
104  int blocksize = 1;
105  int numrhs = 1;
106  std::string filename("mhd1280b.cua");
107  MT tol = 1.0e-5; // relative residual tolerance
108 
109  CommandLineProcessor cmdp(false,true);
110  cmdp.setOption("verbose","quiet",&verbose,"Print messages and results.");
111  cmdp.setOption("frequency",&frequency,"Solvers frequency for printing residuals (#iters).");
112  cmdp.setOption("filename",&filename,"Filename for Harwell-Boeing test matrix.");
113  cmdp.setOption("tol",&tol,"Relative residual tolerance used by CG solver.");
114  cmdp.setOption("num-rhs",&numrhs,"Number of right-hand sides to be solved for.");
115  cmdp.setOption("blocksize",&blocksize,"Block size used by CG .");
116  if (cmdp.parse(argc,argv) != CommandLineProcessor::PARSE_SUCCESSFUL) {
117  return -1;
118  }
119 
120  proc_verbose = verbose && (MyPID==0); /* Only print on the zero processor */
121  if (proc_verbose) {
122  std::cout << Belos::Belos_Version() << std::endl << std::endl;
123  }
124  if (!verbose)
125  frequency = -1; // reset frequency if test is not verbose
126 
127 
128 #ifndef HAVE_BELOS_TRIUTILS
129  std::cout << "This test requires Triutils. Please configure with --enable-triutils." << std::endl;
130  if (MyPID==0) {
131  std::cout << "End Result: TEST FAILED" << std::endl;
132  }
133  return -1;
134 #endif
135 
136  // Get the data from the HB file
137  int dim,dim2,nnz;
138  MT *dvals;
139  int *colptr,*rowind;
140  ST *cvals;
141  nnz = -1;
142  info = readHB_newmat_double(filename.c_str(),&dim,&dim2,&nnz,
143  &colptr,&rowind,&dvals);
144  if (info == 0 || nnz < 0) {
145  if (MyPID==0) {
146  std::cout << "Error reading '" << filename << "'" << std::endl;
147  std::cout << "End Result: TEST FAILED" << std::endl;
148  }
149  return -1;
150  }
151  // Convert interleaved doubles to std::complex values
152  cvals = new ST[nnz];
153  for (int ii=0; ii<nnz; ii++) {
154  cvals[ii] = ST(dvals[ii*2],dvals[ii*2+1]);
155  }
156  // Build the problem matrix
157  RCP< MyBetterOperator<ST> > A
158  = rcp( new MyBetterOperator<ST>(dim,colptr,nnz,rowind,cvals) );
159  //
160  // ********Other information used by block solver***********
161  // *****************(can be user specified)******************
162  //
163  int maxits = dim/blocksize; // maximum number of iterations to run
164  //
165  ParameterList belosList;
166  belosList.set( "Block Size", blocksize ); // Blocksize to be used by iterative solver
167  belosList.set( "Maximum Iterations", maxits ); // Maximum number of iterations allowed
168  belosList.set( "Convergence Tolerance", tol ); // Relative convergence tolerance requested
169  if (verbose) {
170  belosList.set( "Verbosity", Belos::Errors + Belos::Warnings +
172  if (frequency > 0)
173  belosList.set( "Output Frequency", frequency );
174  }
175  else
176  belosList.set( "Verbosity", Belos::Errors + Belos::Warnings );
177  //
178  // Construct the right-hand side and solution multivectors.
179  // NOTE: The right-hand side will be constructed such that the solution is
180  // a vectors of one.
181  //
182  RCP<MyMultiVec<ST> > soln = rcp( new MyMultiVec<ST>(dim,numrhs) );
183  RCP<MyMultiVec<ST> > rhs = rcp( new MyMultiVec<ST>(dim,numrhs) );
184  MVT::MvRandom( *soln );
185  OPT::Apply( *A, *soln, *rhs );
186  MVT::MvInit( *soln, zero );
187  //
188  // Construct an unpreconditioned linear problem instance.
189  //
190  RCP<Belos::LinearProblem<ST,MV,OP> > problem =
191  rcp( new Belos::LinearProblem<ST,MV,OP>( A, soln, rhs ) );
192  bool set = problem->setProblem();
193  if (set == false) {
194  if (proc_verbose)
195  std::cout << std::endl << "ERROR: Belos::LinearProblem failed to set up correctly!" << std::endl;
196  return -1;
197  }
198  //
199  // *******************************************************************
200  // *************Start the block CG iteration***********************
201  // *******************************************************************
202  //
203  if (proc_verbose) {
204  std::cout << "Attempt to create Belos::BlockCGSolMgr" << std::endl;
205  }
206  Belos::BlockCGSolMgr<ST,MV,OP> solver( problem, rcp(&belosList,false) );
207 
208  //
209  // **********Print out information about problem*******************
210  //
211  if (proc_verbose) {
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;
219  }
220  //
221  // Perform solve
222  //
223  Belos::ReturnType ret = solver.solve();
224  //
225  // Compute actual residuals.
226  //
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) {
234  if (proc_verbose)
235  std::cout << "Relative residual "<<i<<" : " << norm_num[i] / norm_denom[i] << std::endl;
236  if ( norm_num[i] / norm_denom[i] > tol ) {
237  norm_failure = true;
238  }
239  }
240 
241  // Test achievedTol output
242  MT ach_tol = solver.achievedTol();
243  if (proc_verbose)
244  std::cout << "Achieved tol : "<<ach_tol<<std::endl;
245 
246  // Clean up.
247  delete [] dvals;
248  delete [] colptr;
249  delete [] rowind;
250  delete [] cvals;
251 
252  success = ret==Belos::Converged && !norm_failure;
253 
254  if (success) {
255  if (proc_verbose)
256  std::cout << "End Result: TEST PASSED" << std::endl;
257  } else {
258  if (proc_verbose)
259  std::cout << "End Result: TEST FAILED" << std::endl;
260  }
261  }
262  TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success);
263 
264  return ( success ? EXIT_SUCCESS : EXIT_FAILURE );
265 } // end test_bl_cg_complex_hb.cpp
std::string Belos_Version()
int main(int argc, char *argv[])
The Belos::BlockCGSolMgr provides a powerful and fully-featured solver manager over the CG and BlockC...
Traits class which defines basic operations on multivectors.
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.