Belos Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
test_pseudo_tfqmr_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"
51#include "Teuchos_CommandLineProcessor.hpp"
52#include "Teuchos_ParameterList.hpp"
53#include "Teuchos_StandardCatchMacros.hpp"
54
55#ifdef HAVE_MPI
56#include <mpi.h>
57#endif
58
59// I/O for Harwell-Boeing files
60#ifdef HAVE_BELOS_TRIUTILS
61#include "Trilinos_Util_iohb.h"
62#endif
63
64#include "MyMultiVec.hpp"
65#include "MyBetterOperator.hpp"
66#include "MyOperator.hpp"
67
68using namespace Teuchos;
69
70int main(int argc, char *argv[]) {
71 //
72 typedef std::complex<double> ST;
73 typedef ScalarTraits<ST> SCT;
74 typedef SCT::magnitudeType MT;
75 typedef Belos::MultiVec<ST> MV;
76 typedef Belos::Operator<ST> OP;
79 ST one = SCT::one();
80 ST zero = SCT::zero();
81
82 int info = 0;
83 bool norm_failure = false;
84
85 Teuchos::GlobalMPISession session(&argc, &argv, NULL);
86 //
87 int MyPID = session.getRank();
88
89 using Teuchos::RCP;
90 using Teuchos::rcp;
91
92 bool success = false;
93 bool verbose = false;
94 try {
95 bool proc_verbose = false;
96 int frequency = -1; // how often residuals are printed by solver
97 int blocksize = 1;
98 int numrhs = 1;
99 std::string filename("mhd1280b.cua");
100 MT tol = 1.0e-5; // relative residual tolerance
101
102 CommandLineProcessor cmdp(false, true);
103 cmdp.setOption("verbose","quiet",&verbose,"Print messages and results.");
104 cmdp.setOption("frequency",&frequency,"Solvers frequency for printing residuals (#iters).");
105 cmdp.setOption("filename",&filename,"Filename for Harwell-Boeing test matrix.");
106 cmdp.setOption("tol",&tol,"Relative residual tolerance used by pseudo-block TFQMR solver.");
107 cmdp.setOption("num-rhs",&numrhs,"Number of right-hand sides to be solved for.");
108 if (cmdp.parse(argc,argv) != CommandLineProcessor::PARSE_SUCCESSFUL) {
109 return EXIT_FAILURE;
110 }
111
112 proc_verbose = verbose && (MyPID==0); /* Only print on the zero processor */
113 if (proc_verbose) {
114 std::cout << Belos::Belos_Version() << std::endl << std::endl;
115 }
116 if (!verbose)
117 frequency = -1; // reset frequency if test is not verbose
118
119
120#ifndef HAVE_BELOS_TRIUTILS
121 std::cout << "This test requires Triutils. Please configure with --enable-triutils." << std::endl;
122 if (MyPID==0) {
123 std::cout << "End Result: TEST FAILED" << std::endl;
124 }
125 return EXIT_FAILURE;
126#endif
127
128 // Get the data from the HB file
129 int dim,dim2,nnz;
130 MT *dvals;
131 int *colptr,*rowind;
132 ST *cvals;
133 nnz = -1;
134 info = readHB_newmat_double(filename.c_str(),&dim,&dim2,&nnz,
135 &colptr,&rowind,&dvals);
136 if (info == 0 || nnz < 0) {
137 if (MyPID==0) {
138 std::cout << "Error reading '" << filename << "'" << std::endl;
139 std::cout << "End Result: TEST FAILED" << std::endl;
140 }
141 return EXIT_FAILURE;
142 }
143 // Convert interleaved doubles to std::complex values
144 cvals = new ST[nnz];
145 for (int ii=0; ii<nnz; ii++) {
146 cvals[ii] = ST(dvals[ii*2],dvals[ii*2+1]);
147 }
148 // Build the problem matrix
149 RCP< MyBetterOperator<ST> > A
150 = rcp( new MyBetterOperator<ST>(dim,colptr,nnz,rowind,cvals) );
151 //
152 // ********Other information used by block solver***********
153 // *****************(can be user specified)******************
154 //
155 int maxits = dim/blocksize; // maximum number of iterations to run
156 //
157 ParameterList belosList;
158 belosList.set( "Maximum Iterations", maxits ); // Maximum number of iterations allowed
159 belosList.set( "Convergence Tolerance", tol ); // Relative convergence tolerance requested
160 if (verbose) {
161 belosList.set( "Verbosity", Belos::Errors + Belos::Warnings +
163 if (frequency > 0)
164 belosList.set( "Output Frequency", frequency );
165 }
166 else
167 belosList.set( "Verbosity", Belos::Errors + Belos::Warnings );
168 //
169 // Construct the right-hand side and solution multivectors.
170 // NOTE: The right-hand side will be constructed such that the solution is
171 // a vectors of one.
172 //
173 RCP<MyMultiVec<ST> > soln = rcp( new MyMultiVec<ST>(dim,numrhs) );
174 RCP<MyMultiVec<ST> > rhs = rcp( new MyMultiVec<ST>(dim,numrhs) );
175 MVT::MvRandom( *soln );
176 OPT::Apply( *A, *soln, *rhs );
177 MVT::MvInit( *soln, zero );
178 //
179 // Construct an unpreconditioned linear problem instance.
180 //
181 RCP<Belos::LinearProblem<ST,MV,OP> > problem =
182 rcp( new Belos::LinearProblem<ST,MV,OP>( A, soln, 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 // *******************************************************************
191 // *************Start the TFQMR iteration***********************
192 // *******************************************************************
193 //
194 Belos::PseudoBlockTFQMRSolMgr<ST,MV,OP> solver( problem, rcp(&belosList,false) );
195
196 //
197 // **********Print out information about problem*******************
198 //
199 if (proc_verbose) {
200 std::cout << std::endl << std::endl;
201 std::cout << "Dimension of matrix: " << dim << std::endl;
202 std::cout << "Number of right-hand sides: " << numrhs << std::endl;
203 std::cout << "Block size used by solver: " << blocksize << std::endl;
204 std::cout << "Max number of pseudo-block TFQMR iterations: " << maxits << std::endl;
205 std::cout << "Relative residual tolerance: " << tol << std::endl;
206 std::cout << std::endl;
207 }
208 //
209 // Perform solve
210 //
211 Belos::ReturnType ret = solver.solve();
212 //
213 // Compute actual residuals.
214 //
215 RCP<MyMultiVec<ST> > temp = rcp( new MyMultiVec<ST>(dim,numrhs) );
216 OPT::Apply( *A, *soln, *temp );
217 MVT::MvAddMv( one, *rhs, -one, *temp, *temp );
218 std::vector<MT> norm_num(numrhs), norm_denom(numrhs);
219 MVT::MvNorm( *temp, norm_num );
220 MVT::MvNorm( *rhs, norm_denom );
221 for (int i=0; i<numrhs; ++i) {
222 if (proc_verbose)
223 std::cout << "Relative residual "<<i<<" : " << norm_num[i] / norm_denom[i] << std::endl;
224 if ( norm_num[i] / norm_denom[i] > tol ) {
225 norm_failure = true;
226 }
227 }
228
229 // Clean up.
230 delete [] dvals;
231 delete [] colptr;
232 delete [] rowind;
233 delete [] cvals;
234
235 success = ret==Belos::Converged && !norm_failure;
236 if (success) {
237 if (proc_verbose)
238 std::cout << "End Result: TEST PASSED" << std::endl;
239 } else {
240 if (proc_verbose)
241 std::cout << "End Result: TEST FAILED" << std::endl;
242 }
243 }
244 TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success);
245
246 return ( success ? EXIT_SUCCESS : EXIT_FAILURE );
247} // end test_tfqmr_complex_hb.cpp
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::PseudoBlockTFQMRSolMgr provides a solver manager for the pseudo-block TFQMR 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.
The Belos::PseudoBlockTFQMRSolMgr provides a powerful and fully-featured solver manager over the pseu...
ReturnType solve() override
This method performs possibly repeated calls to the underlying linear solver's iterate() routine unti...
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[])