Belos Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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()
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
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.
#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
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...
ReturnType solve() override
Iterate until the status test tells us to stop.
Simple example of a user&#39;s defined Belos::Operator class.