Belos Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
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"
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
67using namespace Teuchos;
68
69int 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;
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
Belos header file which uses auto-configuration information to include necessary C++ headers.
Declaration and definition of Belos::GCRODRSolMgr, which implements the GCRODR (recycling GMRES) solv...
Class which describes the linear problem to be solved by the iterative solver.
Implementation of the GCRODR (Recycling GMRES) iterative linear solver.
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.
Simple example of a user's defined Belos::Operator class.
Simple example of a user's defined Belos::MultiVec class.
@ StatusTestDetails
@ TimingDetails
ReturnType
Whether the Belos solve converged for all linear systems.
std::string Belos_Version()
int main(int argc, char *argv[])