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"
53 #include "Teuchos_CommandLineProcessor.hpp"
54 #include "Teuchos_ParameterList.hpp"
55 #include "Teuchos_StandardCatchMacros.hpp"
56 #include "Teuchos_StandardCatchMacros.hpp"
57 
58 #ifdef HAVE_MPI
59 #include <mpi.h>
60 #endif
61 
62 // I/O for Harwell-Boeing files
63 #ifdef HAVE_BELOS_TRIUTILS
64 #include "Trilinos_Util_iohb.h"
65 #endif
66 
67 #include "MyMultiVec.hpp"
68 #include "MyBetterOperator.hpp"
69 #include "MyOperator.hpp"
70 
71 using namespace Teuchos;
72 
73 int main(int argc, char *argv[]) {
74  //
75  typedef std::complex<double> ST;
76  typedef ScalarTraits<ST> SCT;
77  typedef SCT::magnitudeType MT;
78  typedef Belos::MultiVec<ST> MV;
79  typedef Belos::Operator<ST> OP;
80  typedef Belos::MultiVecTraits<ST,MV> MVT;
82  ST one = SCT::one();
83  ST zero = SCT::zero();
84 
85  Teuchos::GlobalMPISession session(&argc, &argv, NULL);
86  int MyPID = session.getRank();
87  //
88  using Teuchos::RCP;
89  using Teuchos::rcp;
90 
91  bool success = false;
92  bool verbose = false;
93  try {
94  int info = 0;
95  bool norm_failure = false;
96  bool proc_verbose = false;
97  bool pseudo = false; // use pseudo block GMRES to solve this linear system.
98  int frequency = -1; // how often residuals are printed by solver
99  int blocksize = 1;
100  int numrhs = 1;
101  int maxrestarts = 15;
102  int length = 50;
103  std::string filename("mhd1280b.cua");
104  MT tol = 1.0e-5; // relative residual tolerance
105 
106  CommandLineProcessor cmdp(false,true);
107  cmdp.setOption("verbose","quiet",&verbose,"Print messages and results.");
108  cmdp.setOption("pseudo","regular",&pseudo,"Use pseudo-block GMRES to solve the linear systems.");
109  cmdp.setOption("frequency",&frequency,"Solvers frequency for printing residuals (#iters).");
110  cmdp.setOption("filename",&filename,"Filename for Harwell-Boeing test matrix.");
111  cmdp.setOption("tol",&tol,"Relative residual tolerance used by GMRES solver.");
112  cmdp.setOption("num-rhs",&numrhs,"Number of right-hand sides to be solved for.");
113  cmdp.setOption("num-restarts",&maxrestarts,"Maximum number of restarts allowed for the GMRES solver.");
114  cmdp.setOption("blocksize",&blocksize,"Block size used by GMRES.");
115  cmdp.setOption("subspace-length",&length,"Maximum dimension of block-subspace used by GMRES solver.");
116  if (cmdp.parse(argc,argv) != CommandLineProcessor::PARSE_SUCCESSFUL) {
117  return EXIT_FAILURE;
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 #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 EXIT_FAILURE;
133 #endif
134 
135  // Get the data from the HB file
136  int dim,dim2,nnz;
137  MT *dvals;
138  int *colptr,*rowind;
139  ST *cvals;
140  nnz = -1;
141  info = readHB_newmat_double(filename.c_str(),&dim,&dim2,&nnz,
142  &colptr,&rowind,&dvals);
143  if (info == 0 || nnz < 0) {
144  if (MyPID==0) {
145  std::cout << "Error reading '" << filename << "'" << std::endl;
146  std::cout << "End Result: TEST FAILED" << std::endl;
147  }
148  return EXIT_FAILURE;
149  }
150  // Convert interleaved doubles to std::complex values
151  cvals = new ST[nnz];
152  for (int ii=0; ii<nnz; ii++) {
153  cvals[ii] = ST(dvals[ii*2],dvals[ii*2+1]);
154  }
155  // Build the problem matrix
156  RCP< MyBetterOperator<ST> > A
157  = rcp( new MyBetterOperator<ST>(dim,colptr,nnz,rowind,cvals) );
158  //
159  // ********Other information used by block solver***********
160  // *****************(can be user specified)******************
161  //
162  int maxits = dim/blocksize; // maximum number of iterations to run
163  //
164  ParameterList belosList;
165  belosList.set( "Num Blocks", length ); // Maximum number of blocks in Krylov factorization
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( "Maximum Restarts", maxrestarts ); // Maximum number of restarts allowed
169  belosList.set( "Convergence Tolerance", tol ); // Relative convergence tolerance requested
170  if (verbose) {
171  belosList.set( "Verbosity", Belos::Errors + Belos::Warnings +
173  if (frequency > 0)
174  belosList.set( "Output Frequency", frequency );
175  }
176  else
177  belosList.set( "Verbosity", Belos::Errors + Belos::Warnings );
178  //
179  // Construct the right-hand side and solution multivectors.
180  // NOTE: The right-hand side will be constructed such that the solution is
181  // a vectors of one.
182  //
183  RCP<MyMultiVec<ST> > soln = rcp( new MyMultiVec<ST>(dim,numrhs) );
184  RCP<MyMultiVec<ST> > rhs = rcp( new MyMultiVec<ST>(dim,numrhs) );
185  MVT::MvRandom( *soln );
186  OPT::Apply( *A, *soln, *rhs );
187  MVT::MvInit( *soln, zero );
188  //
189  // Construct an unpreconditioned linear problem instance.
190  //
191  RCP<Belos::LinearProblem<ST,MV,OP> > problem =
192  rcp( new Belos::LinearProblem<ST,MV,OP>( A, soln, rhs ) );
193  bool set = problem->setProblem();
194  if (set == false) {
195  if (proc_verbose)
196  std::cout << std::endl << "ERROR: Belos::LinearProblem failed to set up correctly!" << std::endl;
197  return EXIT_FAILURE;
198  }
199 
200  // Use a debugging status test to save absolute residual history.
201  // Debugging status tests are peer to the native status tests that are called whenever convergence is checked.
203 
204  //
205  // *******************************************************************
206  // *************Start the block Gmres iteration***********************
207  // *******************************************************************
208  //
209  Teuchos::RCP< Belos::SolverManager<ST,MV,OP> > solver;
210  if (pseudo)
211  solver = Teuchos::rcp( new Belos::PseudoBlockGmresSolMgr<ST,MV,OP>( problem, Teuchos::rcp(&belosList,false) ) );
212  else
213  solver = Teuchos::rcp( new Belos::BlockGmresSolMgr<ST,MV,OP>( problem, Teuchos::rcp(&belosList,false) ) );
214 
215  solver->setDebugStatusTest( Teuchos::rcp(&debugTest, 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 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  // Print absolute residual norm logging.
251  const std::vector<MT> residualLog = debugTest.getLogResNorm();
252  if (numrhs==1 && proc_verbose && residualLog.size())
253  {
254  std::cout << "Absolute residual 2-norm [ " << residualLog.size() << " ] : ";
255  for (unsigned int i=0; i<residualLog.size(); i++)
256  std::cout << residualLog[i] << " ";
257  std::cout << std::endl;
258  std::cout << "Final abs 2-norm / rhs 2-norm : " << residualLog[residualLog.size()-1] / norm_denom[0] << std::endl;
259  }
260 
261  // Clean up.
262  delete [] dvals;
263  delete [] colptr;
264  delete [] rowind;
265  delete [] cvals;
266 
267  success = ret==Belos::Converged && !norm_failure;
268  if (success) {
269  if (proc_verbose)
270  std::cout << "End Result: TEST PASSED" << std::endl;
271  } else {
272  if (proc_verbose)
273  std::cout << "End Result: TEST FAILED" << std::endl;
274  }
275  }
276  TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success);
277 
278  return ( success ? EXIT_SUCCESS : EXIT_FAILURE );
279 } // end test_bl_gmres_complex_hb.cpp
std::string Belos_Version()
A Belos::StatusTest debugging class for storing the absolute residual norms generated during a solve...
Interface to Block GMRES and Flexible GMRES.
Belos::StatusTest debugging class for storing the absolute residual norms generated during a solve...
Traits class which defines basic operations on multivectors.
Simple example of a user&#39;s defined Belos::MultiVec class.
Definition: MyMultiVec.hpp:65
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.
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:155
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...
const std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > & getLogResNorm() const
Returns the log of the absolute residual norm from the iteration.
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.