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