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"
51 #include "Teuchos_CommandLineProcessor.hpp"
52 #include "Teuchos_ParameterList.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  typedef std::complex<double> ST;
72  typedef ScalarTraits<ST> SCT;
73  typedef SCT::magnitudeType MT;
74  typedef Belos::MultiVec<ST> MV;
75  typedef Belos::Operator<ST> OP;
76  typedef Belos::MultiVecTraits<ST,MV> MVT;
78  ST one = SCT::one();
79  ST zero = SCT::zero();
80 
81  int info = 0;
82  int MyPID = 0;
83  bool norm_failure = false;
84 
85  Teuchos::GlobalMPISession session(&argc, &argv, NULL);
86  //
87  using Teuchos::RCP;
88  using Teuchos::rcp;
89 
90 bool success = false;
91 bool verbose = false;
92 try {
93 bool proc_verbose = false;
94  int frequency = -1; // how often residuals are printed by solver
95  int blocksize = 1;
96  int numrhs = 1;
97  std::string filename("mhd1280b.cua");
98  MT tol = 1.0e-5; // relative residual tolerance
99 
100  CommandLineProcessor cmdp(false,true);
101  cmdp.setOption("verbose","quiet",&verbose,"Print messages and results.");
102  cmdp.setOption("frequency",&frequency,"Solvers frequency for printing residuals (#iters).");
103  cmdp.setOption("filename",&filename,"Filename for Harwell-Boeing test matrix.");
104  cmdp.setOption("tol",&tol,"Relative residual tolerance used by MINRES solver.");
105  cmdp.setOption("num-rhs",&numrhs,"Number of right-hand sides to be solved for.");
106  if (cmdp.parse(argc,argv) != CommandLineProcessor::PARSE_SUCCESSFUL) {
107  return -1;
108  }
109 
110  proc_verbose = verbose && (MyPID==0); /* Only print on the zero processor */
111  if (proc_verbose) {
112  std::cout << Belos::Belos_Version() << std::endl << std::endl;
113  }
114  if (!verbose)
115  frequency = -1; // reset frequency if test is not verbose
116 
117 
118 #ifndef HAVE_BELOS_TRIUTILS
119  std::cout << "This test requires Triutils. Please configure with --enable-triutils." << std::endl;
120  if (MyPID==0) {
121  std::cout << "End Result: TEST FAILED" << std::endl;
122  }
123  return -1;
124 #endif
125 
126  // Get the data from the HB file
127  int dim,dim2,nnz;
128  MT *dvals;
129  int *colptr,*rowind;
130  ST *cvals;
131  nnz = -1;
132  info = readHB_newmat_double(filename.c_str(),&dim,&dim2,&nnz,
133  &colptr,&rowind,&dvals);
134  if (info == 0 || nnz < 0) {
135  if (MyPID==0) {
136  std::cout << "Error reading '" << filename << "'" << std::endl;
137  std::cout << "End Result: TEST FAILED" << std::endl;
138  }
139  return -1;
140  }
141  // Convert interleaved doubles to std::complex values
142  cvals = new ST[nnz];
143  for (int ii=0; ii<nnz; ii++) {
144  cvals[ii] = ST(dvals[ii*2],dvals[ii*2+1]);
145  }
146  // Build the problem matrix
147  RCP< MyBetterOperator<ST> > A
148  = rcp( new MyBetterOperator<ST>(dim,colptr,nnz,rowind,cvals) );
149  //
150  // ********Other information used by block solver***********
151  // *****************(can be user specified)******************
152  //
153  int maxits = dim/blocksize; // maximum number of iterations to run
154  //
155  ParameterList belosList;
156  belosList.set( "Maximum Iterations", maxits ); // Maximum number of iterations allowed
157  belosList.set( "Convergence Tolerance", tol ); // Relative convergence tolerance requested
158  if (verbose) {
159  belosList.set( "Verbosity", Belos::Errors + Belos::Warnings +
161  if (frequency > 0)
162  belosList.set( "Output Frequency", frequency );
163  }
164  else
165  belosList.set( "Verbosity", Belos::Errors + Belos::Warnings );
166  //
167  // Construct the right-hand side and solution multivectors.
168  // NOTE: The right-hand side will be constructed such that the solution is
169  // a vectors of one.
170  //
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  //
177  // Construct an unpreconditioned linear problem instance.
178  //
179  RCP<Belos::LinearProblem<ST,MV,OP> > problem =
180  rcp( new Belos::LinearProblem<ST,MV,OP>( A, soln, rhs ) );
181  bool set = problem->setProblem();
182  if (set == false) {
183  if (proc_verbose)
184  std::cout << std::endl << "ERROR: Belos::LinearProblem failed to set up correctly!" << std::endl;
185  return -1;
186  }
187  //
188  // *******************************************************************
189  // *************Start the MINRES iteration***********************
190  // *******************************************************************
191  //
192  Belos::MinresSolMgr<ST,MV,OP> solver( problem, rcp(&belosList,false) );
193 
194  //
195  // **********Print out information about problem*******************
196  //
197  if (proc_verbose) {
198  std::cout << std::endl << std::endl;
199  std::cout << "Dimension of matrix: " << dim << std::endl;
200  std::cout << "Number of right-hand sides: " << numrhs << std::endl;
201  std::cout << "Block size used by solver: " << blocksize << std::endl;
202  std::cout << "Max number of MINRES iterations: " << maxits << std::endl;
203  std::cout << "Relative residual tolerance: " << tol << std::endl;
204  std::cout << std::endl;
205  }
206  //
207  // Perform solve
208  //
209  Belos::ReturnType ret = solver.solve();
210  //
211  // Compute actual residuals.
212  //
213  RCP<MyMultiVec<ST> > temp = rcp( new MyMultiVec<ST>(dim,numrhs) );
214  OPT::Apply( *A, *soln, *temp );
215  MVT::MvAddMv( one, *rhs, -one, *temp, *temp );
216  std::vector<MT> norm_num(numrhs), norm_denom(numrhs);
217  MVT::MvNorm( *temp, norm_num );
218  MVT::MvNorm( *rhs, norm_denom );
219  for (int i=0; i<numrhs; ++i) {
220  if (proc_verbose)
221  std::cout << "Relative residual "<<i<<" : " << norm_num[i] / norm_denom[i] << std::endl;
222  if ( norm_num[i] / norm_denom[i] > tol ) {
223  norm_failure = true;
224  }
225  }
226 
227  // Clean up.
228  delete [] dvals;
229  delete [] colptr;
230  delete [] rowind;
231  delete [] cvals;
232 
233 success = ret==Belos::Converged && !norm_failure;
234 
235 if (success) {
236  if (proc_verbose)
237  std::cout << "End Result: TEST PASSED" << std::endl;
238 } else {
239 if (proc_verbose)
240  std::cout << "End Result: TEST FAILED" << std::endl;
241 }
242 }
243 TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success);
244 
245 return ( success ? EXIT_SUCCESS : EXIT_FAILURE );
246 } // end test_minres_complex_hb.cpp
Solver manager for the MINRES linear solver.
std::string Belos_Version()
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:65
Alternative run-time polymorphic interface for operators.
MINRES linear solver solution manager.
A linear system to solve, and its associated information.
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: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.