Teuchos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
default_blas_rot.cpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Teuchos: Common Tools Package
5// Copyright (2004) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
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
53
54#include <Teuchos_as.hpp>
55#include <Teuchos_BLAS.hpp>
59
60// Anonymous namespace to avoid name collisions.
61namespace {
62
66 template<class OrdinalType, class ScalarType>
67 class GivensTester {
68 public:
69 typedef ScalarType scalar_type;
70 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
71 typedef MagnitudeType magnitude_type;
72
73 private:
75 std::ostream& out_;
76
81 MagnitudeType trigErrorBound_;
82
85
87 MagnitudeType norm2 (const ScalarType a, const ScalarType& b)
88 {
89 const MagnitudeType scale = STS::magnitude(a) + STS::magnitude(b);
90
91 if (scale == STM::zero()) {
92 return STM::zero();
93 } else {
94 const ScalarType a_scaled = a / scale;
95 const ScalarType b_scaled = b / scale;
96 return scale * STM::squareroot (a_scaled*a_scaled + b_scaled*b_scaled);
97 }
98 }
99
109 static MagnitudeType trigErrorBound ()
110 {
111 // NOTE (mfh 12 Oct 2011) I'm not sure if this error bound holds
112 // for complex arithmetic. Does STS report the right machine
113 // precision for complex numbers?
114 const MagnitudeType u = STS::eps();
115
116 // In Higham's notation, $\gamma_k = ku / (1 - ku)$.
117 return 2 * (4*u) / (1 - 4*u);
118 }
119
120 public:
124 GivensTester (std::ostream& out) :
125 out_ (out), trigErrorBound_ (trigErrorBound())
126 {}
127
129 bool compare (ScalarType a, ScalarType b)
130 {
131 using std::endl;
132 typedef Teuchos::DefaultBLASImpl<int, ScalarType> generic_blas_type;
133 typedef Teuchos::BLAS<int, ScalarType> library_blas_type;
134
135 out_ << "Comparing Givens rotations for [a; b] = ["
136 << a << "; " << b << "]" << endl;
137
138 generic_blas_type genericBlas;
139 library_blas_type libraryBlas;
140
141 // ROTG overwrites its input arguments.
142 ScalarType a_generic = a;
143 ScalarType b_generic = b;
144 MagnitudeType c_generic;
145 ScalarType s_generic;
146 genericBlas.ROTG (&a_generic, &b_generic, &c_generic, &s_generic);
147
148 ScalarType a_library = a;
149 ScalarType b_library = b;
150 MagnitudeType c_library;
151 ScalarType s_library;
152 libraryBlas.ROTG (&a_library, &b_library, &c_library, &s_library);
153
154 out_ << "-- DefaultBLASImpl results: a,b,c,s = "
155 << a_generic << ", " << b_generic << ", "
156 << c_generic << ", " << s_generic << endl;
157 out_ << "-- (Library) BLAS results: a,b,c,s = "
158 << a_library << ", " << b_library << ", "
159 << c_library << ", " << s_library << endl;
160
161 bool success = true; // Innocent until proven guilty.
162
163 // Test the difference between the computed cosines.
164 out_ << "-- |c_generic - c_library| = "
165 << STS::magnitude(c_generic - c_library) << endl;
166 if (STS::magnitude(c_generic - c_library) > trigErrorBound_) {
167 success = false;
168 out_ << "---- Difference exceeded error bound " << trigErrorBound_ << endl;
169 }
170
171 // Test the difference between the computed sines.
172 out_ << "-- |s_generic - s_library| = "
173 << STS::magnitude(s_generic - s_library) << endl;
174 if (STS::magnitude(s_generic - s_library) > trigErrorBound_) {
175 success = false;
176 out_ << "---- Difference exceeded error bound " << trigErrorBound_ << endl;
177 }
178
179 // Test the forward error of applying the Givens rotation.
180 // Remember that ROTG applies the rotation to its input
181 // arguments [a; b], overwriting them with the resulting [r; z].
182 //
183 // See Higham's Lemma 19.8.
184 const MagnitudeType inputNorm = norm2 (a, b);
185 const MagnitudeType outputDiffNorm =
186 norm2 (a_generic - a_library, b_generic - b_library);
187
188 out_ << "-- ||[a; b]||_2 = " << inputNorm << endl;
189 out_ << "-- ||[a_generic - a_library; b_generic - b_library]||_2 = "
190 << outputDiffNorm << endl;
191
192 // Multiply by a fudge factor of the base, just in case the
193 // forward error bound wasn't computed accurately. Also
194 // multiply by 2, since we don't know the exact result. The
195 // latter is because the two computed results could be on either
196 // side of the exact result: sqrt((2 * x_diff)^2 + (2 *
197 // y_diff)^2) = sqrt(4) * sqrt(x_diff^2 + y_diff^2).
198 const MagnitudeType two = STM::one() + STM::one();
199 const MagnitudeType fwdErrorBound =
200 2 * STS::base() * STM::squareroot(two) * (6*STS::eps() / (1 - 6*STS::eps()));
201
202 if (outputDiffNorm > fwdErrorBound * inputNorm) {
203 success = false;
204 out_ << "---- Forward error exceeded relative error bound "
205 << fwdErrorBound << endl;
206 }
207 return success;
208 }
209
213 bool test ()
214 {
215 using Teuchos::as;
216
217 const ScalarType zero = STS::zero();
218 const ScalarType one = STS::one();
219 //const ScalarType two = one + one;
220 //const ScalarType four = two + two;
221
222 bool success = true; // Innocent until proven guilty.
223
224 // First test the corner cases: [\pm 1, 0] and [0, \pm 1].
225 success = success && compare (one, zero);
226 success = success && compare (zero, one);
227 success = success && compare (-one, zero);
228 success = success && compare (zero, -one);
229
230 // Test a range of other values.
231 {
232 const ScalarType incr = one / as<ScalarType> (10);
233 for (int k = -30; k < 30; ++k) {
234 const ScalarType a = as<ScalarType> (k) * incr;
235 const ScalarType b = one - as<ScalarType> (k) * incr;
236 success = success && compare (a, b);
237 }
238 }
239
240 //
241 // Try some big values just to see whether ROTG correctly
242 // scales its inputs to avoid overflow.
243 //
244 //success = success && compare (STS::rmax() / four, STS::rmax() / four);
245 //success = success && compare (-STS::rmax() / four, STS::rmax() / four);
246 //success = success && compare (STS::rmax() / four, -STS::rmax() / four);
247
248 //success = success && compare (STS::rmax() / two, STS::rmax() / two);
249 //success = success && compare (-STS::rmax() / two, STS::rmax() / two);
250 //success = success && compare (-STS::rmax() / two, -STS::rmax() / two);
251
252 //success = success && compare (STS::rmax() / two, zero);
253 //success = success && compare (zero, STS::rmax() / two);
254 //success = success && compare (-STS::rmax() / two, zero);
255 //success = success && compare (zero, -STS::rmax() / two);
256
257 //success = success && compare (STS::rmax() / two, one);
258 //success = success && compare (one, STS::rmax() / two);
259 //success = success && compare (-STS::rmax() / two, one);
260 //success = success && compare (one, -STS::rmax() / two);
261
262 return success;
263 }
264 };
265} // namespace (anonymous)
266
267
268int
269main (int argc, char *argv[])
270{
271 using std::endl;
273
274 Teuchos::GlobalMPISession mpiSession (&argc, &argv, NULL);
275 const int myRank = mpiSession.getRank();
277 std::ostream& out = (myRank == 0) ? std::cout : blackHole;
278
279 bool verbose = true;
280 CommandLineProcessor cmdp(false,true);
281 cmdp.setOption ("verbose", "quiet", &verbose, "Print messages and results.");
282
283 // Parse the command-line arguments.
284 {
285 const CommandLineProcessor::EParseCommandLineReturn parseResult =
286 cmdp.parse (argc,argv);
287 // If the caller asks us to print the documentation, let the
288 // "test" pass trivially.
289 if (parseResult == CommandLineProcessor::PARSE_HELP_PRINTED) {
290 out << "End Result: TEST PASSED" << endl;
291 return EXIT_SUCCESS;
292 }
293 TEUCHOS_TEST_FOR_EXCEPTION(parseResult != CommandLineProcessor::PARSE_SUCCESSFUL,
294 std::invalid_argument,
295 "Failed to parse command-line arguments");
296 }
297
298 // Only let the tester print if in verbose mode.
299 GivensTester<int, double> tester (verbose ? out : blackHole);
300 const bool success = tester.test ();
301
302 if (success) {
303 out << "End Result: TEST PASSED" << endl;
304 return EXIT_SUCCESS;
305 } else {
306 out << "End Result: TEST FAILED" << endl;
307 return EXIT_FAILURE;
308 }
309}
310
Templated interface class to BLAS routines.
Basic command line parser for input from (argc,argv[])
A MPI utilities class, providing methods for initializing, finalizing, and querying the global MPI se...
Definition of Teuchos::as, for conversions between types.
Templated BLAS wrapper.
Class that helps parse command line input arguments from (argc,argv[]) and set options.
Default implementation for BLAS routines.
void ROTG(ScalarType *da, ScalarType *db, rotg_c_type *c, ScalarType *s) const
Computes a Givens plane rotation.
Initialize, finalize, and query the global MPI session.
static int getRank()
The rank of the calling process in MPI_COMM_WORLD.
basic_ostream<> subclass that does nothing but discard output.
int main()
Definition: evilMain.cpp:75
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Macro for throwing an exception with breakpointing to ease debugging.
TypeTo as(const TypeFrom &t)
Convert from one value type to another.
This structure defines some basic traits for a scalar field type.
static magnitudeType magnitude(ScalarType a)
Returns the magnitudeType of the scalar type a.
static T one()
Returns representation of one for this scalar type.
T magnitudeType
Mandatory typedef for result of magnitude.
static T zero()
Returns representation of zero for this scalar type.
static magnitudeType base()
Returns the base of the machine.
static magnitudeType eps()
Returns relative machine precision.
static T squareroot(T x)
Returns a number of magnitudeType that is the square root of this scalar type x.