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"
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  bool norm_failure = false;
83 
84  Teuchos::GlobalMPISession session(&argc, &argv, NULL);
85  int MyPID = session.getRank();
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 GCRODR 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 EXIT_FAILURE;
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 #ifndef HAVE_BELOS_TRIUTILS
118  std::cout << "This test requires Triutils. Please configure with --enable-triutils." << std::endl;
119  if (MyPID==0) {
120  std::cout << "End Result: TEST FAILED" << std::endl;
121  }
122  return EXIT_FAILURE;
123 #endif
124  // Get the data from the HB file
125  int dim,dim2,nnz;
126  MT *dvals;
127  int *colptr,*rowind;
128  ST *cvals;
129  nnz = -1;
130  info = readHB_newmat_double(filename.c_str(),&dim,&dim2,&nnz,
131  &colptr,&rowind,&dvals);
132  if (info == 0 || nnz < 0) {
133  if (MyPID==0) {
134  std::cout << "Error reading '" << filename << "'" << std::endl;
135  std::cout << "End Result: TEST FAILED" << std::endl;
136  }
137  return EXIT_FAILURE;
138  }
139  // Convert interleaved doubles to std::complex values
140  cvals = new ST[nnz];
141  for (int ii=0; ii<nnz; ii++) {
142  cvals[ii] = ST(dvals[ii*2],dvals[ii*2+1]);
143  }
144  // Build the problem matrix
145  RCP< MyBetterOperator<ST> > A
146  = rcp( new MyBetterOperator<ST>(dim,colptr,nnz,rowind,cvals) );
147  // ********Other information used by block solver***********
148  // *****************(can be user specified)******************
149  int maxits = dim/blocksize; // maximum number of iterations to run
150  int numBlocks = 100;
151  int numRecycledBlocks = 20;
152  int numIters1, numIters2, numIters3;
153  ParameterList belosList;
154  belosList.set( "Maximum Iterations", maxits ); // Maximum number of iterations allowed
155  belosList.set( "Convergence Tolerance", tol ); // Relative convergence tolerance requested
157  belosList.set( "Num Blocks", numBlocks );
158  belosList.set( "Num Recycled Blocks", numRecycledBlocks );
159  // Construct the right-hand side and solution multivectors.
160  // NOTE: The right-hand side will be constructed such that the solution is
161  // a vectors of one.
162  RCP<MyMultiVec<ST> > soln = rcp( new MyMultiVec<ST>(dim,numrhs) );
163  RCP<MyMultiVec<ST> > rhs = rcp( new MyMultiVec<ST>(dim,numrhs) );
164  MVT::MvRandom( *soln );
165  OPT::Apply( *A, *soln, *rhs );
166  MVT::MvInit( *soln, zero );
167  // Construct an unpreconditioned linear problem instance.
168  RCP<Belos::LinearProblem<ST,MV,OP> > problem =
169  rcp( new Belos::LinearProblem<ST,MV,OP>( A, soln, rhs ) );
170  bool set = problem->setProblem();
171  if (set == false) {
172  if (proc_verbose)
173  std::cout << std::endl << "ERROR: Belos::LinearProblem failed to set up correctly!" << std::endl;
174  return EXIT_FAILURE;
175  }
176  // *******************************************************************
177  // *************Start the GCRODR iteration***********************
178  // *******************************************************************
179  Belos::GCRODRSolMgr<ST,MV,OP> solver( problem, rcp(&belosList,false) );
180  // **********Print out information about problem*******************
181  if (proc_verbose) {
182  std::cout << std::endl << std::endl;
183  std::cout << "Dimension of matrix: " << dim << std::endl;
184  std::cout << "Number of right-hand sides: " << numrhs << std::endl;
185  std::cout << "Block size used by solver: " << blocksize << std::endl;
186  std::cout << "Max number of GCRODR iterations: " << maxits << std::endl;
187  std::cout << "Relative residual tolerance: " << tol << std::endl;
188  std::cout << std::endl;
189  }
190  // Perform solve
191  Belos::ReturnType ret = solver.solve();
192  // Get number of iterations
193  numIters1=solver.getNumIters();
194  // Compute actual residuals.
195  RCP<MyMultiVec<ST> > temp = rcp( new MyMultiVec<ST>(dim,numrhs) );
196  std::vector<MT> norm_num(numrhs), norm_denom(numrhs);
197  OPT::Apply( *A, *soln, *temp );
198  MVT::MvAddMv( one, *rhs, -one, *temp, *temp );
199  MVT::MvNorm( *temp, norm_num );
200  MVT::MvNorm( *rhs, norm_denom );
201  for (int i=0; i<numrhs; ++i) {
202  if (proc_verbose)
203  std::cout << "Relative residual "<<i<<" : " << norm_num[i] / norm_denom[i] << std::endl;
204  if ( norm_num[i] / norm_denom[i] > tol ) {
205  norm_failure = true;
206  }
207  }
208  // Resolve linear system with same rhs and recycled space
209  MVT::MvInit( *soln, zero );
210  solver.reset(Belos::Problem);
211  ret = solver.solve();
212  numIters2=solver.getNumIters();
213 
214  // Resolve linear system (again) with same rhs and recycled space
215  MVT::MvInit( *soln, zero );
216  solver.reset(Belos::Problem);
217  ret = solver.solve();
218  numIters3=solver.getNumIters();
219  // Clean up.
220  delete [] dvals;
221  delete [] colptr;
222  delete [] rowind;
223  delete [] cvals;
224  // Test for failures
225  if ( ret!=Belos::Converged || norm_failure || numIters1 < numIters2 || numIters2 < numIters3 ) {
226  success = false;
227  if (proc_verbose)
228  std::cout << "End Result: TEST FAILED" << std::endl;
229  } else {
230  success = true;
231  if (proc_verbose)
232  std::cout << "End Result: TEST PASSED" << std::endl;
233  }
234  }
235  TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success);
236 
237  return ( success ? EXIT_SUCCESS : EXIT_FAILURE );
238 } // end test_gcrodr_complex_hb.cpp
int main(int argc, char *argv[])
std::string Belos_Version()
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...
Alternative run-time polymorphic interface for operators.
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: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.