Belos Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
test_bl_gmres_complex_diag.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
57#ifdef HAVE_MPI
58#include <mpi.h>
59#endif
60
61#include "MyMultiVec.hpp"
62#include "MyOperator.hpp"
63
64using namespace Teuchos;
65
66int main(int argc, char *argv[]) {
67 //
68 typedef std::complex<double> ST;
69 typedef ScalarTraits<ST> SCT;
70 typedef SCT::magnitudeType MT;
71 typedef Belos::MultiVec<ST> MV;
72 typedef Belos::Operator<ST> OP;
75 ST one = SCT::one();
76 ST zero = SCT::zero();
77
78 Teuchos::GlobalMPISession session(&argc, &argv, NULL);
79 int MyPID = session.getRank();
80 //
81 using Teuchos::RCP;
82 using Teuchos::rcp;
83
84 bool success = false;
85 bool verbose = false;
86 try {
87 bool norm_failure = false;
88 bool proc_verbose = false;
89 bool pseudo = false; // use pseudo block GMRES to solve this linear system.
90 int frequency = -1; // how often residuals are printed by solver
91 int blocksize = 1;
92 int numrhs = 1;
93 int maxrestarts = 15;
94 int length = 50;
95 std::string ortho("DGKS"); // The Belos documentation obscures the fact that
96 MT tol = 1.0e-5; // relative residual tolerance
97
98 CommandLineProcessor cmdp(false,true);
99 cmdp.setOption("verbose","quiet",&verbose,"Print messages and results.");
100 cmdp.setOption("pseudo","regular",&pseudo,"Use pseudo-block GMRES to solve the linear systems.");
101 cmdp.setOption("frequency",&frequency,"Solvers frequency for printing residuals (#iters).");
102 cmdp.setOption("tol",&tol,"Relative residual tolerance used by GMRES solver.");
103 cmdp.setOption("num-rhs",&numrhs,"Number of right-hand sides to be solved for.");
104 cmdp.setOption("num-restarts",&maxrestarts,"Maximum number of restarts allowed for the GMRES solver.");
105 cmdp.setOption("blocksize",&blocksize,"Block size used by GMRES.");
106 cmdp.setOption("subspace-length",&length,"Maximum dimension of block-subspace used by GMRES solver.");
107 cmdp.setOption("ortho-type",&ortho,"Orthogonalization type, either DGKS, ICGS or IMGS");
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#ifndef HAVE_BELOS_TRIUTILS
120 std::cout << "This test requires Triutils. Please configure with --enable-triutils." << std::endl;
121 if (MyPID==0) {
122 std::cout << "End Result: TEST FAILED" << std::endl;
123 }
124 return EXIT_FAILURE;
125#endif
126
127 // Get the data from the HB file
128 int dim=100;
129
130 // Build the problem matrix
131 std::vector<ST> diag( dim, (ST)4.0 );
132 RCP< MyOperator<ST> > A
133 = rcp( new MyOperator<ST>( diag ) );
134 //
135 // ********Other information used by block solver***********
136 // *****************(can be user specified)******************
137 //
138 int maxits = dim/blocksize; // maximum number of iterations to run
139 //
140 ParameterList belosList;
141 belosList.set( "Num Blocks", length ); // Maximum number of blocks in Krylov factorization
142 belosList.set( "Block Size", blocksize ); // Blocksize to be used by iterative solver
143 belosList.set( "Maximum Iterations", maxits ); // Maximum number of iterations allowed
144 belosList.set( "Maximum Restarts", maxrestarts ); // Maximum number of restarts allowed
145 belosList.set( "Convergence Tolerance", tol ); // Relative convergence tolerance requested
146 belosList.set( "Orthogonalization", ortho ); // Orthogonalization type
147 if (verbose) {
148 belosList.set( "Verbosity", Belos::Errors + Belos::Warnings +
150 if (frequency > 0)
151 belosList.set( "Output Frequency", frequency );
152 }
153 else
154 belosList.set( "Verbosity", Belos::Errors + Belos::Warnings );
155 //
156 // Construct the right-hand side and solution multivectors.
157 // NOTE: The right-hand side will be constructed such that the solution is
158 // a vectors of one.
159 //
160 RCP<MyMultiVec<ST> > soln = rcp( new MyMultiVec<ST>(dim,numrhs) );
161 RCP<MyMultiVec<ST> > rhs = rcp( new MyMultiVec<ST>(dim,numrhs) );
162 MVT::MvInit( *rhs, 1.0 );
163 MVT::MvInit( *soln, zero );
164 //
165 // Construct an unpreconditioned linear problem instance.
166 //
167 RCP<Belos::LinearProblem<ST,MV,OP> > problem =
168 rcp( new Belos::LinearProblem<ST,MV,OP>( A, soln, rhs ) );
169 bool set = problem->setProblem();
170 if (set == false) {
171 if (proc_verbose)
172 std::cout << std::endl << "ERROR: Belos::LinearProblem failed to set up correctly!" << std::endl;
173 return EXIT_FAILURE;
174 }
175
176 // Use a debugging status test to save absolute residual history.
177 // Debugging status tests are peer to the native status tests that are called whenever convergence is checked.
179
180 //
181 // *******************************************************************
182 // *************Start the block Gmres iteration***********************
183 // *******************************************************************
184 //
185 Teuchos::RCP< Belos::SolverManager<ST,MV,OP> > solver;
186 if (pseudo)
187 solver = Teuchos::rcp( new Belos::PseudoBlockGmresSolMgr<ST,MV,OP>( problem, Teuchos::rcp(&belosList,false) ) );
188 else
189 solver = Teuchos::rcp( new Belos::BlockGmresSolMgr<ST,MV,OP>( problem, Teuchos::rcp(&belosList,false) ) );
190
191 solver->setDebugStatusTest( Teuchos::rcp(&debugTest, false) );
192
193 //
194 // **********Print out information about problem*******************
195 //
196 if (proc_verbose) {
197 std::cout << std::endl << std::endl;
198 std::cout << "Dimension of matrix: " << dim << std::endl;
199 std::cout << "Number of right-hand sides: " << numrhs << std::endl;
200 std::cout << "Block size used by solver: " << blocksize << std::endl;
201 std::cout << "Max number of Gmres iterations: " << maxits << std::endl;
202 std::cout << "Relative residual tolerance: " << tol << std::endl;
203 std::cout << std::endl;
204 }
205 //
206 // Perform solve
207 //
208 Belos::ReturnType ret = solver->solve();
209 //
210 // Compute actual residuals.
211 //
212 RCP<MyMultiVec<ST> > temp = rcp( new MyMultiVec<ST>(dim,numrhs) );
213 OPT::Apply( *A, *soln, *temp );
214 MVT::MvAddMv( one, *rhs, -one, *temp, *temp );
215 std::vector<MT> norm_num(numrhs), norm_denom(numrhs);
216 MVT::MvNorm( *temp, norm_num );
217 MVT::MvNorm( *rhs, norm_denom );
218 for (int i=0; i<numrhs; ++i) {
219 if (proc_verbose)
220 std::cout << "Relative residual "<<i<<" : " << norm_num[i] / norm_denom[i] << std::endl;
221 if ( norm_num[i] / norm_denom[i] > tol ) {
222 norm_failure = true;
223 }
224 }
225
226 // Print absolute residual norm logging.
227 const std::vector<MT> residualLog = debugTest.getLogResNorm();
228 if (numrhs==1 && proc_verbose && residualLog.size())
229 {
230 std::cout << "Absolute residual 2-norm [ " << residualLog.size() << " ] : ";
231 for (unsigned int i=0; i<residualLog.size(); i++)
232 std::cout << residualLog[i] << " ";
233 std::cout << std::endl;
234 std::cout << "Final abs 2-norm / rhs 2-norm : " << residualLog[residualLog.size()-1] / norm_denom[0] << std::endl;
235 }
236
237 success = ret==Belos::Converged && !norm_failure;
238 if (success) {
239 if (proc_verbose)
240 std::cout << "End Result: TEST PASSED" << std::endl;
241 } else {
242 if (proc_verbose)
243 std::cout << "End Result: TEST FAILED" << std::endl;
244 }
245 }
246 TEUCHOS_STANDARD_CATCH_STATEMENTS(verbose, std::cerr, success);
247
248 return ( success ? EXIT_SUCCESS : EXIT_FAILURE );
249} // 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::MultiVec class.
Definition: MyMultiVec.hpp:66
Simple example of a user's defined Belos::Operator class.
Definition: MyOperator.hpp:66
@ StatusTestDetails
Definition: BelosTypes.hpp:261
@ 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[])