Belos Package Browser (Single Doxygen Collection)  Development
test_bl_gmres_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"
52 #include "Teuchos_CommandLineProcessor.hpp"
53 #include "Teuchos_ParameterList.hpp"
54 #include "Teuchos_StandardCatchMacros.hpp"
55 #include "Teuchos_StandardCatchMacros.hpp"
56 
57 #ifdef HAVE_MPI
58 #include <mpi.h>
59 #endif
60 
61 // I/O for Harwell-Boeing files
62 #ifdef HAVE_BELOS_TRIUTILS
63 #include "Trilinos_Util_iohb.h"
64 #endif
65 
66 #include "MyMultiVec.hpp"
67 #include "MyBetterOperator.hpp"
68 #include "MyOperator.hpp"
69 
70 using namespace Teuchos;
71 
72 int main(int argc, char *argv[]) {
73  //
74 #ifdef HAVE_COMPLEX
75  typedef std::complex<double> ST;
76 #elif HAVE_COMPLEX_H
77  typedef std::complex<double> ST;
78 #else
79  std::cout << "Not compiled with std::complex support." << std::endl;
80  std::cout << "End Result: TEST FAILED" << std::endl;
81  return EXIT_FAILURE;
82 #endif
83 
84  typedef ScalarTraits<ST> SCT;
85  typedef SCT::magnitudeType MT;
86  typedef Belos::MultiVec<ST> MV;
87  typedef Belos::Operator<ST> OP;
88  typedef Belos::MultiVecTraits<ST,MV> MVT;
90  ST one = SCT::one();
91  ST zero = SCT::zero();
92 
93  Teuchos::GlobalMPISession session(&argc, &argv, NULL);
94  int MyPID = session.getRank();
95  //
96  using Teuchos::RCP;
97  using Teuchos::rcp;
98 
99  bool success = false;
100  bool verbose = false;
101  try {
102  int info = 0;
103  bool norm_failure = false;
104  bool proc_verbose = false;
105  bool pseudo = false; // use pseudo block GMRES to solve this linear system.
106  int frequency = -1; // how often residuals are printed by solver
107  int blocksize = 1;
108  int numrhs = 1;
109  int maxrestarts = 15;
110  int length = 50;
111  std::string filename("mhd1280b.cua");
112  MT tol = 1.0e-5; // relative residual tolerance
113 
114  CommandLineProcessor cmdp(false,true);
115  cmdp.setOption("verbose","quiet",&verbose,"Print messages and results.");
116  cmdp.setOption("pseudo","regular",&pseudo,"Use pseudo-block GMRES to solve the linear systems.");
117  cmdp.setOption("frequency",&frequency,"Solvers frequency for printing residuals (#iters).");
118  cmdp.setOption("filename",&filename,"Filename for Harwell-Boeing test matrix.");
119  cmdp.setOption("tol",&tol,"Relative residual tolerance used by GMRES solver.");
120  cmdp.setOption("num-rhs",&numrhs,"Number of right-hand sides to be solved for.");
121  cmdp.setOption("num-restarts",&maxrestarts,"Maximum number of restarts allowed for the GMRES solver.");
122  cmdp.setOption("blocksize",&blocksize,"Block size used by GMRES.");
123  cmdp.setOption("subspace-length",&length,"Maximum dimension of block-subspace used by GMRES solver.");
124  if (cmdp.parse(argc,argv) != CommandLineProcessor::PARSE_SUCCESSFUL) {
125  return EXIT_FAILURE;
126  }
127 
128  proc_verbose = verbose && (MyPID==0); /* Only print on the zero processor */
129  if (proc_verbose) {
130  std::cout << Belos::Belos_Version() << std::endl << std::endl;
131  }
132  if (!verbose)
133  frequency = -1; // reset frequency if test is not verbose
134 
135 #ifndef HAVE_BELOS_TRIUTILS
136  std::cout << "This test requires Triutils. Please configure with --enable-triutils." << std::endl;
137  if (MyPID==0) {
138  std::cout << "End Result: TEST FAILED" << std::endl;
139  }
140  return EXIT_FAILURE;
141 #endif
142 
143  // Get the data from the HB file
144  int dim,dim2,nnz;
145  MT *dvals;
146  int *colptr,*rowind;
147  ST *cvals;
148  nnz = -1;
149  info = readHB_newmat_double(filename.c_str(),&dim,&dim2,&nnz,
150  &colptr,&rowind,&dvals);
151  if (info == 0 || nnz < 0) {
152  if (MyPID==0) {
153  std::cout << "Error reading '" << filename << "'" << std::endl;
154  std::cout << "End Result: TEST FAILED" << std::endl;
155  }
156  return EXIT_FAILURE;
157  }
158  // Convert interleaved doubles to std::complex values
159  cvals = new ST[nnz];
160  for (int ii=0; ii<nnz; ii++) {
161  cvals[ii] = ST(dvals[ii*2],dvals[ii*2+1]);
162  }
163  // Build the problem matrix
164  RCP< MyBetterOperator<ST> > A
165  = rcp( new MyBetterOperator<ST>(dim,colptr,nnz,rowind,cvals) );
166  //
167  // ********Other information used by block solver***********
168  // *****************(can be user specified)******************
169  //
170  int maxits = dim/blocksize; // maximum number of iterations to run
171  //
172  ParameterList belosList;
173  belosList.set( "Num Blocks", length ); // Maximum number of blocks in Krylov factorization
174  belosList.set( "Block Size", blocksize ); // Blocksize to be used by iterative solver
175  belosList.set( "Maximum Iterations", maxits ); // Maximum number of iterations allowed
176  belosList.set( "Maximum Restarts", maxrestarts ); // Maximum number of restarts allowed
177  belosList.set( "Convergence Tolerance", tol ); // Relative convergence tolerance requested
178  if (verbose) {
179  belosList.set( "Verbosity", Belos::Errors + Belos::Warnings +
181  if (frequency > 0)
182  belosList.set( "Output Frequency", frequency );
183  }
184  else
185  belosList.set( "Verbosity", Belos::Errors + Belos::Warnings );
186  //
187  // Construct the right-hand side and solution multivectors.
188  // NOTE: The right-hand side will be constructed such that the solution is
189  // a vectors of one.
190  //
191  RCP<MyMultiVec<ST> > soln = rcp( new MyMultiVec<ST>(dim,numrhs) );
192  RCP<MyMultiVec<ST> > rhs = rcp( new MyMultiVec<ST>(dim,numrhs) );
193  MVT::MvRandom( *soln );
194  OPT::Apply( *A, *soln, *rhs );
195  MVT::MvInit( *soln, zero );
196  //
197  // Construct an unpreconditioned linear problem instance.
198  //
199  RCP<Belos::LinearProblem<ST,MV,OP> > problem =
200  rcp( new Belos::LinearProblem<ST,MV,OP>( A, soln, rhs ) );
201  bool set = problem->setProblem();
202  if (set == false) {
203  if (proc_verbose)
204  std::cout << std::endl << "ERROR: Belos::LinearProblem failed to set up correctly!" << std::endl;
205  return EXIT_FAILURE;
206  }
207  //
208  // *******************************************************************
209  // *************Start the block Gmres iteration***********************
210  // *******************************************************************
211  //
212  Teuchos::RCP< Belos::SolverManager<ST,MV,OP> > solver;
213  if (pseudo)
214  solver = Teuchos::rcp( new Belos::PseudoBlockGmresSolMgr<ST,MV,OP>( problem, Teuchos::rcp(&belosList,false) ) );
215  else
216  solver = Teuchos::rcp( new Belos::BlockGmresSolMgr<ST,MV,OP>( problem, Teuchos::rcp(&belosList,false) ) );
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 Gmres 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  // Clean up.
251  delete [] dvals;
252  delete [] colptr;
253  delete [] rowind;
254  delete [] cvals;
255 
256  success = ret==Belos::Converged && !norm_failure;
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_gmres_complex_hb.cpp
std::string Belos_Version()
Interface to Block GMRES and Flexible GMRES.
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.
The Belos::BlockGmresSolMgr provides a solver manager for the BlockGmres linear solver.
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.
Interface to standard and "pseudoblock" GMRES.
ReturnType
Whether the Belos solve converged for all linear systems.
Definition: BelosTypes.hpp:154
Interface for multivectors used by Belos&#39; linear solvers.
Class which defines basic traits for the operator type.
The Belos::PseudoBlockGmresSolMgr provides a solver manager for the BlockGmres linear solver...
Belos header file which uses auto-configuration information to include necessary C++ headers...
int main(int argc, char *argv[])
Simple example of a user&#39;s defined Belos::Operator class.