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

Generated on Tue Jul 13 09:23:00 2010 for Teuchos - Trilinos Tools Package by  doxygen 1.4.7