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"
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 #ifdef HAVE_COMPLEX
76  typedef std::complex<double> ST;
77 #elif HAVE_COMPLEX_H
78  typedef std::complex<double> ST;
79 #else
80  std::cout << "Not compiled with std::complex support." << std::endl;
81  std::cout << "End Result: TEST FAILED" << std::endl;
82  return EXIT_FAILURE;
83 #endif
84 
85  typedef ScalarTraits<ST> SCT;
86  typedef SCT::magnitudeType MT;
87  typedef Belos::MultiVec<ST> MV;
88  typedef Belos::Operator<ST> OP;
89  typedef Belos::MultiVecTraits<ST,MV> MVT;
91  ST one = SCT::one();
92  ST zero = SCT::zero();
93 
94  Teuchos::GlobalMPISession session(&argc, &argv, NULL);
95  int MyPID = session.getRank();
96  //
97  using Teuchos::RCP;
98  using Teuchos::rcp;
99 
100  bool success = false;
101  bool verbose = false;
102  try {
103  int info = 0;
104  bool norm_failure = false;
105  bool proc_verbose = false;
106  bool pseudo = false; // use pseudo block GMRES to solve this linear system.
107  int frequency = -1; // how often residuals are printed by solver
108  int blocksize = 1;
109  int numrhs = 1;
110  int maxrestarts = 15;
111  int length = 50;
112  std::string filename("mhd1280b.cua");
113  MT tol = 1.0e-5; // relative residual tolerance
114 
115  CommandLineProcessor cmdp(false,true);
116  cmdp.setOption("verbose","quiet",&verbose,"Print messages and results.");
117  cmdp.setOption("pseudo","regular",&pseudo,"Use pseudo-block GMRES to solve the linear systems.");
118  cmdp.setOption("frequency",&frequency,"Solvers frequency for printing residuals (#iters).");
119  cmdp.setOption("filename",&filename,"Filename for Harwell-Boeing test matrix.");
120  cmdp.setOption("tol",&tol,"Relative residual tolerance used by GMRES solver.");
121  cmdp.setOption("num-rhs",&numrhs,"Number of right-hand sides to be solved for.");
122  cmdp.setOption("num-restarts",&maxrestarts,"Maximum number of restarts allowed for the GMRES solver.");
123  cmdp.setOption("blocksize",&blocksize,"Block size used by GMRES.");
124  cmdp.setOption("subspace-length",&length,"Maximum dimension of block-subspace used by GMRES solver.");
125  if (cmdp.parse(argc,argv) != CommandLineProcessor::PARSE_SUCCESSFUL) {
126  return EXIT_FAILURE;
127  }
128 
129  proc_verbose = verbose && (MyPID==0); /* Only print on the zero processor */
130  if (proc_verbose) {
131  std::cout << Belos::Belos_Version() << std::endl << std::endl;
132  }
133  if (!verbose)
134  frequency = -1; // reset frequency if test is not verbose
135 
136 #ifndef HAVE_BELOS_TRIUTILS
137  std::cout << "This test requires Triutils. Please configure with --enable-triutils." << std::endl;
138  if (MyPID==0) {
139  std::cout << "End Result: TEST FAILED" << std::endl;
140  }
141  return EXIT_FAILURE;
142 #endif
143 
144  // Get the data from the HB file
145  int dim,dim2,nnz;
146  MT *dvals;
147  int *colptr,*rowind;
148  ST *cvals;
149  nnz = -1;
150  info = readHB_newmat_double(filename.c_str(),&dim,&dim2,&nnz,
151  &colptr,&rowind,&dvals);
152  if (info == 0 || nnz < 0) {
153  if (MyPID==0) {
154  std::cout << "Error reading '" << filename << "'" << std::endl;
155  std::cout << "End Result: TEST FAILED" << std::endl;
156  }
157  return EXIT_FAILURE;
158  }
159  // Convert interleaved doubles to std::complex values
160  cvals = new ST[nnz];
161  for (int ii=0; ii<nnz; ii++) {
162  cvals[ii] = ST(dvals[ii*2],dvals[ii*2+1]);
163  }
164  // Build the problem matrix
166  = rcp( new MyBetterOperator<ST>(dim,colptr,nnz,rowind,cvals) );
167  //
168  // ********Other information used by block solver***********
169  // *****************(can be user specified)******************
170  //
171  int maxits = dim/blocksize; // maximum number of iterations to run
172  //
173  ParameterList belosList;
174  belosList.set( "Num Blocks", length ); // Maximum number of blocks in Krylov factorization
175  belosList.set( "Block Size", blocksize ); // Blocksize to be used by iterative solver
176  belosList.set( "Maximum Iterations", maxits ); // Maximum number of iterations allowed
177  belosList.set( "Maximum Restarts", maxrestarts ); // Maximum number of restarts allowed
178  belosList.set( "Convergence Tolerance", tol ); // Relative convergence tolerance requested
179  if (verbose) {
180  belosList.set( "Verbosity", Belos::Errors + Belos::Warnings +
182  if (frequency > 0)
183  belosList.set( "Output Frequency", frequency );
184  }
185  else
186  belosList.set( "Verbosity", Belos::Errors + Belos::Warnings );
187  //
188  // Construct the right-hand side and solution multivectors.
189  // NOTE: The right-hand side will be constructed such that the solution is
190  // a vectors of one.
191  //
192  RCP<MyMultiVec<ST> > soln = rcp( new MyMultiVec<ST>(dim,numrhs) );
193  RCP<MyMultiVec<ST> > rhs = rcp( new MyMultiVec<ST>(dim,numrhs) );
194  MVT::MvRandom( *soln );
195  OPT::Apply( *A, *soln, *rhs );
196  MVT::MvInit( *soln, zero );
197  //
198  // Construct an unpreconditioned linear problem instance.
199  //
201  rcp( new Belos::LinearProblem<ST,MV,OP>( A, soln, rhs ) );
202  bool set = problem->setProblem();
203  if (set == false) {
204  if (proc_verbose)
205  std::cout << std::endl << "ERROR: Belos::LinearProblem failed to set up correctly!" << std::endl;
206  return EXIT_FAILURE;
207  }
208 
209  // Use a debugging status test to save absolute residual history.
210  // Debugging status tests are peer to the native status tests that are called whenever convergence is checked.
212 
213  //
214  // *******************************************************************
215  // *************Start the block Gmres iteration***********************
216  // *******************************************************************
217  //
219  if (pseudo)
220  solver = Teuchos::rcp( new Belos::PseudoBlockGmresSolMgr<ST,MV,OP>( problem, Teuchos::rcp(&belosList,false) ) );
221  else
222  solver = Teuchos::rcp( new Belos::BlockGmresSolMgr<ST,MV,OP>( problem, Teuchos::rcp(&belosList,false) ) );
223 
224  solver->setDebugStatusTest( Teuchos::rcp(&debugTest, false) );
225 
226  //
227  // **********Print out information about problem*******************
228  //
229  if (proc_verbose) {
230  std::cout << std::endl << std::endl;
231  std::cout << "Dimension of matrix: " << dim << std::endl;
232  std::cout << "Number of right-hand sides: " << numrhs << std::endl;
233  std::cout << "Block size used by solver: " << blocksize << std::endl;
234  std::cout << "Max number of Gmres iterations: " << maxits << std::endl;
235  std::cout << "Relative residual tolerance: " << tol << std::endl;
236  std::cout << std::endl;
237  }
238  //
239  // Perform solve
240  //
241  Belos::ReturnType ret = solver->solve();
242  //
243  // Compute actual residuals.
244  //
245  RCP<MyMultiVec<ST> > temp = rcp( new MyMultiVec<ST>(dim,numrhs) );
246  OPT::Apply( *A, *soln, *temp );
247  MVT::MvAddMv( one, *rhs, -one, *temp, *temp );
248  std::vector<MT> norm_num(numrhs), norm_denom(numrhs);
249  MVT::MvNorm( *temp, norm_num );
250  MVT::MvNorm( *rhs, norm_denom );
251  for (int i=0; i<numrhs; ++i) {
252  if (proc_verbose)
253  std::cout << "Relative residual "<<i<<" : " << norm_num[i] / norm_denom[i] << std::endl;
254  if ( norm_num[i] / norm_denom[i] > tol ) {
255  norm_failure = true;
256  }
257  }
258 
259  // Print absolute residual norm logging.
260  const std::vector<MT> residualLog = debugTest.getLogResNorm();
261  if (numrhs==1 && proc_verbose && residualLog.size())
262  {
263  std::cout << "Absolute residual 2-norm [ " << residualLog.size() << " ] : ";
264  for (unsigned int i=0; i<residualLog.size(); i++)
265  std::cout << residualLog[i] << " ";
266  std::cout << std::endl;
267  std::cout << "Final abs 2-norm / rhs 2-norm : " << residualLog[residualLog.size()-1] / norm_denom[0] << std::endl;
268  }
269 
270  // Clean up.
271  delete [] dvals;
272  delete [] colptr;
273  delete [] rowind;
274  delete [] cvals;
275 
276  success = ret==Belos::Converged && !norm_failure;
277  if (success) {
278  if (proc_verbose)
279  std::cout << "End Result: TEST PASSED" << std::endl;
280  } else {
281  if (proc_verbose)
282  std::cout << "End Result: TEST FAILED" << std::endl;
283  }
284  }
285  TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success);
286 
287  return ( success ? EXIT_SUCCESS : EXIT_FAILURE );
288 } // end test_bl_gmres_complex_hb.cpp
std::string Belos_Version()
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
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:62
std::string filename
Alternative run-time polymorphic interface for operators.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
The Belos::BlockGmresSolMgr provides a solver manager for the BlockGmres linear solver.
void setOption(const char option_true[], const char option_false[], bool *option_val, const char documentation[]=NULL)
EParseCommandLineReturn parse(int argc, char *argv[], std::ostream *errout=&std::cerr) const
#define TEUCHOS_STANDARD_CATCH_STATEMENTS(VERBOSE, ERR_STREAM, SUCCESS_FLAG)
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: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...
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.