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