Belos Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
test_hybrid_gmres_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"
52#include "Teuchos_CommandLineProcessor.hpp"
53#include "Teuchos_ParameterList.hpp"
54#include "Teuchos_StandardCatchMacros.hpp"
55#include "Teuchos_StandardCatchMacros.hpp"
56
57#ifdef HAVE_MPI
58#include <mpi.h>
59#endif
60
61// I/O for Harwell-Boeing files
62#ifdef HAVE_BELOS_TRIUTILS
63#include "Trilinos_Util_iohb.h"
64#endif
65
66#include "MyMultiVec.hpp"
67#include "MyBetterOperator.hpp"
68#include "MyOperator.hpp"
69
70using namespace Teuchos;
71
72int main(int argc, char *argv[]) {
73 //
74 typedef std::complex<double> ST;
75 typedef ScalarTraits<ST> SCT;
76 typedef SCT::magnitudeType MT;
77 typedef Belos::MultiVec<ST> MV;
78 typedef Belos::Operator<ST> OP;
81 ST one = SCT::one();
82 ST zero = SCT::zero();
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 int info = 0;
94 bool norm_failure = false;
95 bool proc_verbose = false;
96 bool userandomrhs = true;
97 int frequency = -1; // frequency of status test output.
98 int blocksize = 1; // blocksize
99 int numrhs = 1; // number of right-hand sides to solve for
100 int maxiters = -1; // maximum number of iterations allowed per linear system
101 int maxdegree = 25; // maximum degree of polynomial
102 int maxsubspace = 50; // maximum number of blocks the solver can use for the subspace
103 int maxrestarts = 15; // number of restarts allowed
104 std::string outersolver("Block Gmres");
105 std::string filename("mhd1280b.cua");
106 std::string polyType("Arnoldi");
107 MT tol = 1.0e-5; // relative residual tolerance
108 MT polytol = tol/10; // relative residual tolerance for polynomial construction
109
110 Teuchos::CommandLineProcessor cmdp(false,true);
111 cmdp.setOption("verbose","quiet",&verbose,"Print messages and results.");
112 cmdp.setOption("use-random-rhs","use-rhs",&userandomrhs,"Use linear system RHS or random RHS to generate polynomial.");
113 cmdp.setOption("frequency",&frequency,"Solvers frequency for printing residuals (#iters).");
114 cmdp.setOption("filename",&filename,"Filename for test matrix. Acceptable file extensions: *.hb,*.mtx,*.triU,*.triS");
115 cmdp.setOption("outersolver",&outersolver,"Name of outer solver to be used with GMRES poly");
116 cmdp.setOption("tol",&tol,"Relative residual tolerance used by GMRES solver.");
117 cmdp.setOption("poly-tol",&polytol,"Relative residual tolerance used to construct the GMRES polynomial.");
118 cmdp.setOption("poly-type",&polyType,"Polynomial type (Roots, Arnoldi, or Gmres).");
119 cmdp.setOption("num-rhs",&numrhs,"Number of right-hand sides to be solved for.");
120 cmdp.setOption("block-size",&blocksize,"Block size used by GMRES.");
121 cmdp.setOption("max-iters",&maxiters,"Maximum number of iterations per linear system (-1 = adapted to problem/block size).");
122 cmdp.setOption("max-degree",&maxdegree,"Maximum degree of the GMRES polynomial.");
123 cmdp.setOption("max-subspace",&maxsubspace,"Maximum number of blocks the solver can use for the subspace.");
124 cmdp.setOption("max-restarts",&maxrestarts,"Maximum number of restarts allowed for GMRES solver.");
125 if (cmdp.parse(argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
126 return EXIT_FAILURE;
127 }
128
129 proc_verbose = verbose && (MyPID==0); /* Only print on the zero processor */
130 if (proc_verbose) {
131 std::cout << Belos::Belos_Version() << std::endl << std::endl;
132 }
133 if (!verbose)
134 frequency = -1; // reset frequency if test is not verbose
135
136#ifndef HAVE_BELOS_TRIUTILS
137 std::cout << "This test requires Triutils. Please configure with --enable-triutils." << std::endl;
138 if (MyPID==0) {
139 std::cout << "End Result: TEST FAILED" << std::endl;
140 }
141 return EXIT_FAILURE;
142#endif
143
144 // Get the data from the HB file
145 int dim,dim2,nnz;
146 MT *dvals;
147 int *colptr,*rowind;
148 ST *cvals;
149 nnz = -1;
150 info = readHB_newmat_double(filename.c_str(),&dim,&dim2,&nnz,
151 &colptr,&rowind,&dvals);
152 if (info == 0 || nnz < 0) {
153 if (MyPID==0) {
154 std::cout << "Error reading '" << filename << "'" << std::endl;
155 std::cout << "End Result: TEST FAILED" << std::endl;
156 }
157 return EXIT_FAILURE;
158 }
159 // Convert interleaved doubles to std::complex values
160 cvals = new ST[nnz];
161 for (int ii=0; ii<nnz; ii++) {
162 cvals[ii] = ST(dvals[ii*2],dvals[ii*2+1]);
163 }
164 // Build the problem matrix
165 RCP< MyBetterOperator<ST> > A
166 = rcp( new MyBetterOperator<ST>(dim,colptr,nnz,rowind,cvals) );
167 //
168 // Construct the right-hand side and solution multivectors.
169 // NOTE: The right-hand side will be constructed such that the solution is
170 // a vectors of one.
171 //
172 RCP<MyMultiVec<ST> > soln = rcp( new MyMultiVec<ST>(dim,numrhs) );
173 RCP<MyMultiVec<ST> > rhs = rcp( new MyMultiVec<ST>(dim,numrhs) );
174 MVT::MvRandom( *soln );
175 OPT::Apply( *A, *soln, *rhs );
176 MVT::MvInit( *soln, zero );
177 //
178 // Construct an unpreconditioned linear problem instance.
179 //
180 RCP<Belos::LinearProblem<ST,MV,OP> > problem =
181 rcp( new Belos::LinearProblem<ST,MV,OP>( A, soln, rhs ) );
182 problem->setInitResVec( rhs );
183 bool set = problem->setProblem();
184 if (set == false) {
185 if (proc_verbose)
186 std::cout << std::endl << "ERROR: Belos::LinearProblem failed to set up correctly!" << std::endl;
187 return EXIT_FAILURE;
188 }
189 //
190 // ********Other information used by block solver***********
191 // *****************(can be user specified)******************
192 //
193 if (maxiters == -1)
194 maxiters = dim/blocksize - 1; // maximum number of iterations to run
195
196 ParameterList belosList;
197 belosList.set( "Num Blocks", maxsubspace); // Maximum number of blocks in Krylov factorization
198 belosList.set( "Block Size", blocksize ); // Blocksize to be used by iterative solver
199 belosList.set( "Maximum Iterations", maxiters ); // Maximum number of iterations allowed
200 belosList.set( "Maximum Restarts", maxrestarts ); // Maximum number of restarts allowed
201 belosList.set( "Convergence Tolerance", tol ); // Relative convergence tolerance requested
202 int verbosity = Belos::Errors + Belos::Warnings;
203 if (verbose) {
205 if (frequency > 0)
206 belosList.set( "Output Frequency", frequency );
207 }
208 belosList.set( "Verbosity", verbosity );
209
210 ParameterList polyList;
211 polyList.set( "Maximum Degree", maxdegree ); // Maximum degree of the GMRES polynomial
212 polyList.set( "Polynomial Tolerance", polytol ); // Polynomial convergence tolerance requested
213 polyList.set( "Polynomial Type", polyType ); // Type of polynomial to construct
214 polyList.set( "Verbosity", verbosity ); // Verbosity for polynomial construction
215 polyList.set( "Random RHS", userandomrhs ); // Use RHS from linear system or random vector
216 if ( outersolver != "" ) {
217 polyList.set( "Outer Solver", outersolver );
218 polyList.set( "Outer Solver Params", belosList );
219 }
220
221 // Use a debugging status test to save absolute residual history.
222 // Debugging status tests are peer to the native status tests that are called whenever convergence is checked.
224
225 //
226 // *******************************************************************
227 // *************Start the block Gmres iteration***********************
228 // *******************************************************************
229 //
230 RCP< Belos::SolverManager<ST,MV,OP> > solver = rcp( new Belos::GmresPolySolMgr<ST,MV,OP>( problem, rcp(&polyList,false) ) );
231
232 // The debug status test does not work for the GmresPolySolMgr right now.
233 // solver->setDebugStatusTest( Teuchos::rcp(&debugTest, false) );
234
235 //
236 // **********Print out information about problem*******************
237 //
238 if (proc_verbose) {
239 std::cout << std::endl << std::endl;
240 std::cout << "Dimension of matrix: " << dim << std::endl;
241 std::cout << "Number of right-hand sides: " << numrhs << std::endl;
242 std::cout << "Block size used by solver: " << blocksize << std::endl;
243 std::cout << "Max number of Gmres iterations: " << maxiters << std::endl;
244 std::cout << "Relative residual tolerance: " << tol << std::endl;
245 std::cout << std::endl;
246 }
247 //
248 // Perform solve
249 //
250 Belos::ReturnType ret = solver->solve();
251 //
252 std::cout << "Belos::GmresPolySolMgr returned achievedTol: " << solver->achievedTol() << std::endl << std::endl;
253 //
254 // Compute actual residuals.
255 RCP<MyMultiVec<ST> > temp = rcp( new MyMultiVec<ST>(dim,numrhs) );
256 OPT::Apply( *A, *soln, *temp );
257 MVT::MvAddMv( one, *rhs, -one, *temp, *temp );
258 std::vector<MT> norm_num(numrhs), norm_denom(numrhs);
259 MVT::MvNorm( *temp, norm_num );
260 MVT::MvNorm( *rhs, norm_denom );
261 for (int i=0; i<numrhs; ++i) {
262 if (proc_verbose)
263 std::cout << "Relative residual "<<i<<" : " << norm_num[i] / norm_denom[i] << std::endl;
264 if ( norm_num[i] / norm_denom[i] > tol ) {
265 norm_failure = true;
266 }
267 }
268
269 // Print absolute residual norm logging.
270 const std::vector<MT> residualLog = debugTest.getLogResNorm();
271 if (numrhs==1 && proc_verbose && residualLog.size())
272 {
273 std::cout << "Absolute residual 2-norm [ " << residualLog.size() << " ] : ";
274 for (unsigned int i=0; i<residualLog.size(); i++)
275 std::cout << residualLog[i] << " ";
276 std::cout << std::endl;
277 std::cout << "Final abs 2-norm / rhs 2-norm : " << residualLog[residualLog.size()-1] / norm_denom[0] << std::endl;
278 }
279
280 // Clean up.
281 delete [] dvals;
282 delete [] colptr;
283 delete [] rowind;
284 delete [] cvals;
285
286 success = ret==Belos::Converged && !norm_failure;
287 if (success) {
288 if (proc_verbose)
289 std::cout << "End Result: TEST PASSED" << std::endl;
290 } else {
291 if (proc_verbose)
292 std::cout << "End Result: TEST FAILED" << std::endl;
293 }
294 }
295 TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success);
296
297 return ( success ? EXIT_SUCCESS : EXIT_FAILURE );
298} // end test_bl_gmres_complex_hb.cpp
Belos header file which uses auto-configuration information to include necessary C++ headers.
Declaration and definition of Belos::GmresPolySolMgr (hybrid block GMRES linear solver).
Class which describes the linear problem to be solved by the iterative solver.
Belos::StatusTest debugging class for storing the absolute residual norms generated during a solve.
The GMRES polynomial can be created in conjunction with any standard preconditioner.
A linear system to solve, and its associated information.
Traits class which defines basic operations on multivectors.
Interface for multivectors used by Belos' linear solvers.
Class which defines basic traits for the operator type.
Alternative run-time polymorphic interface for operators.
A Belos::StatusTest debugging class for storing the absolute residual norms generated during a solve.
const std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > & getLogResNorm() const
Returns the log of the absolute residual norm from the iteration.
Simple example of a user's defined Belos::Operator class.
Simple example of a user's defined Belos::MultiVec class.
Definition: MyMultiVec.hpp:66
@ StatusTestDetails
Definition: BelosTypes.hpp:261
@ FinalSummary
Definition: BelosTypes.hpp:259
@ TimingDetails
Definition: BelosTypes.hpp:260
@ Warnings
Definition: BelosTypes.hpp:256
ReturnType
Whether the Belos solve converged for all linear systems.
Definition: BelosTypes.hpp:155
@ Converged
Definition: BelosTypes.hpp:156
std::string Belos_Version()
int main(int argc, char *argv[])