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