Belos Package Browser (Single Doxygen Collection)  Development
test_minres_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 "BelosMinresSolMgr.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 -1;
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  int MyPID = 0;
92  bool norm_failure = false;
93 
94  Teuchos::GlobalMPISession session(&argc, &argv, NULL);
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 MINRES 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 -1;
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 
127 #ifndef HAVE_BELOS_TRIUTILS
128  std::cout << "This test requires Triutils. Please configure with --enable-triutils." << std::endl;
129  if (MyPID==0) {
130  std::cout << "End Result: TEST FAILED" << std::endl;
131  }
132  return -1;
133 #endif
134 
135  // Get the data from the HB file
136  int dim,dim2,nnz;
137  MT *dvals;
138  int *colptr,*rowind;
139  ST *cvals;
140  nnz = -1;
141  info = readHB_newmat_double(filename.c_str(),&dim,&dim2,&nnz,
142  &colptr,&rowind,&dvals);
143  if (info == 0 || nnz < 0) {
144  if (MyPID==0) {
145  std::cout << "Error reading '" << filename << "'" << std::endl;
146  std::cout << "End Result: TEST FAILED" << std::endl;
147  }
148  return -1;
149  }
150  // Convert interleaved doubles to std::complex values
151  cvals = new ST[nnz];
152  for (int ii=0; ii<nnz; ii++) {
153  cvals[ii] = ST(dvals[ii*2],dvals[ii*2+1]);
154  }
155  // Build the problem matrix
157  = rcp( new MyBetterOperator<ST>(dim,colptr,nnz,rowind,cvals) );
158  //
159  // ********Other information used by block solver***********
160  // *****************(can be user specified)******************
161  //
162  int maxits = dim/blocksize; // maximum number of iterations to run
163  //
164  ParameterList belosList;
165  belosList.set( "Maximum Iterations", maxits ); // Maximum number of iterations allowed
166  belosList.set( "Convergence Tolerance", tol ); // Relative convergence tolerance requested
167  if (verbose) {
168  belosList.set( "Verbosity", Belos::Errors + Belos::Warnings +
170  if (frequency > 0)
171  belosList.set( "Output Frequency", frequency );
172  }
173  else
174  belosList.set( "Verbosity", Belos::Errors + Belos::Warnings );
175  //
176  // Construct the right-hand side and solution multivectors.
177  // NOTE: The right-hand side will be constructed such that the solution is
178  // a vectors of one.
179  //
180  RCP<MyMultiVec<ST> > soln = rcp( new MyMultiVec<ST>(dim,numrhs) );
181  RCP<MyMultiVec<ST> > rhs = rcp( new MyMultiVec<ST>(dim,numrhs) );
182  MVT::MvRandom( *soln );
183  OPT::Apply( *A, *soln, *rhs );
184  MVT::MvInit( *soln, zero );
185  //
186  // Construct an unpreconditioned linear problem instance.
187  //
189  rcp( new Belos::LinearProblem<ST,MV,OP>( A, soln, 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 -1;
195  }
196  //
197  // *******************************************************************
198  // *************Start the MINRES iteration***********************
199  // *******************************************************************
200  //
201  Belos::MinresSolMgr<ST,MV,OP> solver( problem, rcp(&belosList,false) );
202 
203  //
204  // **********Print out information about problem*******************
205  //
206  if (proc_verbose) {
207  std::cout << std::endl << std::endl;
208  std::cout << "Dimension of matrix: " << dim << std::endl;
209  std::cout << "Number of right-hand sides: " << numrhs << std::endl;
210  std::cout << "Block size used by solver: " << blocksize << std::endl;
211  std::cout << "Max number of MINRES iterations: " << maxits << std::endl;
212  std::cout << "Relative residual tolerance: " << tol << std::endl;
213  std::cout << std::endl;
214  }
215  //
216  // Perform solve
217  //
218  Belos::ReturnType ret = solver.solve();
219  //
220  // Compute actual residuals.
221  //
222  RCP<MyMultiVec<ST> > temp = rcp( new MyMultiVec<ST>(dim,numrhs) );
223  OPT::Apply( *A, *soln, *temp );
224  MVT::MvAddMv( one, *rhs, -one, *temp, *temp );
225  std::vector<MT> norm_num(numrhs), norm_denom(numrhs);
226  MVT::MvNorm( *temp, norm_num );
227  MVT::MvNorm( *rhs, norm_denom );
228  for (int i=0; i<numrhs; ++i) {
229  if (proc_verbose)
230  std::cout << "Relative residual "<<i<<" : " << norm_num[i] / norm_denom[i] << std::endl;
231  if ( norm_num[i] / norm_denom[i] > tol ) {
232  norm_failure = true;
233  }
234  }
235 
236  // Clean up.
237  delete [] dvals;
238  delete [] colptr;
239  delete [] rowind;
240  delete [] cvals;
241 
242 success = ret==Belos::Converged && !norm_failure;
243 
244 if (success) {
245  if (proc_verbose)
246  std::cout << "End Result: TEST PASSED" << std::endl;
247 } else {
248 if (proc_verbose)
249  std::cout << "End Result: TEST FAILED" << std::endl;
250 }
251 }
252 TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success);
253 
254 return ( success ? EXIT_SUCCESS : EXIT_FAILURE );
255 } // end test_minres_complex_hb.cpp
Solver manager for the MINRES linear solver.
std::string Belos_Version()
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
int main(int argc, char *argv[])
Traits class which defines basic operations on multivectors.
Simple example of a user&#39;s defined Belos::MultiVec class.
Definition: MyMultiVec.hpp:62
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)
MINRES linear solver solution manager.
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.
const double tol
Class which describes the linear problem to be solved by the iterative 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...
ReturnType solve() override
Iterate until the status test tells us to stop.
Simple example of a user&#39;s defined Belos::Operator class.