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

Generated on Wed May 12 21:24:41 2010 for Teuchos - Trilinos Tools Package by  doxygen 1.4.7