Belos Package Browser (Single Doxygen Collection)  Development
test_hybrid_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"
50 #include "BelosGmresPolySolMgr.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 userandomrhs = true;
106  int frequency = -1; // frequency of status test output.
107  int blocksize = 1; // blocksize
108  int numrhs = 1; // number of right-hand sides to solve for
109  int maxiters = -1; // maximum number of iterations allowed per linear system
110  int maxdegree = 25; // maximum degree of polynomial
111  int maxsubspace = 50; // maximum number of blocks the solver can use for the subspace
112  int maxrestarts = 15; // number of restarts allowed
113  std::string outersolver("Block Gmres");
114  std::string filename("mhd1280b.cua");
115  std::string polyType("Arnoldi");
116  MT tol = 1.0e-5; // relative residual tolerance
117  MT polytol = tol/10; // relative residual tolerance for polynomial construction
118 
119  Teuchos::CommandLineProcessor cmdp(false,true);
120  cmdp.setOption("verbose","quiet",&verbose,"Print messages and results.");
121  cmdp.setOption("use-random-rhs","use-rhs",&userandomrhs,"Use linear system RHS or random RHS to generate polynomial.");
122  cmdp.setOption("frequency",&frequency,"Solvers frequency for printing residuals (#iters).");
123  cmdp.setOption("filename",&filename,"Filename for test matrix. Acceptable file extensions: *.hb,*.mtx,*.triU,*.triS");
124  cmdp.setOption("outersolver",&outersolver,"Name of outer solver to be used with GMRES poly");
125  cmdp.setOption("tol",&tol,"Relative residual tolerance used by GMRES solver.");
126  cmdp.setOption("poly-tol",&polytol,"Relative residual tolerance used to construct the GMRES polynomial.");
127  cmdp.setOption("poly-type",&polyType,"Polynomial type (Roots, Arnoldi, or Gmres).");
128  cmdp.setOption("num-rhs",&numrhs,"Number of right-hand sides to be solved for.");
129  cmdp.setOption("block-size",&blocksize,"Block size used by GMRES.");
130  cmdp.setOption("max-iters",&maxiters,"Maximum number of iterations per linear system (-1 = adapted to problem/block size).");
131  cmdp.setOption("max-degree",&maxdegree,"Maximum degree of the GMRES polynomial.");
132  cmdp.setOption("max-subspace",&maxsubspace,"Maximum number of blocks the solver can use for the subspace.");
133  cmdp.setOption("max-restarts",&maxrestarts,"Maximum number of restarts allowed for GMRES solver.");
134  if (cmdp.parse(argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
135  return EXIT_FAILURE;
136  }
137 
138  proc_verbose = verbose && (MyPID==0); /* Only print on the zero processor */
139  if (proc_verbose) {
140  std::cout << Belos::Belos_Version() << std::endl << std::endl;
141  }
142  if (!verbose)
143  frequency = -1; // reset frequency if test is not verbose
144 
145 #ifndef HAVE_BELOS_TRIUTILS
146  std::cout << "This test requires Triutils. Please configure with --enable-triutils." << std::endl;
147  if (MyPID==0) {
148  std::cout << "End Result: TEST FAILED" << std::endl;
149  }
150  return EXIT_FAILURE;
151 #endif
152 
153  // Get the data from the HB file
154  int dim,dim2,nnz;
155  MT *dvals;
156  int *colptr,*rowind;
157  ST *cvals;
158  nnz = -1;
159  info = readHB_newmat_double(filename.c_str(),&dim,&dim2,&nnz,
160  &colptr,&rowind,&dvals);
161  if (info == 0 || nnz < 0) {
162  if (MyPID==0) {
163  std::cout << "Error reading '" << filename << "'" << std::endl;
164  std::cout << "End Result: TEST FAILED" << std::endl;
165  }
166  return EXIT_FAILURE;
167  }
168  // Convert interleaved doubles to std::complex values
169  cvals = new ST[nnz];
170  for (int ii=0; ii<nnz; ii++) {
171  cvals[ii] = ST(dvals[ii*2],dvals[ii*2+1]);
172  }
173  // Build the problem matrix
174  RCP< MyBetterOperator<ST> > A
175  = rcp( new MyBetterOperator<ST>(dim,colptr,nnz,rowind,cvals) );
176  //
177  // Construct the right-hand side and solution multivectors.
178  // NOTE: The right-hand side will be constructed such that the solution is
179  // a vectors of one.
180  //
181  RCP<MyMultiVec<ST> > soln = rcp( new MyMultiVec<ST>(dim,numrhs) );
182  RCP<MyMultiVec<ST> > rhs = rcp( new MyMultiVec<ST>(dim,numrhs) );
183  MVT::MvRandom( *soln );
184  OPT::Apply( *A, *soln, *rhs );
185  MVT::MvInit( *soln, zero );
186  //
187  // Construct an unpreconditioned linear problem instance.
188  //
189  RCP<Belos::LinearProblem<ST,MV,OP> > problem =
190  rcp( new Belos::LinearProblem<ST,MV,OP>( A, soln, rhs ) );
191  problem->setInitResVec( 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 EXIT_FAILURE;
197  }
198  //
199  // ********Other information used by block solver***********
200  // *****************(can be user specified)******************
201  //
202  if (maxiters == -1)
203  maxiters = dim/blocksize - 1; // maximum number of iterations to run
204 
205  ParameterList belosList;
206  belosList.set( "Num Blocks", maxsubspace); // Maximum number of blocks in Krylov factorization
207  belosList.set( "Block Size", blocksize ); // Blocksize to be used by iterative solver
208  belosList.set( "Maximum Iterations", maxiters ); // Maximum number of iterations allowed
209  belosList.set( "Maximum Restarts", maxrestarts ); // Maximum number of restarts allowed
210  belosList.set( "Convergence Tolerance", tol ); // Relative convergence tolerance requested
211  int verbosity = Belos::Errors + Belos::Warnings;
212  if (verbose) {
214  if (frequency > 0)
215  belosList.set( "Output Frequency", frequency );
216  }
217  belosList.set( "Verbosity", verbosity );
218 
219  ParameterList polyList;
220  polyList.set( "Maximum Degree", maxdegree ); // Maximum degree of the GMRES polynomial
221  polyList.set( "Polynomial Tolerance", polytol ); // Polynomial convergence tolerance requested
222  polyList.set( "Polynomial Type", polyType ); // Type of polynomial to construct
223  polyList.set( "Verbosity", verbosity ); // Verbosity for polynomial construction
224  polyList.set( "Random RHS", userandomrhs ); // Use RHS from linear system or random vector
225  if ( outersolver != "" ) {
226  polyList.set( "Outer Solver", outersolver );
227  polyList.set( "Outer Solver Params", belosList );
228  }
229 
230  // Use a debugging status test to save absolute residual history.
231  // Debugging status tests are peer to the native status tests that are called whenever convergence is checked.
233 
234  //
235  // *******************************************************************
236  // *************Start the block Gmres iteration***********************
237  // *******************************************************************
238  //
239  RCP< Belos::SolverManager<ST,MV,OP> > solver = rcp( new Belos::GmresPolySolMgr<ST,MV,OP>( problem, rcp(&polyList,false) ) );
240 
241  // The debug status test does not work for the GmresPolySolMgr right now.
242  // solver->setDebugStatusTest( Teuchos::rcp(&debugTest, false) );
243 
244  //
245  // **********Print out information about problem*******************
246  //
247  if (proc_verbose) {
248  std::cout << std::endl << std::endl;
249  std::cout << "Dimension of matrix: " << dim << std::endl;
250  std::cout << "Number of right-hand sides: " << numrhs << std::endl;
251  std::cout << "Block size used by solver: " << blocksize << std::endl;
252  std::cout << "Max number of Gmres iterations: " << maxiters << std::endl;
253  std::cout << "Relative residual tolerance: " << tol << std::endl;
254  std::cout << std::endl;
255  }
256  //
257  // Perform solve
258  //
259  Belos::ReturnType ret = solver->solve();
260  //
261  // Compute actual residuals.
262  //
263  RCP<MyMultiVec<ST> > temp = rcp( new MyMultiVec<ST>(dim,numrhs) );
264  OPT::Apply( *A, *soln, *temp );
265  MVT::MvAddMv( one, *rhs, -one, *temp, *temp );
266  std::vector<MT> norm_num(numrhs), norm_denom(numrhs);
267  MVT::MvNorm( *temp, norm_num );
268  MVT::MvNorm( *rhs, norm_denom );
269  for (int i=0; i<numrhs; ++i) {
270  if (proc_verbose)
271  std::cout << "Relative residual "<<i<<" : " << norm_num[i] / norm_denom[i] << std::endl;
272  if ( norm_num[i] / norm_denom[i] > tol ) {
273  norm_failure = true;
274  }
275  }
276 
277  // Print absolute residual norm logging.
278  const std::vector<MT> residualLog = debugTest.getLogResNorm();
279  if (numrhs==1 && proc_verbose && residualLog.size())
280  {
281  std::cout << "Absolute residual 2-norm [ " << residualLog.size() << " ] : ";
282  for (unsigned int i=0; i<residualLog.size(); i++)
283  std::cout << residualLog[i] << " ";
284  std::cout << std::endl;
285  std::cout << "Final abs 2-norm / rhs 2-norm : " << residualLog[residualLog.size()-1] / norm_denom[0] << std::endl;
286  }
287 
288  // Clean up.
289  delete [] dvals;
290  delete [] colptr;
291  delete [] rowind;
292  delete [] cvals;
293 
294  success = ret==Belos::Converged && !norm_failure;
295  if (success) {
296  if (proc_verbose)
297  std::cout << "End Result: TEST PASSED" << std::endl;
298  } else {
299  if (proc_verbose)
300  std::cout << "End Result: TEST FAILED" << std::endl;
301  }
302  }
303  TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success);
304 
305  return ( success ? EXIT_SUCCESS : EXIT_FAILURE );
306 } // 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...
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
int main(int argc, char *argv[])
Alternative run-time polymorphic interface for operators.
Declaration and definition of Belos::GmresPolySolMgr (hybrid block GMRES linear solver).
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
The GMRES polynomial can be created in conjunction with any standard preconditioner.
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.
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...
Simple example of a user&#39;s defined Belos::Operator class.