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 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 // Kris
00030 // 06.18.03 -- Minor formatting changes
00031 //          -- Changed calls to LAPACK objects to use new <OType, SType> templates
00032 // 07.08.03 -- Move into Teuchos package/namespace
00033 // 07.11.03 -- Added ScalarTraits for ARPREC::mp_real
00034 // 07.14.03 -- Fixed int rand() function (was set up to return a floating-point style random number)
00035 // 07.17.03 -- Added squareroot() function
00036 
00037 #ifndef _TEUCHOS_SCALARTRAITS_HPP_
00038 #define _TEUCHOS_SCALARTRAITS_HPP_
00039 
00044 #include "Teuchos_LAPACK.hpp"
00045 
00046 #ifdef HAVE_TEUCHOS_ARPREC
00047 #include "mp/mpreal.h"
00048 #endif
00049 
00050 #ifdef HAVE_TEUCHOS_GNU_MP
00051 #include "gmp.h"
00052 #include "gmpxx.h"
00053 #endif
00054 
00082 namespace Teuchos {
00083   template<class T>
00084   struct UndefinedScalarTraits
00085   {
00087     static inline T notDefined() { return T::this_type_is_missing_a_specialization(); };
00088   };
00089 
00090   template<class T>
00091   struct ScalarTraits
00092   {
00094     typedef T magnitudeType;
00096     static inline bool haveMachineParameters() { return false; };
00098     static inline magnitudeType eps()   { return UndefinedScalarTraits<T>::notDefined(); };
00100     static inline magnitudeType sfmin() { return UndefinedScalarTraits<T>::notDefined(); };
00102     static inline magnitudeType base()  { return UndefinedScalarTraits<T>::notDefined(); };
00104     static inline magnitudeType prec()  { return UndefinedScalarTraits<T>::notDefined(); };
00106     static inline magnitudeType t()     { return UndefinedScalarTraits<T>::notDefined(); };
00108     static inline magnitudeType rnd()   { return UndefinedScalarTraits<T>::notDefined(); };
00110     static inline magnitudeType emin()  { return UndefinedScalarTraits<T>::notDefined(); };
00112     static inline magnitudeType rmin()  { return UndefinedScalarTraits<T>::notDefined(); };
00114     static inline magnitudeType emax()  { return UndefinedScalarTraits<T>::notDefined(); };
00116     static inline magnitudeType rmax()  { return UndefinedScalarTraits<T>::notDefined(); };
00118     static inline magnitudeType magnitude(T a) { return UndefinedScalarTraits<T>::notDefined(); };
00120     static inline T zero()                     { return UndefinedScalarTraits<T>::notDefined(); };
00122     static inline T one()                      { return UndefinedScalarTraits<T>::notDefined(); };
00124     static inline T conjugate(T a) { return UndefinedScalarTraits<T>::notDefined(); };
00126     static inline T nan()                      { return UndefinedScalarTraits<T>::notDefined(); };
00128     static inline bool isnaninf(const T& x)     { return UndefinedScalarTraits<T>::notDefined(); };
00130     static inline void seedrandom(unsigned int s) { int i; T t = &i; };
00132     static inline T random()                   { return UndefinedScalarTraits<T>::notDefined(); };
00134     static inline std::string name()           { (void)UndefinedScalarTraits<T>::notDefined(); return 0; };
00136     static inline magnitudeType squareroot(T x) { return UndefinedScalarTraits<T>::notDefined(); };
00137   };
00138   
00139 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00140 
00141   template<>
00142   struct ScalarTraits<int>
00143   {
00144     typedef int magnitudeType;
00145     static inline bool haveMachineParameters() { return false; };
00146     // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
00147     static inline magnitudeType magnitude(int a) { return ::abs(a); };
00148     static inline int zero()  { return 0; };
00149     static inline int one()   { return 1; };
00150     static inline int conjugate(int x) { return x; };
00151     static inline void seedrandom(unsigned int s) { srand(s); };
00152     //static inline int random() { return (-1 + 2*rand()); };  // RAB: This version should be used to be consistent with others
00153     static inline int random() { return rand(); };             // RAB: This version should be used for an unsigned int, not int
00154     static inline std::string name() { return "int"; };
00155     static inline int squareroot(int x) { return (int) ::sqrt((double) x); };
00156   };
00157 
00158   extern float flt_nan;
00159  
00160   template<>
00161   struct ScalarTraits<float>
00162   {
00163     typedef float magnitudeType;
00164     static inline bool haveMachineParameters() { return true; };
00165     static inline float eps()   {
00166 #ifdef HAVE_NUMERIC_LIMITS
00167       return std::numeric_limits<float>::epsilon();
00168 #else
00169       LAPACK<int, float> lp; return lp.LAMCH('E');
00170 #endif
00171     }
00172     static inline float sfmin() {
00173 #ifdef HAVE_NUMERIC_LIMITS
00174       return std::numeric_limits<float>::min();
00175 #else
00176       LAPACK<int, float> lp; return lp.LAMCH('S');
00177 #endif
00178     };
00179     static inline float base()  {
00180 #ifdef HAVE_NUMERIC_LIMITS
00181       return std::numeric_limits<float>::radix;
00182 #else
00183       LAPACK<int, float> lp; return lp.LAMCH('B');
00184 #endif
00185     };
00186     static inline float prec()  {
00187 #ifdef HAVE_NUMERIC_LIMITS
00188       return eps()*base();
00189 #else
00190       LAPACK<int, float> lp; return lp.LAMCH('P');
00191 #endif
00192     };
00193     static inline float t()     {
00194 #ifdef HAVE_NUMERIC_LIMITS
00195       return std::numeric_limits<float>::digits;
00196 #else
00197       LAPACK<int, float> lp; return lp.LAMCH('N');
00198 #endif
00199     };
00200     static inline float rnd()   {
00201 #ifdef HAVE_NUMERIC_LIMITS
00202       return ( std::numeric_limits<float>::round_style == std::round_to_nearest ? float(1.0) : float(0.0) );
00203 #else
00204       LAPACK<int, float> lp; return lp.LAMCH('R');
00205 #endif
00206     };
00207     static inline float emin()  {
00208 #ifdef HAVE_NUMERIC_LIMITS
00209       return std::numeric_limits<float>::min_exponent;
00210 #else
00211       LAPACK<int, float> lp; return lp.LAMCH('M');
00212 #endif
00213     };
00214     static inline float rmin()  {
00215 #ifdef HAVE_NUMERIC_LIMITS
00216       return std::numeric_limits<float>::min();
00217 #else
00218       LAPACK<int, float> lp; return lp.LAMCH('U');
00219 #endif
00220     };
00221     static inline float emax()  {
00222 #ifdef HAVE_NUMERIC_LIMITS
00223       return std::numeric_limits<float>::max_exponent;
00224 #else
00225       LAPACK<int, float> lp; return lp.LAMCH('L');
00226 #endif
00227     };
00228     static inline float rmax()  {
00229 #ifdef HAVE_NUMERIC_LIMITS
00230       return std::numeric_limits<float>::max();
00231 #else
00232       LAPACK<int, float> lp; return lp.LAMCH('O');
00233 #endif
00234     };
00235     static inline magnitudeType magnitude(float a) { return fabs(a); };    
00236     static inline float zero()  { return(0.0); };
00237     static inline float one()   { return(1.0); };    
00238     static inline float conjugate(float x)   { return(x); };    
00239     static inline float nan()   { return flt_nan; };
00240     static inline bool isnaninf(float x) { // RAB: 2004/05/28: Taken from NOX_StatusTest_FiniteValue.C
00241       const float tol = 1e-6; // Any (bounded) number should do!
00242       if( !(x <= tol) && !(x > tol) ) return true;                 // IEEE says this should fail for NaN
00243       float z=0.0*x; if( !(z <= tol) && !(z > tol) ) return true;  // Use fact that Inf*0 = NaN
00244       return false;
00245     }
00246     static inline void seedrandom(unsigned int s) { srand(s); };
00247     static inline float random() { float rnd = (float) rand() / RAND_MAX; return (float)(-1.0 + 2.0 * rnd); };
00248     static inline std::string name() { return "float"; };
00249     static inline float squareroot(float x) { return sqrt(x); };
00250   };
00251 
00252   extern double dbl_nan;
00253  
00254   template<>
00255   struct ScalarTraits<double>
00256   {
00257     typedef double magnitudeType;
00258     static inline bool haveMachineParameters() { return true; };
00259     static inline double eps()   {
00260 #ifdef HAVE_NUMERIC_LIMITS
00261       return std::numeric_limits<double>::epsilon();
00262 #else
00263       LAPACK<int, double> lp; return lp.LAMCH('E');
00264 #endif
00265     }
00266     static inline double sfmin() {
00267 #ifdef HAVE_NUMERIC_LIMITS
00268       return std::numeric_limits<double>::min();
00269 #else
00270       LAPACK<int, double> lp; return lp.LAMCH('S');
00271 #endif
00272     };
00273     static inline double base()  {
00274 #ifdef HAVE_NUMERIC_LIMITS
00275       return std::numeric_limits<double>::radix;
00276 #else
00277       LAPACK<int, double> lp; return lp.LAMCH('B');
00278 #endif
00279     };
00280     static inline double prec()  {
00281 #ifdef HAVE_NUMERIC_LIMITS
00282       return eps()*base();
00283 #else
00284       LAPACK<int, double> lp; return lp.LAMCH('P');
00285 #endif
00286     };
00287     static inline double t()     {
00288 #ifdef HAVE_NUMERIC_LIMITS
00289       return std::numeric_limits<double>::digits;
00290 #else
00291       LAPACK<int, double> lp; return lp.LAMCH('N');
00292 #endif
00293     };
00294     static inline double rnd()   {
00295 #ifdef HAVE_NUMERIC_LIMITS
00296       return ( std::numeric_limits<double>::round_style == std::round_to_nearest ? double(1.0) : double(0.0) );
00297 #else
00298       LAPACK<int, double> lp; return lp.LAMCH('R');
00299 #endif
00300     };
00301     static inline double emin()  {
00302 #ifdef HAVE_NUMERIC_LIMITS
00303       return std::numeric_limits<double>::min_exponent;
00304 #else
00305       LAPACK<int, double> lp; return lp.LAMCH('M');
00306 #endif
00307     };
00308     static inline double rmin()  {
00309 #ifdef HAVE_NUMERIC_LIMITS
00310       return std::numeric_limits<double>::min();
00311 #else
00312       LAPACK<int, double> lp; return lp.LAMCH('U');
00313 #endif
00314     };
00315     static inline double emax()  {
00316 #ifdef HAVE_NUMERIC_LIMITS
00317       return std::numeric_limits<double>::max_exponent;
00318 #else
00319       LAPACK<int, double> lp; return lp.LAMCH('L');
00320 #endif
00321     };
00322     static inline double rmax()  {
00323 #ifdef HAVE_NUMERIC_LIMITS
00324       return std::numeric_limits<double>::max();
00325 #else
00326       LAPACK<int, double> lp; return lp.LAMCH('O');
00327 #endif
00328     };
00329     static inline magnitudeType magnitude(double a) { return fabs(a); };
00330     static inline double zero()  { return 0.0; };
00331     static inline double one()   { return 1.0; };
00332     static inline double conjugate(double x)   { return(x); };    
00333     static inline double nan() { return dbl_nan; };
00334     static inline bool isnaninf(double x) { // RAB: 2004/05/28: Taken from NOX_StatusTest_FiniteValue.C
00335       const double tol = 1e-6; // Any (bounded) number should do!
00336       if( !(x <= tol) && !(x > tol) ) return true;                  // IEEE says this should fail for NaN
00337       double z=0.0*x; if( !(z <= tol) && !(z > tol) ) return true;  // Use fact that Inf*0 = NaN
00338       return false;
00339     }
00340     static inline void seedrandom(unsigned int s) { srand(s); };
00341     static inline double random() { double rnd = (double) rand() / RAND_MAX; return (double)(-1.0 + 2.0 * rnd); };
00342     static inline std::string name() { return "double"; };
00343     static inline double squareroot(double x) { return ::sqrt(x); };
00344   };
00345 
00346 #ifdef HAVE_TEUCHOS_GNU_MP
00347 
00348   extern gmp_randclass gmp_rng; 
00349 
00350   template<>
00351   struct ScalarTraits<mpf_class>
00352   {
00353     typedef mpf_class magnitudeType;
00354     static inline bool haveMachineParameters() { return false; };
00355     // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
00356     static magnitudeType magnitude(mpf_class a) { return abs(a); };
00357     static inline mpf_class zero() { mpf_class zero = 0.0; return zero; };
00358     static inline mpf_class one() { mpf_class one = 1.0; return one; };    
00359     static inline mpf_class conjugate(mpf_class x) { return x; };
00360     static inline void seedrandom(unsigned int s) { 
00361       unsigned long int seedVal = static_cast<unsigned long int>(s);
00362       gmp_rng.seed( seedVal );  
00363     };
00364     static inline mpf_class random() { 
00365       return gmp_rng.get_f(); 
00366     };
00367     static inline std::string name() { return "mpf_class"; };
00368     static inline mpf_class squareroot(mpf_class x) { return sqrt(x); };
00369     // Todo: RAB: 2004/05/28: Add nan() and isnaninf() functions when needed!
00370   };
00371 
00372 #endif  
00373 
00374 #ifdef HAVE_TEUCHOS_ARPREC
00375 
00376   template<>
00377   struct ScalarTraits<mp_real>
00378   {
00379     typedef mp_real magnitudeType;
00380     static inline bool haveMachineParameters() { return false; };
00381     // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
00382     static magnitudeType magnitude(mp_real a) { return abs(a); };
00383     static inline mp_real zero() { mp_real zero = 0.0; return zero; };
00384     static inline mp_real one() { mp_real one = 1.0; return one; };    
00385     static inline mp_real conjugate(mp_real x) { return x; };
00386     static inline void seedrandom(unsigned int s) { 
00387       long int seedVal = static_cast<long int>(s);
00388       srand48(seedVal);
00389     };
00390     static inline mp_real random() { return mp_rand(); };
00391     static inline std::string name() { return "mp_real"; };
00392     static inline mp_real squareroot(mp_real x) { return sqrt(x); };
00393     // Todo: RAB: 2004/05/28: Add nan() and isnaninf() functions when needed!
00394   };
00395   
00396 #endif // HAVE_TEUCHOS_ARPREC
00397  
00398 #if ( defined(HAVE_COMPLEX) || defined(HAVE_COMPLEX_H) ) && defined(HAVE_TEUCHOS_COMPLEX)
00399 
00400   // Partial specialization for complex numbers templated on real type T
00401   template<class T> 
00402   struct ScalarTraits<
00403 #if defined(HAVE_COMPLEX)
00404     std::complex<T>
00405 #elif  defined(HAVE_COMPLEX_H)
00406     ::complex<T>
00407 #endif
00408     >
00409   {
00410 #if defined(HAVE_COMPLEX)
00411     typedef std::complex<T>  ComplexT;
00412 #elif  defined(HAVE_COMPLEX_H)
00413     typedef ::complex<T>     ComplexT;
00414 #endif
00415     typedef typename ScalarTraits<T>::magnitudeType magnitudeType;
00416     static inline bool haveMachineParameters() { return true; };
00417     static inline magnitudeType eps()          { return ScalarTraits<magnitudeType>::eps(); };
00418     static inline magnitudeType sfmin()        { return ScalarTraits<magnitudeType>::sfmin(); };
00419     static inline magnitudeType base()         { return ScalarTraits<magnitudeType>::base(); };
00420     static inline magnitudeType prec()         { return ScalarTraits<magnitudeType>::prec(); };
00421     static inline magnitudeType t()            { return ScalarTraits<magnitudeType>::t(); };
00422     static inline magnitudeType rnd()          { return ScalarTraits<magnitudeType>::rnd(); };
00423     static inline magnitudeType emin()         { return ScalarTraits<magnitudeType>::emin(); };
00424     static inline magnitudeType rmin()         { return ScalarTraits<magnitudeType>::rmin(); };
00425     static inline magnitudeType emax()         { return ScalarTraits<magnitudeType>::emax(); };
00426     static inline magnitudeType rmax()         { return ScalarTraits<magnitudeType>::rmax(); };
00427     static magnitudeType magnitude(ComplexT a) { return std::abs(a); };
00428     static inline ComplexT zero()              { return ComplexT(ScalarTraits<magnitudeType>::zero(),ScalarTraits<magnitudeType>::zero()); };
00429     static inline ComplexT one()               { return ComplexT(ScalarTraits<magnitudeType>::one(),ScalarTraits<magnitudeType>::zero()); };
00430     static inline ComplexT conjugate(ComplexT a){ return ComplexT(a.real(),-a.imag()); };
00431     static inline ComplexT nan()               { return ComplexT(ScalarTraits<magnitudeType>::nan(),ScalarTraits<magnitudeType>::nan()); };
00432     static inline bool isnaninf(ComplexT x)    { return ScalarTraits<magnitudeType>::isnaninf(x.real()) || ScalarTraits<magnitudeType>::isnaninf(x.imag()); };
00433     static inline void seedrandom(unsigned int s) { ScalarTraits<magnitudeType>::seedrandom(s); };
00434     static inline ComplexT random()
00435     {
00436       const T rnd1 = ScalarTraits<magnitudeType>::random();
00437       const T rnd2 = ScalarTraits<magnitudeType>::random();
00438       return ComplexT(rnd1,rnd2);
00439     };
00440     static inline std::string name() { return std::string("std::complex<")+std::string(ScalarTraits<magnitudeType>::name())+std::string(">"); };
00441     // This will only return one of the square roots of x, the other can be obtained by taking its conjugate
00442     static inline ComplexT squareroot(ComplexT x)
00443     {
00444       typedef ScalarTraits<magnitudeType>  STMT;
00445       const T r  = x.real(), i = x.imag();
00446       const T a  = STMT::squareroot((r*r)+(i*i));
00447       const T nr = STMT::squareroot((a+r)/2);
00448       const T ni = STMT::squareroot((a-r)/2);
00449       return ComplexT(nr,ni);
00450     };
00451   };
00452 
00453 #endif //  HAVE_COMPLEX || HAVE_COMPLEX_H
00454 
00455 #endif // DOXYGEN_SHOULD_SKIP_THIS
00456 
00457 } // Teuchos namespace
00458 
00459 #endif // _TEUCHOS_SCALARTRAITS_HPP_

Generated on Thu Sep 18 12:41:17 2008 for Teuchos - Trilinos Tools Package by doxygen 1.3.9.1