Teuchos Package Browser (Single Doxygen Collection) Version of the Day
default_blas_rot.cpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //                    Teuchos: Common Tools Package
00005 //                 Copyright (2004) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ***********************************************************************
00040 // @HEADER
00041 
00053 
00054 #include <Teuchos_as.hpp>
00055 #include <Teuchos_BLAS.hpp>
00056 #include <Teuchos_CommandLineProcessor.hpp>
00057 #include <Teuchos_GlobalMPISession.hpp>
00058 #include <Teuchos_oblackholestream.hpp>
00059 
00060 // Anonymous namespace to avoid name collisions.
00061 namespace {
00062 
00066   template<class OrdinalType, class ScalarType>
00067   class GivensTester {
00068   public:
00069     typedef ScalarType scalar_type;
00070     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00071     typedef MagnitudeType magnitude_type;
00072 
00073   private:
00075     std::ostream& out_;
00076 
00081     MagnitudeType trigErrorBound_;
00082 
00083     typedef Teuchos::ScalarTraits<ScalarType> STS;
00084     typedef Teuchos::ScalarTraits<MagnitudeType> STM;
00085 
00087     MagnitudeType norm2 (const ScalarType a, const ScalarType& b)
00088     {
00089       const MagnitudeType scale = STS::magnitude(a) + STS::magnitude(b);
00090 
00091       if (scale == STM::zero()) {
00092   return STM::zero();
00093       } else {
00094   const ScalarType a_scaled = a / scale;
00095   const ScalarType b_scaled = b / scale;
00096   return scale * STM::squareroot (a_scaled*a_scaled + b_scaled*b_scaled);
00097       }
00098     }
00099 
00109     static MagnitudeType trigErrorBound ()
00110     {
00111       // NOTE (mfh 12 Oct 2011) I'm not sure if this error bound holds
00112       // for complex arithmetic.  Does STS report the right machine
00113       // precision for complex numbers?
00114       const MagnitudeType u = STS::eps();
00115 
00116       // In Higham's notation, $\gamma_k = ku / (1 - ku)$.
00117       return 2 * (4*u) / (1 - 4*u);
00118     }
00119 
00120   public:
00124     GivensTester (std::ostream& out) :
00125       out_ (out), trigErrorBound_ (trigErrorBound())
00126     {}
00127 
00129     bool compare (ScalarType a, ScalarType b) 
00130     {
00131       using std::endl;
00132       typedef Teuchos::DefaultBLASImpl<int, ScalarType> generic_blas_type;
00133       typedef Teuchos::BLAS<int, ScalarType> library_blas_type;
00134 
00135       out_ << "Comparing Givens rotations for [a; b] = [" 
00136      << a << "; " << b << "]" << endl;
00137 
00138       generic_blas_type genericBlas;
00139       library_blas_type libraryBlas;
00140 
00141       // ROTG overwrites its input arguments.
00142       ScalarType a_generic = a;
00143       ScalarType b_generic = b;
00144       MagnitudeType c_generic;
00145       ScalarType s_generic;
00146       genericBlas.ROTG (&a_generic, &b_generic, &c_generic, &s_generic);
00147 
00148       ScalarType a_library = a;
00149       ScalarType b_library = b;
00150       MagnitudeType c_library;
00151       ScalarType s_library;
00152       libraryBlas.ROTG (&a_library, &b_library, &c_library, &s_library);
00153 
00154       out_ << "-- DefaultBLASImpl results: a,b,c,s = " 
00155      << a_generic << ", " << b_generic << ", " 
00156      << c_generic << ", " << s_generic << endl;
00157       out_ << "-- (Library) BLAS results: a,b,c,s = " 
00158      << a_library << ", " << b_library << ", " 
00159      << c_library << ", " << s_library << endl;
00160 
00161       bool success = true; // Innocent until proven guilty.
00162 
00163       // Test the difference between the computed cosines.
00164       out_ << "-- |c_generic - c_library| = " 
00165      << STS::magnitude(c_generic - c_library) << endl;
00166       if (STS::magnitude(c_generic - c_library) > trigErrorBound_) {
00167   success = false;
00168   out_ << "---- Difference exceeded error bound " << trigErrorBound_ << endl;
00169       }
00170 
00171       // Test the difference between the computed sines.
00172       out_ << "-- |s_generic - s_library| = " 
00173      << STS::magnitude(s_generic - s_library) << endl;
00174       if (STS::magnitude(s_generic - s_library) > trigErrorBound_) {
00175   success = false;
00176   out_ << "---- Difference exceeded error bound " << trigErrorBound_ << endl;
00177       } 
00178       
00179       // Test the forward error of applying the Givens rotation.
00180       // Remember that ROTG applies the rotation to its input
00181       // arguments [a; b], overwriting them with the resulting [r; z].
00182       //
00183       // See Higham's Lemma 19.8.
00184       const MagnitudeType inputNorm = norm2 (a, b);
00185       const MagnitudeType outputDiffNorm = 
00186   norm2 (a_generic - a_library, b_generic - b_library);
00187 
00188       out_ << "-- ||[a; b]||_2 = " << inputNorm << endl;
00189       out_ << "-- ||[a_generic - a_library; b_generic - b_library]||_2 = " 
00190      << outputDiffNorm << endl;
00191       
00192       // Multiply by a fudge factor of the base, just in case the
00193       // forward error bound wasn't computed accurately.  Also
00194       // multiply by 2, since we don't know the exact result.  The
00195       // latter is because the two computed results could be on either
00196       // side of the exact result: sqrt((2 * x_diff)^2 + (2 *
00197       // y_diff)^2) = sqrt(4) * sqrt(x_diff^2 + y_diff^2).
00198       const MagnitudeType two = STM::one() + STM::one();
00199       const MagnitudeType fwdErrorBound = 
00200   2 * STS::base() * STM::squareroot(two) * (6*STS::eps() / (1 - 6*STS::eps()));
00201 
00202       if (outputDiffNorm > fwdErrorBound * inputNorm) {
00203   success = false;
00204   out_ << "---- Forward error exceeded relative error bound "
00205        << fwdErrorBound << endl;
00206       }
00207       return success;
00208     }
00209 
00213     bool test ()
00214     {
00215       using Teuchos::as;
00216 
00217       const ScalarType zero = STS::zero();
00218       const ScalarType one = STS::one();
00219       const ScalarType two = one + one;
00220       const ScalarType four = two + two;
00221 
00222       bool success = true; // Innocent until proven guilty.
00223 
00224       // First test the corner cases: [\pm 1, 0] and [0, \pm 1].
00225       success = success && compare (one, zero);
00226       success = success && compare (zero, one);
00227       success = success && compare (-one, zero);
00228       success = success && compare (zero, -one);
00229 
00230       // Test a range of other values.
00231       {
00232   const ScalarType incr = one / as<ScalarType> (10);
00233   for (int k = -30; k < 30; ++k) {
00234     const ScalarType a = as<ScalarType> (k) * incr;
00235     const ScalarType b = one - as<ScalarType> (k) * incr;
00236     success = success && compare (a, b);
00237   }
00238       }
00239 
00240       //
00241       // Try some big values just to see whether ROTG correctly
00242       // scales its inputs to avoid overflow.
00243       //
00244       success = success && compare (STS::rmax() / four, STS::rmax() / four);
00245       success = success && compare (-STS::rmax() / four, STS::rmax() / four);
00246       success = success && compare (STS::rmax() / four, -STS::rmax() / four);
00247 
00248       success = success && compare (STS::rmax() / two, STS::rmax() / two);
00249       success = success && compare (-STS::rmax() / two, STS::rmax() / two);
00250       success = success && compare (-STS::rmax() / two, -STS::rmax() / two);
00251 
00252       success = success && compare (STS::rmax() / two, zero);
00253       success = success && compare (zero, STS::rmax() / two);
00254       success = success && compare (-STS::rmax() / two, zero);
00255       success = success && compare (zero, -STS::rmax() / two);
00256 
00257       success = success && compare (STS::rmax() / two, one);
00258       success = success && compare (one, STS::rmax() / two);
00259       success = success && compare (-STS::rmax() / two, one);
00260       success = success && compare (one, -STS::rmax() / two);
00261 
00262       return success;
00263     }
00264   };
00265 } // namespace (anonymous)
00266 
00267 
00268 int 
00269 main (int argc, char *argv[])
00270 {
00271   using std::endl;
00272   using Teuchos::CommandLineProcessor;
00273 
00274   Teuchos::GlobalMPISession mpiSession (&argc, &argv, NULL);
00275   const int myRank = mpiSession.getRank();
00276   Teuchos::oblackholestream blackHole;
00277   std::ostream& out = (myRank == 0) ? std::cout : blackHole;
00278 
00279   bool verbose = true;
00280   CommandLineProcessor cmdp(false,true);
00281   cmdp.setOption ("verbose", "quiet", &verbose, "Print messages and results.");
00282 
00283   // Parse the command-line arguments.
00284   {
00285     const CommandLineProcessor::EParseCommandLineReturn parseResult = 
00286       cmdp.parse (argc,argv);
00287     // If the caller asks us to print the documentation, let the
00288     // "test" pass trivially.
00289     if (parseResult == CommandLineProcessor::PARSE_HELP_PRINTED) {
00290       out << "End Result: TEST PASSED" << endl;
00291       return EXIT_SUCCESS;
00292     }
00293     TEUCHOS_TEST_FOR_EXCEPTION(parseResult != CommandLineProcessor::PARSE_SUCCESSFUL, 
00294            std::invalid_argument, 
00295            "Failed to parse command-line arguments");
00296   }
00297 
00298   // Only let the tester print if in verbose mode.
00299   GivensTester<int, double> tester (verbose ? out : blackHole);
00300   const bool success = tester.test ();
00301 
00302   if (success) {
00303     out << "End Result: TEST PASSED" << endl;
00304     return EXIT_SUCCESS;
00305   } else {
00306     out << "End Result: TEST FAILED" << endl;
00307     return EXIT_FAILURE;
00308   }
00309 }
00310 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines