Belos Package Browser (Single Doxygen Collection)  Development
test_gcrodr_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 "BelosGCRODRSolMgr.hpp"
53 
54 #ifdef HAVE_MPI
55 #include <mpi.h>
56 #endif
57 
58 // I/O for Harwell-Boeing files
59 #ifdef HAVE_BELOS_TRIUTILS
60 #include "Trilinos_Util_iohb.h"
61 #endif
62 
63 #include "MyMultiVec.hpp"
64 #include "MyBetterOperator.hpp"
65 #include "MyOperator.hpp"
66 
67 using namespace Teuchos;
68 
69 int main(int argc, char *argv[]) {
70  //
71 #ifdef HAVE_COMPLEX
72  typedef std::complex<double> ST;
73 #elif HAVE_COMPLEX_H
74  typedef std::complex<double> ST;
75 #else
76  std::cout << "Not compiled with std::complex support." << std::endl;
77  std::cout << "End Result: TEST FAILED" << std::endl;
78  return EXIT_FAILURE;
79 #endif
80 
81  typedef ScalarTraits<ST> SCT;
82  typedef SCT::magnitudeType MT;
83  typedef Belos::MultiVec<ST> MV;
84  typedef Belos::Operator<ST> OP;
85  typedef Belos::MultiVecTraits<ST,MV> MVT;
87  ST one = SCT::one();
88  ST zero = SCT::zero();
89 
90  int info = 0;
91  bool norm_failure = false;
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  bool proc_verbose = false;
103  int frequency = -1; // how often residuals are printed by solver
104  int blocksize = 1;
105  int numrhs = 1;
106  std::string filename("mhd1280b.cua");
107  MT tol = 1.0e-5; // relative residual tolerance
108 
109  CommandLineProcessor cmdp(false,true);
110  cmdp.setOption("verbose","quiet",&verbose,"Print messages and results.");
111  cmdp.setOption("frequency",&frequency,"Solvers frequency for printing residuals (#iters).");
112  cmdp.setOption("filename",&filename,"Filename for Harwell-Boeing test matrix.");
113  cmdp.setOption("tol",&tol,"Relative residual tolerance used by GCRODR solver.");
114  cmdp.setOption("num-rhs",&numrhs,"Number of right-hand sides to be solved for.");
115  if (cmdp.parse(argc,argv) != CommandLineProcessor::PARSE_SUCCESSFUL) {
116  return EXIT_FAILURE;
117  }
118 
119  proc_verbose = verbose && (MyPID==0); /* Only print on the zero processor */
120  if (proc_verbose) {
121  std::cout << Belos::Belos_Version() << std::endl << std::endl;
122  }
123  if (!verbose)
124  frequency = -1; // reset frequency if test is not verbose
125 
126 #ifndef HAVE_BELOS_TRIUTILS
127  std::cout << "This test requires Triutils. Please configure with --enable-triutils." << std::endl;
128  if (MyPID==0) {
129  std::cout << "End Result: TEST FAILED" << std::endl;
130  }
131  return EXIT_FAILURE;
132 #endif
133  // Get the data from the HB file
134  int dim,dim2,nnz;
135  MT *dvals;
136  int *colptr,*rowind;
137  ST *cvals;
138  nnz = -1;
139  info = readHB_newmat_double(filename.c_str(),&dim,&dim2,&nnz,
140  &colptr,&rowind,&dvals);
141  if (info == 0 || nnz < 0) {
142  if (MyPID==0) {
143  std::cout << "Error reading '" << filename << "'" << std::endl;
144  std::cout << "End Result: TEST FAILED" << std::endl;
145  }
146  return EXIT_FAILURE;
147  }
148  // Convert interleaved doubles to std::complex values
149  cvals = new ST[nnz];
150  for (int ii=0; ii<nnz; ii++) {
151  cvals[ii] = ST(dvals[ii*2],dvals[ii*2+1]);
152  }
153  // Build the problem matrix
155  = rcp( new MyBetterOperator<ST>(dim,colptr,nnz,rowind,cvals) );
156  // ********Other information used by block solver***********
157  // *****************(can be user specified)******************
158  int maxits = dim/blocksize; // maximum number of iterations to run
159  int numBlocks = 100;
160  int numRecycledBlocks = 80;
161  int numIters1, numIters2, numIters3;
162  ParameterList belosList;
163  belosList.set( "Maximum Iterations", maxits ); // Maximum number of iterations allowed
164  belosList.set( "Convergence Tolerance", tol ); // Relative convergence tolerance requested
165  belosList.set( "Verbosity", Belos::Errors + Belos::Warnings + Belos::StatusTestDetails);
166  belosList.set( "Num Blocks", numBlocks );
167  belosList.set( "Num Recycled Blocks", numRecycledBlocks );
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  RCP<MyMultiVec<ST> > soln = rcp( new MyMultiVec<ST>(dim,numrhs) );
172  RCP<MyMultiVec<ST> > rhs = rcp( new MyMultiVec<ST>(dim,numrhs) );
173  MVT::MvRandom( *soln );
174  OPT::Apply( *A, *soln, *rhs );
175  MVT::MvInit( *soln, zero );
176  // Construct an unpreconditioned linear problem instance.
178  rcp( new Belos::LinearProblem<ST,MV,OP>( A, soln, rhs ) );
179  bool set = problem->setProblem();
180  if (set == false) {
181  if (proc_verbose)
182  std::cout << std::endl << "ERROR: Belos::LinearProblem failed to set up correctly!" << std::endl;
183  return EXIT_FAILURE;
184  }
185  // *******************************************************************
186  // *************Start the GCRODR iteration***********************
187  // *******************************************************************
188  Belos::GCRODRSolMgr<ST,MV,OP> solver( problem, rcp(&belosList,false) );
189  // **********Print out information about problem*******************
190  if (proc_verbose) {
191  std::cout << std::endl << std::endl;
192  std::cout << "Dimension of matrix: " << dim << std::endl;
193  std::cout << "Number of right-hand sides: " << numrhs << std::endl;
194  std::cout << "Block size used by solver: " << blocksize << std::endl;
195  std::cout << "Max number of GCRODR iterations: " << maxits << std::endl;
196  std::cout << "Relative residual tolerance: " << tol << std::endl;
197  std::cout << std::endl;
198  }
199  // Perform solve
200  Belos::ReturnType ret = solver.solve();
201  // Get number of iterations
202  numIters1=solver.getNumIters();
203  // Compute actual residuals.
204  RCP<MyMultiVec<ST> > temp = rcp( new MyMultiVec<ST>(dim,numrhs) );
205  std::vector<MT> norm_num(numrhs), norm_denom(numrhs);
206  OPT::Apply( *A, *soln, *temp );
207  MVT::MvAddMv( one, *rhs, -one, *temp, *temp );
208  MVT::MvNorm( *temp, norm_num );
209  MVT::MvNorm( *rhs, norm_denom );
210  for (int i=0; i<numrhs; ++i) {
211  if (proc_verbose)
212  std::cout << "Relative residual "<<i<<" : " << norm_num[i] / norm_denom[i] << std::endl;
213  if ( norm_num[i] / norm_denom[i] > tol ) {
214  norm_failure = true;
215  }
216  }
217  // Setup new rhs, perform solve
218  MVT::MvRandom( *soln );
219  OPT::Apply( *A, *soln, *rhs );
220  MVT::MvInit( *soln, zero );
221  solver.reset(Belos::Problem);
222  ret = solver.solve();
223  numIters2=solver.getNumIters();
224  // Setup new rhs, perform solve
225  MVT::MvRandom( *soln );
226  OPT::Apply( *A, *soln, *rhs );
227  MVT::MvInit( *soln, zero );
228  solver.reset(Belos::Problem);
229  ret = solver.solve();
230  numIters3=solver.getNumIters();
231  // Clean up.
232  delete [] dvals;
233  delete [] colptr;
234  delete [] rowind;
235  delete [] cvals;
236  // Test for failures
237  if ( ret!=Belos::Converged || norm_failure || numIters1 < numIters2 || numIters1 < numIters3 ) {
238  success = false;
239  if (proc_verbose)
240  std::cout << "End Result: TEST FAILED" << std::endl;
241  } else {
242  success = true;
243  if (proc_verbose)
244  std::cout << "End Result: TEST PASSED" << std::endl;
245  }
246  }
247  TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success);
248 
249  return ( success ? EXIT_SUCCESS : EXIT_FAILURE );
250 } // end test_gcrodr_complex_hb.cpp
int main(int argc, char *argv[])
std::string Belos_Version()
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Traits class which defines basic operations on multivectors.
Simple example of a user&#39;s defined Belos::MultiVec class.
Definition: MyMultiVec.hpp:62
Declaration and definition of Belos::GCRODRSolMgr, which implements the GCRODR (recycling GMRES) solv...
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)
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.
Implementation of the GCRODR (Recycling GMRES) iterative linear solver.
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.
Belos header file which uses auto-configuration information to include necessary C++ headers...
Simple example of a user&#39;s defined Belos::Operator class.