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  typedef std::complex<double> ST;
75  typedef ScalarTraits<ST> SCT;
76  typedef SCT::magnitudeType MT;
77  typedef Belos::MultiVec<ST> MV;
78  typedef Belos::Operator<ST> OP;
79  typedef Belos::MultiVecTraits<ST,MV> MVT;
81  ST one = SCT::one();
82  ST zero = SCT::zero();
83 
84  Teuchos::GlobalMPISession session(&argc, &argv, NULL);
85  int MyPID = session.getRank();
86  //
87  using Teuchos::RCP;
88  using Teuchos::rcp;
89 
90  bool success = false;
91  bool verbose = false;
92  try {
93  int info = 0;
94  bool norm_failure = false;
95  bool proc_verbose = false;
96  bool userandomrhs = true;
97  int frequency = -1; // frequency of status test output.
98  int blocksize = 1; // blocksize
99  int numrhs = 1; // number of right-hand sides to solve for
100  int maxiters = -1; // maximum number of iterations allowed per linear system
101  int maxdegree = 25; // maximum degree of polynomial
102  int maxsubspace = 50; // maximum number of blocks the solver can use for the subspace
103  int maxrestarts = 15; // number of restarts allowed
104  std::string outersolver("Block Gmres");
105  std::string filename("mhd1280b.cua");
106  std::string polyType("Arnoldi");
107  MT tol = 1.0e-5; // relative residual tolerance
108  MT polytol = tol/10; // relative residual tolerance for polynomial construction
109 
110  Teuchos::CommandLineProcessor cmdp(false,true);
111  cmdp.setOption("verbose","quiet",&verbose,"Print messages and results.");
112  cmdp.setOption("use-random-rhs","use-rhs",&userandomrhs,"Use linear system RHS or random RHS to generate polynomial.");
113  cmdp.setOption("frequency",&frequency,"Solvers frequency for printing residuals (#iters).");
114  cmdp.setOption("filename",&filename,"Filename for test matrix. Acceptable file extensions: *.hb,*.mtx,*.triU,*.triS");
115  cmdp.setOption("outersolver",&outersolver,"Name of outer solver to be used with GMRES poly");
116  cmdp.setOption("tol",&tol,"Relative residual tolerance used by GMRES solver.");
117  cmdp.setOption("poly-tol",&polytol,"Relative residual tolerance used to construct the GMRES polynomial.");
118  cmdp.setOption("poly-type",&polyType,"Polynomial type (Roots, Arnoldi, or Gmres).");
119  cmdp.setOption("num-rhs",&numrhs,"Number of right-hand sides to be solved for.");
120  cmdp.setOption("block-size",&blocksize,"Block size used by GMRES.");
121  cmdp.setOption("max-iters",&maxiters,"Maximum number of iterations per linear system (-1 = adapted to problem/block size).");
122  cmdp.setOption("max-degree",&maxdegree,"Maximum degree of the GMRES polynomial.");
123  cmdp.setOption("max-subspace",&maxsubspace,"Maximum number of blocks the solver can use for the subspace.");
124  cmdp.setOption("max-restarts",&maxrestarts,"Maximum number of restarts allowed for GMRES solver.");
125  if (cmdp.parse(argc,argv) != Teuchos::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
165  RCP< MyBetterOperator<ST> > A
166  = rcp( new MyBetterOperator<ST>(dim,colptr,nnz,rowind,cvals) );
167  //
168  // Construct the right-hand side and solution multivectors.
169  // NOTE: The right-hand side will be constructed such that the solution is
170  // a vectors of one.
171  //
172  RCP<MyMultiVec<ST> > soln = rcp( new MyMultiVec<ST>(dim,numrhs) );
173  RCP<MyMultiVec<ST> > rhs = rcp( new MyMultiVec<ST>(dim,numrhs) );
174  MVT::MvRandom( *soln );
175  OPT::Apply( *A, *soln, *rhs );
176  MVT::MvInit( *soln, zero );
177  //
178  // Construct an unpreconditioned linear problem instance.
179  //
180  RCP<Belos::LinearProblem<ST,MV,OP> > problem =
181  rcp( new Belos::LinearProblem<ST,MV,OP>( A, soln, rhs ) );
182  problem->setInitResVec( rhs );
183  bool set = problem->setProblem();
184  if (set == false) {
185  if (proc_verbose)
186  std::cout << std::endl << "ERROR: Belos::LinearProblem failed to set up correctly!" << std::endl;
187  return EXIT_FAILURE;
188  }
189  //
190  // ********Other information used by block solver***********
191  // *****************(can be user specified)******************
192  //
193  if (maxiters == -1)
194  maxiters = dim/blocksize - 1; // maximum number of iterations to run
195 
196  ParameterList belosList;
197  belosList.set( "Num Blocks", maxsubspace); // Maximum number of blocks in Krylov factorization
198  belosList.set( "Block Size", blocksize ); // Blocksize to be used by iterative solver
199  belosList.set( "Maximum Iterations", maxiters ); // Maximum number of iterations allowed
200  belosList.set( "Maximum Restarts", maxrestarts ); // Maximum number of restarts allowed
201  belosList.set( "Convergence Tolerance", tol ); // Relative convergence tolerance requested
202  int verbosity = Belos::Errors + Belos::Warnings;
203  if (verbose) {
205  if (frequency > 0)
206  belosList.set( "Output Frequency", frequency );
207  }
208  belosList.set( "Verbosity", verbosity );
209 
210  ParameterList polyList;
211  polyList.set( "Maximum Degree", maxdegree ); // Maximum degree of the GMRES polynomial
212  polyList.set( "Polynomial Tolerance", polytol ); // Polynomial convergence tolerance requested
213  polyList.set( "Polynomial Type", polyType ); // Type of polynomial to construct
214  polyList.set( "Verbosity", verbosity ); // Verbosity for polynomial construction
215  polyList.set( "Random RHS", userandomrhs ); // Use RHS from linear system or random vector
216  if ( outersolver != "" ) {
217  polyList.set( "Outer Solver", outersolver );
218  polyList.set( "Outer Solver Params", belosList );
219  }
220 
221  // Use a debugging status test to save absolute residual history.
222  // Debugging status tests are peer to the native status tests that are called whenever convergence is checked.
224 
225  //
226  // *******************************************************************
227  // *************Start the block Gmres iteration***********************
228  // *******************************************************************
229  //
230  RCP< Belos::SolverManager<ST,MV,OP> > solver = rcp( new Belos::GmresPolySolMgr<ST,MV,OP>( problem, rcp(&polyList,false) ) );
231 
232  // The debug status test does not work for the GmresPolySolMgr right now.
233  // solver->setDebugStatusTest( Teuchos::rcp(&debugTest, false) );
234 
235  //
236  // **********Print out information about problem*******************
237  //
238  if (proc_verbose) {
239  std::cout << std::endl << std::endl;
240  std::cout << "Dimension of matrix: " << dim << std::endl;
241  std::cout << "Number of right-hand sides: " << numrhs << std::endl;
242  std::cout << "Block size used by solver: " << blocksize << std::endl;
243  std::cout << "Max number of Gmres iterations: " << maxiters << std::endl;
244  std::cout << "Relative residual tolerance: " << tol << std::endl;
245  std::cout << std::endl;
246  }
247  //
248  // Perform solve
249  //
250  Belos::ReturnType ret = solver->solve();
251  //
252  // Compute actual residuals.
253  //
254  RCP<MyMultiVec<ST> > temp = rcp( new MyMultiVec<ST>(dim,numrhs) );
255  OPT::Apply( *A, *soln, *temp );
256  MVT::MvAddMv( one, *rhs, -one, *temp, *temp );
257  std::vector<MT> norm_num(numrhs), norm_denom(numrhs);
258  MVT::MvNorm( *temp, norm_num );
259  MVT::MvNorm( *rhs, norm_denom );
260  for (int i=0; i<numrhs; ++i) {
261  if (proc_verbose)
262  std::cout << "Relative residual "<<i<<" : " << norm_num[i] / norm_denom[i] << std::endl;
263  if ( norm_num[i] / norm_denom[i] > tol ) {
264  norm_failure = true;
265  }
266  }
267 
268  // Print absolute residual norm logging.
269  const std::vector<MT> residualLog = debugTest.getLogResNorm();
270  if (numrhs==1 && proc_verbose && residualLog.size())
271  {
272  std::cout << "Absolute residual 2-norm [ " << residualLog.size() << " ] : ";
273  for (unsigned int i=0; i<residualLog.size(); i++)
274  std::cout << residualLog[i] << " ";
275  std::cout << std::endl;
276  std::cout << "Final abs 2-norm / rhs 2-norm : " << residualLog[residualLog.size()-1] / norm_denom[0] << std::endl;
277  }
278 
279  // Clean up.
280  delete [] dvals;
281  delete [] colptr;
282  delete [] rowind;
283  delete [] cvals;
284 
285  success = ret==Belos::Converged && !norm_failure;
286  if (success) {
287  if (proc_verbose)
288  std::cout << "End Result: TEST PASSED" << std::endl;
289  } else {
290  if (proc_verbose)
291  std::cout << "End Result: TEST FAILED" << std::endl;
292  }
293  }
294  TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success);
295 
296  return ( success ? EXIT_SUCCESS : EXIT_FAILURE );
297 } // 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.