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.
Class that helps parse command line input arguments from (argc,argv[]) and set options.
Initialize, finalize, and query the global MPI session.
static int getRank()
The rank of the calling process in MPI_COMM_WORLD.
Concrete serial communicator subclass.
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.
static magnitudeType magnitude(ScalarType a)
Returns the magnitudeType of the scalar type a.
static ScalarType one()
Returns representation of one for this scalar type.
T magnitudeType
Mandatory typedef for result of magnitude.
static ScalarType 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.