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 #ifndef _TEUCHOS_SCALARTRAITS_HPP_
00041 #define _TEUCHOS_SCALARTRAITS_HPP_
00042 
00047 #include "Teuchos_ConfigDefs.hpp"
00048 
00049 #ifndef HAVE_NUMERIC_LIMITS
00050 #include "Teuchos_LAPACK.hpp"
00051 #endif
00052 
00053 #ifdef HAVE_TEUCHOS_ARPREC
00054 #include "mp/mpreal.h"
00055 #endif
00056 
00057 #ifdef HAVE_TEUCHOS_GNU_MP
00058 #include "gmp.h"
00059 #include "gmpxx.h"
00060 #endif
00061 
00084 /* This is the default structure used by ScalarTraits<T> to produce a compile time
00085   error when the specialization does not exist for type <tt>T</tt>.
00086 */
00087 namespace Teuchos {
00088 
00089 template<class T>
00090 struct UndefinedScalarTraits
00091 {
00093   static inline T notDefined() { return T::this_type_is_missing_a_specialization(); }
00094 };
00095 
00096 template<class T>
00097 struct ScalarTraits
00098 {
00100   typedef T magnitudeType;
00102   static const bool isComplex = false;
00104   static const bool isComparable = false;
00106   static const bool hasMachineParameters = false;
00108   static inline magnitudeType eps()   { return UndefinedScalarTraits<T>::notDefined(); }
00110   static inline magnitudeType sfmin() { return UndefinedScalarTraits<T>::notDefined(); }
00112   static inline magnitudeType base()  { return UndefinedScalarTraits<T>::notDefined(); }
00114   static inline magnitudeType prec()  { return UndefinedScalarTraits<T>::notDefined(); }
00116   static inline magnitudeType t()     { return UndefinedScalarTraits<T>::notDefined(); }
00118   static inline magnitudeType rnd()   { return UndefinedScalarTraits<T>::notDefined(); }
00120   static inline magnitudeType emin()  { return UndefinedScalarTraits<T>::notDefined(); }
00122   static inline magnitudeType rmin()  { return UndefinedScalarTraits<T>::notDefined(); }
00124   static inline magnitudeType emax()  { return UndefinedScalarTraits<T>::notDefined(); }
00126   static inline magnitudeType rmax()  { return UndefinedScalarTraits<T>::notDefined(); }
00128   static inline magnitudeType magnitude(T a) { return UndefinedScalarTraits<T>::notDefined(); }
00130   static inline T zero()                     { return UndefinedScalarTraits<T>::notDefined(); }
00132   static inline T one()                      { return UndefinedScalarTraits<T>::notDefined(); }
00134   static inline magnitudeType real(T a) { return UndefinedScalarTraits<T>::notDefined(); }
00136   static inline magnitudeType imag(T a) { return UndefinedScalarTraits<T>::notDefined(); }
00138   static inline T conjugate(T a) { return UndefinedScalarTraits<T>::notDefined(); }
00140   static inline T nan()                      { return UndefinedScalarTraits<T>::notDefined(); }
00142   static inline bool isnaninf(const T& x)     { return UndefinedScalarTraits<T>::notDefined(); }
00144   static inline void seedrandom(unsigned int s) { int i; T t = &i; }
00146   static inline T random()                   { return UndefinedScalarTraits<T>::notDefined(); }
00148   static inline std::string name()           { (void)UndefinedScalarTraits<T>::notDefined(); return 0; }
00150   static inline T squareroot(T x) { return UndefinedScalarTraits<T>::notDefined(); }
00152   static inline T pow(T x, T y) { return UndefinedScalarTraits<T>::notDefined(); }
00153 };
00154   
00155 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00156 
00157 
00158 void throwScalarTraitsNanInfError( const std::string &errMsg );
00159 
00160 
00161 #define TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR( VALUE, MSG ) \
00162   if (isnaninf(VALUE)) { \
00163     std::ostringstream omsg; \
00164     omsg << MSG; \
00165     throwScalarTraitsNanInfError(omsg.str()); \
00166   }
00167 
00168 
00169 template<>
00170 struct ScalarTraits<char>
00171 {
00172   typedef char magnitudeType;
00173   static const bool isComplex = false;
00174   static const bool isComparable = true;
00175   static const bool hasMachineParameters = false;
00176   // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
00177   static inline magnitudeType magnitude(char a) { return static_cast<char>(std::fabs(static_cast<double>(a))); }
00178   static inline char zero()  { return 0; }
00179   static inline char one()   { return 1; }
00180   static inline char conjugate(char x) { return x; }
00181   static inline char real(char x) { return x; }
00182   static inline char imag(char) { return 0; }
00183   static inline bool isnaninf(char ) { return false; }
00184   static inline void seedrandom(unsigned int s) { 
00185     std::srand(s); 
00186 #ifdef __APPLE__
00187     // throw away first random number to address bug 3655
00188     // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
00189     random();
00190 #endif
00191   }
00192   //static inline char random() { return (-1 + 2*rand()); } // RAB: This version should be used to be consistent with others
00193   static inline char random() { return std::rand(); } // RAB: This version should be used for an unsigned char, not char
00194   static inline std::string name() { return "char"; }
00195   static inline char squareroot(char x) { return (char) std::sqrt((double) x); }
00196   static inline char pow(char x, char y) { return (char) std::pow((double)x,(double)y); }
00197 };
00198 
00199 template<>
00200 struct ScalarTraits<short int>
00201 {
00202   typedef short int magnitudeType;
00203   static const bool isComplex = false;
00204   static const bool isComparable = true;
00205   static const bool hasMachineParameters = false;
00206   // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
00207   static inline magnitudeType magnitude(short int a) { return static_cast<short int>(std::fabs(static_cast<double>(a))); }
00208   static inline short int zero()  { return 0; }
00209   static inline short int one()   { return 1; }
00210   static inline short int conjugate(short int x) { return x; }
00211   static inline short int real(short int x) { return x; }
00212   static inline short int imag(short int) { return 0; }
00213   static inline void seedrandom(unsigned int s) { 
00214     std::srand(s); 
00215 #ifdef __APPLE__
00216     // throw away first random number to address bug 3655
00217     // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
00218     random();
00219 #endif
00220   }
00221   //static inline int random() { return (-1 + 2*rand()); }  // RAB: This version should be used to be consistent with others
00222   static inline short int random() { return std::rand(); }             // RAB: This version should be used for an unsigned int, not int
00223   static inline std::string name() { return "short int"; }
00224   static inline short int squareroot(short int x) { return (short int) std::sqrt((double) x); }
00225   static inline short int pow(short int x, short int y) { return (short int) std::pow((double)x,(double)y); }
00226 };
00227 
00228 template<>
00229 struct ScalarTraits<int>
00230 {
00231   typedef int magnitudeType;
00232   static const bool isComplex = false;
00233   static const bool isComparable = true;
00234   static const bool hasMachineParameters = false;
00235   // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
00236   static inline magnitudeType magnitude(int a) { return static_cast<int>(std::fabs(static_cast<double>(a))); }
00237   static inline int zero()  { return 0; }
00238   static inline int one()   { return 1; }
00239   static inline int conjugate(int x) { return x; }
00240   static inline int real(int x) { return x; }
00241   static inline int imag(int) { return 0; }
00242   static inline bool isnaninf(int) { return false; }
00243   static inline void seedrandom(unsigned int s) { 
00244     std::srand(s); 
00245 #ifdef __APPLE__
00246     // throw away first random number to address bug 3655
00247     // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
00248     random();
00249 #endif
00250   }
00251   //static inline int random() { return (-1 + 2*rand()); }  // RAB: This version should be used to be consistent with others
00252   static inline int random() { return std::rand(); }             // RAB: This version should be used for an unsigned int, not int
00253   static inline std::string name() { return "int"; }
00254   static inline int squareroot(int x) { return (int) std::sqrt((double) x); }
00255   static inline int pow(int x, int y) { return (int) std::pow((double)x,(double)y); }
00256 };
00257 
00258 template<>
00259 struct ScalarTraits<long int>
00260 {
00261   typedef long int magnitudeType;
00262   static const bool isComplex = false;
00263   static const bool isComparable = true;
00264   static const bool hasMachineParameters = false;
00265   // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
00266   static inline magnitudeType magnitude(long int a) { return static_cast<long int>(std::fabs(static_cast<double>(a))); }
00267   static inline long int zero()  { return 0; }
00268   static inline long int one()   { return 1; }
00269   static inline long int conjugate(long int x) { return x; }
00270   static inline long int real(long int x) { return x; }
00271   static inline long int imag(long int) { return 0; }
00272   static inline void seedrandom(unsigned int s) { 
00273     std::srand(s); 
00274 #ifdef __APPLE__
00275     // throw away first random number to address bug 3655
00276     // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
00277     random();
00278 #endif
00279   }
00280   //static inline int random() { return (-1 + 2*rand()); }  // RAB: This version should be used to be consistent with others
00281   static inline long int random() { return std::rand(); }             // RAB: This version should be used for an unsigned int, not int
00282   static inline std::string name() { return "long int"; }
00283   static inline long int squareroot(long int x) { return (long int) std::sqrt((double) x); }
00284   static inline long int pow(long int x, long int y) { return (long int) std::pow((double)x,(double)y); }
00285 };
00286 
00287 #ifndef __sun
00288 extern const float flt_nan;
00289 #endif
00290  
00291 template<>
00292 struct ScalarTraits<float>
00293 {
00294   typedef float magnitudeType;
00295   static const bool isComplex = false;
00296   static const bool isComparable = true;
00297   static const bool hasMachineParameters = true;
00298   static inline float eps()   {
00299 #ifdef HAVE_NUMERIC_LIMITS
00300     return std::numeric_limits<float>::epsilon();
00301 #else
00302     LAPACK<int, float> lp; return lp.LAMCH('E');
00303 #endif
00304   }
00305   static inline float sfmin() {
00306 #ifdef HAVE_NUMERIC_LIMITS
00307     return std::numeric_limits<float>::min();
00308 #else
00309     LAPACK<int, float> lp; return lp.LAMCH('S');
00310 #endif
00311   }
00312   static inline float base()  {
00313 #ifdef HAVE_NUMERIC_LIMITS
00314     return std::numeric_limits<float>::radix;
00315 #else
00316     LAPACK<int, float> lp; return lp.LAMCH('B');
00317 #endif
00318   }
00319   static inline float prec()  {
00320 #ifdef HAVE_NUMERIC_LIMITS
00321     return eps()*base();
00322 #else
00323     LAPACK<int, float> lp; return lp.LAMCH('P');
00324 #endif
00325   }
00326   static inline float t()     {
00327 #ifdef HAVE_NUMERIC_LIMITS
00328     return std::numeric_limits<float>::digits;
00329 #else
00330     LAPACK<int, float> lp; return lp.LAMCH('N');
00331 #endif
00332   }
00333   static inline float rnd()   {
00334 #ifdef HAVE_NUMERIC_LIMITS
00335     return ( std::numeric_limits<float>::round_style == std::round_to_nearest ? float(1.0) : float(0.0) );
00336 #else
00337     LAPACK<int, float> lp; return lp.LAMCH('R');
00338 #endif
00339   }
00340   static inline float emin()  {
00341 #ifdef HAVE_NUMERIC_LIMITS
00342     return std::numeric_limits<float>::min_exponent;
00343 #else
00344     LAPACK<int, float> lp; return lp.LAMCH('M');
00345 #endif
00346   }
00347   static inline float rmin()  {
00348 #ifdef HAVE_NUMERIC_LIMITS
00349     return std::numeric_limits<float>::min();
00350 #else
00351     LAPACK<int, float> lp; return lp.LAMCH('U');
00352 #endif
00353   }
00354   static inline float emax()  {
00355 #ifdef HAVE_NUMERIC_LIMITS
00356     return std::numeric_limits<float>::max_exponent;
00357 #else
00358     LAPACK<int, float> lp; return lp.LAMCH('L');
00359 #endif
00360   }
00361   static inline float rmax()  {
00362 #ifdef HAVE_NUMERIC_LIMITS
00363     return std::numeric_limits<float>::max();
00364 #else
00365     LAPACK<int, float> lp; return lp.LAMCH('O');
00366 #endif
00367   }
00368   static inline magnitudeType magnitude(float a)
00369     {
00370 #ifdef TEUCHOS_DEBUG
00371       TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
00372         a, "Error, the input value to magnitude(...) a = " << a << " can not be NaN!" );
00373 #endif      
00374       return std::fabs(a);
00375     }    
00376   static inline float zero()  { return(0.0); }
00377   static inline float one()   { return(1.0); }    
00378   static inline float conjugate(float x)   { return(x); }    
00379   static inline float real(float x) { return x; }
00380   static inline float imag(float) { return 0; }
00381   static inline float nan() {
00382 #ifdef __sun
00383     return 0.0/std::sin(0.0);
00384 #else
00385     return flt_nan;
00386 #endif
00387   }
00388   static inline bool isnaninf(float x) { // RAB: 2004/05/28: Taken from NOX_StatusTest_FiniteValue.C
00389     const float tol = 1e-6; // Any (bounded) number should do!
00390     if( !(x <= tol) && !(x > tol) ) return true;                 // IEEE says this should fail for NaN
00391     float z=0.0*x; if( !(z <= tol) && !(z > tol) ) return true;  // Use fact that Inf*0 = NaN
00392     return false;
00393   }
00394   static inline void seedrandom(unsigned int s) { 
00395     std::srand(s); 
00396 #ifdef __APPLE__
00397     // throw away first random number to address bug 3655
00398     // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
00399     random();
00400 #endif
00401   }
00402   static inline float random() { float rnd = (float) std::rand() / RAND_MAX; return (float)(-1.0 + 2.0 * rnd); }
00403   static inline std::string name() { return "float"; }
00404   static inline float squareroot(float x)
00405     {
00406 #ifdef TEUCHOS_DEBUG
00407       TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
00408         x, "Error, the input value to squareroot(...) x = " << x << " can not be NaN!" );
00409 #endif
00410       errno = 0;
00411       const float rtn = std::sqrt(x);
00412       if (errno)
00413         return nan();
00414       return rtn;
00415     }
00416   static inline float pow(float x, float y) { return std::pow(x,y); }
00417 };
00418 
00419 #ifndef __sun
00420 extern const double dbl_nan;
00421 #endif
00422  
00423 template<>
00424 struct ScalarTraits<double>
00425 {
00426   typedef double magnitudeType;
00427   static const bool isComplex = false;
00428   static const bool isComparable = true;
00429   static const bool hasMachineParameters = true;
00430   static inline double eps()   {
00431 #ifdef HAVE_NUMERIC_LIMITS
00432     return std::numeric_limits<double>::epsilon();
00433 #else
00434     LAPACK<int, double> lp; return lp.LAMCH('E');
00435 #endif
00436   }
00437   static inline double sfmin() {
00438 #ifdef HAVE_NUMERIC_LIMITS
00439     return std::numeric_limits<double>::min();
00440 #else
00441     LAPACK<int, double> lp; return lp.LAMCH('S');
00442 #endif
00443   }
00444   static inline double base()  {
00445 #ifdef HAVE_NUMERIC_LIMITS
00446     return std::numeric_limits<double>::radix;
00447 #else
00448     LAPACK<int, double> lp; return lp.LAMCH('B');
00449 #endif
00450   }
00451   static inline double prec()  {
00452 #ifdef HAVE_NUMERIC_LIMITS
00453     return eps()*base();
00454 #else
00455     LAPACK<int, double> lp; return lp.LAMCH('P');
00456 #endif
00457   }
00458   static inline double t()     {
00459 #ifdef HAVE_NUMERIC_LIMITS
00460     return std::numeric_limits<double>::digits;
00461 #else
00462     LAPACK<int, double> lp; return lp.LAMCH('N');
00463 #endif
00464   }
00465   static inline double rnd()   {
00466 #ifdef HAVE_NUMERIC_LIMITS
00467     return ( std::numeric_limits<double>::round_style == std::round_to_nearest ? double(1.0) : double(0.0) );
00468 #else
00469     LAPACK<int, double> lp; return lp.LAMCH('R');
00470 #endif
00471   }
00472   static inline double emin()  {
00473 #ifdef HAVE_NUMERIC_LIMITS
00474     return std::numeric_limits<double>::min_exponent;
00475 #else
00476     LAPACK<int, double> lp; return lp.LAMCH('M');
00477 #endif
00478   }
00479   static inline double rmin()  {
00480 #ifdef HAVE_NUMERIC_LIMITS
00481     return std::numeric_limits<double>::min();
00482 #else
00483     LAPACK<int, double> lp; return lp.LAMCH('U');
00484 #endif
00485   }
00486   static inline double emax()  {
00487 #ifdef HAVE_NUMERIC_LIMITS
00488     return std::numeric_limits<double>::max_exponent;
00489 #else
00490     LAPACK<int, double> lp; return lp.LAMCH('L');
00491 #endif
00492   }
00493   static inline double rmax()  {
00494 #ifdef HAVE_NUMERIC_LIMITS
00495     return std::numeric_limits<double>::max();
00496 #else
00497     LAPACK<int, double> lp; return lp.LAMCH('O');
00498 #endif
00499   }
00500   static inline magnitudeType magnitude(double a)
00501     {
00502 #ifdef TEUCHOS_DEBUG
00503       TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
00504         a, "Error, the input value to magnitude(...) a = " << a << " can not be NaN!" );
00505 #endif      
00506       return std::fabs(a);
00507     }
00508   static inline double zero()  { return 0.0; }
00509   static inline double one()   { return 1.0; }
00510   static inline double conjugate(double x)   { return(x); }    
00511   static inline double real(double x) { return(x); }
00512   static inline double imag(double) { return(0); }
00513   static inline double nan() {
00514 #ifdef __sun
00515     return 0.0/std::sin(0.0);
00516 #else
00517     return dbl_nan;
00518 #endif
00519   }
00520   static inline bool isnaninf(double x) { // RAB: 2004/05/28: Taken from NOX_StatusTest_FiniteValue.C
00521     const double tol = 1e-6; // Any (bounded) number should do!
00522     if( !(x <= tol) && !(x > tol) ) return true;                  // IEEE says this should fail for NaN
00523     double z=0.0*x; if( !(z <= tol) && !(z > tol) ) return true;  // Use fact that Inf*0 = NaN
00524     return false;
00525   }
00526   static inline void seedrandom(unsigned int s) { 
00527     std::srand(s); 
00528 #ifdef __APPLE__
00529     // throw away first random number to address bug 3655
00530     // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
00531     random();
00532 #endif
00533   }
00534   static inline double random() { double rnd = (double) std::rand() / RAND_MAX; return (double)(-1.0 + 2.0 * rnd); }
00535   static inline std::string name() { return "double"; }
00536   static inline double squareroot(double x)
00537     {
00538 #ifdef TEUCHOS_DEBUG
00539       TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
00540         x, "Error, the input value to squareroot(...) x = " << x << " can not be NaN!" );
00541 #endif      
00542       errno = 0;
00543       const double rtn = std::sqrt(x);
00544       if (errno)
00545         return nan();
00546       return rtn;
00547     }
00548   static inline double pow(double x, double y) { return std::pow(x,y); }
00549 };
00550 
00551 #ifdef HAVE_TEUCHOS_GNU_MP
00552 
00553 extern gmp_randclass gmp_rng; 
00554 
00555 template<>
00556 struct ScalarTraits<mpf_class>
00557 {
00558   typedef mpf_class magnitudeType;
00559   static const bool isComplex = false;
00560   static const bool isComparable = true;
00561   static const bool hasMachineParameters = false;
00562   // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
00563   static magnitudeType magnitude(mpf_class a) { return std::abs(a); }
00564   static inline mpf_class zero() { mpf_class zero = 0.0; return zero; }
00565   static inline mpf_class one() { mpf_class one = 1.0; return one; }    
00566   static inline mpf_class conjugate(mpf_class x) { return x; }
00567   static inline mpf_class real(mpf_class x) { return(x); }
00568   static inline mpf_class imag(mpf_class x) { return(0); }
00569   static inline bool isnaninf(mpf_class x) { return false; } // mpf_class currently can't handle nan or inf!
00570   static inline void seedrandom(unsigned int s) { 
00571     unsigned long int seedVal = static_cast<unsigned long int>(s);
00572     gmp_rng.seed( seedVal );  
00573   }
00574   static inline mpf_class random() { 
00575     return gmp_rng.get_f(); 
00576   }
00577   static inline std::string name() { return "mpf_class"; }
00578   static inline mpf_class squareroot(mpf_class x) { return std::sqrt(x); }
00579   static inline mpf_class pow(mpf_class x, mpf_class y) { return pow(x,y); }
00580   // Todo: RAB: 2004/05/28: Add nan() and isnaninf() functions when needed!
00581 };
00582 
00583 #endif  
00584 
00585 #ifdef HAVE_TEUCHOS_ARPREC
00586 
00587 template<>
00588 struct ScalarTraits<mp_real>
00589 {
00590   typedef mp_real magnitudeType;
00591   static const bool isComplex = false;
00592   static const bool isComparable = true;
00593   static const bool hasMachineParameters = false;
00594   // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
00595   static magnitudeType magnitude(mp_real a) { return std::abs(a); }
00596   static inline mp_real zero() { mp_real zero = 0.0; return zero; }
00597   static inline mp_real one() { mp_real one = 1.0; return one; }    
00598   static inline mp_real conjugate(mp_real x) { return x; }
00599   static inline mp_real real(mp_real x) { return(x); }
00600   static inline mp_real imag(mp_real x) { return(0); }
00601   static inline bool isnaninf(mp_real x) { return false; } // ToDo: Change this?
00602   static inline void seedrandom(unsigned int s) { 
00603     long int seedVal = static_cast<long int>(s);
00604     srand48(seedVal);
00605   }
00606   static inline mp_real random() { return mp_rand(); }
00607   static inline std::string name() { return "mp_real"; }
00608   static inline mp_real squareroot(mp_real x) { return std::sqrt(x); }
00609   static inline mp_real pow(mp_real x, mp_real y) { return pow(x,y); }
00610   // Todo: RAB: 2004/05/28: Add nan() and isnaninf() functions when needed!
00611 };
00612   
00613 #endif // HAVE_TEUCHOS_ARPREC
00614  
00615 #if ( defined(HAVE_COMPLEX) || defined(HAVE_COMPLEX_H) ) && defined(HAVE_TEUCHOS_COMPLEX)
00616 
00617 // Partial specialization for std::complex numbers templated on real type T
00618 template<class T> 
00619 struct ScalarTraits<
00620 #if defined(HAVE_COMPLEX)
00621   std::complex<T>
00622 #elif  defined(HAVE_COMPLEX_H)
00623 std::complex<T>
00624 #endif
00625 >
00626 {
00627 #if defined(HAVE_COMPLEX)
00628   typedef std::complex<T>  ComplexT;
00629 #elif  defined(HAVE_COMPLEX_H)
00630   typedef std::complex<T>     ComplexT;
00631 #endif
00632   typedef typename ScalarTraits<T>::magnitudeType magnitudeType;
00633   static const bool isComplex = true;
00634   static const bool isComparable = false;
00635   static const bool hasMachineParameters = true;
00636   static inline magnitudeType eps()          { return ScalarTraits<magnitudeType>::eps(); }
00637   static inline magnitudeType sfmin()        { return ScalarTraits<magnitudeType>::sfmin(); }
00638   static inline magnitudeType base()         { return ScalarTraits<magnitudeType>::base(); }
00639   static inline magnitudeType prec()         { return ScalarTraits<magnitudeType>::prec(); }
00640   static inline magnitudeType t()            { return ScalarTraits<magnitudeType>::t(); }
00641   static inline magnitudeType rnd()          { return ScalarTraits<magnitudeType>::rnd(); }
00642   static inline magnitudeType emin()         { return ScalarTraits<magnitudeType>::emin(); }
00643   static inline magnitudeType rmin()         { return ScalarTraits<magnitudeType>::rmin(); }
00644   static inline magnitudeType emax()         { return ScalarTraits<magnitudeType>::emax(); }
00645   static inline magnitudeType rmax()         { return ScalarTraits<magnitudeType>::rmax(); }
00646   static magnitudeType magnitude(ComplexT a)
00647     {
00648 #ifdef TEUCHOS_DEBUG
00649       TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
00650         a, "Error, the input value to magnitude(...) a = " << a << " can not be NaN!" );
00651 #endif      
00652       return std::abs(a);
00653     }
00654   static inline ComplexT zero()              { return ComplexT(ScalarTraits<magnitudeType>::zero(),ScalarTraits<magnitudeType>::zero()); }
00655   static inline ComplexT one()               { return ComplexT(ScalarTraits<magnitudeType>::one(),ScalarTraits<magnitudeType>::zero()); }
00656   static inline ComplexT conjugate(ComplexT a){ return ComplexT(a.real(),-a.imag()); }
00657   static inline magnitudeType real(ComplexT a) { return a.real(); }
00658   static inline magnitudeType imag(ComplexT a) { return a.imag(); }
00659   static inline ComplexT nan()               { return ComplexT(ScalarTraits<magnitudeType>::nan(),ScalarTraits<magnitudeType>::nan()); }
00660   static inline bool isnaninf(ComplexT x)    { return ScalarTraits<magnitudeType>::isnaninf(x.real()) || ScalarTraits<magnitudeType>::isnaninf(x.imag()); }
00661   static inline void seedrandom(unsigned int s) { ScalarTraits<magnitudeType>::seedrandom(s); }
00662   static inline ComplexT random()
00663     {
00664       const T rnd1 = ScalarTraits<magnitudeType>::random();
00665       const T rnd2 = ScalarTraits<magnitudeType>::random();
00666       return ComplexT(rnd1,rnd2);
00667     }
00668   static inline std::string name() { return std::string("std::complex<")+std::string(ScalarTraits<magnitudeType>::name())+std::string(">"); }
00669   // This will only return one of the square roots of x, the other can be obtained by taking its conjugate
00670   static inline ComplexT squareroot(ComplexT x)
00671     {
00672 #ifdef TEUCHOS_DEBUG
00673       TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
00674         x, "Error, the input value to squareroot(...) x = " << x << " can not be NaN!" );
00675 #endif
00676       typedef ScalarTraits<magnitudeType>  STMT;
00677       const T r  = x.real(), i = x.imag();
00678       const T a  = STMT::squareroot((r*r)+(i*i));
00679       const T nr = STMT::squareroot((a+r)/2);
00680       const T ni = STMT::squareroot((a-r)/2);
00681       return ComplexT(nr,ni);
00682     }
00683   static inline ComplexT pow(ComplexT x, ComplexT y) { return pow(x,y); }
00684 };
00685 
00686 #endif //  HAVE_COMPLEX || HAVE_COMPLEX_H
00687 
00688 #endif // DOXYGEN_SHOULD_SKIP_THIS
00689 
00690 } // Teuchos namespace
00691 
00692 #endif // _TEUCHOS_SCALARTRAITS_HPP_

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