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