Teuchos - Trilinos Tools Package: Teuchos_ScalarTraits.hpp Source 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_ConfigDefs.hpp"
00045 
00046 #ifndef HAVE_NUMERIC_LIMITS
00047 #include "Teuchos_LAPACK.hpp"
00048 #endif
00049 
00050 #ifdef HAVE_TEUCHOS_ARPREC
00051 #include "mp/mpreal.h"
00052 #endif
00053 
00054 #ifdef HAVE_TEUCHOS_GNU_MP
00055 #include "gmp.h"
00056 #include "gmpxx.h"
00057 #endif
00058 
00081 /* This is the default structure used by ScalarTraits<T> to produce a compile time
00082   error when the specialization does not exist for type <tt>T</tt>.
00083 */
00084 namespace Teuchos {
00085   template<class T>
00086   struct UndefinedScalarTraits
00087   {
00089     static inline T notDefined() { return T::this_type_is_missing_a_specialization(); };
00090   };
00091 
00092   template<class T>
00093   struct ScalarTraits
00094   {
00096     typedef T magnitudeType;
00098     static const bool isComplex = false;
00100     static const bool isComparable = false;
00102     static const bool hasMachineParameters = false;
00104     static inline magnitudeType eps()   { return UndefinedScalarTraits<T>::notDefined(); };
00106     static inline magnitudeType sfmin() { return UndefinedScalarTraits<T>::notDefined(); };
00108     static inline magnitudeType base()  { return UndefinedScalarTraits<T>::notDefined(); };
00110     static inline magnitudeType prec()  { return UndefinedScalarTraits<T>::notDefined(); };
00112     static inline magnitudeType t()     { return UndefinedScalarTraits<T>::notDefined(); };
00114     static inline magnitudeType rnd()   { return UndefinedScalarTraits<T>::notDefined(); };
00116     static inline magnitudeType emin()  { return UndefinedScalarTraits<T>::notDefined(); };
00118     static inline magnitudeType rmin()  { return UndefinedScalarTraits<T>::notDefined(); };
00120     static inline magnitudeType emax()  { return UndefinedScalarTraits<T>::notDefined(); };
00122     static inline magnitudeType rmax()  { return UndefinedScalarTraits<T>::notDefined(); };
00124     static inline magnitudeType magnitude(T a) { return UndefinedScalarTraits<T>::notDefined(); };
00126     static inline T zero()                     { return UndefinedScalarTraits<T>::notDefined(); };
00128     static inline T one()                      { return UndefinedScalarTraits<T>::notDefined(); };
00130     static inline magnitudeType real(T a) { return UndefinedScalarTraits<T>::notDefined(); };
00132     static inline magnitudeType imag(T a) { return UndefinedScalarTraits<T>::notDefined(); };
00134     static inline T conjugate(T a) { return UndefinedScalarTraits<T>::notDefined(); };
00136     static inline T nan()                      { return UndefinedScalarTraits<T>::notDefined(); };
00138     static inline bool isnaninf(const T& x)     { return UndefinedScalarTraits<T>::notDefined(); };
00140     static inline void seedrandom(unsigned int s) { int i; T t = &i; };
00142     static inline T random()                   { return UndefinedScalarTraits<T>::notDefined(); };
00144     static inline std::string name()           { (void)UndefinedScalarTraits<T>::notDefined(); return 0; };
00146     static inline T squareroot(T x) { return UndefinedScalarTraits<T>::notDefined(); };
00147   };
00148   
00149 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00150 
00151   template<>
00152   struct ScalarTraits<char>
00153   {
00154     typedef char magnitudeType;
00155     static const bool isComplex = false;
00156     static const bool isComparable = true;
00157     static const bool hasMachineParameters = false;
00158     // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
00159     static inline magnitudeType magnitude(char a) { return static_cast<char>(::fabs(static_cast<double>(a))); };
00160     static inline char zero()  { return 0; };
00161     static inline char one()   { return 1; };
00162     static inline char conjugate(char x) { return x; };
00163     static inline char real(char x) { return x; };
00164     static inline char imag(char x) { return 0; };
00165     static inline void seedrandom(unsigned int s) { srand(s); };
00166     //static inline char random() { return (-1 + 2*rand()); };  // RAB: This version should be used to be consistent with others
00167     static inline char random() { return rand(); };             // RAB: This version should be used for an unsigned char, not char
00168     static inline std::string name() { return "char"; };
00169     static inline char squareroot(char x) { return (char) ::sqrt((double) x); };
00170   };
00171 
00172   template<>
00173   struct ScalarTraits<int>
00174   {
00175     typedef int magnitudeType;
00176     static const bool isComplex = false;
00177     static const bool isComparable = true;
00178     static const bool hasMachineParameters = false;
00179     // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
00180     static inline magnitudeType magnitude(int a) { return static_cast<int>(::fabs(static_cast<double>(a))); };
00181     static inline int zero()  { return 0; };
00182     static inline int one()   { return 1; };
00183     static inline int conjugate(int x) { return x; };
00184     static inline int real(int x) { return x; };
00185     static inline int imag(int x) { return 0; };
00186     static inline void seedrandom(unsigned int s) { srand(s); };
00187     //static inline int random() { return (-1 + 2*rand()); };  // RAB: This version should be used to be consistent with others
00188     static inline int random() { return rand(); };             // RAB: This version should be used for an unsigned int, not int
00189     static inline std::string name() { return "int"; };
00190     static inline int squareroot(int x) { return (int) ::sqrt((double) x); };
00191   };
00192 
00193 #ifndef __sun
00194   extern const float flt_nan;
00195 #endif
00196  
00197   template<>
00198   struct ScalarTraits<float>
00199   {
00200     typedef float magnitudeType;
00201     static const bool isComplex = false;
00202     static const bool isComparable = true;
00203     static const bool hasMachineParameters = true;
00204     static inline float eps()   {
00205 #ifdef HAVE_NUMERIC_LIMITS
00206       return std::numeric_limits<float>::epsilon();
00207 #else
00208       LAPACK<int, float> lp; return lp.LAMCH('E');
00209 #endif
00210     }
00211     static inline float sfmin() {
00212 #ifdef HAVE_NUMERIC_LIMITS
00213       return std::numeric_limits<float>::min();
00214 #else
00215       LAPACK<int, float> lp; return lp.LAMCH('S');
00216 #endif
00217     };
00218     static inline float base()  {
00219 #ifdef HAVE_NUMERIC_LIMITS
00220       return std::numeric_limits<float>::radix;
00221 #else
00222       LAPACK<int, float> lp; return lp.LAMCH('B');
00223 #endif
00224     };
00225     static inline float prec()  {
00226 #ifdef HAVE_NUMERIC_LIMITS
00227       return eps()*base();
00228 #else
00229       LAPACK<int, float> lp; return lp.LAMCH('P');
00230 #endif
00231     };
00232     static inline float t()     {
00233 #ifdef HAVE_NUMERIC_LIMITS
00234       return std::numeric_limits<float>::digits;
00235 #else
00236       LAPACK<int, float> lp; return lp.LAMCH('N');
00237 #endif
00238     };
00239     static inline float rnd()   {
00240 #ifdef HAVE_NUMERIC_LIMITS
00241       return ( std::numeric_limits<float>::round_style == std::round_to_nearest ? float(1.0) : float(0.0) );
00242 #else
00243       LAPACK<int, float> lp; return lp.LAMCH('R');
00244 #endif
00245     };
00246     static inline float emin()  {
00247 #ifdef HAVE_NUMERIC_LIMITS
00248       return std::numeric_limits<float>::min_exponent;
00249 #else
00250       LAPACK<int, float> lp; return lp.LAMCH('M');
00251 #endif
00252     };
00253     static inline float rmin()  {
00254 #ifdef HAVE_NUMERIC_LIMITS
00255       return std::numeric_limits<float>::min();
00256 #else
00257       LAPACK<int, float> lp; return lp.LAMCH('U');
00258 #endif
00259     };
00260     static inline float emax()  {
00261 #ifdef HAVE_NUMERIC_LIMITS
00262       return std::numeric_limits<float>::max_exponent;
00263 #else
00264       LAPACK<int, float> lp; return lp.LAMCH('L');
00265 #endif
00266     };
00267     static inline float rmax()  {
00268 #ifdef HAVE_NUMERIC_LIMITS
00269       return std::numeric_limits<float>::max();
00270 #else
00271       LAPACK<int, float> lp; return lp.LAMCH('O');
00272 #endif
00273     };
00274     static inline magnitudeType magnitude(float a) { return fabs(a); };    
00275     static inline float zero()  { return(0.0); };
00276     static inline float one()   { return(1.0); };    
00277     static inline float conjugate(float x)   { return(x); };    
00278     static inline float real(float x) { return x; };
00279     static inline float imag(float x) { return 0; };
00280     static inline float nan() {
00281 #ifdef __sun
00282       return 0.0/sin(0.0);
00283 #else
00284       return flt_nan;
00285 #endif
00286     };
00287     static inline bool isnaninf(float x) { // RAB: 2004/05/28: Taken from NOX_StatusTest_FiniteValue.C
00288       const float tol = 1e-6; // Any (bounded) number should do!
00289       if( !(x <= tol) && !(x > tol) ) return true;                 // IEEE says this should fail for NaN
00290       float z=0.0*x; if( !(z <= tol) && !(z > tol) ) return true;  // Use fact that Inf*0 = NaN
00291       return false;
00292     }
00293     static inline void seedrandom(unsigned int s) { srand(s); };
00294     static inline float random() { float rnd = (float) rand() / RAND_MAX; return (float)(-1.0 + 2.0 * rnd); };
00295     static inline std::string name() { return "float"; };
00296     static inline float squareroot(float x) { return sqrt(x); };
00297   };
00298 
00299 #ifndef __sun
00300   extern const double dbl_nan;
00301 #endif
00302  
00303   template<>
00304   struct ScalarTraits<double>
00305   {
00306     typedef double magnitudeType;
00307     static const bool isComplex = false;
00308     static const bool isComparable = true;
00309     static const bool hasMachineParameters = true;
00310     static inline double eps()   {
00311 #ifdef HAVE_NUMERIC_LIMITS
00312       return std::numeric_limits<double>::epsilon();
00313 #else
00314       LAPACK<int, double> lp; return lp.LAMCH('E');
00315 #endif
00316     }
00317     static inline double sfmin() {
00318 #ifdef HAVE_NUMERIC_LIMITS
00319       return std::numeric_limits<double>::min();
00320 #else
00321       LAPACK<int, double> lp; return lp.LAMCH('S');
00322 #endif
00323     };
00324     static inline double base()  {
00325 #ifdef HAVE_NUMERIC_LIMITS
00326       return std::numeric_limits<double>::radix;
00327 #else
00328       LAPACK<int, double> lp; return lp.LAMCH('B');
00329 #endif
00330     };
00331     static inline double prec()  {
00332 #ifdef HAVE_NUMERIC_LIMITS
00333       return eps()*base();
00334 #else
00335       LAPACK<int, double> lp; return lp.LAMCH('P');
00336 #endif
00337     };
00338     static inline double t()     {
00339 #ifdef HAVE_NUMERIC_LIMITS
00340       return std::numeric_limits<double>::digits;
00341 #else
00342       LAPACK<int, double> lp; return lp.LAMCH('N');
00343 #endif
00344     };
00345     static inline double rnd()   {
00346 #ifdef HAVE_NUMERIC_LIMITS
00347       return ( std::numeric_limits<double>::round_style == std::round_to_nearest ? double(1.0) : double(0.0) );
00348 #else
00349       LAPACK<int, double> lp; return lp.LAMCH('R');
00350 #endif
00351     };
00352     static inline double emin()  {
00353 #ifdef HAVE_NUMERIC_LIMITS
00354       return std::numeric_limits<double>::min_exponent;
00355 #else
00356       LAPACK<int, double> lp; return lp.LAMCH('M');
00357 #endif
00358     };
00359     static inline double rmin()  {
00360 #ifdef HAVE_NUMERIC_LIMITS
00361       return std::numeric_limits<double>::min();
00362 #else
00363       LAPACK<int, double> lp; return lp.LAMCH('U');
00364 #endif
00365     };
00366     static inline double emax()  {
00367 #ifdef HAVE_NUMERIC_LIMITS
00368       return std::numeric_limits<double>::max_exponent;
00369 #else
00370       LAPACK<int, double> lp; return lp.LAMCH('L');
00371 #endif
00372     };
00373     static inline double rmax()  {
00374 #ifdef HAVE_NUMERIC_LIMITS
00375       return std::numeric_limits<double>::max();
00376 #else
00377       LAPACK<int, double> lp; return lp.LAMCH('O');
00378 #endif
00379     };
00380     static inline magnitudeType magnitude(double a) { return fabs(a); };
00381     static inline double zero()  { return 0.0; };
00382     static inline double one()   { return 1.0; };
00383     static inline double conjugate(double x)   { return(x); };    
00384     static inline double real(double x) { return(x); };
00385     static inline double imag(double x) { return(0); };
00386     static inline double nan() {
00387 #ifdef __sun
00388       return 0.0/sin(0.0);
00389 #else
00390       return dbl_nan;
00391 #endif
00392     };
00393     static inline bool isnaninf(double x) { // RAB: 2004/05/28: Taken from NOX_StatusTest_FiniteValue.C
00394       const double tol = 1e-6; // Any (bounded) number should do!
00395       if( !(x <= tol) && !(x > tol) ) return true;                  // IEEE says this should fail for NaN
00396       double z=0.0*x; if( !(z <= tol) && !(z > tol) ) return true;  // Use fact that Inf*0 = NaN
00397       return false;
00398     }
00399     static inline void seedrandom(unsigned int s) { srand(s); };
00400     static inline double random() { double rnd = (double) rand() / RAND_MAX; return (double)(-1.0 + 2.0 * rnd); };
00401     static inline std::string name() { return "double"; };
00402     static inline double squareroot(double x) { return ::sqrt(x); };
00403   };
00404 
00405 #ifdef HAVE_TEUCHOS_GNU_MP
00406 
00407   extern gmp_randclass gmp_rng; 
00408 
00409   template<>
00410   struct ScalarTraits<mpf_class>
00411   {
00412     typedef mpf_class magnitudeType;
00413     static const bool isComplex = false;
00414     static const bool isComparable = true;
00415     static const bool hasMachineParameters = false;
00416     // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
00417     static magnitudeType magnitude(mpf_class a) { return abs(a); };
00418     static inline mpf_class zero() { mpf_class zero = 0.0; return zero; };
00419     static inline mpf_class one() { mpf_class one = 1.0; return one; };    
00420     static inline mpf_class conjugate(mpf_class x) { return x; };
00421     static inline mpf_class real(mpf_class x) { return(x); };
00422     static inline mpf_class imag(mpf_class x) { return(0); };
00423     static inline bool isnaninf(mpf_class x) { return false; } // mpf_class currently can't handle nan or inf!
00424     static inline void seedrandom(unsigned int s) { 
00425       unsigned long int seedVal = static_cast<unsigned long int>(s);
00426       gmp_rng.seed( seedVal );  
00427     };
00428     static inline mpf_class random() { 
00429       return gmp_rng.get_f(); 
00430     };
00431     static inline std::string name() { return "mpf_class"; };
00432     static inline mpf_class squareroot(mpf_class x) { return sqrt(x); };
00433     // Todo: RAB: 2004/05/28: Add nan() and isnaninf() functions when needed!
00434   };
00435 
00436 #endif  
00437 
00438 #ifdef HAVE_TEUCHOS_ARPREC
00439 
00440   template<>
00441   struct ScalarTraits<mp_real>
00442   {
00443     typedef mp_real magnitudeType;
00444     static const bool isComplex = false;
00445     static const bool isComparable = true;
00446     static const bool hasMachineParameters = false;
00447     // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
00448     static magnitudeType magnitude(mp_real a) { return abs(a); };
00449     static inline mp_real zero() { mp_real zero = 0.0; return zero; };
00450     static inline mp_real one() { mp_real one = 1.0; return one; };    
00451     static inline mp_real conjugate(mp_real x) { return x; };
00452     static inline mp_real real(mp_real x) { return(x); };
00453     static inline mp_real imag(mp_real x) { return(0); };
00454     static inline bool isnaninf(mp_real x) { return false; } // ToDo: Change this?
00455     static inline void seedrandom(unsigned int s) { 
00456       long int seedVal = static_cast<long int>(s);
00457       srand48(seedVal);
00458     };
00459     static inline mp_real random() { return mp_rand(); };
00460     static inline std::string name() { return "mp_real"; };
00461     static inline mp_real squareroot(mp_real x) { return sqrt(x); };
00462     // Todo: RAB: 2004/05/28: Add nan() and isnaninf() functions when needed!
00463   };
00464   
00465 #endif // HAVE_TEUCHOS_ARPREC
00466  
00467 #if ( defined(HAVE_COMPLEX) || defined(HAVE_COMPLEX_H) ) && defined(HAVE_TEUCHOS_COMPLEX)
00468 
00469   // Partial specialization for complex numbers templated on real type T
00470   template<class T> 
00471   struct ScalarTraits<
00472 #if defined(HAVE_COMPLEX)
00473     std::complex<T>
00474 #elif  defined(HAVE_COMPLEX_H)
00475     ::complex<T>
00476 #endif
00477     >
00478   {
00479 #if defined(HAVE_COMPLEX)
00480     typedef std::complex<T>  ComplexT;
00481 #elif  defined(HAVE_COMPLEX_H)
00482     typedef ::complex<T>     ComplexT;
00483 #endif
00484     typedef typename ScalarTraits<T>::magnitudeType magnitudeType;
00485     static const bool isComplex = true;
00486     static const bool isComparable = false;
00487     static const bool hasMachineParameters = true;
00488     static inline magnitudeType eps()          { return ScalarTraits<magnitudeType>::eps(); };
00489     static inline magnitudeType sfmin()        { return ScalarTraits<magnitudeType>::sfmin(); };
00490     static inline magnitudeType base()         { return ScalarTraits<magnitudeType>::base(); };
00491     static inline magnitudeType prec()         { return ScalarTraits<magnitudeType>::prec(); };
00492     static inline magnitudeType t()            { return ScalarTraits<magnitudeType>::t(); };
00493     static inline magnitudeType rnd()          { return ScalarTraits<magnitudeType>::rnd(); };
00494     static inline magnitudeType emin()         { return ScalarTraits<magnitudeType>::emin(); };
00495     static inline magnitudeType rmin()         { return ScalarTraits<magnitudeType>::rmin(); };
00496     static inline magnitudeType emax()         { return ScalarTraits<magnitudeType>::emax(); };
00497     static inline magnitudeType rmax()         { return ScalarTraits<magnitudeType>::rmax(); };
00498     static magnitudeType magnitude(ComplexT a) { return std::abs(a); };
00499     static inline ComplexT zero()              { return ComplexT(ScalarTraits<magnitudeType>::zero(),ScalarTraits<magnitudeType>::zero()); };
00500     static inline ComplexT one()               { return ComplexT(ScalarTraits<magnitudeType>::one(),ScalarTraits<magnitudeType>::zero()); };
00501     static inline ComplexT conjugate(ComplexT a){ return ComplexT(a.real(),-a.imag()); };
00502     static inline magnitudeType real(ComplexT a) { return a.real(); };
00503     static inline magnitudeType imag(ComplexT a) { return a.imag(); };
00504     static inline ComplexT nan()               { return ComplexT(ScalarTraits<magnitudeType>::nan(),ScalarTraits<magnitudeType>::nan()); };
00505     static inline bool isnaninf(ComplexT x)    { return ScalarTraits<magnitudeType>::isnaninf(x.real()) || ScalarTraits<magnitudeType>::isnaninf(x.imag()); };
00506     static inline void seedrandom(unsigned int s) { ScalarTraits<magnitudeType>::seedrandom(s); };
00507     static inline ComplexT random()
00508     {
00509       const T rnd1 = ScalarTraits<magnitudeType>::random();
00510       const T rnd2 = ScalarTraits<magnitudeType>::random();
00511       return ComplexT(rnd1,rnd2);
00512     };
00513     static inline std::string name() { return std::string("std::complex<")+std::string(ScalarTraits<magnitudeType>::name())+std::string(">"); };
00514     // This will only return one of the square roots of x, the other can be obtained by taking its conjugate
00515     static inline ComplexT squareroot(ComplexT x)
00516     {
00517       typedef ScalarTraits<magnitudeType>  STMT;
00518       const T r  = x.real(), i = x.imag();
00519       const T a  = STMT::squareroot((r*r)+(i*i));
00520       const T nr = STMT::squareroot((a+r)/2);
00521       const T ni = STMT::squareroot((a-r)/2);
00522       return ComplexT(nr,ni);
00523     };
00524   };
00525 
00526 #endif //  HAVE_COMPLEX || HAVE_COMPLEX_H
00527 
00528 #endif // DOXYGEN_SHOULD_SKIP_THIS
00529 
00530 } // Teuchos namespace
00531 
00532 #endif // _TEUCHOS_SCALARTRAITS_HPP_