66 template<
class OrdinalType,
class ScalarType>
69 typedef ScalarType scalar_type;
71 typedef MagnitudeType magnitude_type;
81 MagnitudeType trigErrorBound_;
87 MagnitudeType norm2 (
const ScalarType a,
const ScalarType& b)
91 if (scale == STM::zero()) {
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);
109 static MagnitudeType trigErrorBound ()
117 return 2 * (4*u) / (1 - 4*u);
124 GivensTester (std::ostream& out) :
125 out_ (out), trigErrorBound_ (trigErrorBound())
129 bool compare (ScalarType a, ScalarType b)
135 out_ <<
"Comparing Givens rotations for [a; b] = ["
136 << a <<
"; " << b <<
"]" << endl;
138 generic_blas_type genericBlas;
139 library_blas_type libraryBlas;
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);
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);
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;
164 out_ <<
"-- |c_generic - c_library| = "
168 out_ <<
"---- Difference exceeded error bound " << trigErrorBound_ << endl;
172 out_ <<
"-- |s_generic - s_library| = "
176 out_ <<
"---- Difference exceeded error bound " << trigErrorBound_ << endl;
184 const MagnitudeType inputNorm = norm2 (a, b);
185 const MagnitudeType outputDiffNorm =
186 norm2 (a_generic - a_library, b_generic - b_library);
188 out_ <<
"-- ||[a; b]||_2 = " << inputNorm << endl;
189 out_ <<
"-- ||[a_generic - a_library; b_generic - b_library]||_2 = "
190 << outputDiffNorm << endl;
198 const MagnitudeType two = STM::one() + STM::one();
199 const MagnitudeType fwdErrorBound =
202 if (outputDiffNorm > fwdErrorBound * inputNorm) {
204 out_ <<
"---- Forward error exceeded relative error bound "
205 << fwdErrorBound << endl;
225 success = success && compare (one, zero);
226 success = success && compare (zero, one);
227 success = success && compare (-one, zero);
228 success = success && compare (zero, -one);
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);
275 const int myRank = mpiSession.
getRank();
277 std::ostream& out = (myRank == 0) ? std::cout : blackHole;
280 CommandLineProcessor cmdp(
false,
true);
281 cmdp.setOption (
"verbose",
"quiet", &verbose,
"Print messages and results.");
285 const CommandLineProcessor::EParseCommandLineReturn parseResult =
286 cmdp.parse (argc,argv);
289 if (parseResult == CommandLineProcessor::PARSE_HELP_PRINTED) {
290 out <<
"End Result: TEST PASSED" << endl;
294 std::invalid_argument,
295 "Failed to parse command-line arguments");
299 GivensTester<int, double> tester (verbose ? out : blackHole);
300 const bool success = tester.test ();
303 out <<
"End Result: TEST PASSED" << endl;
306 out <<
"End Result: TEST FAILED" << endl;
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.
#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.