Teuchos - Trilinos Tools Package Version of the Day
Teuchos_ScalarTraits.hpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //                    Teuchos: Common Tools Package
00005 //                 Copyright (2004) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ***********************************************************************
00040 // @HEADER
00041 
00042 // NOTE: Before adding specializations of ScalarTraits, make sure that they do
00043 // not duplicate specializations already present in PyTrilinos (see
00044 // packages/PyTrilinos/src/Teuchos_Traits.i)
00045 
00046 // NOTE: halfPrecision and doublePrecision are not currently implemented for
00047 // ARPREC, GMP or the ordinal types (e.g., int, char)
00048 
00049 #ifndef _TEUCHOS_SCALARTRAITS_HPP_
00050 #define _TEUCHOS_SCALARTRAITS_HPP_
00051 
00056 #include "Teuchos_ConfigDefs.hpp"
00057 
00058 #ifdef HAVE_TEUCHOS_ARPREC
00059 #include <arprec/mp_real.h>
00060 #endif
00061 
00062 #ifdef HAVE_TEUCHOS_QD
00063 #include <qd/qd_real.h>
00064 #include <qd/dd_real.h>
00065 #endif
00066 
00067 #ifdef HAVE_TEUCHOS_GNU_MP
00068 #include <gmp.h>
00069 #include <gmpxx.h>
00070 #endif
00071 
00072 
00073 #include "Teuchos_ScalarTraitsDecl.hpp"
00074 
00075 
00076 namespace Teuchos {
00077 
00078 
00079 #ifndef DOXYGEN_SHOULD_SKIP_THIS
00080 
00081 
00082 TEUCHOSCORE_LIB_DLL_EXPORT
00083 void throwScalarTraitsNanInfError( const std::string &errMsg );
00084 
00085 
00086 template<class Scalar>
00087 bool generic_real_isnaninf(const Scalar &x)
00088 {
00089   typedef std::numeric_limits<Scalar> NL;
00090   // IEEE says this should fail for NaN (not all compilers do not obey IEEE)
00091   const Scalar tol = 1.0; // Any (bounded) number should do!
00092   if (!(x <= tol) && !(x > tol)) return true;
00093   // Use fact that Inf*0 = NaN (again, all compilers do not obey IEEE)
00094   Scalar z = static_cast<Scalar>(0.0) * x;
00095   if (!(z <= tol) && !(z > tol)) return true;
00096   // As a last result use comparisons
00097   if (x == NL::infinity() || x == -NL::infinity()) return true;
00098   // We give up and assume the number is finite
00099   return false;
00100 }
00101 
00102 
00103 #define TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR( VALUE, MSG ) \
00104   if (isnaninf(VALUE)) { \
00105     std::ostringstream omsg; \
00106     omsg << MSG; \
00107     Teuchos::throwScalarTraitsNanInfError(omsg.str());  \
00108   }
00109 
00110 
00111 template<>
00112 struct ScalarTraits<char>
00113 {
00114   typedef char magnitudeType;
00115   typedef char halfPrecision;
00116   typedef char doublePrecision;
00117   static const bool isComplex = false;
00118   static const bool isOrdinal = true;
00119   static const bool isComparable = true;
00120   static const bool hasMachineParameters = false;
00121   // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
00122   static inline magnitudeType magnitude(char a) { return static_cast<char>(std::fabs(static_cast<double>(a))); }
00123   static inline char zero()  { return 0; }
00124   static inline char one()   { return 1; }
00125   static inline char conjugate(char x) { return x; }
00126   static inline char real(char x) { return x; }
00127   static inline char imag(char) { return 0; }
00128   static inline bool isnaninf(char ) { return false; }
00129   static inline void seedrandom(unsigned int s) { 
00130     std::srand(s); 
00131 #ifdef __APPLE__
00132     // throw away first random number to address bug 3655
00133     // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
00134     random();
00135 #endif
00136   }
00137   //static inline char random() { return (-1 + 2*rand()); } // RAB: This version should be used to be consistent with others
00138   static inline char random() { return std::rand(); } // RAB: This version should be used for an unsigned char, not char
00139   static inline std::string name() { return "char"; }
00140   static inline char squareroot(char x) { return (char) std::sqrt((double) x); }
00141   static inline char pow(char x, char y) { return (char) std::pow((double)x,(double)y); }
00142   static inline char log(char x) { return static_cast<char> (std::log (static_cast<double> (x))); }
00143   static inline char log10(char x) { return static_cast<char> (std::log10 (static_cast<double> (x))); }
00144 };
00145 
00146 
00147 template<>
00148 struct ScalarTraits<short int>
00149 {
00150   typedef short int magnitudeType;
00151   typedef short int halfPrecision;
00152   typedef short int doublePrecision;
00153   static const bool isComplex = false;
00154   static const bool isOrdinal = true;
00155   static const bool isComparable = true;
00156   static const bool hasMachineParameters = false;
00157   // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
00158   static inline magnitudeType magnitude(short int a) { return static_cast<short int>(std::fabs(static_cast<double>(a))); }
00159   static inline short int zero()  { return 0; }
00160   static inline short int one()   { return 1; }
00161   static inline short int conjugate(short int x) { return x; }
00162   static inline short int real(short int x) { return x; }
00163   static inline short int imag(short int) { return 0; }
00164   static inline bool isnaninf(short int) { return false; }
00165   static inline void seedrandom(unsigned int s) { 
00166     std::srand(s); 
00167 #ifdef __APPLE__
00168     // throw away first random number to address bug 3655
00169     // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
00170     random();
00171 #endif
00172   }
00173   //static inline int random() { return (-1 + 2*rand()); }  // RAB: This version should be used to be consistent with others
00174   static inline short int random() { return std::rand(); }             // RAB: This version should be used for an unsigned int, not int
00175   static inline std::string name() { return "short int"; }
00176   static inline short int squareroot(short int x) { return (short int) std::sqrt((double) x); }
00177   static inline short int pow(short int x, short int y) { return (short int) std::pow((double)x,(double)y); }
00178   static inline short int log(short int x) { return static_cast<short int> (std::log (static_cast<double> (x))); }
00179   static inline short int log10(short int x) { return static_cast<short int> (std::log10 (static_cast<double> (x))); }
00180 };
00181 
00182 template<>
00183 struct ScalarTraits<unsigned short int>
00184 {
00185   typedef unsigned short int magnitudeType;
00186   typedef unsigned short int halfPrecision;
00187   typedef unsigned short int doublePrecision;
00188   static const bool isComplex = false;
00189   static const bool isOrdinal = true;
00190   static const bool isComparable = true;
00191   static const bool hasMachineParameters = false;
00192   // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
00193   static inline magnitudeType magnitude(unsigned short int a) { return static_cast<unsigned short int>(std::fabs(static_cast<double>(a))); }
00194   static inline unsigned short int zero()  { return 0; }
00195   static inline unsigned short int one()   { return 1; }
00196   static inline unsigned short int conjugate(unsigned short int x) { return x; }
00197   static inline unsigned short int real(unsigned short int x) { return x; }
00198   static inline unsigned short int imag(unsigned short int) { return 0; }
00199   static inline bool isnaninf(unsigned short int) { return false; }
00200   static inline void seedrandom(unsigned int s) { 
00201     std::srand(s); 
00202 #ifdef __APPLE__
00203     // throw away first random number to address bug 3655
00204     // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
00205     random();
00206 #endif
00207   }
00208   //static inline int random() { return (-1 + 2*rand()); }  // RAB: This version should be used to be consistent with others
00209   static inline unsigned short int random() { return std::rand(); }             // RAB: This version should be used for an unsigned int, not int
00210   static inline std::string name() { return "unsigned short int"; }
00211   static inline unsigned short int squareroot(unsigned short int x) { return (unsigned short int) std::sqrt((double) x); }
00212   static inline unsigned short int pow(unsigned short int x, unsigned short int y) { return (unsigned short int) std::pow((double)x,(double)y); }
00213   static inline unsigned short int log(unsigned short int x) { return static_cast<unsigned short int> (std::log (static_cast<double> (x))); }
00214   static inline unsigned short int log10(unsigned short int x) { return static_cast<unsigned short int> (std::log10 (static_cast<double> (x))); }
00215 };
00216 
00217 
00218 template<>
00219 struct ScalarTraits<int>
00220 {
00221   typedef int magnitudeType;
00222   typedef int halfPrecision;
00223   typedef int doublePrecision;
00224   static const bool isComplex = false;
00225   static const bool isOrdinal = true;
00226   static const bool isComparable = true;
00227   static const bool hasMachineParameters = false;
00228   // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
00229   static inline magnitudeType magnitude(int a) { return static_cast<int>(std::fabs(static_cast<double>(a))); }
00230   static inline int zero()  { return 0; }
00231   static inline int one()   { return 1; }
00232   static inline int conjugate(int x) { return x; }
00233   static inline int real(int x) { return x; }
00234   static inline int imag(int) { return 0; }
00235   static inline bool isnaninf(int) { return false; }
00236   static inline void seedrandom(unsigned int s) { 
00237     std::srand(s); 
00238 #ifdef __APPLE__
00239     // throw away first random number to address bug 3655
00240     // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
00241     random();
00242 #endif
00243   }
00244   //static inline int random() { return (-1 + 2*rand()); }  // RAB: This version should be used to be consistent with others
00245   static inline int random() { return std::rand(); }             // RAB: This version should be used for an unsigned int, not int
00246   static inline std::string name() { return "int"; }
00247   static inline int squareroot(int x) { return (int) std::sqrt((double) x); }
00248   static inline int pow(int x, int y) { return (int) std::pow((double)x,(double)y); }
00249   static inline int log(int x) { return static_cast<int> (std::log (static_cast<double> (x))); }
00250   static inline int log10(int x) { return static_cast<int> (std::log10 (static_cast<double> (x))); }
00251 };
00252 
00253 
00254 template<>
00255 struct ScalarTraits<unsigned int>
00256 {
00257   typedef unsigned int magnitudeType;
00258   typedef unsigned int halfPrecision;
00259   typedef unsigned int doublePrecision;
00260   static const bool isComplex = false;
00261   static const bool isOrdinal = true;
00262   static const bool isComparable = true;
00263   static const bool hasMachineParameters = false;
00264   // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
00265   static inline magnitudeType magnitude(unsigned int a) { return static_cast<unsigned int>(std::fabs(static_cast<double>(a))); }
00266   static inline unsigned int zero()  { return 0; }
00267   static inline unsigned int one()   { return 1; }
00268   static inline unsigned int conjugate(unsigned int x) { return x; }
00269   static inline unsigned int real(unsigned int x) { return x; }
00270   static inline unsigned int imag(unsigned int) { return 0; }
00271   static inline bool isnaninf(unsigned int) { return false; }
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 unsigned int random() { return std::rand(); }             // RAB: This version should be used for an unsigned int, not int
00282   static inline std::string name() { return "unsigned int"; }
00283   static inline unsigned int squareroot(unsigned int x) { return (unsigned int) std::sqrt((double) x); }
00284   static inline unsigned int pow(unsigned int x, unsigned int y) { return (unsigned int) std::pow((double)x,(double)y); }
00285   static inline unsigned int log(unsigned int x) { return static_cast<unsigned int> (std::log (static_cast<double> (x))); }
00286   static inline unsigned int log10(unsigned int x) { return static_cast<unsigned int> (std::log10 (static_cast<double> (x))); }
00287 };
00288 
00289 
00290 template<>
00291 struct ScalarTraits<long int>
00292 {
00293   typedef long int magnitudeType;
00294   typedef long int halfPrecision;
00295   typedef long int doublePrecision;
00296   static const bool isComplex = false;
00297   static const bool isOrdinal = true;
00298   static const bool isComparable = true;
00299   static const bool hasMachineParameters = false;
00300   // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
00301   static inline magnitudeType magnitude(long int a) { return static_cast<long int>(std::fabs(static_cast<double>(a))); }
00302   static inline long int zero()  { return 0; }
00303   static inline long int one()   { return 1; }
00304   static inline long int conjugate(long int x) { return x; }
00305   static inline long int real(long int x) { return x; }
00306   static inline long int imag(long int) { return 0; }
00307   static inline bool isnaninf(long int) { return false; }
00308   static inline void seedrandom(unsigned int s) { 
00309     std::srand(s); 
00310 #ifdef __APPLE__
00311     // throw away first random number to address bug 3655
00312     // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
00313     random();
00314 #endif
00315   }
00316   //static inline int random() { return (-1 + 2*rand()); }  // RAB: This version should be used to be consistent with others
00317   static inline long int random() { return std::rand(); }             // RAB: This version should be used for an unsigned int, not int
00318   static inline std::string name() { return "long int"; }
00319   static inline long int squareroot(long int x) { return (long int) std::sqrt((double) x); }
00320   static inline long int pow(long int x, long int y) { return (long int) std::pow((double)x,(double)y); }
00321   // Note: Depending on the number of bits in long int, the cast from
00322   // long int to double may not be exact.
00323   static inline long int log(long int x) { return static_cast<long int> (std::log (static_cast<double> (x))); }
00324   static inline long int log10(long int x) { return static_cast<long int> (std::log10 (static_cast<double> (x))); }
00325 };
00326 
00327 
00328 template<>
00329 struct ScalarTraits<long unsigned int>
00330 {
00331   typedef long unsigned int magnitudeType;
00332   typedef long unsigned int halfPrecision;
00333   typedef long unsigned int doublePrecision;
00334   static const bool isComplex = false;
00335   static const bool isOrdinal = true;
00336   static const bool isComparable = true;
00337   static const bool hasMachineParameters = false;
00338   // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
00339   static inline magnitudeType magnitude(long unsigned int a) { return static_cast<long unsigned int>(std::fabs(static_cast<double>(a))); }
00340   static inline long unsigned int zero()  { return 0; }
00341   static inline long unsigned int one()   { return 1; }
00342   static inline long unsigned int conjugate(long unsigned int x) { return x; }
00343   static inline long unsigned int real(long unsigned int x) { return x; }
00344   static inline long unsigned int imag(long unsigned int) { return 0; }
00345   static inline bool isnaninf(long unsigned int) { return false; }
00346   static inline void seedrandom(unsigned int s) { 
00347     std::srand(s); 
00348 #ifdef __APPLE__
00349     // throw away first random number to address bug 3655
00350     // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
00351     random();
00352 #endif
00353   }
00354   //static inline int random() { return (-1 + 2*rand()); }  // RAB: This version should be used to be consistent with others
00355   static inline long unsigned int random() { return std::rand(); }             // RAB: This version should be used for an unsigned int, not int
00356   static inline std::string name() { return "long unsigned int"; }
00357   static inline long unsigned int squareroot(long unsigned int x) { return (long unsigned int) std::sqrt((double) x); }
00358   static inline long unsigned int pow(long unsigned int x, long unsigned int y) { return (long unsigned int) std::pow((double)x,(double)y); }
00359   // Note: Depending on the number of bits in long unsigned int, the
00360   // cast from long unsigned int to double may not be exact.
00361   static inline long unsigned int log(long unsigned int x) { return static_cast<long unsigned int> (std::log (static_cast<double> (x))); }
00362   static inline long unsigned int log10(long unsigned int x) { return static_cast<long unsigned int> (std::log10 (static_cast<double> (x))); }
00363 };
00364 
00365 
00366 #ifdef HAVE_TEUCHOS_LONG_LONG_INT
00367 template<>
00368 struct ScalarTraits<long long int>
00369 {
00370   typedef long long int magnitudeType;
00371   typedef long long int halfPrecision;
00372   typedef long long int doublePrecision;
00373   static const bool isComplex = false;
00374   static const bool isOrdinal = true;
00375   static const bool isComparable = true;
00376   static const bool hasMachineParameters = false;
00377   // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
00378   static inline magnitudeType magnitude(long long int a) { return static_cast<long long int>(std::fabs(static_cast<double>(a))); }
00379   static inline long long int zero()  { return 0; }
00380   static inline long long int one()   { return 1; }
00381   static inline long long int conjugate(long long int x) { return x; }
00382   static inline long long int real(long long int x) { return x; }
00383   static inline long long int imag(long long int) { return 0; }
00384   static inline bool isnaninf(long long int) { return false; }
00385   static inline void seedrandom(unsigned int s) { 
00386     std::srand(s); 
00387 #ifdef __APPLE__
00388     // throw away first random number to address bug 3655
00389     // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
00390     random();
00391 #endif
00392   }
00393   //static inline int random() { return (-1 + 2*rand()); }  // RAB: This version should be used to be consistent with others
00394   static inline long long int random() { return std::rand(); }             // RAB: This version should be used for an unsigned int, not int
00395   static inline std::string name() { return "long long int"; }
00396   static inline long long int squareroot(long long int x) { return (long long int) std::sqrt((double) x); }
00397   static inline long long int pow(long long int x, long long int y) { return (long long int) std::pow((double)x,(double)y); }
00398   // Note: Depending on the number of bits in long long int, the cast
00399   // from long long int to double may not be exact.
00400   static inline long long int log(long long int x) { return static_cast<long long int> (std::log (static_cast<double> (x))); }
00401   static inline long long int log10(long long int x) { return static_cast<long long int> (std::log10 (static_cast<double> (x))); }
00402 };
00403 
00404 template<>
00405 struct ScalarTraits<unsigned long long int>
00406 {
00407   typedef unsigned long long int magnitudeType;
00408   typedef unsigned long long int halfPrecision;
00409   typedef unsigned long long int doublePrecision;
00410   static const bool isComplex = false;
00411   static const bool isOrdinal = true;
00412   static const bool isComparable = true;
00413   static const bool hasMachineParameters = false;
00414   // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
00415   static inline magnitudeType magnitude(unsigned long long int a) { return static_cast<unsigned long long int>(std::fabs(static_cast<double>(a))); }
00416   static inline unsigned long long int zero()  { return 0; }
00417   static inline unsigned long long int one()   { return 1; }
00418   static inline unsigned long long int conjugate(unsigned long long int x) { return x; }
00419   static inline unsigned long long int real(unsigned long long int x) { return x; }
00420   static inline unsigned long long int imag(unsigned long long int) { return 0; }
00421   static inline bool isnaninf(unsigned long long int) { return false; }
00422   static inline void seedrandom(unsigned int s) { 
00423     std::srand(s); 
00424 #ifdef __APPLE__
00425     // throw away first random number to address bug 3655
00426     // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
00427     random();
00428 #endif
00429   }
00430   //static inline int random() { return (-1 + 2*rand()); }  // RAB: This version should be used to be consistent with others
00431   static inline unsigned long long int random() { return std::rand(); }             // RAB: This version should be used for an unsigned int, not int
00432   static inline std::string name() { return "unsigned long long int"; }
00433   static inline unsigned long long int squareroot(unsigned long long int x) { return (unsigned long long int) std::sqrt((double) x); }
00434   static inline unsigned long long int pow(unsigned long long int x, unsigned long long int y) { return (unsigned long long int) std::pow((double)x,(double)y); }
00435   // Note: Depending on the number of bits in unsigned long long int,
00436   // the cast from unsigned long long int to double may not be exact.
00437   static inline unsigned long long int log(unsigned long long int x) { return static_cast<unsigned long long int> (std::log (static_cast<double> (x))); }
00438   static inline unsigned long long int log10(unsigned long long int x) { return static_cast<unsigned long long int> (std::log10 (static_cast<double> (x))); }
00439 };
00440 #endif // HAVE_TEUCHOS_LONG_LONG_INT
00441 
00442 
00443 #ifdef HAVE_TEUCHOS___INT64
00444 
00445 template<>
00446 struct ScalarTraits<__int64>
00447 {
00448   typedef __int64 magnitudeType;
00449   typedef __int64 halfPrecision;
00450   typedef __int64 doublePrecision;
00451   static const bool isComplex = false;
00452   static const bool isOrdinal = true;
00453   static const bool isComparable = true;
00454   static const bool hasMachineParameters = false;
00455   // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
00456   static inline magnitudeType magnitude(__int64 a) { return static_cast<__int64>(std::fabs(static_cast<double>(a))); }
00457   static inline __int64 zero()  { return 0; }
00458   static inline __int64 one()   { return 1; }
00459   static inline __int64 conjugate(__int64 x) { return x; }
00460   static inline __int64 real(__int64 x) { return x; }
00461   static inline __int64 imag(__int64) { return 0; }
00462   static inline void seedrandom(unsigned int s) { 
00463     std::srand(s); 
00464 #ifdef __APPLE__
00465     // throw away first random number to address bug 3655
00466     // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
00467     random();
00468 #endif
00469   }
00470   //static inline int random() { return (-1 + 2*rand()); }  // RAB: This version should be used to be consistent with others
00471   static inline __int64 random() { return std::rand(); }             // RAB: This version should be used for an unsigned int, not int
00472   static inline std::string name() { return "__int64"; }
00473   static inline __int64 squareroot(__int64 x) { return (__int64) std::sqrt((double) x); }
00474   static inline __int64 pow(__int64 x, __int64 y) { return (__int64) std::pow((double)x,(double)y); }
00475   // Note: Depending on the number of bits in __int64, the cast
00476   // from __int64 to double may not be exact.
00477   static inline __int64 log(__int64 x) { return static_cast<__int64> (std::log (static_cast<double> (x))); }
00478   static inline __int64 log10(__int64 x) { return static_cast<__int64> (std::log10 (static_cast<double> (x))); }
00479 };
00480 
00481 template<>
00482 struct ScalarTraits<unsigned __int64>
00483 {
00484   typedef unsigned __int64 magnitudeType;
00485   typedef unsigned __int64 halfPrecision;
00486   typedef unsigned __int64 doublePrecision;
00487   static const bool isComplex = false;
00488   static const bool isOrdinal = true;
00489   static const bool isComparable = true;
00490   static const bool hasMachineParameters = false;
00491   // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
00492   static inline magnitudeType magnitude(unsigned __int64 a) { return static_cast<unsigned __int64>(std::fabs(static_cast<double>(a))); }
00493   static inline unsigned __int64 zero()  { return 0; }
00494   static inline unsigned __int64 one()   { return 1; }
00495   static inline unsigned __int64 conjugate(unsigned __int64 x) { return x; }
00496   static inline unsigned __int64 real(unsigned __int64 x) { return x; }
00497   static inline unsigned __int64 imag(unsigned __int64) { return 0; }
00498   static inline void seedrandom(unsigned int s) { 
00499     std::srand(s); 
00500 #ifdef __APPLE__
00501     // throw away first random number to address bug 3655
00502     // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
00503     random();
00504 #endif
00505   }
00506   //static inline int random() { return (-1 + 2*rand()); }  // RAB: This version should be used to be consistent with others
00507   static inline unsigned __int64 random() { return std::rand(); }             // RAB: This version should be used for an unsigned int, not int
00508   static inline std::string name() { return "unsigned __int64"; }
00509   static inline unsigned __int64 squareroot(unsigned __int64 x) { return (unsigned __int64) std::sqrt((double) x); }
00510   static inline unsigned __int64 pow(unsigned __int64 x, unsigned __int64 y) { return (unsigned __int64) std::pow((double)x,(double)y); }
00511   // Note: Depending on the number of bits in unsigned __int64,
00512   // the cast from unsigned __int64 to double may not be exact.
00513   static inline unsigned __int64 log(unsigned __int64 x) { return static_cast<unsigned __int64> (std::log (static_cast<double> (x))); }
00514   static inline unsigned __int64 log10(unsigned __int64 x) { return static_cast<unsigned __int64> (std::log10 (static_cast<double> (x))); }
00515 };
00516 
00517 #endif // HAVE_TEUCHOS___INT64
00518 
00519 
00520 #ifndef __sun
00521 extern TEUCHOSCORE_LIB_DLL_EXPORT const float flt_nan;
00522 #endif
00523 
00524  
00525 template<>
00526 struct ScalarTraits<float>
00527 {
00528   typedef float magnitudeType;
00529   typedef float halfPrecision; // should become IEEE754-2008 binary16 or fp16 later, perhaps specified at configure according to architectural support
00530   typedef double doublePrecision;
00531   static const bool isComplex = false;
00532   static const bool isOrdinal = false;
00533   static const bool isComparable = true;
00534   static const bool hasMachineParameters = true;
00535   static inline float eps()   {
00536     return std::numeric_limits<float>::epsilon();
00537   }
00538   static inline float sfmin() {
00539     return std::numeric_limits<float>::min();
00540   }
00541   static inline float base()  {
00542     return static_cast<float>(std::numeric_limits<float>::radix);
00543   }
00544   static inline float prec()  {
00545     return eps()*base();
00546   }
00547   static inline float t()     {
00548     return static_cast<float>(std::numeric_limits<float>::digits);
00549   }
00550   static inline float rnd()   {
00551     return ( std::numeric_limits<float>::round_style == std::round_to_nearest ? one() : zero() );
00552   }
00553   static inline float emin()  {
00554     return static_cast<float>(std::numeric_limits<float>::min_exponent);
00555   }
00556   static inline float rmin()  {
00557     return std::numeric_limits<float>::min();
00558   }
00559   static inline float emax()  {
00560     return static_cast<float>(std::numeric_limits<float>::max_exponent);
00561   }
00562   static inline float rmax()  {
00563     return std::numeric_limits<float>::max();
00564   }
00565   static inline magnitudeType magnitude(float a)
00566     {
00567 #ifdef TEUCHOS_DEBUG
00568       TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
00569         a, "Error, the input value to magnitude(...) a = " << a << " can not be NaN!" );
00570 #endif      
00571       return std::fabs(a);
00572     }    
00573   static inline float zero()  { return(0.0f); }
00574   static inline float one()   { return(1.0f); }    
00575   static inline float conjugate(float x)   { return(x); }    
00576   static inline float real(float x) { return x; }
00577   static inline float imag(float) { return zero(); }
00578   static inline float nan() {
00579 #ifdef __sun
00580     return 0.0f/std::sin(0.0f);
00581 #else
00582     return flt_nan;
00583 #endif
00584   }
00585   static inline bool isnaninf(float x) {
00586     return generic_real_isnaninf<float>(x);
00587   }
00588   static inline void seedrandom(unsigned int s) { 
00589     std::srand(s); 
00590 #ifdef __APPLE__
00591     // throw away first random number to address bug 3655
00592     // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
00593     random();
00594 #endif
00595   }
00596   static inline float random() { float rnd = (float) std::rand() / RAND_MAX; return (-1.0f + 2.0f * rnd); }
00597   static inline std::string name() { return "float"; }
00598   static inline float squareroot(float x)
00599     {
00600 #ifdef TEUCHOS_DEBUG
00601       TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
00602         x, "Error, the input value to squareroot(...) x = " << x << " can not be NaN!" );
00603 #endif
00604       errno = 0;
00605       const float rtn = std::sqrt(x);
00606       if (errno)
00607         return nan();
00608       return rtn;
00609     }
00610   static inline float pow(float x, float y) { return std::pow(x,y); }
00611   static inline float log(float x) { return std::log(x); }
00612   static inline float log10(float x) { return std::log10(x); }
00613 };
00614 
00615 
00616 #ifndef __sun
00617 extern TEUCHOSCORE_LIB_DLL_EXPORT const double dbl_nan;
00618 #endif
00619 
00620  
00621 template<>
00622 struct ScalarTraits<double>
00623 {
00624   typedef double magnitudeType;
00625   typedef float halfPrecision;
00626   /* there are different options as to how to double "double"
00627      - QD's DD (if available)
00628      - ARPREC
00629      - GNU MP
00630      - a true hardware quad
00631 
00632      in the shortterm, this should be specified at configure time. I have inserted a configure-time option (--enable-teuchos-double-to-dd) 
00633      which uses QD's DD when available. This must be used alongside --enable-teuchos-qd.
00634    */
00635 #if defined(HAVE_TEUCHOS_DOUBLE_TO_QD)
00636   typedef dd_real doublePrecision;
00637 #elif defined(HAVE_TEUCHOS_DOUBLE_TO_ARPREC)
00638   typedef mp_real doublePrecision;
00639 #else
00640   typedef double doublePrecision;     // don't double "double" in this case
00641 #endif
00642   static const bool isComplex = false;
00643   static const bool isOrdinal = false;
00644   static const bool isComparable = true;
00645   static const bool hasMachineParameters = true;
00646   static inline double eps()   {
00647     return std::numeric_limits<double>::epsilon();
00648   }
00649   static inline double sfmin() {
00650     return std::numeric_limits<double>::min();
00651   }
00652   static inline double base()  {
00653     return std::numeric_limits<double>::radix;
00654   }
00655   static inline double prec()  {
00656     return eps()*base();
00657   }
00658   static inline double t()     {
00659     return std::numeric_limits<double>::digits;
00660   }
00661   static inline double rnd()   {
00662     return ( std::numeric_limits<double>::round_style == std::round_to_nearest ? double(1.0) : double(0.0) );
00663   }
00664   static inline double emin()  {
00665     return std::numeric_limits<double>::min_exponent;
00666   }
00667   static inline double rmin()  {
00668     return std::numeric_limits<double>::min();
00669   }
00670   static inline double emax()  {
00671     return std::numeric_limits<double>::max_exponent;
00672   }
00673   static inline double rmax()  {
00674     return std::numeric_limits<double>::max();
00675   }
00676   static inline magnitudeType magnitude(double a)
00677     {
00678 #ifdef TEUCHOS_DEBUG
00679       TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
00680         a, "Error, the input value to magnitude(...) a = " << a << " can not be NaN!" );
00681 #endif      
00682       return std::fabs(a);
00683     }
00684   static inline double zero()  { return 0.0; }
00685   static inline double one()   { return 1.0; }
00686   static inline double conjugate(double x)   { return(x); }    
00687   static inline double real(double x) { return(x); }
00688   static inline double imag(double) { return(0); }
00689   static inline double nan() {
00690 #ifdef __sun
00691     return 0.0/std::sin(0.0);
00692 #else
00693     return dbl_nan;
00694 #endif
00695   }
00696   static inline bool isnaninf(double x) {
00697     return generic_real_isnaninf<double>(x);
00698   }
00699   static inline void seedrandom(unsigned int s) { 
00700     std::srand(s); 
00701 #ifdef __APPLE__
00702     // throw away first random number to address bug 3655
00703     // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
00704     random();
00705 #endif
00706   }
00707   static inline double random() { double rnd = (double) std::rand() / RAND_MAX; return (double)(-1.0 + 2.0 * rnd); }
00708   static inline std::string name() { return "double"; }
00709   static inline double squareroot(double x)
00710     {
00711 #ifdef TEUCHOS_DEBUG
00712       TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
00713         x, "Error, the input value to squareroot(...) x = " << x << " can not be NaN!" );
00714 #endif      
00715       errno = 0;
00716       const double rtn = std::sqrt(x);
00717       if (errno)
00718         return nan();
00719       return rtn;
00720     }
00721   static inline double pow(double x, double y) { return std::pow(x,y); }
00722   static inline double log(double x) { return std::log(x); }
00723   static inline double log10(double x) { return std::log10(x); }
00724 };
00725 
00726 
00727 #ifdef HAVE_TEUCHOS_QD
00728 
00729 
00730 bool operator&&(const dd_real &a, const dd_real &b);
00731 bool operator&&(const qd_real &a, const qd_real &b);
00732 
00733 
00734 template<>
00735 struct ScalarTraits<dd_real>
00736 {
00737   typedef dd_real magnitudeType;
00738   typedef double halfPrecision;
00739   typedef qd_real doublePrecision;
00740   static const bool isComplex = false;
00741   static const bool isOrdinal = false;
00742   static const bool isComparable = true;
00743   static const bool hasMachineParameters = true;
00744   static inline dd_real eps()   { return std::numeric_limits<dd_real>::epsilon(); }
00745   static inline dd_real sfmin() { return std::numeric_limits<dd_real>::min(); }
00746   static inline dd_real base()  { return std::numeric_limits<dd_real>::radix; }
00747   static inline dd_real prec()  { return eps()*base(); }
00748   static inline dd_real t()     { return std::numeric_limits<dd_real>::digits; }
00749   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) ); }
00750   static inline dd_real emin()  { return std::numeric_limits<dd_real>::min_exponent; }
00751   static inline dd_real rmin()  { return std::numeric_limits<dd_real>::min(); }
00752   static inline dd_real emax()  { return std::numeric_limits<dd_real>::max_exponent; }
00753   static inline dd_real rmax()  { return std::numeric_limits<dd_real>::max(); }
00754   static inline magnitudeType magnitude(dd_real a)
00755   {
00756 #ifdef TEUCHOS_DEBUG
00757       TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
00758         a, "Error, the input value to magnitude(...) a = " << a << " can not be NaN!" );
00759 #endif      
00760       return ::abs(a);
00761   }
00762   static inline dd_real zero()  { return dd_real(0.0); }
00763   static inline dd_real one()   { return dd_real(1.0); }
00764   static inline dd_real conjugate(dd_real x)   { return(x); }    
00765   static inline dd_real real(dd_real x) { return x ; }
00766   static inline dd_real imag(dd_real) { return zero(); }
00767   static inline dd_real nan() { return dd_real::_nan; }
00768   static inline bool isnaninf(dd_real x) { return isnan(x) || isinf(x); }
00769   static inline void seedrandom(unsigned int s) {
00770     // ddrand() uses std::rand(), so the std::srand() is our seed
00771     std::srand(s); 
00772 #ifdef __APPLE__
00773     // throw away first random number to address bug 3655
00774     // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
00775     random();
00776 #endif
00777   }
00778   static inline dd_real random() { return ddrand(); }
00779   static inline std::string name() { return "dd_real"; }
00780   static inline dd_real squareroot(dd_real x)
00781   {
00782 #ifdef TEUCHOS_DEBUG
00783       TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
00784         x, "Error, the input value to squareroot(...) x = " << x << " can not be NaN!" );
00785 #endif      
00786       return ::sqrt(x);
00787   }
00788   static inline dd_real pow(dd_real x, dd_real y) { return ::pow(x,y); }
00789   // dd_real puts its transcendental functions in the global namespace.
00790   static inline dd_real log(dd_real x) { return ::log(x); }
00791   static inline dd_real log10(dd_real x) { return ::log10(x); }
00792 };
00793 
00794 
00795 template<>
00796 struct ScalarTraits<qd_real>
00797 {
00798   typedef qd_real magnitudeType;
00799   typedef dd_real halfPrecision;
00800   typedef qd_real doublePrecision;
00801   static const bool isComplex = false;
00802   static const bool isOrdinal = false;
00803   static const bool isComparable = true;
00804   static const bool hasMachineParameters = true;
00805   static inline qd_real eps()   { return std::numeric_limits<qd_real>::epsilon(); }
00806   static inline qd_real sfmin() { return std::numeric_limits<qd_real>::min(); }
00807   static inline qd_real base()  { return std::numeric_limits<qd_real>::radix; }
00808   static inline qd_real prec()  { return eps()*base(); }
00809   static inline qd_real t()     { return std::numeric_limits<qd_real>::digits; }
00810   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) ); }
00811   static inline qd_real emin()  { return std::numeric_limits<qd_real>::min_exponent; }
00812   static inline qd_real rmin()  { return std::numeric_limits<qd_real>::min(); }
00813   static inline qd_real emax()  { return std::numeric_limits<qd_real>::max_exponent; }
00814   static inline qd_real rmax()  { return std::numeric_limits<qd_real>::max(); }
00815   static inline magnitudeType magnitude(qd_real a)
00816   {
00817 #ifdef TEUCHOS_DEBUG
00818       TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
00819         a, "Error, the input value to magnitude(...) a = " << a << " can not be NaN!" );
00820 #endif      
00821       return ::abs(a);
00822   }
00823   static inline qd_real zero()  { return qd_real(0.0); }
00824   static inline qd_real one()   { return qd_real(1.0); }
00825   static inline qd_real conjugate(qd_real x)   { return(x); }    
00826   static inline qd_real real(qd_real x) { return x ; }
00827   static inline qd_real imag(qd_real) { return zero(); }
00828   static inline qd_real nan() { return qd_real::_nan; }
00829   static inline bool isnaninf(qd_real x) { return isnan(x) || isinf(x); }
00830   static inline void seedrandom(unsigned int s) {
00831     // qdrand() uses std::rand(), so the std::srand() is our seed
00832     std::srand(s); 
00833 #ifdef __APPLE__
00834     // throw away first random number to address bug 3655
00835     // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
00836     random();
00837 #endif
00838   }
00839   static inline qd_real random() { return qdrand(); }
00840   static inline std::string name() { return "qd_real"; }
00841   static inline qd_real squareroot(qd_real x)
00842   {
00843 #ifdef TEUCHOS_DEBUG
00844       TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
00845         x, "Error, the input value to squareroot(...) x = " << x << " can not be NaN!" );
00846 #endif      
00847       return ::sqrt(x);
00848   }
00849   static inline qd_real pow(qd_real x, qd_real y) { return ::pow(x,y); }
00850   // qd_real puts its transcendental functions in the global namespace.
00851   static inline qd_real log(qd_real x) { return ::log(x); }
00852   static inline qd_real log10(qd_real x) { return ::log10(x); }
00853 };
00854 
00855 
00856 #endif  // HAVE_TEUCHOS_QD
00857 
00858 
00859 #ifdef HAVE_TEUCHOS_GNU_MP
00860 
00861 
00862 extern gmp_randclass gmp_rng; 
00863 
00864 
00865 /* Regarding halfPrecision, doublePrecision and mpf_class: 
00866    Because the precision of an mpf_class float is not determined by the data type, 
00867    there is no way to fill the typedefs for this object. 
00868 
00869    Instead, we could create new data classes (e.g., Teuchos::MPF128, Teuchos::MPF256) for 
00870    commonly used levels of precision, and fill out ScalarTraits for these. This would allow us
00871    to typedef the promotions and demotions in the appropriate way. These classes would serve to 
00872    wrap an mpf_class object, calling the constructor for the appropriate precision, exposing the 
00873    arithmetic routines but hiding the precision-altering routines.
00874    
00875    Alternatively (perhaps, preferably), would could create a single class templated on the precision (e.g., Teuchos::MPF<N>). 
00876    Then we have a single (partially-specialized) implementation of ScalarTraits. This class, as above, must expose all of the 
00877    operations expected of a scalar type; however, most of these can be trivially stolen from the gmpcxx.h header file
00878 
00879    CGB/RAB, 01/05/2009
00880 */
00881 template<>
00882 struct ScalarTraits<mpf_class>
00883 {
00884   typedef mpf_class magnitudeType;
00885   typedef mpf_class halfPrecision;
00886   typedef mpf_class doublePrecision;
00887   static const bool isComplex = false;
00888   static const bool hasMachineParameters = false;
00889   // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
00890   static magnitudeType magnitude(mpf_class a) { return std::abs(a); }
00891   static inline mpf_class zero() { mpf_class zero = 0.0; return zero; }
00892   static inline mpf_class one() { mpf_class one = 1.0; return one; }    
00893   static inline mpf_class conjugate(mpf_class x) { return x; }
00894   static inline mpf_class real(mpf_class x) { return(x); }
00895   static inline mpf_class imag(mpf_class x) { return(0); }
00896   static inline bool isnaninf(mpf_class x) { return false; } // mpf_class currently can't handle nan or inf!
00897   static inline void seedrandom(unsigned int s) { 
00898     unsigned long int seedVal = static_cast<unsigned long int>(s);
00899     gmp_rng.seed( seedVal );  
00900   }
00901   static inline mpf_class random() { 
00902     return gmp_rng.get_f(); 
00903   }
00904   static inline std::string name() { return "mpf_class"; }
00905   static inline mpf_class squareroot(mpf_class x) { return std::sqrt(x); }
00906   static inline mpf_class pow(mpf_class x, mpf_class y) { return pow(x,y); }
00907   // Todo: RAB: 2004/05/28: Add nan() and isnaninf() functions when needed!
00908 };
00909 
00910 #endif  // HAVE_TEUCHOS_GNU_MP
00911 
00912 #ifdef HAVE_TEUCHOS_ARPREC
00913 
00914 /* See discussion above for mpf_class, regarding halfPrecision and doublePrecision. Something similar will need to be done
00915    for ARPREC. */
00916 template<>
00917 struct ScalarTraits<mp_real>
00918 {
00919   typedef mp_real magnitudeType;
00920   typedef double halfPrecision;
00921   typedef mp_real doublePrecision;
00922   static const bool isComplex = false;
00923   static const bool isComparable = true;
00924   static const bool isOrdinal = false;
00925   static const bool hasMachineParameters = false;
00926   // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
00927   static magnitudeType magnitude(mp_real a) { return abs(a); }
00928   static inline mp_real zero() { mp_real zero = 0.0; return zero; }
00929   static inline mp_real one() { mp_real one = 1.0; return one; }    
00930   static inline mp_real conjugate(mp_real x) { return x; }
00931   static inline mp_real real(mp_real x) { return(x); }
00932   static inline mp_real imag(mp_real x) { return zero(); }
00933   static inline bool isnaninf(mp_real x) { return false; } // ToDo: Change this?
00934   static inline void seedrandom(unsigned int s) { 
00935     long int seedVal = static_cast<long int>(s);
00936     srand48(seedVal);
00937   }
00938   static inline mp_real random() { return mp_rand(); }
00939   static inline std::string name() { return "mp_real"; }
00940   static inline mp_real squareroot(mp_real x) { return sqrt(x); }
00941   static inline mp_real pow(mp_real x, mp_real y) { return pow(x,y); }
00942   // Todo: RAB: 2004/05/28: Add nan() and isnaninf() functions when needed!
00943 };
00944 
00945   
00946 #endif // HAVE_TEUCHOS_ARPREC
00947 
00948  
00949 #ifdef HAVE_TEUCHOS_COMPLEX
00950 
00951 
00952 // Partial specialization for std::complex numbers templated on real type T
00953 template<class T> 
00954 struct ScalarTraits<
00955   std::complex<T>
00956 >
00957 {
00958   typedef std::complex<T>  ComplexT;
00959   typedef std::complex<typename ScalarTraits<T>::halfPrecision> halfPrecision;
00960   typedef std::complex<typename ScalarTraits<T>::doublePrecision> doublePrecision;
00961   typedef typename ScalarTraits<T>::magnitudeType magnitudeType;
00962   static const bool isComplex = true;
00963   static const bool isOrdinal = ScalarTraits<T>::isOrdinal;
00964   static const bool isComparable = false;
00965   static const bool hasMachineParameters = true;
00966   static inline magnitudeType eps()          { return ScalarTraits<magnitudeType>::eps(); }
00967   static inline magnitudeType sfmin()        { return ScalarTraits<magnitudeType>::sfmin(); }
00968   static inline magnitudeType base()         { return ScalarTraits<magnitudeType>::base(); }
00969   static inline magnitudeType prec()         { return ScalarTraits<magnitudeType>::prec(); }
00970   static inline magnitudeType t()            { return ScalarTraits<magnitudeType>::t(); }
00971   static inline magnitudeType rnd()          { return ScalarTraits<magnitudeType>::rnd(); }
00972   static inline magnitudeType emin()         { return ScalarTraits<magnitudeType>::emin(); }
00973   static inline magnitudeType rmin()         { return ScalarTraits<magnitudeType>::rmin(); }
00974   static inline magnitudeType emax()         { return ScalarTraits<magnitudeType>::emax(); }
00975   static inline magnitudeType rmax()         { return ScalarTraits<magnitudeType>::rmax(); }
00976   static magnitudeType magnitude(ComplexT a)
00977     {
00978 #ifdef TEUCHOS_DEBUG
00979       TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
00980         a, "Error, the input value to magnitude(...) a = " << a << " can not be NaN!" );
00981 #endif      
00982       return std::abs(a);
00983     }
00984   static inline ComplexT zero()              { return ComplexT(ScalarTraits<magnitudeType>::zero(),ScalarTraits<magnitudeType>::zero()); }
00985   static inline ComplexT one()               { return ComplexT(ScalarTraits<magnitudeType>::one(),ScalarTraits<magnitudeType>::zero()); }
00986   static inline ComplexT conjugate(ComplexT a){ return ComplexT(a.real(),-a.imag()); }
00987   static inline magnitudeType real(ComplexT a) { return a.real(); }
00988   static inline magnitudeType imag(ComplexT a) { return a.imag(); }
00989   static inline ComplexT nan()               { return ComplexT(ScalarTraits<magnitudeType>::nan(),ScalarTraits<magnitudeType>::nan()); }
00990   static inline bool isnaninf(ComplexT x)    { return ScalarTraits<magnitudeType>::isnaninf(x.real()) || ScalarTraits<magnitudeType>::isnaninf(x.imag()); }
00991   static inline void seedrandom(unsigned int s) { ScalarTraits<magnitudeType>::seedrandom(s); }
00992   static inline ComplexT random()
00993     {
00994       const T rnd1 = ScalarTraits<magnitudeType>::random();
00995       const T rnd2 = ScalarTraits<magnitudeType>::random();
00996       return ComplexT(rnd1,rnd2);
00997     }
00998   static inline std::string name() { return std::string("std::complex<")+std::string(ScalarTraits<magnitudeType>::name())+std::string(">"); }
00999   // This will only return one of the square roots of x, the other can be obtained by taking its conjugate
01000   static inline ComplexT squareroot(ComplexT x)
01001     {
01002 #ifdef TEUCHOS_DEBUG
01003       TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
01004         x, "Error, the input value to squareroot(...) x = " << x << " can not be NaN!" );
01005 #endif
01006       typedef ScalarTraits<magnitudeType>  STMT;
01007       const T r  = x.real(), i = x.imag(), zero = STMT::zero(), two = 2.0;
01008       const T a  = STMT::squareroot((r*r)+(i*i));
01009       const T nr = STMT::squareroot((a+r)/two);
01010       const T ni = ( i == zero ? zero : STMT::squareroot((a-r)/two) );
01011       return ComplexT(nr,ni);
01012     }
01013     // 2010/03/19: rabartl: Above, I had to add the check for i == zero
01014     // to avoid a returned NaN on the Intel 10.1 compiler.  For some
01015     // reason, having these two squareroot calls in a row produce a NaN
01016     // in an optimized build with this compiler.  Amazingly, when I put
01017     // in print statements (i.e. std::cout << ...) the NaN went away and
01018     // the second squareroot((a-r)/two) returned zero correctly.  I
01019     // guess that calling the output routine flushed the registers or
01020     // something and restarted the squareroot rountine for this compiler
01021     // or something similar.  Actually, due to roundoff, it is possible that a-r
01022     // might be slightly less than zero (i.e. -1e-16) and that would cause
01023     // a possbile NaN return.  THe above if test is the right thing to do
01024     // I think and is very cheap.
01025   static inline ComplexT pow(ComplexT x, ComplexT y) { return pow(x,y); }
01026 };
01027 
01028 
01029 #endif //  HAVE_TEUCHOS_COMPLEX
01030 
01031 
01032 #endif // DOXYGEN_SHOULD_SKIP_THIS
01033 
01034 
01035 } // Teuchos namespace
01036 
01037 
01038 #endif // _TEUCHOS_SCALARTRAITS_HPP_
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines