Belos Package Browser (Single Doxygen Collection)  Development
test_bl_gmres_complex_diag.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"
53 #include "Teuchos_CommandLineProcessor.hpp"
54 #include "Teuchos_ParameterList.hpp"
55 #include "Teuchos_StandardCatchMacros.hpp"
56 
57 #ifdef HAVE_MPI
58 #include <mpi.h>
59 #endif
60 
61 #include "MyMultiVec.hpp"
62 #include "MyOperator.hpp"
63 
64 using namespace Teuchos;
65 
66 int main(int argc, char *argv[]) {
67  //
68 #ifdef HAVE_COMPLEX
69  typedef std::complex<double> ST;
70 #elif HAVE_COMPLEX_H
71  typedef std::complex<double> ST;
72 #else
73  std::cout << "Not compiled with std::complex support." << std::endl;
74  std::cout << "End Result: TEST FAILED" << std::endl;
75  return EXIT_FAILURE;
76 #endif
77 
78  typedef ScalarTraits<ST> SCT;
79  typedef SCT::magnitudeType MT;
80  typedef Belos::MultiVec<ST> MV;
81  typedef Belos::Operator<ST> OP;
82  typedef Belos::MultiVecTraits<ST,MV> MVT;
84  ST one = SCT::one();
85  ST zero = SCT::zero();
86 
87  Teuchos::GlobalMPISession session(&argc, &argv, NULL);
88  int MyPID = session.getRank();
89  //
90  using Teuchos::RCP;
91  using Teuchos::rcp;
92 
93  bool success = false;
94  bool verbose = false;
95  try {
96  bool norm_failure = false;
97  bool proc_verbose = false;
98  bool pseudo = false; // use pseudo block GMRES to solve this linear system.
99  int frequency = -1; // how often residuals are printed by solver
100  int blocksize = 1;
101  int numrhs = 1;
102  int maxrestarts = 15;
103  int length = 50;
104  std::string ortho("DGKS"); // The Belos documentation obscures the fact that
105  MT tol = 1.0e-5; // relative residual tolerance
106 
107  CommandLineProcessor cmdp(false,true);
108  cmdp.setOption("verbose","quiet",&verbose,"Print messages and results.");
109  cmdp.setOption("pseudo","regular",&pseudo,"Use pseudo-block GMRES to solve the linear systems.");
110  cmdp.setOption("frequency",&frequency,"Solvers frequency for printing residuals (#iters).");
111  cmdp.setOption("tol",&tol,"Relative residual tolerance used by GMRES solver.");
112  cmdp.setOption("num-rhs",&numrhs,"Number of right-hand sides to be solved for.");
113  cmdp.setOption("num-restarts",&maxrestarts,"Maximum number of restarts allowed for the GMRES solver.");
114  cmdp.setOption("blocksize",&blocksize,"Block size used by GMRES.");
115  cmdp.setOption("subspace-length",&length,"Maximum dimension of block-subspace used by GMRES solver.");
116  cmdp.setOption("ortho-type",&ortho,"Orthogonalization type, either DGKS, ICGS or IMGS");
117  if (cmdp.parse(argc,argv) != CommandLineProcessor::PARSE_SUCCESSFUL) {
118  return EXIT_FAILURE;
119  }
120 
121  proc_verbose = verbose && (MyPID==0); /* Only print on the zero processor */
122  if (proc_verbose) {
123  std::cout << Belos::Belos_Version() << std::endl << std::endl;
124  }
125  if (!verbose)
126  frequency = -1; // reset frequency if test is not verbose
127 
128 #ifndef HAVE_BELOS_TRIUTILS
129  std::cout << "This test requires Triutils. Please configure with --enable-triutils." << std::endl;
130  if (MyPID==0) {
131  std::cout << "End Result: TEST FAILED" << std::endl;
132  }
133  return EXIT_FAILURE;
134 #endif
135 
136  // Get the data from the HB file
137  int dim=100;
138 
139  // Build the problem matrix
140  std::vector<ST> diag( dim, (ST)4.0 );
141  RCP< MyOperator<ST> > A
142  = rcp( new MyOperator<ST>( diag ) );
143  //
144  // ********Other information used by block solver***********
145  // *****************(can be user specified)******************
146  //
147  int maxits = dim/blocksize; // maximum number of iterations to run
148  //
149  ParameterList belosList;
150  belosList.set( "Num Blocks", length ); // Maximum number of blocks in Krylov factorization
151  belosList.set( "Block Size", blocksize ); // Blocksize to be used by iterative solver
152  belosList.set( "Maximum Iterations", maxits ); // Maximum number of iterations allowed
153  belosList.set( "Maximum Restarts", maxrestarts ); // Maximum number of restarts allowed
154  belosList.set( "Convergence Tolerance", tol ); // Relative convergence tolerance requested
155  belosList.set( "Orthogonalization", ortho ); // Orthogonalization type
156  if (verbose) {
157  belosList.set( "Verbosity", Belos::Errors + Belos::Warnings +
159  if (frequency > 0)
160  belosList.set( "Output Frequency", frequency );
161  }
162  else
163  belosList.set( "Verbosity", Belos::Errors + Belos::Warnings );
164  //
165  // Construct the right-hand side and solution multivectors.
166  // NOTE: The right-hand side will be constructed such that the solution is
167  // a vectors of one.
168  //
169  RCP<MyMultiVec<ST> > soln = rcp( new MyMultiVec<ST>(dim,numrhs) );
170  RCP<MyMultiVec<ST> > rhs = rcp( new MyMultiVec<ST>(dim,numrhs) );
171  MVT::MvInit( *rhs, 1.0 );
172  MVT::MvInit( *soln, zero );
173  //
174  // Construct an unpreconditioned linear problem instance.
175  //
176  RCP<Belos::LinearProblem<ST,MV,OP> > problem =
177  rcp( new Belos::LinearProblem<ST,MV,OP>( A, soln, rhs ) );
178  bool set = problem->setProblem();
179  if (set == false) {
180  if (proc_verbose)
181  std::cout << std::endl << "ERROR: Belos::LinearProblem failed to set up correctly!" << std::endl;
182  return EXIT_FAILURE;
183  }
184 
185  // Use a debugging status test to save absolute residual history.
186  // Debugging status tests are peer to the native status tests that are called whenever convergence is checked.
188 
189  //
190  // *******************************************************************
191  // *************Start the block Gmres iteration***********************
192  // *******************************************************************
193  //
194  Teuchos::RCP< Belos::SolverManager<ST,MV,OP> > solver;
195  if (pseudo)
196  solver = Teuchos::rcp( new Belos::PseudoBlockGmresSolMgr<ST,MV,OP>( problem, Teuchos::rcp(&belosList,false) ) );
197  else
198  solver = Teuchos::rcp( new Belos::BlockGmresSolMgr<ST,MV,OP>( problem, Teuchos::rcp(&belosList,false) ) );
199 
200  solver->setDebugStatusTest( Teuchos::rcp(&debugTest, false) );
201 
202  //
203  // **********Print out information about problem*******************
204  //
205  if (proc_verbose) {
206  std::cout << std::endl << std::endl;
207  std::cout << "Dimension of matrix: " << dim << std::endl;
208  std::cout << "Number of right-hand sides: " << numrhs << std::endl;
209  std::cout << "Block size used by solver: " << blocksize << std::endl;
210  std::cout << "Max number of Gmres iterations: " << maxits << std::endl;
211  std::cout << "Relative residual tolerance: " << tol << std::endl;
212  std::cout << std::endl;
213  }
214  //
215  // Perform solve
216  //
217  Belos::ReturnType ret = solver->solve();
218  //
219  // Compute actual residuals.
220  //
221  RCP<MyMultiVec<ST> > temp = rcp( new MyMultiVec<ST>(dim,numrhs) );
222  OPT::Apply( *A, *soln, *temp );
223  MVT::MvAddMv( one, *rhs, -one, *temp, *temp );
224  std::vector<MT> norm_num(numrhs), norm_denom(numrhs);
225  MVT::MvNorm( *temp, norm_num );
226  MVT::MvNorm( *rhs, norm_denom );
227  for (int i=0; i<numrhs; ++i) {
228  if (proc_verbose)
229  std::cout << "Relative residual "<<i<<" : " << norm_num[i] / norm_denom[i] << std::endl;
230  if ( norm_num[i] / norm_denom[i] > tol ) {
231  norm_failure = true;
232  }
233  }
234 
235  // Print absolute residual norm logging.
236  const std::vector<MT> residualLog = debugTest.getLogResNorm();
237  if (numrhs==1 && proc_verbose && residualLog.size())
238  {
239  std::cout << "Absolute residual 2-norm [ " << residualLog.size() << " ] : ";
240  for (unsigned int i=0; i<residualLog.size(); i++)
241  std::cout << residualLog[i] << " ";
242  std::cout << std::endl;
243  std::cout << "Final abs 2-norm / rhs 2-norm : " << residualLog[residualLog.size()-1] / norm_denom[0] << std::endl;
244  }
245 
246  success = ret==Belos::Converged && !norm_failure;
247  if (success) {
248  if (proc_verbose)
249  std::cout << "End Result: TEST PASSED" << std::endl;
250  } else {
251  if (proc_verbose)
252  std::cout << "End Result: TEST FAILED" << std::endl;
253  }
254  }
255  TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success);
256 
257  return ( success ? EXIT_SUCCESS : EXIT_FAILURE );
258 } // end test_bl_gmres_complex_hb.cpp
std::string Belos_Version()
A Belos::StatusTest debugging class for storing the absolute residual norms generated during a solve...
Interface to Block GMRES and Flexible GMRES.
Belos::StatusTest debugging class for storing the absolute residual norms generated during a solve...
Traits class which defines basic operations on multivectors.
Simple example of a user&#39;s defined Belos::MultiVec class.
Definition: MyMultiVec.hpp:62
Alternative run-time polymorphic interface for operators.
The Belos::BlockGmresSolMgr provides a solver manager for the BlockGmres linear solver.
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
Interface to standard and "pseudoblock" GMRES.
ReturnType
Whether the Belos solve converged for all linear systems.
Definition: BelosTypes.hpp:154
Simple example of a user&#39;s defined Belos::Operator class.
Definition: MyOperator.hpp:65
Interface for multivectors used by Belos&#39; linear solvers.
Class which defines basic traits for the operator type.
The Belos::PseudoBlockGmresSolMgr provides a solver manager for the BlockGmres linear solver...
int main(int argc, char *argv[])
const std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > & getLogResNorm() const
Returns the log of the absolute residual norm from the iteration.
Belos header file which uses auto-configuration information to include necessary C++ headers...