Belos Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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"
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.");
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
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  //
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()
int main(int argc, char *argv[])
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
const std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > & getLogResNorm() const
Returns the log of the absolute residual norm from the iteration.
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
std::string filename
Alternative run-time polymorphic interface for operators.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void setOption(const char option_true[], const char option_false[], bool *option_val, const char documentation[]=NULL)
Declaration and definition of Belos::GmresPolySolMgr (hybrid block GMRES linear solver).
#define TEUCHOS_STANDARD_CATCH_STATEMENTS(VERBOSE, ERR_STREAM, SUCCESS_FLAG)
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.
EParseCommandLineReturn parse(int argc, char *argv[], std::ostream *errout=&std::cerr) const
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.
Belos header file which uses auto-configuration information to include necessary C++ headers...
Simple example of a user&#39;s defined Belos::Operator class.