Belos Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
test_bl_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"
53#include "Teuchos_CommandLineProcessor.hpp"
54#include "Teuchos_ParameterList.hpp"
55#include "Teuchos_StandardCatchMacros.hpp"
56#include "Teuchos_StandardCatchMacros.hpp"
57
58#ifdef HAVE_MPI
59#include <mpi.h>
60#endif
61
62// I/O for Harwell-Boeing files
63#ifdef HAVE_BELOS_TRIUTILS
64#include "Trilinos_Util_iohb.h"
65#endif
66
67#include "MyMultiVec.hpp"
68#include "MyBetterOperator.hpp"
69#include "MyOperator.hpp"
70
71using namespace Teuchos;
72
73int main(int argc, char *argv[]) {
74 //
75 typedef std::complex<double> ST;
76 typedef ScalarTraits<ST> SCT;
77 typedef SCT::magnitudeType MT;
78 typedef Belos::MultiVec<ST> MV;
79 typedef Belos::Operator<ST> OP;
82 ST one = SCT::one();
83 ST zero = SCT::zero();
84
85 Teuchos::GlobalMPISession session(&argc, &argv, NULL);
86 int MyPID = session.getRank();
87 //
88 using Teuchos::RCP;
89 using Teuchos::rcp;
90
91 bool success = false;
92 bool verbose = false;
93 bool debug = false;
94 try {
95 int info = 0;
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 filename("mhd1280b.cua");
105 std::string ortho("ICGS");
106 MT tol = 1.0e-5; // relative residual tolerance
107
108 CommandLineProcessor cmdp(false,true);
109 cmdp.setOption("verbose","quiet",&verbose,"Print messages and results.");
110 cmdp.setOption("debug","no-debug",&debug,"Print debug messages.");
111 cmdp.setOption("pseudo","regular",&pseudo,"Use pseudo-block GMRES to solve the linear systems.");
112 cmdp.setOption("frequency",&frequency,"Solvers frequency for printing residuals (#iters).");
113 cmdp.setOption("filename",&filename,"Filename for Harwell-Boeing test matrix.");
114 cmdp.setOption("tol",&tol,"Relative residual tolerance used by GMRES solver.");
115 cmdp.setOption("num-rhs",&numrhs,"Number of right-hand sides to be solved for.");
116 cmdp.setOption("num-restarts",&maxrestarts,"Maximum number of restarts allowed for the GMRES solver.");
117 cmdp.setOption("blocksize",&blocksize,"Block size used by GMRES.");
118 cmdp.setOption("subspace-length",&length,"Maximum dimension of block-subspace used by GMRES solver.");
119 cmdp.setOption("ortho-type",&ortho,"Orthogonalization type, either DGKS, ICGS or IMGS");
120 if (cmdp.parse(argc,argv) != CommandLineProcessor::PARSE_SUCCESSFUL) {
121 return EXIT_FAILURE;
122 }
123
124 proc_verbose = verbose && (MyPID==0); /* Only print on the zero processor */
125 if (proc_verbose) {
126 std::cout << Belos::Belos_Version() << std::endl << std::endl;
127 }
128 if (!verbose)
129 frequency = -1; // reset frequency if test is not verbose
130
131#ifndef HAVE_BELOS_TRIUTILS
132 std::cout << "This test requires Triutils. Please configure with --enable-triutils." << std::endl;
133 if (MyPID==0) {
134 std::cout << "End Result: TEST FAILED" << std::endl;
135 }
136 return EXIT_FAILURE;
137#endif
138
139 // Get the data from the HB file
140 int dim,dim2,nnz;
141 MT *dvals;
142 int *colptr,*rowind;
143 ST *cvals;
144 nnz = -1;
145 info = readHB_newmat_double(filename.c_str(),&dim,&dim2,&nnz,
146 &colptr,&rowind,&dvals);
147 if (info == 0 || nnz < 0) {
148 if (MyPID==0) {
149 std::cout << "Error reading '" << filename << "'" << std::endl;
150 std::cout << "End Result: TEST FAILED" << std::endl;
151 }
152 return EXIT_FAILURE;
153 }
154 // Convert interleaved doubles to std::complex values
155 cvals = new ST[nnz];
156 for (int ii=0; ii<nnz; ii++) {
157 cvals[ii] = ST(dvals[ii*2],dvals[ii*2+1]);
158 }
159 // Build the problem matrix
160 RCP< MyBetterOperator<ST> > A
161 = rcp( new MyBetterOperator<ST>(dim,colptr,nnz,rowind,cvals) );
162 //
163 // ********Other information used by block solver***********
164 // *****************(can be user specified)******************
165 //
166 int maxits = dim/blocksize; // maximum number of iterations to run
167 //
168 ParameterList belosList;
169 belosList.set( "Num Blocks", length ); // Maximum number of blocks in Krylov factorization
170 belosList.set( "Block Size", blocksize ); // Blocksize to be used by iterative solver
171 belosList.set( "Maximum Iterations", maxits ); // Maximum number of iterations allowed
172 belosList.set( "Maximum Restarts", maxrestarts ); // Maximum number of restarts allowed
173 belosList.set( "Convergence Tolerance", tol ); // Relative convergence tolerance requested
174 belosList.set( "Orthogonalization", ortho ); // Orthogonalization type
175 if (verbose) {
176 int verbosity = Belos::Errors + Belos::Warnings +
178 if (debug)
179 verbosity += Belos::OrthoDetails + Belos::Debug;
180 belosList.set( "Verbosity", verbosity );
181 if (frequency > 0)
182 belosList.set( "Output Frequency", frequency );
183 }
184 else
185 belosList.set( "Verbosity", Belos::Errors + Belos::Warnings );
186 //
187 // Construct the right-hand side and solution multivectors.
188 // NOTE: The right-hand side will be constructed such that the solution is
189 // a vectors of one.
190 //
191 RCP<MyMultiVec<ST> > soln = rcp( new MyMultiVec<ST>(dim,numrhs) );
192 RCP<MyMultiVec<ST> > rhs = rcp( new MyMultiVec<ST>(dim,numrhs) );
193 MVT::MvRandom( *soln );
194 OPT::Apply( *A, *soln, *rhs );
195 MVT::MvInit( *soln, zero );
196 //
197 // Construct an unpreconditioned linear problem instance.
198 //
199 RCP<Belos::LinearProblem<ST,MV,OP> > problem =
200 rcp( new Belos::LinearProblem<ST,MV,OP>( A, soln, rhs ) );
201 bool set = problem->setProblem();
202 if (set == false) {
203 if (proc_verbose)
204 std::cout << std::endl << "ERROR: Belos::LinearProblem failed to set up correctly!" << std::endl;
205 return EXIT_FAILURE;
206 }
207
208 // Use a debugging status test to save absolute residual history.
209 // Debugging status tests are peer to the native status tests that are called whenever convergence is checked.
211
212 //
213 // *******************************************************************
214 // *************Start the block Gmres iteration***********************
215 // *******************************************************************
216 //
217 Teuchos::RCP< Belos::SolverManager<ST,MV,OP> > solver;
218 if (pseudo)
219 solver = Teuchos::rcp( new Belos::PseudoBlockGmresSolMgr<ST,MV,OP>( problem, Teuchos::rcp(&belosList,false) ) );
220 else
221 solver = Teuchos::rcp( new Belos::BlockGmresSolMgr<ST,MV,OP>( problem, Teuchos::rcp(&belosList,false) ) );
222
223 solver->setDebugStatusTest( Teuchos::rcp(&debugTest, false) );
224
225 //
226 // **********Print out information about problem*******************
227 //
228 if (proc_verbose) {
229 std::cout << std::endl << std::endl;
230 std::cout << "Dimension of matrix: " << dim << std::endl;
231 std::cout << "Number of right-hand sides: " << numrhs << std::endl;
232 std::cout << "Block size used by solver: " << blocksize << std::endl;
233 std::cout << "Max number of Gmres iterations: " << maxits << std::endl;
234 std::cout << "Relative residual tolerance: " << tol << std::endl;
235 std::cout << std::endl;
236 }
237 //
238 // Perform solve
239 //
240 Belos::ReturnType ret = solver->solve();
241 //
242 // Compute actual residuals.
243 //
244 RCP<MyMultiVec<ST> > temp = rcp( new MyMultiVec<ST>(dim,numrhs) );
245 OPT::Apply( *A, *soln, *temp );
246 MVT::MvAddMv( one, *rhs, -one, *temp, *temp );
247 std::vector<MT> norm_num(numrhs), norm_denom(numrhs);
248 MVT::MvNorm( *temp, norm_num );
249 MVT::MvNorm( *rhs, norm_denom );
250 for (int i=0; i<numrhs; ++i) {
251 if (proc_verbose)
252 std::cout << "Relative residual "<<i<<" : " << norm_num[i] / norm_denom[i] << std::endl;
253 if ( norm_num[i] / norm_denom[i] > tol ) {
254 norm_failure = true;
255 }
256 }
257
258 // Print absolute residual norm logging.
259 const std::vector<MT> residualLog = debugTest.getLogResNorm();
260 if (numrhs==1 && proc_verbose && residualLog.size())
261 {
262 std::cout << "Absolute residual 2-norm [ " << residualLog.size() << " ] : ";
263 for (unsigned int i=0; i<residualLog.size(); i++)
264 std::cout << residualLog[i] << " ";
265 std::cout << std::endl;
266 std::cout << "Final abs 2-norm / rhs 2-norm : " << residualLog[residualLog.size()-1] / norm_denom[0] << std::endl;
267 }
268
269 // Clean up.
270 delete [] dvals;
271 delete [] colptr;
272 delete [] rowind;
273 delete [] cvals;
274
275 success = ret==Belos::Converged && !norm_failure;
276 if (success) {
277 if (proc_verbose)
278 std::cout << "End Result: TEST PASSED" << std::endl;
279 } else {
280 if (proc_verbose)
281 std::cout << "End Result: TEST FAILED" << std::endl;
282 }
283 }
284 TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success);
285
286 return ( success ? EXIT_SUCCESS : EXIT_FAILURE );
287} // end test_bl_gmres_complex_hb.cpp
The Belos::BlockGmresSolMgr provides a solver manager for the BlockGmres linear solver.
Belos header file which uses auto-configuration information to include necessary C++ headers.
Class which describes the linear problem to be solved by the iterative solver.
The Belos::PseudoBlockGmresSolMgr provides a solver manager for the BlockGmres linear solver.
Belos::StatusTest debugging class for storing the absolute residual norms generated during a solve.
Interface to Block GMRES and Flexible GMRES.
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.
Interface to standard and "pseudoblock" GMRES.
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.
@ OrthoDetails
@ StatusTestDetails
@ TimingDetails
ReturnType
Whether the Belos solve converged for all linear systems.
std::string Belos_Version()
int main(int argc, char *argv[])