Belos Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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
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)
Traits class which defines basic operations on multivectors.
Simple example of a user&#39;s defined Belos::MultiVec class.
Definition: MyMultiVec.hpp:65
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)
#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
Implementation of the GCRODR (Recycling 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.