Teuchos - Trilinos Tools Package Version of the Day
Teuchos_ScalarTraits.hpp
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 
00042 // Kris
00043 // 06.18.03 -- Minor formatting changes
00044 //          -- Changed calls to LAPACK objects to use new <OType, SType> templates
00045 // 07.08.03 -- Move into Teuchos package/namespace
00046 // 07.11.03 -- Added ScalarTraits for ARPREC mp_real
00047 // 07.14.03 -- Fixed int rand() function (was set up to return a floating-point style random number)
00048 // 07.17.03 -- Added squareroot() function
00049 
00050 // NOTE: Before adding specializations of ScalarTraits, make sure that they do not duplicate 
00051 // specializations already present in PyTrilinos (see packages/PyTrilinos/src/Teuchos_Traits.i)
00052 
00053 // NOTE: halfPrecision and doublePrecision are not currently implemented for ARPREC, GMP or the ordinal types (e.g., int, char)
00054 
00055 #ifndef _TEUCHOS_SCALARTRAITS_HPP_
00056 #define _TEUCHOS_SCALARTRAITS_HPP_
00057 
00062 #include "Teuchos_ConfigDefs.hpp"
00063 
00064 #ifdef HAVE_TEUCHOS_ARPREC
00065 #include <arprec/mp_real.h>
00066 #endif
00067 
00068 #ifdef HAVE_TEUCHOS_QD
00069 #include <qd/qd_real.h>
00070 #include <qd/dd_real.h>
00071 #endif
00072 
00073 #ifdef HAVE_TEUCHOS_GNU_MP
00074 #include <gmp.h>
00075 #include <gmpxx.h>
00076 #endif
00077 
00078 
00079 #include "Teuchos_ScalarTraitsDecl.hpp"
00080 
00081 
00082 namespace Teuchos {
00083 
00084 
00085 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00086 
00087 
00088 void throwScalarTraitsNanInfError( const std::string &errMsg );
00089 
00090 
00091 template<class Scalar>
00092 bool generic_real_isnaninf(const Scalar &x)
00093 {
00094   typedef std::numeric_limits<Scalar> NL;
00095   // IEEE says this should fail for NaN (not all compilers do not obey IEEE)
00096   const Scalar tol = 1.0; // Any (bounded) number should do!
00097   if (!(x <= tol) && !(x > tol)) return true;
00098   // Use fact that Inf*0 = NaN (again, all compilers do not obey IEEE)
00099   Scalar z = static_cast<Scalar>(0.0) * x;
00100   if (!(z <= tol) && !(z > tol)) return true;
00101   // As a last result use comparisons
00102   if (x == NL::infinity() || x == -NL::infinity()) return true;
00103   // We give up and assume the number is finite
00104   return false;
00105 }
00106 
00107 
00108 #define TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR( VALUE, MSG ) \
00109   if (isnaninf(VALUE)) { \
00110     std::ostringstream omsg; \
00111     omsg << MSG; \
00112     Teuchos::throwScalarTraitsNanInfError(omsg.str());  \
00113   }
00114 
00115 
00116 template<>
00117 struct ScalarTraits<char>
00118 {
00119   typedef char magnitudeType;
00120   typedef char halfPrecision;
00121   typedef char doublePrecision;
00122   static const bool isComplex = false;
00123   static const bool isOrdinal = true;
00124   static const bool isComparable = true;
00125   static const bool hasMachineParameters = false;
00126   // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
00127   static inline magnitudeType magnitude(char a) { return static_cast<char>(std::fabs(static_cast<double>(a))); }
00128   static inline char zero()  { return 0; }
00129   static inline char one()   { return 1; }
00130   static inline char conjugate(char x) { return x; }
00131   static inline char real(char x) { return x; }
00132   static inline char imag(char) { return 0; }
00133   static inline bool isnaninf(char ) { return false; }
00134   static inline void seedrandom(unsigned int s) { 
00135     std::srand(s); 
00136 #ifdef __APPLE__
00137     // throw away first random number to address bug 3655
00138     // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
00139     random();
00140 #endif
00141   }
00142   //static inline char random() { return (-1 + 2*rand()); } // RAB: This version should be used to be consistent with others
00143   static inline char random() { return std::rand(); } // RAB: This version should be used for an unsigned char, not char
00144   static inline std::string name() { return "char"; }
00145   static inline char squareroot(char x) { return (char) std::sqrt((double) x); }
00146   static inline char pow(char x, char y) { return (char) std::pow((double)x,(double)y); }
00147   static inline char log(char x) { return static_cast<char> (std::log (static_cast<double> (x))); }
00148   static inline char log10(char x) { return static_cast<char> (std::log10 (static_cast<double> (x))); }
00149 };
00150 
00151 
00152 template<>
00153 struct ScalarTraits<short int>
00154 {
00155   typedef short int magnitudeType;
00156   typedef short int halfPrecision;
00157   typedef short int doublePrecision;
00158   static const bool isComplex = false;
00159   static const bool isOrdinal = true;
00160   static const bool isComparable = true;
00161   static const bool hasMachineParameters = false;
00162   // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
00163   static inline magnitudeType magnitude(short int a) { return static_cast<short int>(std::fabs(static_cast<double>(a))); }
00164   static inline short int zero()  { return 0; }
00165   static inline short int one()   { return 1; }
00166   static inline short int conjugate(short int x) { return x; }
00167   static inline short int real(short int x) { return x; }
00168   static inline short int imag(short int) { return 0; }
00169   static inline bool isnaninf(short int) { return false; }
00170   static inline void seedrandom(unsigned int s) { 
00171     std::srand(s); 
00172 #ifdef __APPLE__
00173     // throw away first random number to address bug 3655
00174     // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
00175     random();
00176 #endif
00177   }
00178   //static inline int random() { return (-1 + 2*rand()); }  // RAB: This version should be used to be consistent with others
00179   static inline short int random() { return std::rand(); }             // RAB: This version should be used for an unsigned int, not int
00180   static inline std::string name() { return "short int"; }
00181   static inline short int squareroot(short int x) { return (short int) std::sqrt((double) x); }
00182   static inline short int pow(short int x, short int y) { return (short int) std::pow((double)x,(double)y); }
00183   static inline short int log(short int x) { return static_cast<short int> (std::log (static_cast<double> (x))); }
00184   static inline short int log10(short int x) { return static_cast<short int> (std::log10 (static_cast<double> (x))); }
00185 };
00186 
00187 template<>
00188 struct ScalarTraits<unsigned short int>
00189 {
00190   typedef unsigned short int magnitudeType;
00191   typedef unsigned short int halfPrecision;
00192   typedef unsigned short int doublePrecision;
00193   static const bool isComplex = false;
00194   static const bool isOrdinal = true;
00195   static const bool isComparable = true;
00196   static const bool hasMachineParameters = false;
00197   // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
00198   static inline magnitudeType magnitude(unsigned short int a) { return static_cast<unsigned short int>(std::fabs(static_cast<double>(a))); }
00199   static inline unsigned short int zero()  { return 0; }
00200   static inline unsigned short int one()   { return 1; }
00201   static inline unsigned short int conjugate(unsigned short int x) { return x; }
00202   static inline unsigned short int real(unsigned short int x) { return x; }
00203   static inline unsigned short int imag(unsigned short int) { return 0; }
00204   static inline bool isnaninf(unsigned short int) { return false; }
00205   static inline void seedrandom(unsigned int s) { 
00206     std::srand(s); 
00207 #ifdef __APPLE__
00208     // throw away first random number to address bug 3655
00209     // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
00210     random();
00211 #endif
00212   }
00213   //static inline int random() { return (-1 + 2*rand()); }  // RAB: This version should be used to be consistent with others
00214   static inline unsigned short int random() { return std::rand(); }             // RAB: This version should be used for an unsigned int, not int
00215   static inline std::string name() { return "unsigned short int"; }
00216   static inline unsigned short int squareroot(unsigned short int x) { return (unsigned short int) std::sqrt((double) x); }
00217   static inline unsigned short int pow(unsigned short int x, unsigned short int y) { return (unsigned short int) std::pow((double)x,(double)y); }
00218   static inline unsigned short int log(unsigned short int x) { return static_cast<unsigned short int> (std::log (static_cast<double> (x))); }
00219   static inline unsigned short int log10(unsigned short int x) { return static_cast<unsigned short int> (std::log10 (static_cast<double> (x))); }
00220 };
00221 
00222 
00223 template<>
00224 struct ScalarTraits<int>
00225 {
00226   typedef int magnitudeType;
00227   typedef int halfPrecision;
00228   typedef int doublePrecision;
00229   static const bool isComplex = false;
00230   static const bool isOrdinal = true;
00231   static const bool isComparable = true;
00232   static const bool hasMachineParameters = false;
00233   // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
00234   static inline magnitudeType magnitude(int a) { return static_cast<int>(std::fabs(static_cast<double>(a))); }
00235   static inline int zero()  { return 0; }
00236   static inline int one()   { return 1; }
00237   static inline int conjugate(int x) { return x; }
00238   static inline int real(int x) { return x; }
00239   static inline int imag(int) { return 0; }
00240   static inline bool isnaninf(int) { return false; }
00241   static inline void seedrandom(unsigned int s) { 
00242     std::srand(s); 
00243 #ifdef __APPLE__
00244     // throw away first random number to address bug 3655
00245     // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
00246     random();
00247 #endif
00248   }
00249   //static inline int random() { return (-1 + 2*rand()); }  // RAB: This version should be used to be consistent with others
00250   static inline int random() { return std::rand(); }             // RAB: This version should be used for an unsigned int, not int
00251   static inline std::string name() { return "int"; }
00252   static inline int squareroot(int x) { return (int) std::sqrt((double) x); }
00253   static inline int pow(int x, int y) { return (int) std::pow((double)x,(double)y); }
00254   static inline int log(int x) { return static_cast<int> (std::log (static_cast<double> (x))); }
00255   static inline int log10(int x) { return static_cast<int> (std::log10 (static_cast<double> (x))); }
00256 };
00257 
00258 
00259 template<>
00260 struct ScalarTraits<unsigned int>
00261 {
00262   typedef unsigned int magnitudeType;
00263   typedef unsigned int halfPrecision;
00264   typedef unsigned int doublePrecision;
00265   static const bool isComplex = false;
00266   static const bool isOrdinal = true;
00267   static const bool isComparable = true;
00268   static const bool hasMachineParameters = false;
00269   // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
00270   static inline magnitudeType magnitude(unsigned int a) { return static_cast<unsigned int>(std::fabs(static_cast<double>(a))); }
00271   static inline unsigned int zero()  { return 0; }
00272   static inline unsigned int one()   { return 1; }
00273   static inline unsigned int conjugate(unsigned int x) { return x; }
00274   static inline unsigned int real(unsigned int x) { return x; }
00275   static inline unsigned int imag(unsigned int) { return 0; }
00276   static inline void seedrandom(unsigned int s) { 
00277     std::srand(s); 
00278 #ifdef __APPLE__
00279     // throw away first random number to address bug 3655
00280     // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
00281     random();
00282 #endif
00283   }
00284   //static inline int random() { return (-1 + 2*rand()); }  // RAB: This version should be used to be consistent with others
00285   static inline unsigned int random() { return std::rand(); }             // RAB: This version should be used for an unsigned int, not int
00286   static inline std::string name() { return "unsigned int"; }
00287   static inline unsigned int squareroot(unsigned int x) { return (unsigned int) std::sqrt((double) x); }
00288   static inline unsigned int pow(unsigned int x, unsigned int y) { return (unsigned int) std::pow((double)x,(double)y); }
00289   static inline unsigned int log(unsigned int x) { return static_cast<unsigned int> (std::log (static_cast<double> (x))); }
00290   static inline unsigned int log10(unsigned int x) { return static_cast<unsigned int> (std::log10 (static_cast<double> (x))); }
00291 };
00292 
00293 
00294 template<>
00295 struct ScalarTraits<long int>
00296 {
00297   typedef long int magnitudeType;
00298   typedef long int halfPrecision;
00299   typedef long int doublePrecision;
00300   static const bool isComplex = false;
00301   static const bool isOrdinal = true;
00302   static const bool isComparable = true;
00303   static const bool hasMachineParameters = false;
00304   // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
00305   static inline magnitudeType magnitude(long int a) { return static_cast<long int>(std::fabs(static_cast<double>(a))); }
00306   static inline long int zero()  { return 0; }
00307   static inline long int one()   { return 1; }
00308   static inline long int conjugate(long int x) { return x; }
00309   static inline long int real(long int x) { return x; }
00310   static inline long int imag(long int) { return 0; }
00311   static inline void seedrandom(unsigned int s) { 
00312     std::srand(s); 
00313 #ifdef __APPLE__
00314     // throw away first random number to address bug 3655
00315     // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
00316     random();
00317 #endif
00318   }
00319   //static inline int random() { return (-1 + 2*rand()); }  // RAB: This version should be used to be consistent with others
00320   static inline long int random() { return std::rand(); }             // RAB: This version should be used for an unsigned int, not int
00321   static inline std::string name() { return "long int"; }
00322   static inline long int squareroot(long int x) { return (long int) std::sqrt((double) x); }
00323   static inline long int pow(long int x, long int y) { return (long int) std::pow((double)x,(double)y); }
00324   // Note: Depending on the number of bits in long int, the cast from
00325   // long int to double may not be exact.
00326   static inline long int log(long int x) { return static_cast<long int> (std::log (static_cast<double> (x))); }
00327   static inline long int log10(long int x) { return static_cast<long int> (std::log10 (static_cast<double> (x))); }
00328 };
00329 
00330 
00331 template<>
00332 struct ScalarTraits<long unsigned int>
00333 {
00334   typedef long unsigned int magnitudeType;
00335   typedef long unsigned int halfPrecision;
00336   typedef long unsigned int doublePrecision;
00337   static const bool isComplex = false;
00338   static const bool isOrdinal = true;
00339   static const bool isComparable = true;
00340   static const bool hasMachineParameters = false;
00341   // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
00342   static inline magnitudeType magnitude(long unsigned int a) { return static_cast<long unsigned int>(std::fabs(static_cast<double>(a))); }
00343   static inline long unsigned int zero()  { return 0; }
00344   static inline long unsigned int one()   { return 1; }
00345   static inline long unsigned int conjugate(long unsigned int x) { return x; }
00346   static inline long unsigned int real(long unsigned int x) { return x; }
00347   static inline long unsigned int imag(long unsigned int) { return 0; }
00348   static inline void seedrandom(unsigned int s) { 
00349     std::srand(s); 
00350 #ifdef __APPLE__
00351     // throw away first random number to address bug 3655
00352     // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
00353     random();
00354 #endif
00355   }
00356   //static inline int random() { return (-1 + 2*rand()); }  // RAB: This version should be used to be consistent with others
00357   static inline long unsigned int random() { return std::rand(); }             // RAB: This version should be used for an unsigned int, not int
00358   static inline std::string name() { return "long unsigned int"; }
00359   static inline long unsigned int squareroot(long unsigned int x) { return (long unsigned int) std::sqrt((double) x); }
00360   static inline long unsigned int pow(long unsigned int x, long unsigned int y) { return (long unsigned int) std::pow((double)x,(double)y); }
00361   // Note: Depending on the number of bits in long unsigned int, the
00362   // cast from long unsigned int to double may not be exact.
00363   static inline long unsigned int log(long unsigned int x) { return static_cast<long unsigned int> (std::log (static_cast<double> (x))); }
00364   static inline long unsigned int log10(long unsigned int x) { return static_cast<long unsigned int> (std::log10 (static_cast<double> (x))); }
00365 };
00366 
00367 
00368 #ifdef HAVE_TEUCHOS_LONG_LONG_INT
00369 template<>
00370 struct ScalarTraits<long long int>
00371 {
00372   typedef long long int magnitudeType;
00373   typedef long long int halfPrecision;
00374   typedef long long int doublePrecision;
00375   static const bool isComplex = false;
00376   static const bool isOrdinal = true;
00377   static const bool isComparable = true;
00378   static const bool hasMachineParameters = false;
00379   // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
00380   static inline magnitudeType magnitude(long long int a) { return static_cast<long long int>(std::fabs(static_cast<double>(a))); }
00381   static inline long long int zero()  { return 0; }
00382   static inline long long int one()   { return 1; }
00383   static inline long long int conjugate(long long int x) { return x; }
00384   static inline long long int real(long long int x) { return x; }
00385   static inline long long int imag(long long int) { return 0; }
00386   static inline void seedrandom(unsigned int s) { 
00387     std::srand(s); 
00388 #ifdef __APPLE__
00389     // throw away first random number to address bug 3655
00390     // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
00391     random();
00392 #endif
00393   }
00394   //static inline int random() { return (-1 + 2*rand()); }  // RAB: This version should be used to be consistent with others
00395   static inline long long int random() { return std::rand(); }             // RAB: This version should be used for an unsigned int, not int
00396   static inline std::string name() { return "long long int"; }
00397   static inline long long int squareroot(long long int x) { return (long long int) std::sqrt((double) x); }
00398   static inline long long int pow(long long int x, long long int y) { return (long long int) std::pow((double)x,(double)y); }
00399   // Note: Depending on the number of bits in long long int, the cast
00400   // from long long int to double may not be exact.
00401   static inline long long int log(long long int x) { return static_cast<long long int> (std::log (static_cast<double> (x))); }
00402   static inline long long int log10(long long int x) { return static_cast<long long int> (std::log10 (static_cast<double> (x))); }
00403 };
00404 
00405 template<>
00406 struct ScalarTraits<unsigned long long int>
00407 {
00408   typedef unsigned long long int magnitudeType;
00409   typedef unsigned long long int halfPrecision;
00410   typedef unsigned long long int doublePrecision;
00411   static const bool isComplex = false;
00412   static const bool isOrdinal = true;
00413   static const bool isComparable = true;
00414   static const bool hasMachineParameters = false;
00415   // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
00416   static inline magnitudeType magnitude(unsigned long long int a) { return static_cast<unsigned long long int>(std::fabs(static_cast<double>(a))); }
00417   static inline unsigned long long int zero()  { return 0; }
00418   static inline unsigned long long int one()   { return 1; }
00419   static inline unsigned long long int conjugate(unsigned long long int x) { return x; }
00420   static inline unsigned long long int real(unsigned long long int x) { return x; }
00421   static inline unsigned long long int imag(unsigned long long int) { return 0; }
00422   static inline void seedrandom(unsigned int s) { 
00423     std::srand(s); 
00424 #ifdef __APPLE__
00425     // throw away first random number to address bug 3655
00426     // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
00427     random();
00428 #endif
00429   }
00430   //static inline int random() { return (-1 + 2*rand()); }  // RAB: This version should be used to be consistent with others
00431   static inline unsigned long long int random() { return std::rand(); }             // RAB: This version should be used for an unsigned int, not int
00432   static inline std::string name() { return "unsigned long long int"; }
00433   static inline unsigned long long int squareroot(unsigned long long int x) { return (unsigned long long int) std::sqrt((double) x); }
00434   static inline unsigned long long int pow(unsigned long long int x, unsigned long long int y) { return (unsigned long long int) std::pow((double)x,(double)y); }
00435   // Note: Depending on the number of bits in unsigned long long int,
00436   // the cast from unsigned long long int to double may not be exact.
00437   static inline unsigned long long int log(unsigned long long int x) { return static_cast<unsigned long long int> (std::log (static_cast<double> (x))); }
00438   static inline unsigned long long int log10(unsigned long long int x) { return static_cast<unsigned long long int> (std::log10 (static_cast<double> (x))); }
00439 };
00440 #endif // HAVE_TEUCHOS_LONG_LONG_INT
00441 
00442 
00443 #ifndef __sun
00444 extern TEUCHOS_LIB_DLL_EXPORT const float flt_nan;
00445 #endif
00446 
00447  
00448 template<>
00449 struct ScalarTraits<float>
00450 {
00451   typedef float magnitudeType;
00452   typedef float halfPrecision; // should become IEEE754-2008 binary16 or fp16 later, perhaps specified at configure according to architectural support
00453   typedef double doublePrecision;
00454   static const bool isComplex = false;
00455   static const bool isOrdinal = false;
00456   static const bool isComparable = true;
00457   static const bool hasMachineParameters = true;
00458   static inline float eps()   {
00459     return std::numeric_limits<float>::epsilon();
00460   }
00461   static inline float sfmin() {
00462     return std::numeric_limits<float>::min();
00463   }
00464   static inline float base()  {
00465     return static_cast<float>(std::numeric_limits<float>::radix);
00466   }
00467   static inline float prec()  {
00468     return eps()*base();
00469   }
00470   static inline float t()     {
00471     return static_cast<float>(std::numeric_limits<float>::digits);
00472   }
00473   static inline float rnd()   {
00474     return ( std::numeric_limits<float>::round_style == std::round_to_nearest ? one() : zero() );
00475   }
00476   static inline float emin()  {
00477     return static_cast<float>(std::numeric_limits<float>::min_exponent);
00478   }
00479   static inline float rmin()  {
00480     return std::numeric_limits<float>::min();
00481   }
00482   static inline float emax()  {
00483     return static_cast<float>(std::numeric_limits<float>::max_exponent);
00484   }
00485   static inline float rmax()  {
00486     return std::numeric_limits<float>::max();
00487   }
00488   static inline magnitudeType magnitude(float a)
00489     {
00490 #ifdef TEUCHOS_DEBUG
00491       TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
00492         a, "Error, the input value to magnitude(...) a = " << a << " can not be NaN!" );
00493 #endif      
00494       return std::fabs(a);
00495     }    
00496   static inline float zero()  { return(0.0f); }
00497   static inline float one()   { return(1.0f); }    
00498   static inline float conjugate(float x)   { return(x); }    
00499   static inline float real(float x) { return x; }
00500   static inline float imag(float) { return zero(); }
00501   static inline float nan() {
00502 #ifdef __sun
00503     return 0.0f/std::sin(0.0f);
00504 #else
00505     return flt_nan;
00506 #endif
00507   }
00508   static inline bool isnaninf(float x) {
00509     return generic_real_isnaninf<float>(x);
00510   }
00511   static inline void seedrandom(unsigned int s) { 
00512     std::srand(s); 
00513 #ifdef __APPLE__
00514     // throw away first random number to address bug 3655
00515     // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
00516     random();
00517 #endif
00518   }
00519   static inline float random() { float rnd = (float) std::rand() / RAND_MAX; return (-1.0f + 2.0f * rnd); }
00520   static inline std::string name() { return "float"; }
00521   static inline float squareroot(float x)
00522     {
00523 #ifdef TEUCHOS_DEBUG
00524       TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
00525         x, "Error, the input value to squareroot(...) x = " << x << " can not be NaN!" );
00526 #endif
00527       errno = 0;
00528       const float rtn = std::sqrt(x);
00529       if (errno)
00530         return nan();
00531       return rtn;
00532     }
00533   static inline float pow(float x, float y) { return std::pow(x,y); }
00534   static inline float log(float x) { return std::log(x); }
00535   static inline float log10(float x) { return std::log10(x); }
00536 };
00537 
00538 
00539 #ifndef __sun
00540 extern TEUCHOS_LIB_DLL_EXPORT const double dbl_nan;
00541 #endif
00542 
00543  
00544 template<>
00545 struct ScalarTraits<double>
00546 {
00547   typedef double magnitudeType;
00548   typedef float halfPrecision;
00549   /* there are different options as to how to double "double"
00550      - QD's DD (if available)
00551      - ARPREC
00552      - GNU MP
00553      - a true hardware quad
00554 
00555      in the shortterm, this should be specified at configure time. I have inserted a configure-time option (--enable-teuchos-double-to-dd) 
00556      which uses QD's DD when available. This must be used alongside --enable-teuchos-qd.
00557    */
00558 #if defined(HAVE_TEUCHOS_DOUBLE_TO_QD)
00559   typedef dd_real doublePrecision;
00560 #elif defined(HAVE_TEUCHOS_DOUBLE_TO_ARPREC)
00561   typedef mp_real doublePrecision;
00562 #else
00563   typedef double doublePrecision;     // don't double "double" in this case
00564 #endif
00565   static const bool isComplex = false;
00566   static const bool isOrdinal = false;
00567   static const bool isComparable = true;
00568   static const bool hasMachineParameters = true;
00569   static inline double eps()   {
00570     return std::numeric_limits<double>::epsilon();
00571   }
00572   static inline double sfmin() {
00573     return std::numeric_limits<double>::min();
00574   }
00575   static inline double base()  {
00576     return std::numeric_limits<double>::radix;
00577   }
00578   static inline double prec()  {
00579     return eps()*base();
00580   }
00581   static inline double t()     {
00582     return std::numeric_limits<double>::digits;
00583   }
00584   static inline double rnd()   {
00585     return ( std::numeric_limits<double>::round_style == std::round_to_nearest ? double(1.0) : double(0.0) );
00586   }
00587   static inline double emin()  {
00588     return std::numeric_limits<double>::min_exponent;
00589   }
00590   static inline double rmin()  {
00591     return std::numeric_limits<double>::min();
00592   }
00593   static inline double emax()  {
00594     return std::numeric_limits<double>::max_exponent;
00595   }
00596   static inline double rmax()  {
00597     return std::numeric_limits<double>::max();
00598   }
00599   static inline magnitudeType magnitude(double a)
00600     {
00601 #ifdef TEUCHOS_DEBUG
00602       TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
00603         a, "Error, the input value to magnitude(...) a = " << a << " can not be NaN!" );
00604 #endif      
00605       return std::fabs(a);
00606     }
00607   static inline double zero()  { return 0.0; }
00608   static inline double one()   { return 1.0; }
00609   static inline double conjugate(double x)   { return(x); }    
00610   static inline double real(double x) { return(x); }
00611   static inline double imag(double) { return(0); }
00612   static inline double nan() {
00613 #ifdef __sun
00614     return 0.0/std::sin(0.0);
00615 #else
00616     return dbl_nan;
00617 #endif
00618   }
00619   static inline bool isnaninf(double x) {
00620     return generic_real_isnaninf<double>(x);
00621   }
00622   static inline void seedrandom(unsigned int s) { 
00623     std::srand(s); 
00624 #ifdef __APPLE__
00625     // throw away first random number to address bug 3655
00626     // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
00627     random();
00628 #endif
00629   }
00630   static inline double random() { double rnd = (double) std::rand() / RAND_MAX; return (double)(-1.0 + 2.0 * rnd); }
00631   static inline std::string name() { return "double"; }
00632   static inline double squareroot(double x)
00633     {
00634 #ifdef TEUCHOS_DEBUG
00635       TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
00636         x, "Error, the input value to squareroot(...) x = " << x << " can not be NaN!" );
00637 #endif      
00638       errno = 0;
00639       const double rtn = std::sqrt(x);
00640       if (errno)
00641         return nan();
00642       return rtn;
00643     }
00644   static inline double pow(double x, double y) { return std::pow(x,y); }
00645   static inline double log(double x) { return std::log(x); }
00646   static inline double log10(double x) { return std::log10(x); }
00647 };
00648 
00649 
00650 #ifdef HAVE_TEUCHOS_QD
00651 
00652 
00653 bool operator&&(const dd_real &a, const dd_real &b);
00654 bool operator&&(const qd_real &a, const qd_real &b);
00655 
00656 
00657 template<>
00658 struct ScalarTraits<dd_real>
00659 {
00660   typedef dd_real magnitudeType;
00661   typedef double halfPrecision;
00662   typedef qd_real doublePrecision;
00663   static const bool isComplex = false;
00664   static const bool isOrdinal = false;
00665   static const bool isComparable = true;
00666   static const bool hasMachineParameters = true;
00667   static inline dd_real eps()   { return std::numeric_limits<dd_real>::epsilon(); }
00668   static inline dd_real sfmin() { return std::numeric_limits<dd_real>::min(); }
00669   static inline dd_real base()  { return std::numeric_limits<dd_real>::radix; }
00670   static inline dd_real prec()  { return eps()*base(); }
00671   static inline dd_real t()     { return std::numeric_limits<dd_real>::digits; }
00672   static inline dd_real rnd()   { return ( std::numeric_limits<dd_real>::round_style == std::round_to_nearest ? dd_real(1.0) : dd_real(0.0) ); }
00673   static inline dd_real emin()  { return std::numeric_limits<dd_real>::min_exponent; }
00674   static inline dd_real rmin()  { return std::numeric_limits<dd_real>::min(); }
00675   static inline dd_real emax()  { return std::numeric_limits<dd_real>::max_exponent; }
00676   static inline dd_real rmax()  { return std::numeric_limits<dd_real>::max(); }
00677   static inline magnitudeType magnitude(dd_real a)
00678   {
00679 #ifdef TEUCHOS_DEBUG
00680       TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
00681         a, "Error, the input value to magnitude(...) a = " << a << " can not be NaN!" );
00682 #endif      
00683       return ::abs(a);
00684   }
00685   static inline dd_real zero()  { return dd_real(0.0); }
00686   static inline dd_real one()   { return dd_real(1.0); }
00687   static inline dd_real conjugate(dd_real x)   { return(x); }    
00688   static inline dd_real real(dd_real x) { return x ; }
00689   static inline dd_real imag(dd_real) { return zero(); }
00690   static inline dd_real nan() { return dd_real::_nan; }
00691   static inline bool isnaninf(dd_real x) { return isnan(x) || isinf(x); }
00692   static inline void seedrandom(unsigned int s) {
00693     // ddrand() uses std::rand(), so the std::srand() is our seed
00694     std::srand(s); 
00695 #ifdef __APPLE__
00696     // throw away first random number to address bug 3655
00697     // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
00698     random();
00699 #endif
00700   }
00701   static inline dd_real random() { return ddrand(); }
00702   static inline std::string name() { return "dd_real"; }
00703   static inline dd_real squareroot(dd_real x)
00704   {
00705 #ifdef TEUCHOS_DEBUG
00706       TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
00707         x, "Error, the input value to squareroot(...) x = " << x << " can not be NaN!" );
00708 #endif      
00709       return ::sqrt(x);
00710   }
00711   static inline dd_real pow(dd_real x, dd_real y) { return ::pow(x,y); }
00712   // dd_real puts its transcendental functions in the global namespace.
00713   static inline dd_real log(dd_real x) { return ::log(x); }
00714   static inline dd_real log10(dd_real x) { return ::log10(x); }
00715 };
00716 
00717 
00718 template<>
00719 struct ScalarTraits<qd_real>
00720 {
00721   typedef qd_real magnitudeType;
00722   typedef dd_real halfPrecision;
00723   typedef qd_real doublePrecision;
00724   static const bool isComplex = false;
00725   static const bool isOrdinal = false;
00726   static const bool isComparable = true;
00727   static const bool hasMachineParameters = true;
00728   static inline qd_real eps()   { return std::numeric_limits<qd_real>::epsilon(); }
00729   static inline qd_real sfmin() { return std::numeric_limits<qd_real>::min(); }
00730   static inline qd_real base()  { return std::numeric_limits<qd_real>::radix; }
00731   static inline qd_real prec()  { return eps()*base(); }
00732   static inline qd_real t()     { return std::numeric_limits<qd_real>::digits; }
00733   static inline qd_real rnd()   { return ( std::numeric_limits<qd_real>::round_style == std::round_to_nearest ? qd_real(1.0) : qd_real(0.0) ); }
00734   static inline qd_real emin()  { return std::numeric_limits<qd_real>::min_exponent; }
00735   static inline qd_real rmin()  { return std::numeric_limits<qd_real>::min(); }
00736   static inline qd_real emax()  { return std::numeric_limits<qd_real>::max_exponent; }
00737   static inline qd_real rmax()  { return std::numeric_limits<qd_real>::max(); }
00738   static inline magnitudeType magnitude(qd_real a)
00739   {
00740 #ifdef TEUCHOS_DEBUG
00741       TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
00742         a, "Error, the input value to magnitude(...) a = " << a << " can not be NaN!" );
00743 #endif      
00744       return ::abs(a);
00745   }
00746   static inline qd_real zero()  { return qd_real(0.0); }
00747   static inline qd_real one()   { return qd_real(1.0); }
00748   static inline qd_real conjugate(qd_real x)   { return(x); }    
00749   static inline qd_real real(qd_real x) { return x ; }
00750   static inline qd_real imag(qd_real) { return zero(); }
00751   static inline qd_real nan() { return qd_real::_nan; }
00752   static inline bool isnaninf(qd_real x) { return isnan(x) || isinf(x); }
00753   static inline void seedrandom(unsigned int s) {
00754     // qdrand() uses std::rand(), so the std::srand() is our seed
00755     std::srand(s); 
00756 #ifdef __APPLE__
00757     // throw away first random number to address bug 3655
00758     // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
00759     random();
00760 #endif
00761   }
00762   static inline qd_real random() { return qdrand(); }
00763   static inline std::string name() { return "qd_real"; }
00764   static inline qd_real squareroot(qd_real x)
00765   {
00766 #ifdef TEUCHOS_DEBUG
00767       TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
00768         x, "Error, the input value to squareroot(...) x = " << x << " can not be NaN!" );
00769 #endif      
00770       return ::sqrt(x);
00771   }
00772   static inline qd_real pow(qd_real x, qd_real y) { return ::pow(x,y); }
00773   // qd_real puts its transcendental functions in the global namespace.
00774   static inline qd_real log(qd_real x) { return ::log(x); }
00775   static inline qd_real log10(qd_real x) { return ::log10(x); }
00776 };
00777 
00778 
00779 #endif  // HAVE_TEUCHOS_QD
00780 
00781 
00782 #ifdef HAVE_TEUCHOS_GNU_MP
00783 
00784 
00785 extern gmp_randclass gmp_rng; 
00786 
00787 
00788 /* Regarding halfPrecision, doublePrecision and mpf_class: 
00789    Because the precision of an mpf_class float is not determined by the data type, 
00790    there is no way to fill the typedefs for this object. 
00791 
00792    Instead, we could create new data classes (e.g., Teuchos::MPF128, Teuchos::MPF256) for 
00793    commonly used levels of precision, and fill out ScalarTraits for these. This would allow us
00794    to typedef the promotions and demotions in the appropriate way. These classes would serve to 
00795    wrap an mpf_class object, calling the constructor for the appropriate precision, exposing the 
00796    arithmetic routines but hiding the precision-altering routines.
00797    
00798    Alternatively (perhaps, preferably), would could create a single class templated on the precision (e.g., Teuchos::MPF<N>). 
00799    Then we have a single (partially-specialized) implementation of ScalarTraits. This class, as above, must expose all of the 
00800    operations expected of a scalar type; however, most of these can be trivially stolen from the gmpcxx.h header file
00801 
00802    CGB/RAB, 01/05/2009
00803 */
00804 template<>
00805 struct ScalarTraits<mpf_class>
00806 {
00807   typedef mpf_class magnitudeType;
00808   typedef mpf_class halfPrecision;
00809   typedef mpf_class doublePrecision;
00810   static const bool isComplex = false;
00811   static const bool hasMachineParameters = false;
00812   // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
00813   static magnitudeType magnitude(mpf_class a) { return std::abs(a); }
00814   static inline mpf_class zero() { mpf_class zero = 0.0; return zero; }
00815   static inline mpf_class one() { mpf_class one = 1.0; return one; }    
00816   static inline mpf_class conjugate(mpf_class x) { return x; }
00817   static inline mpf_class real(mpf_class x) { return(x); }
00818   static inline mpf_class imag(mpf_class x) { return(0); }
00819   static inline bool isnaninf(mpf_class x) { return false; } // mpf_class currently can't handle nan or inf!
00820   static inline void seedrandom(unsigned int s) { 
00821     unsigned long int seedVal = static_cast<unsigned long int>(s);
00822     gmp_rng.seed( seedVal );  
00823   }
00824   static inline mpf_class random() { 
00825     return gmp_rng.get_f(); 
00826   }
00827   static inline std::string name() { return "mpf_class"; }
00828   static inline mpf_class squareroot(mpf_class x) { return std::sqrt(x); }
00829   static inline mpf_class pow(mpf_class x, mpf_class y) { return pow(x,y); }
00830   // Todo: RAB: 2004/05/28: Add nan() and isnaninf() functions when needed!
00831 };
00832 
00833 #endif  // HAVE_TEUCHOS_GNU_MP
00834 
00835 #ifdef HAVE_TEUCHOS_ARPREC
00836 
00837 /* See discussion above for mpf_class, regarding halfPrecision and doublePrecision. Something similar will need to be done
00838    for ARPREC. */
00839 template<>
00840 struct ScalarTraits<mp_real>
00841 {
00842   typedef mp_real magnitudeType;
00843   typedef double halfPrecision;
00844   typedef mp_real doublePrecision;
00845   static const bool isComplex = false;
00846   static const bool isComparable = true;
00847   static const bool isOrdinal = false;
00848   static const bool hasMachineParameters = false;
00849   // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
00850   static magnitudeType magnitude(mp_real a) { return abs(a); }
00851   static inline mp_real zero() { mp_real zero = 0.0; return zero; }
00852   static inline mp_real one() { mp_real one = 1.0; return one; }    
00853   static inline mp_real conjugate(mp_real x) { return x; }
00854   static inline mp_real real(mp_real x) { return(x); }
00855   static inline mp_real imag(mp_real x) { return zero(); }
00856   static inline bool isnaninf(mp_real x) { return false; } // ToDo: Change this?
00857   static inline void seedrandom(unsigned int s) { 
00858     long int seedVal = static_cast<long int>(s);
00859     srand48(seedVal);
00860   }
00861   static inline mp_real random() { return mp_rand(); }
00862   static inline std::string name() { return "mp_real"; }
00863   static inline mp_real squareroot(mp_real x) { return sqrt(x); }
00864   static inline mp_real pow(mp_real x, mp_real y) { return pow(x,y); }
00865   // Todo: RAB: 2004/05/28: Add nan() and isnaninf() functions when needed!
00866 };
00867 
00868   
00869 #endif // HAVE_TEUCHOS_ARPREC
00870 
00871  
00872 #ifdef HAVE_TEUCHOS_COMPLEX
00873 
00874 
00875 // Partial specialization for std::complex numbers templated on real type T
00876 template<class T> 
00877 struct ScalarTraits<
00878   std::complex<T>
00879 >
00880 {
00881   typedef std::complex<T>  ComplexT;
00882   typedef std::complex<typename ScalarTraits<T>::halfPrecision> halfPrecision;
00883   typedef std::complex<typename ScalarTraits<T>::doublePrecision> doublePrecision;
00884   typedef typename ScalarTraits<T>::magnitudeType magnitudeType;
00885   static const bool isComplex = true;
00886   static const bool isOrdinal = ScalarTraits<T>::isOrdinal;
00887   static const bool isComparable = false;
00888   static const bool hasMachineParameters = true;
00889   static inline magnitudeType eps()          { return ScalarTraits<magnitudeType>::eps(); }
00890   static inline magnitudeType sfmin()        { return ScalarTraits<magnitudeType>::sfmin(); }
00891   static inline magnitudeType base()         { return ScalarTraits<magnitudeType>::base(); }
00892   static inline magnitudeType prec()         { return ScalarTraits<magnitudeType>::prec(); }
00893   static inline magnitudeType t()            { return ScalarTraits<magnitudeType>::t(); }
00894   static inline magnitudeType rnd()          { return ScalarTraits<magnitudeType>::rnd(); }
00895   static inline magnitudeType emin()         { return ScalarTraits<magnitudeType>::emin(); }
00896   static inline magnitudeType rmin()         { return ScalarTraits<magnitudeType>::rmin(); }
00897   static inline magnitudeType emax()         { return ScalarTraits<magnitudeType>::emax(); }
00898   static inline magnitudeType rmax()         { return ScalarTraits<magnitudeType>::rmax(); }
00899   static magnitudeType magnitude(ComplexT a)
00900     {
00901 #ifdef TEUCHOS_DEBUG
00902       TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
00903         a, "Error, the input value to magnitude(...) a = " << a << " can not be NaN!" );
00904 #endif      
00905       return std::abs(a);
00906     }
00907   static inline ComplexT zero()              { return ComplexT(ScalarTraits<magnitudeType>::zero(),ScalarTraits<magnitudeType>::zero()); }
00908   static inline ComplexT one()               { return ComplexT(ScalarTraits<magnitudeType>::one(),ScalarTraits<magnitudeType>::zero()); }
00909   static inline ComplexT conjugate(ComplexT a){ return ComplexT(a.real(),-a.imag()); }
00910   static inline magnitudeType real(ComplexT a) { return a.real(); }
00911   static inline magnitudeType imag(ComplexT a) { return a.imag(); }
00912   static inline ComplexT nan()               { return ComplexT(ScalarTraits<magnitudeType>::nan(),ScalarTraits<magnitudeType>::nan()); }
00913   static inline bool isnaninf(ComplexT x)    { return ScalarTraits<magnitudeType>::isnaninf(x.real()) || ScalarTraits<magnitudeType>::isnaninf(x.imag()); }
00914   static inline void seedrandom(unsigned int s) { ScalarTraits<magnitudeType>::seedrandom(s); }
00915   static inline ComplexT random()
00916     {
00917       const T rnd1 = ScalarTraits<magnitudeType>::random();
00918       const T rnd2 = ScalarTraits<magnitudeType>::random();
00919       return ComplexT(rnd1,rnd2);
00920     }
00921   static inline std::string name() { return std::string("std::complex<")+std::string(ScalarTraits<magnitudeType>::name())+std::string(">"); }
00922   // This will only return one of the square roots of x, the other can be obtained by taking its conjugate
00923   static inline ComplexT squareroot(ComplexT x)
00924     {
00925 #ifdef TEUCHOS_DEBUG
00926       TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
00927         x, "Error, the input value to squareroot(...) x = " << x << " can not be NaN!" );
00928 #endif
00929       typedef ScalarTraits<magnitudeType>  STMT;
00930       const T r  = x.real(), i = x.imag(), zero = STMT::zero(), two = 2.0;
00931       const T a  = STMT::squareroot((r*r)+(i*i));
00932       const T nr = STMT::squareroot((a+r)/two);
00933       const T ni = ( i == zero ? zero : STMT::squareroot((a-r)/two) );
00934       return ComplexT(nr,ni);
00935     }
00936     // 2010/03/19: rabartl: Above, I had to add the check for i == zero
00937     // to avoid a returned NaN on the Intel 10.1 compiler.  For some
00938     // reason, having these two squareroot calls in a row produce a NaN
00939     // in an optimized build with this compiler.  Amazingly, when I put
00940     // in print statements (i.e. std::cout << ...) the NaN went away and
00941     // the second squareroot((a-r)/two) returned zero correctly.  I
00942     // guess that calling the output routine flushed the registers or
00943     // something and restarted the squareroot rountine for this compiler
00944     // or something similar.  Actually, due to roundoff, it is possible that a-r
00945     // might be slightly less than zero (i.e. -1e-16) and that would cause
00946     // a possbile NaN return.  THe above if test is the right thing to do
00947     // I think and is very cheap.
00948   static inline ComplexT pow(ComplexT x, ComplexT y) { return pow(x,y); }
00949 };
00950 
00951 
00952 #endif //  HAVE_TEUCHOS_COMPLEX
00953 
00954 
00955 #endif // DOXYGEN_SHOULD_SKIP_THIS
00956 
00957 
00958 } // Teuchos namespace
00959 
00960 
00961 #endif // _TEUCHOS_SCALARTRAITS_HPP_
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines