Teuchos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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.
61 namespace {
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 
268 int
269 main (int argc, char *argv[])
270 {
271  using std::endl;
273 
274  Teuchos::GlobalMPISession mpiSession (&argc, &argv, NULL);
275  const int myRank = mpiSession.getRank();
276  Teuchos::oblackholestream blackHole;
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 
static int getRank()
The rank of the calling process in MPI_COMM_WORLD.
T magnitudeType
Mandatory typedef for result of magnitude.
Templated interface class to BLAS routines.
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.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Macro for throwing an exception with breakpointing to ease debugging.
basic_ostream&lt;&gt; subclass that does nothing but discard output.
Templated BLAS wrapper.
static magnitudeType base()
Returns the base of the machine.
Initialize, finalize, and query the global MPI session.
TypeTo as(const TypeFrom &t)
Convert from one value type to another.
int main(int argc, char *argv[])
static magnitudeType magnitude(ScalarTypea)
Returns the magnitudeType of the scalar type a.
Default implementation for BLAS routines.
A MPI utilities class, providing methods for initializing, finalizing, and querying the global MPI se...
void ROTG(ScalarType *da, ScalarType *db, rotg_c_type *c, ScalarType *s) const
Computes a Givens plane rotation.
Basic command line parser for input from (argc,argv[])
static T zero()
Returns representation of zero for this scalar type.
Definition of Teuchos::as, for conversions between types.
static T one()
Returns representation of one for this scalar type.
Class that helps parse command line input arguments from (argc,argv[]) and set options.