Teuchos - Trilinos Tools Package  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Teuchos_ScalarTraits.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Teuchos: Common Tools Package
5 // Copyright (2004) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
42 // NOTE: Before adding specializations of ScalarTraits, make sure that they do
43 // not duplicate specializations already present in PyTrilinos (see
44 // packages/PyTrilinos/src/Teuchos_Traits.i)
45 
46 // NOTE: halfPrecision and doublePrecision are not currently implemented for
47 // ARPREC, GMP or the ordinal types (e.g., int, char)
48 
49 #ifndef _TEUCHOS_SCALARTRAITS_HPP_
50 #define _TEUCHOS_SCALARTRAITS_HPP_
51 
56 #include "Teuchos_ConfigDefs.hpp"
57 
58 #ifdef HAVE_TEUCHOS_ARPREC
59 #include <arprec/mp_real.h>
60 #endif
61 
62 #ifdef HAVE_TEUCHOSCORE_QUADMATH
63 #include <quadmath.h>
64 
65 // Teuchos_ConfigDefs.hpp includes <iostream>, which includes
66 // <ostream>. If this ever changes, include <ostream> here.
67 
68 namespace std {
69 
79 ostream&
80 operator<< (std::ostream& out, const __float128& x);
81 
91 istream&
92 operator>> (std::istream& in, __float128& x);
93 
94 } // namespace std
95 
96 #endif // HAVE_TEUCHOSCORE_QUADMATH
97 
98 #ifdef HAVE_TEUCHOS_QD
99 #include <qd/qd_real.h>
100 #include <qd/dd_real.h>
101 #endif
102 
103 #ifdef HAVE_TEUCHOS_GNU_MP
104 #include <gmp.h>
105 #include <gmpxx.h>
106 #endif
107 
108 
109 #include "Teuchos_ScalarTraitsDecl.hpp"
110 
111 
112 namespace Teuchos {
113 
114 
115 #ifndef DOXYGEN_SHOULD_SKIP_THIS
116 
117 
118 TEUCHOSCORE_LIB_DLL_EXPORT
119 void throwScalarTraitsNanInfError( const std::string &errMsg );
120 
121 
122 template<class Scalar>
123 bool generic_real_isnaninf(const Scalar &x)
124 {
125 #ifdef HAVE_TEUCHOSCORE_CXX11
126  if (std::isnan(x)) return true;
127  if (std::isinf(x)) return true;
128  return false;
129 #else
130  typedef std::numeric_limits<Scalar> STD_NL;
131  // IEEE says this should fail for NaN (not all compilers do not obey IEEE)
132  const Scalar tol = 1.0; // Any (bounded) number should do!
133  if (!(x <= tol) && !(x > tol)) return true;
134  // Use fact that Inf*0 = NaN (again, all compilers do not obey IEEE)
135  Scalar z = static_cast<Scalar>(0.0) * x;
136  if (!(z <= tol) && !(z > tol)) return true;
137  // As a last result use comparisons
138  if (x == STD_NL::infinity() || x == -STD_NL::infinity()) return true;
139  // We give up and assume the number is finite
140  return false;
141 #endif
142 }
143 
144 
145 #define TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR( VALUE, MSG ) \
146  if (isnaninf(VALUE)) { \
147  std::ostringstream omsg; \
148  omsg << MSG; \
149  Teuchos::throwScalarTraitsNanInfError(omsg.str()); \
150  }
151 
152 
153 template<>
154 struct ScalarTraits<char>
155 {
156  typedef char magnitudeType;
157  typedef char halfPrecision;
158  typedef char doublePrecision;
159  typedef char coordinateType;
160  static const bool isComplex = false;
161  static const bool isOrdinal = true;
162  static const bool isComparable = true;
163  static const bool hasMachineParameters = false;
164  // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
165  static inline magnitudeType magnitude(char a) { return static_cast<char>(std::fabs(static_cast<double>(a))); }
166  static inline char zero() { return 0; }
167  static inline char one() { return 1; }
168  static inline char conjugate(char x) { return x; }
169  static inline char real(char x) { return x; }
170  static inline char imag(char) { return 0; }
171  static inline bool isnaninf(char ) { return false; }
172  static inline void seedrandom(unsigned int s) {
173  std::srand(s);
174 #ifdef __APPLE__
175  // throw away first random number to address bug 3655
176  // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
177  random();
178 #endif
179  }
180  //static inline char random() { return (-1 + 2*rand()); } // RAB: This version should be used to be consistent with others
181  static inline char random() { return std::rand(); } // RAB: This version should be used for an unsigned char, not char
182  static inline std::string name() { return "char"; }
183  static inline char squareroot(char x) { return (char) std::sqrt((double) x); }
184  static inline char pow(char x, char y) { return (char) std::pow((double)x,(double)y); }
185  static inline char log(char x) { return static_cast<char> (std::log (static_cast<double> (x))); }
186  static inline char log10(char x) { return static_cast<char> (std::log10 (static_cast<double> (x))); }
187 };
188 
189 
190 template<>
191 struct ScalarTraits<short int>
192 {
193  typedef short int magnitudeType;
194  typedef short int halfPrecision;
195  typedef short int doublePrecision;
196  typedef short int coordinateType;
197  static const bool isComplex = false;
198  static const bool isOrdinal = true;
199  static const bool isComparable = true;
200  static const bool hasMachineParameters = false;
201  // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
202  static inline magnitudeType magnitude(short int a) { return static_cast<short int>(std::fabs(static_cast<double>(a))); }
203  static inline short int zero() { return 0; }
204  static inline short int one() { return 1; }
205  static inline short int conjugate(short int x) { return x; }
206  static inline short int real(short int x) { return x; }
207  static inline short int imag(short int) { return 0; }
208  static inline bool isnaninf(short int) { return false; }
209  static inline void seedrandom(unsigned int s) {
210  std::srand(s);
211 #ifdef __APPLE__
212  // throw away first random number to address bug 3655
213  // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
214  random();
215 #endif
216  }
217  //static inline int random() { return (-1 + 2*rand()); } // RAB: This version should be used to be consistent with others
218  static inline short int random() { return std::rand(); } // RAB: This version should be used for an unsigned int, not int
219  static inline std::string name() { return "short int"; }
220  static inline short int squareroot(short int x) { return (short int) std::sqrt((double) x); }
221  static inline short int pow(short int x, short int y) { return (short int) std::pow((double)x,(double)y); }
222  static inline short int log(short int x) { return static_cast<short int> (std::log (static_cast<double> (x))); }
223  static inline short int log10(short int x) { return static_cast<short int> (std::log10 (static_cast<double> (x))); }
224 };
225 
226 template<>
227 struct ScalarTraits<unsigned short int>
228 {
229  typedef unsigned short int magnitudeType;
230  typedef unsigned short int halfPrecision;
231  typedef unsigned short int doublePrecision;
232  typedef unsigned short int coordinateType;
233  static const bool isComplex = false;
234  static const bool isOrdinal = true;
235  static const bool isComparable = true;
236  static const bool hasMachineParameters = false;
237  // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
238  static inline magnitudeType magnitude(unsigned short int a) { return static_cast<unsigned short int>(std::fabs(static_cast<double>(a))); }
239  static inline unsigned short int zero() { return 0; }
240  static inline unsigned short int one() { return 1; }
241  static inline unsigned short int conjugate(unsigned short int x) { return x; }
242  static inline unsigned short int real(unsigned short int x) { return x; }
243  static inline unsigned short int imag(unsigned short int) { return 0; }
244  static inline bool isnaninf(unsigned short int) { return false; }
245  static inline void seedrandom(unsigned int s) {
246  std::srand(s);
247 #ifdef __APPLE__
248  // throw away first random number to address bug 3655
249  // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
250  random();
251 #endif
252  }
253  //static inline int random() { return (-1 + 2*rand()); } // RAB: This version should be used to be consistent with others
254  static inline unsigned short int random() { return std::rand(); } // RAB: This version should be used for an unsigned int, not int
255  static inline std::string name() { return "unsigned short int"; }
256  static inline unsigned short int squareroot(unsigned short int x) { return (unsigned short int) std::sqrt((double) x); }
257  static inline unsigned short int pow(unsigned short int x, unsigned short int y) { return (unsigned short int) std::pow((double)x,(double)y); }
258  static inline unsigned short int log(unsigned short int x) { return static_cast<unsigned short int> (std::log (static_cast<double> (x))); }
259  static inline unsigned short int log10(unsigned short int x) { return static_cast<unsigned short int> (std::log10 (static_cast<double> (x))); }
260 };
261 
262 
263 template<>
264 struct ScalarTraits<int>
265 {
266  typedef int magnitudeType;
267  typedef int halfPrecision;
268  typedef int doublePrecision;
269  typedef int coordinateType;
270  static const bool isComplex = false;
271  static const bool isOrdinal = true;
272  static const bool isComparable = true;
273  static const bool hasMachineParameters = false;
274  // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
275  static inline magnitudeType magnitude(int a) { return static_cast<int>(std::fabs(static_cast<double>(a))); }
276  static inline int zero() { return 0; }
277  static inline int one() { return 1; }
278  static inline int conjugate(int x) { return x; }
279  static inline int real(int x) { return x; }
280  static inline int imag(int) { return 0; }
281  static inline bool isnaninf(int) { return false; }
282  static inline void seedrandom(unsigned int s) {
283  std::srand(s);
284 #ifdef __APPLE__
285  // throw away first random number to address bug 3655
286  // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
287  random();
288 #endif
289  }
290  //static inline int random() { return (-1 + 2*rand()); } // RAB: This version should be used to be consistent with others
291  static inline int random() { return std::rand(); } // RAB: This version should be used for an unsigned int, not int
292  static inline std::string name() { return "int"; }
293  static inline int squareroot(int x) { return (int) std::sqrt((double) x); }
294  static inline int pow(int x, int y) { return (int) std::pow((double)x,(double)y); }
295  static inline int log(int x) { return static_cast<int> (std::log (static_cast<double> (x))); }
296  static inline int log10(int x) { return static_cast<int> (std::log10 (static_cast<double> (x))); }
297 };
298 
299 
300 template<>
301 struct ScalarTraits<unsigned int>
302 {
303  typedef unsigned int magnitudeType;
304  typedef unsigned int halfPrecision;
305  typedef unsigned int doublePrecision;
306  typedef unsigned int coordinateType;
307  static const bool isComplex = false;
308  static const bool isOrdinal = true;
309  static const bool isComparable = true;
310  static const bool hasMachineParameters = false;
311  // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
312  static inline magnitudeType magnitude(unsigned int a) { return static_cast<unsigned int>(std::fabs(static_cast<double>(a))); }
313  static inline unsigned int zero() { return 0; }
314  static inline unsigned int one() { return 1; }
315  static inline unsigned int conjugate(unsigned int x) { return x; }
316  static inline unsigned int real(unsigned int x) { return x; }
317  static inline unsigned int imag(unsigned int) { return 0; }
318  static inline bool isnaninf(unsigned int) { return false; }
319  static inline void seedrandom(unsigned int s) {
320  std::srand(s);
321 #ifdef __APPLE__
322  // throw away first random number to address bug 3655
323  // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
324  random();
325 #endif
326  }
327  //static inline int random() { return (-1 + 2*rand()); } // RAB: This version should be used to be consistent with others
328  static inline unsigned int random() { return std::rand(); } // RAB: This version should be used for an unsigned int, not int
329  static inline std::string name() { return "unsigned int"; }
330  static inline unsigned int squareroot(unsigned int x) { return (unsigned int) std::sqrt((double) x); }
331  static inline unsigned int pow(unsigned int x, unsigned int y) { return (unsigned int) std::pow((double)x,(double)y); }
332  static inline unsigned int log(unsigned int x) { return static_cast<unsigned int> (std::log (static_cast<double> (x))); }
333  static inline unsigned int log10(unsigned int x) { return static_cast<unsigned int> (std::log10 (static_cast<double> (x))); }
334 };
335 
336 
337 template<>
338 struct ScalarTraits<long int>
339 {
340  typedef long int magnitudeType;
341  typedef long int halfPrecision;
342  typedef long int doublePrecision;
343  typedef long int coordinateType;
344  static const bool isComplex = false;
345  static const bool isOrdinal = true;
346  static const bool isComparable = true;
347  static const bool hasMachineParameters = false;
348  // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
349  static inline magnitudeType magnitude(long int a) { return static_cast<long int>(std::fabs(static_cast<double>(a))); }
350  static inline long int zero() { return 0; }
351  static inline long int one() { return 1; }
352  static inline long int conjugate(long int x) { return x; }
353  static inline long int real(long int x) { return x; }
354  static inline long int imag(long int) { return 0; }
355  static inline bool isnaninf(long int) { return false; }
356  static inline void seedrandom(unsigned int s) {
357  std::srand(s);
358 #ifdef __APPLE__
359  // throw away first random number to address bug 3655
360  // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
361  random();
362 #endif
363  }
364  //static inline int random() { return (-1 + 2*rand()); } // RAB: This version should be used to be consistent with others
365  static inline long int random() { return std::rand(); } // RAB: This version should be used for an unsigned int, not int
366  static inline std::string name() { return "long int"; }
367  static inline long int squareroot(long int x) { return (long int) std::sqrt((double) x); }
368  static inline long int pow(long int x, long int y) { return (long int) std::pow((double)x,(double)y); }
369  // Note: Depending on the number of bits in long int, the cast from
370  // long int to double may not be exact.
371  static inline long int log(long int x) { return static_cast<long int> (std::log (static_cast<double> (x))); }
372  static inline long int log10(long int x) { return static_cast<long int> (std::log10 (static_cast<double> (x))); }
373 };
374 
375 
376 template<>
377 struct ScalarTraits<long unsigned int>
378 {
379  typedef long unsigned int magnitudeType;
380  typedef long unsigned int halfPrecision;
381  typedef long unsigned int doublePrecision;
382  typedef long unsigned int coordinateType;
383  static const bool isComplex = false;
384  static const bool isOrdinal = true;
385  static const bool isComparable = true;
386  static const bool hasMachineParameters = false;
387  // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
388  static inline magnitudeType magnitude(long unsigned int a) { return static_cast<long unsigned int>(std::fabs(static_cast<double>(a))); }
389  static inline long unsigned int zero() { return 0; }
390  static inline long unsigned int one() { return 1; }
391  static inline long unsigned int conjugate(long unsigned int x) { return x; }
392  static inline long unsigned int real(long unsigned int x) { return x; }
393  static inline long unsigned int imag(long unsigned int) { return 0; }
394  static inline bool isnaninf(long unsigned int) { return false; }
395  static inline void seedrandom(unsigned int s) {
396  std::srand(s);
397 #ifdef __APPLE__
398  // throw away first random number to address bug 3655
399  // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
400  random();
401 #endif
402  }
403  //static inline int random() { return (-1 + 2*rand()); } // RAB: This version should be used to be consistent with others
404  static inline long unsigned int random() { return std::rand(); } // RAB: This version should be used for an unsigned int, not int
405  static inline std::string name() { return "long unsigned int"; }
406  static inline long unsigned int squareroot(long unsigned int x) { return (long unsigned int) std::sqrt((double) x); }
407  static inline long unsigned int pow(long unsigned int x, long unsigned int y) { return (long unsigned int) std::pow((double)x,(double)y); }
408  // Note: Depending on the number of bits in long unsigned int, the
409  // cast from long unsigned int to double may not be exact.
410  static inline long unsigned int log(long unsigned int x) { return static_cast<long unsigned int> (std::log (static_cast<double> (x))); }
411  static inline long unsigned int log10(long unsigned int x) { return static_cast<long unsigned int> (std::log10 (static_cast<double> (x))); }
412 };
413 
414 
415 template<>
416 struct ScalarTraits<long long int>
417 {
418  typedef long long int magnitudeType;
419  typedef long long int halfPrecision;
420  typedef long long int doublePrecision;
421  typedef long long int coordinateType;
422  static const bool isComplex = false;
423  static const bool isOrdinal = true;
424  static const bool isComparable = true;
425  static const bool hasMachineParameters = false;
426  // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
427  static inline magnitudeType magnitude(long long int a) { return static_cast<long long int>(std::fabs(static_cast<double>(a))); }
428  static inline long long int zero() { return 0; }
429  static inline long long int one() { return 1; }
430  static inline long long int conjugate(long long int x) { return x; }
431  static inline long long int real(long long int x) { return x; }
432  static inline long long int imag(long long int) { return 0; }
433  static inline bool isnaninf(long long int) { return false; }
434  static inline void seedrandom(unsigned int s) {
435  std::srand(s);
436 #ifdef __APPLE__
437  // throw away first random number to address bug 3655
438  // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
439  random();
440 #endif
441  }
442  //static inline int random() { return (-1 + 2*rand()); } // RAB: This version should be used to be consistent with others
443  static inline long long int random() { return std::rand(); } // RAB: This version should be used for an unsigned int, not int
444  static inline std::string name() { return "long long int"; }
445  static inline long long int squareroot(long long int x) { return (long long int) std::sqrt((double) x); }
446  static inline long long int pow(long long int x, long long int y) { return (long long int) std::pow((double)x,(double)y); }
447  // Note: Depending on the number of bits in long long int, the cast
448  // from long long int to double may not be exact.
449  static inline long long int log(long long int x) { return static_cast<long long int> (std::log (static_cast<double> (x))); }
450  static inline long long int log10(long long int x) { return static_cast<long long int> (std::log10 (static_cast<double> (x))); }
451 };
452 
453 template<>
454 struct ScalarTraits<unsigned long long int>
455 {
456  typedef unsigned long long int magnitudeType;
457  typedef unsigned long long int halfPrecision;
458  typedef unsigned long long int doublePrecision;
459  typedef unsigned long long int coordinateType;
460  static const bool isComplex = false;
461  static const bool isOrdinal = true;
462  static const bool isComparable = true;
463  static const bool hasMachineParameters = false;
464  // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
465  static inline magnitudeType magnitude(unsigned long long int a) { return static_cast<unsigned long long int>(std::fabs(static_cast<double>(a))); }
466  static inline unsigned long long int zero() { return 0; }
467  static inline unsigned long long int one() { return 1; }
468  static inline unsigned long long int conjugate(unsigned long long int x) { return x; }
469  static inline unsigned long long int real(unsigned long long int x) { return x; }
470  static inline unsigned long long int imag(unsigned long long int) { return 0; }
471  static inline bool isnaninf(unsigned long long int) { return false; }
472  static inline void seedrandom(unsigned int s) {
473  std::srand(s);
474 #ifdef __APPLE__
475  // throw away first random number to address bug 3655
476  // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
477  random();
478 #endif
479  }
480  //static inline int random() { return (-1 + 2*rand()); } // RAB: This version should be used to be consistent with others
481  static inline unsigned long long int random() { return std::rand(); } // RAB: This version should be used for an unsigned int, not int
482  static inline std::string name() { return "unsigned long long int"; }
483  static inline unsigned long long int squareroot(unsigned long long int x) { return (unsigned long long int) std::sqrt((double) x); }
484  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); }
485  // Note: Depending on the number of bits in unsigned long long int,
486  // the cast from unsigned long long int to double may not be exact.
487  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))); }
488  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))); }
489 };
490 
491 
492 #ifdef HAVE_TEUCHOS___INT64
493 
494 template<>
495 struct ScalarTraits<__int64>
496 {
497  typedef __int64 magnitudeType;
498  typedef __int64 halfPrecision;
499  typedef __int64 doublePrecision;
500  typedef __int64 coordinateType;
501  static const bool isComplex = false;
502  static const bool isOrdinal = true;
503  static const bool isComparable = true;
504  static const bool hasMachineParameters = false;
505  // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
506  static inline magnitudeType magnitude(__int64 a) { return static_cast<__int64>(std::fabs(static_cast<double>(a))); }
507  static inline __int64 zero() { return 0; }
508  static inline __int64 one() { return 1; }
509  static inline __int64 conjugate(__int64 x) { return x; }
510  static inline __int64 real(__int64 x) { return x; }
511  static inline __int64 imag(__int64) { return 0; }
512  static inline void seedrandom(unsigned int s) {
513  std::srand(s);
514 #ifdef __APPLE__
515  // throw away first random number to address bug 3655
516  // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
517  random();
518 #endif
519  }
520  //static inline int random() { return (-1 + 2*rand()); } // RAB: This version should be used to be consistent with others
521  static inline __int64 random() { return std::rand(); } // RAB: This version should be used for an unsigned int, not int
522  static inline std::string name() { return "__int64"; }
523  static inline __int64 squareroot(__int64 x) { return (__int64) std::sqrt((double) x); }
524  static inline __int64 pow(__int64 x, __int64 y) { return (__int64) std::pow((double)x,(double)y); }
525  // Note: Depending on the number of bits in __int64, the cast
526  // from __int64 to double may not be exact.
527  static inline __int64 log(__int64 x) { return static_cast<__int64> (std::log (static_cast<double> (x))); }
528  static inline __int64 log10(__int64 x) { return static_cast<__int64> (std::log10 (static_cast<double> (x))); }
529 };
530 
531 template<>
532 struct ScalarTraits<unsigned __int64>
533 {
534  typedef unsigned __int64 magnitudeType;
535  typedef unsigned __int64 halfPrecision;
536  typedef unsigned __int64 doublePrecision;
537  typedef unsigned __int64 coordinateType;
538  static const bool isComplex = false;
539  static const bool isOrdinal = true;
540  static const bool isComparable = true;
541  static const bool hasMachineParameters = false;
542  // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
543  static inline magnitudeType magnitude(unsigned __int64 a) { return static_cast<unsigned __int64>(std::fabs(static_cast<double>(a))); }
544  static inline unsigned __int64 zero() { return 0; }
545  static inline unsigned __int64 one() { return 1; }
546  static inline unsigned __int64 conjugate(unsigned __int64 x) { return x; }
547  static inline unsigned __int64 real(unsigned __int64 x) { return x; }
548  static inline unsigned __int64 imag(unsigned __int64) { return 0; }
549  static inline void seedrandom(unsigned int s) {
550  std::srand(s);
551 #ifdef __APPLE__
552  // throw away first random number to address bug 3655
553  // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
554  random();
555 #endif
556  }
557  //static inline int random() { return (-1 + 2*rand()); } // RAB: This version should be used to be consistent with others
558  static inline unsigned __int64 random() { return std::rand(); } // RAB: This version should be used for an unsigned int, not int
559  static inline std::string name() { return "unsigned __int64"; }
560  static inline unsigned __int64 squareroot(unsigned __int64 x) { return (unsigned __int64) std::sqrt((double) x); }
561  static inline unsigned __int64 pow(unsigned __int64 x, unsigned __int64 y) { return (unsigned __int64) std::pow((double)x,(double)y); }
562  // Note: Depending on the number of bits in unsigned __int64,
563  // the cast from unsigned __int64 to double may not be exact.
564  static inline unsigned __int64 log(unsigned __int64 x) { return static_cast<unsigned __int64> (std::log (static_cast<double> (x))); }
565  static inline unsigned __int64 log10(unsigned __int64 x) { return static_cast<unsigned __int64> (std::log10 (static_cast<double> (x))); }
566 };
567 
568 #endif // HAVE_TEUCHOS___INT64
569 
570 
571 #ifndef __sun
572 extern TEUCHOSCORE_LIB_DLL_EXPORT const float flt_nan;
573 #endif
574 
575 
576 template<>
577 struct ScalarTraits<float>
578 {
579  typedef float magnitudeType;
580  typedef float halfPrecision; // should become IEEE754-2008 binary16 or fp16 later, perhaps specified at configure according to architectural support
581  typedef double doublePrecision;
582  typedef float coordinateType;
583  static const bool isComplex = false;
584  static const bool isOrdinal = false;
585  static const bool isComparable = true;
586  static const bool hasMachineParameters = true;
587  static inline float eps() {
588  return std::numeric_limits<float>::epsilon();
589  }
590  static inline float sfmin() {
591  return std::numeric_limits<float>::min();
592  }
593  static inline float base() {
594  return static_cast<float>(std::numeric_limits<float>::radix);
595  }
596  static inline float prec() {
597  return eps()*base();
598  }
599  static inline float t() {
600  return static_cast<float>(std::numeric_limits<float>::digits);
601  }
602  static inline float rnd() {
603  return ( std::numeric_limits<float>::round_style == std::round_to_nearest ? one() : zero() );
604  }
605  static inline float emin() {
606  return static_cast<float>(std::numeric_limits<float>::min_exponent);
607  }
608  static inline float rmin() {
609  return std::numeric_limits<float>::min();
610  }
611  static inline float emax() {
612  return static_cast<float>(std::numeric_limits<float>::max_exponent);
613  }
614  static inline float rmax() {
615  return std::numeric_limits<float>::max();
616  }
617  static inline magnitudeType magnitude(float a)
618  {
619 #ifdef TEUCHOS_DEBUG
620  TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
621  a, "Error, the input value to magnitude(...) a = " << a << " can not be NaN!" );
622 #endif
623  return std::fabs(a);
624  }
625  static inline float zero() { return(0.0f); }
626  static inline float one() { return(1.0f); }
627  static inline float conjugate(float x) { return(x); }
628  static inline float real(float x) { return x; }
629  static inline float imag(float) { return zero(); }
630  static inline float nan() {
631 #ifdef __sun
632  return 0.0f/std::sin(0.0f);
633 #else
634  return flt_nan;
635 #endif
636  }
637  static inline bool isnaninf(float x) {
638  return generic_real_isnaninf<float>(x);
639  }
640  static inline void seedrandom(unsigned int s) {
641  std::srand(s);
642 #ifdef __APPLE__
643  // throw away first random number to address bug 3655
644  // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
645  random();
646 #endif
647  }
648  static inline float random() { float rnd = (float) std::rand() / RAND_MAX; return (-1.0f + 2.0f * rnd); }
649  static inline std::string name() { return "float"; }
650  static inline float squareroot(float x)
651  {
652 #ifdef TEUCHOS_DEBUG
653  TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
654  x, "Error, the input value to squareroot(...) x = " << x << " can not be NaN!" );
655 #endif
656  errno = 0;
657  const float rtn = std::sqrt(x);
658  if (errno)
659  return nan();
660  return rtn;
661  }
662  static inline float pow(float x, float y) { return std::pow(x,y); }
663  static inline float pi() { return 3.14159265358979323846f; }
664  static inline float log(float x) { return std::log(x); }
665  static inline float log10(float x) { return std::log10(x); }
666 };
667 
668 
669 #ifndef __sun
670 extern TEUCHOSCORE_LIB_DLL_EXPORT const double dbl_nan;
671 #endif
672 
673 
674 template<>
675 struct ScalarTraits<double>
676 {
677  typedef double magnitudeType;
678  typedef float halfPrecision;
679  /* there are different options as to how to double "double"
680  - QD's DD (if available)
681  - ARPREC
682  - GNU MP
683  - a true hardware quad
684 
685  in the shortterm, this should be specified at configure time. I have inserted a configure-time option (--enable-teuchos-double-to-dd)
686  which uses QD's DD when available. This must be used alongside --enable-teuchos-qd.
687  */
688 #if defined(HAVE_TEUCHOS_DOUBLE_TO_QD)
689  typedef dd_real doublePrecision;
690 #elif defined(HAVE_TEUCHOS_DOUBLE_TO_ARPREC)
691  typedef mp_real doublePrecision;
692 #else
693  typedef double doublePrecision; // don't double "double" in this case
694 #endif
695  typedef double coordinateType;
696  static const bool isComplex = false;
697  static const bool isOrdinal = false;
698  static const bool isComparable = true;
699  static const bool hasMachineParameters = true;
700  static inline double eps() {
701  return std::numeric_limits<double>::epsilon();
702  }
703  static inline double sfmin() {
704  return std::numeric_limits<double>::min();
705  }
706  static inline double base() {
707  return std::numeric_limits<double>::radix;
708  }
709  static inline double prec() {
710  return eps()*base();
711  }
712  static inline double t() {
713  return std::numeric_limits<double>::digits;
714  }
715  static inline double rnd() {
716  return ( std::numeric_limits<double>::round_style == std::round_to_nearest ? double(1.0) : double(0.0) );
717  }
718  static inline double emin() {
719  return std::numeric_limits<double>::min_exponent;
720  }
721  static inline double rmin() {
722  return std::numeric_limits<double>::min();
723  }
724  static inline double emax() {
725  return std::numeric_limits<double>::max_exponent;
726  }
727  static inline double rmax() {
728  return std::numeric_limits<double>::max();
729  }
730  static inline magnitudeType magnitude(double a)
731  {
732 #ifdef TEUCHOS_DEBUG
733  TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
734  a, "Error, the input value to magnitude(...) a = " << a << " can not be NaN!" );
735 #endif
736  return std::fabs(a);
737  }
738  static inline double zero() { return 0.0; }
739  static inline double one() { return 1.0; }
740  static inline double conjugate(double x) { return(x); }
741  static inline double real(double x) { return(x); }
742  static inline double imag(double) { return(0); }
743  static inline double nan() {
744 #ifdef __sun
745  return 0.0/std::sin(0.0);
746 #else
747  return dbl_nan;
748 #endif
749  }
750  static inline bool isnaninf(double x) {
751  return generic_real_isnaninf<double>(x);
752  }
753  static inline void seedrandom(unsigned int s) {
754  std::srand(s);
755 #ifdef __APPLE__
756  // throw away first random number to address bug 3655
757  // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
758  random();
759 #endif
760  }
761  static inline double random() { double rnd = (double) std::rand() / RAND_MAX; return (double)(-1.0 + 2.0 * rnd); }
762  static inline std::string name() { return "double"; }
763  static inline double squareroot(double x)
764  {
765 #ifdef TEUCHOS_DEBUG
766  TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
767  x, "Error, the input value to squareroot(...) x = " << x << " can not be NaN!" );
768 #endif
769  errno = 0;
770  const double rtn = std::sqrt(x);
771  if (errno)
772  return nan();
773  return rtn;
774  }
775  static inline double pow(double x, double y) { return std::pow(x,y); }
776  static inline double pi() { return 3.14159265358979323846; }
777  static inline double log(double x) { return std::log(x); }
778  static inline double log10(double x) { return std::log10(x); }
779 };
780 
781 
782 #ifdef HAVE_TEUCHOSCORE_QUADMATH
783 
784 template<>
785 struct ScalarTraits<__float128> {
786  typedef __float128 magnitudeType;
787  // Unfortunately, we can't rely on a standard __float256 type. It
788  // might make sense to use qd_real, but mixing quadmath and QD might
789  // cause unforeseen issues.
790  typedef __float128 doublePrecision;
791  typedef double halfPrecision;
792  typedef __float128 coordinateType;
793 
794  static const bool isComplex = false;
795  static const bool isOrdinal = false;
796  static const bool isComparable = true;
797  static const bool hasMachineParameters = true;
798 
799  static __float128 eps () {
800  return FLT128_EPSILON;
801  }
802  static __float128 sfmin () {
803  return FLT128_MIN; // ???
804  }
805  static __float128 base () {
806  return 2;
807  }
808  static __float128 prec () {
809  return eps () * base ();
810  }
811  static __float128 t () {
812  return FLT128_MANT_DIG;
813  }
814  static __float128 rnd () {
815  return 1.0;
816  }
817  static __float128 emin () {
818  return FLT128_MIN_EXP;
819  }
820  static __float128 rmin () {
821  return FLT128_MIN; // ??? // should be base^(emin-1)
822  }
823  static __float128 emax () {
824  return FLT128_MAX_EXP;
825  }
826  static __float128 rmax () {
827  return FLT128_MAX; // ??? // should be (base^emax)*(1-eps)
828  }
829  static magnitudeType magnitude (const __float128& x) {
830  return fabsq (x);
831  }
832  static __float128 zero () {
833  return 0.0;
834  }
835  static __float128 one () {
836  return 1.0;
837  }
838  static __float128 conjugate (const __float128& x) {
839  return x;
840  }
841  static __float128 real (const __float128& x) {
842  return x;
843  }
844  static __float128 imag (const __float128& /* x */) {
845  return 0.0;
846  }
847  static __float128 nan () {
848  return strtoflt128 ("NAN()", NULL); // ???
849  }
850  static bool isnaninf (const __float128& x) {
851  return isinfq (x) || isnanq (x);
852  }
853  static inline void seedrandom (unsigned int s) {
854  std::srand (s);
855 #ifdef __APPLE__
856  // throw away first random number to address bug 3655
857  // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
858  random ();
859 #endif
860  }
861  static __float128 random () {
862  // Half the smallest normalized double, is the scaling factor of
863  // the lower-order term in the double-double representation.
864  const __float128 scalingFactor =
865  static_cast<__float128> (std::numeric_limits<double>::min ()) /
866  static_cast<__float128> (2.0);
867  const __float128 higherOrderTerm =
868  static_cast<__float128> (ScalarTraits<double>::random ());
869  const __float128 lowerOrderTerm =
870  static_cast<__float128> (ScalarTraits<double>::random ()) *
871  scalingFactor;
872  return higherOrderTerm + lowerOrderTerm;
873  }
874  static std::string name () {
875  return "__float128";
876  }
877  static __float128 squareroot (const __float128& x) {
878  return sqrtq (x);
879  }
880  static __float128 pow (const __float128& x, const __float128& y) {
881  return powq (x, y);
882  }
883  static __float128 pi() { return 3.14159265358979323846; }
884  static __float128 log (const __float128& x) {
885  return logq (x);
886  }
887  static __float128 log10 (const __float128& x) {
888  return log10q (x);
889  }
890 };
891 #endif // HAVE_TEUCHOSCORE_QUADMATH
892 
893 
894 
895 #ifdef HAVE_TEUCHOS_QD
896 
897 bool operator&&(const dd_real &a, const dd_real &b);
898 bool operator&&(const qd_real &a, const qd_real &b);
899 
900 template<>
901 struct ScalarTraits<dd_real>
902 {
903  typedef dd_real magnitudeType;
904  typedef double halfPrecision;
905  typedef qd_real doublePrecision;
906  typedef dd_real coordinateType;
907  static const bool isComplex = false;
908  static const bool isOrdinal = false;
909  static const bool isComparable = true;
910  static const bool hasMachineParameters = true;
911  static inline dd_real eps() { return std::numeric_limits<dd_real>::epsilon(); }
912  static inline dd_real sfmin() { return std::numeric_limits<dd_real>::min(); }
913  static inline dd_real base() { return std::numeric_limits<dd_real>::radix; }
914  static inline dd_real prec() { return eps()*base(); }
915  static inline dd_real t() { return std::numeric_limits<dd_real>::digits; }
916  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) ); }
917  static inline dd_real emin() { return std::numeric_limits<dd_real>::min_exponent; }
918  static inline dd_real rmin() { return std::numeric_limits<dd_real>::min(); }
919  static inline dd_real emax() { return std::numeric_limits<dd_real>::max_exponent; }
920  static inline dd_real rmax() { return std::numeric_limits<dd_real>::max(); }
921  static inline magnitudeType magnitude(dd_real a)
922  {
923 #ifdef TEUCHOS_DEBUG
924  TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
925  a, "Error, the input value to magnitude(...) a = " << a << " can not be NaN!" );
926 #endif
927  return ::abs(a);
928  }
929  static inline dd_real zero() { return dd_real(0.0); }
930  static inline dd_real one() { return dd_real(1.0); }
931  static inline dd_real conjugate(dd_real x) { return(x); }
932  static inline dd_real real(dd_real x) { return x ; }
933  static inline dd_real imag(dd_real) { return zero(); }
934  static inline dd_real nan() { return dd_real::_nan; }
935  static inline bool isnaninf(dd_real x) { return isnan(x) || isinf(x); }
936  static inline void seedrandom(unsigned int s) {
937  // ddrand() uses std::rand(), so the std::srand() is our seed
938  std::srand(s);
939 #ifdef __APPLE__
940  // throw away first random number to address bug 3655
941  // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
942  random();
943 #endif
944  }
945  static inline dd_real random() { return ddrand(); }
946  static inline std::string name() { return "dd_real"; }
947  static inline dd_real squareroot(dd_real x)
948  {
949 #ifdef TEUCHOS_DEBUG
950  TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
951  x, "Error, the input value to squareroot(...) x = " << x << " can not be NaN!" );
952 #endif
953  return ::sqrt(x);
954  }
955  static inline dd_real pow(dd_real x, dd_real y) { return ::pow(x,y); }
956  static inline dd_real pi() { return 3.14159265358979323846; }
957  // dd_real puts its transcendental functions in the global namespace.
958  static inline dd_real log(dd_real x) { return ::log(x); }
959  static inline dd_real log10(dd_real x) { return ::log10(x); }
960 };
961 
962 
963 template<>
964 struct ScalarTraits<qd_real>
965 {
966  typedef qd_real magnitudeType;
967  typedef dd_real halfPrecision;
968  typedef qd_real doublePrecision;
969  typedef qd_real coordinateType;
970  static const bool isComplex = false;
971  static const bool isOrdinal = false;
972  static const bool isComparable = true;
973  static const bool hasMachineParameters = true;
974  static inline qd_real eps() { return std::numeric_limits<qd_real>::epsilon(); }
975  static inline qd_real sfmin() { return std::numeric_limits<qd_real>::min(); }
976  static inline qd_real base() { return std::numeric_limits<qd_real>::radix; }
977  static inline qd_real prec() { return eps()*base(); }
978  static inline qd_real t() { return std::numeric_limits<qd_real>::digits; }
979  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) ); }
980  static inline qd_real emin() { return std::numeric_limits<qd_real>::min_exponent; }
981  static inline qd_real rmin() { return std::numeric_limits<qd_real>::min(); }
982  static inline qd_real emax() { return std::numeric_limits<qd_real>::max_exponent; }
983  static inline qd_real rmax() { return std::numeric_limits<qd_real>::max(); }
984  static inline magnitudeType magnitude(qd_real a)
985  {
986 #ifdef TEUCHOS_DEBUG
987  TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
988  a, "Error, the input value to magnitude(...) a = " << a << " can not be NaN!" );
989 #endif
990  return ::abs(a);
991  }
992  static inline qd_real zero() { return qd_real(0.0); }
993  static inline qd_real one() { return qd_real(1.0); }
994  static inline qd_real conjugate(qd_real x) { return(x); }
995  static inline qd_real real(qd_real x) { return x ; }
996  static inline qd_real imag(qd_real) { return zero(); }
997  static inline qd_real nan() { return qd_real::_nan; }
998  static inline bool isnaninf(qd_real x) { return isnan(x) || isinf(x); }
999  static inline void seedrandom(unsigned int s) {
1000  // qdrand() uses std::rand(), so the std::srand() is our seed
1001  std::srand(s);
1002 #ifdef __APPLE__
1003  // throw away first random number to address bug 3655
1004  // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
1005  random();
1006 #endif
1007  }
1008  static inline qd_real random() { return qdrand(); }
1009  static inline std::string name() { return "qd_real"; }
1010  static inline qd_real squareroot(qd_real x)
1011  {
1012 #ifdef TEUCHOS_DEBUG
1013  TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
1014  x, "Error, the input value to squareroot(...) x = " << x << " can not be NaN!" );
1015 #endif
1016  return ::sqrt(x);
1017  }
1018  static inline qd_real pow(qd_real x, qd_real y) { return ::pow(x,y); }
1019  static inline qd_real pi() { return 3.14159265358979323846; }
1020  // qd_real puts its transcendental functions in the global namespace.
1021  static inline qd_real log(qd_real x) { return ::log(x); }
1022  static inline qd_real log10(qd_real x) { return ::log10(x); }
1023 };
1024 
1025 
1026 #endif // HAVE_TEUCHOS_QD
1027 
1028 
1029 #ifdef HAVE_TEUCHOS_GNU_MP
1030 
1031 
1032 extern gmp_randclass gmp_rng;
1033 
1034 
1035 /* Regarding halfPrecision, doublePrecision and mpf_class:
1036  Because the precision of an mpf_class float is not determined by the data type,
1037  there is no way to fill the typedefs for this object.
1038 
1039  Instead, we could create new data classes (e.g., Teuchos::MPF128, Teuchos::MPF256) for
1040  commonly used levels of precision, and fill out ScalarTraits for these. This would allow us
1041  to typedef the promotions and demotions in the appropriate way. These classes would serve to
1042  wrap an mpf_class object, calling the constructor for the appropriate precision, exposing the
1043  arithmetic routines but hiding the precision-altering routines.
1044 
1045  Alternatively (perhaps, preferably), would could create a single class templated on the precision (e.g., Teuchos::MPF<N>).
1046  Then we have a single (partially-specialized) implementation of ScalarTraits. This class, as above, must expose all of the
1047  operations expected of a scalar type; however, most of these can be trivially stolen from the gmpcxx.h header file
1048 
1049  CGB/RAB, 01/05/2009
1050 */
1051 template<>
1052 struct ScalarTraits<mpf_class>
1053 {
1054  typedef mpf_class magnitudeType;
1055  typedef mpf_class halfPrecision;
1056  typedef mpf_class doublePrecision;
1057  typedef mpf_class coordinateType;
1058  static const bool isComplex = false;
1059  static const bool hasMachineParameters = false;
1060  // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
1061  static magnitudeType magnitude(mpf_class a) { return std::abs(a); }
1062  static inline mpf_class zero() { mpf_class zero = 0.0; return zero; }
1063  static inline mpf_class one() { mpf_class one = 1.0; return one; }
1064  static inline mpf_class conjugate(mpf_class x) { return x; }
1065  static inline mpf_class real(mpf_class x) { return(x); }
1066  static inline mpf_class imag(mpf_class x) { return(0); }
1067  static inline bool isnaninf(mpf_class x) { return false; } // mpf_class currently can't handle nan or inf!
1068  static inline void seedrandom(unsigned int s) {
1069  unsigned long int seedVal = static_cast<unsigned long int>(s);
1070  gmp_rng.seed( seedVal );
1071  }
1072  static inline mpf_class random() {
1073  return gmp_rng.get_f();
1074  }
1075  static inline std::string name() { return "mpf_class"; }
1076  static inline mpf_class squareroot(mpf_class x) { return std::sqrt(x); }
1077  static inline mpf_class pow(mpf_class x, mpf_class y) { return pow(x,y); }
1078  // Todo: RAB: 2004/05/28: Add nan() and isnaninf() functions when needed!
1079 };
1080 
1081 #endif // HAVE_TEUCHOS_GNU_MP
1082 
1083 #ifdef HAVE_TEUCHOS_ARPREC
1084 
1085 /* See discussion above for mpf_class, regarding halfPrecision and doublePrecision. Something similar will need to be done
1086  for ARPREC. */
1087 template<>
1088 struct ScalarTraits<mp_real>
1089 {
1090  typedef mp_real magnitudeType;
1091  typedef double halfPrecision;
1092  typedef mp_real doublePrecision;
1093  typedef mp_real coordinateType;
1094  static const bool isComplex = false;
1095  static const bool isComparable = true;
1096  static const bool isOrdinal = false;
1097  static const bool hasMachineParameters = false;
1098  // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
1099  static magnitudeType magnitude(mp_real a) { return abs(a); }
1100  static inline mp_real zero() { mp_real zero = 0.0; return zero; }
1101  static inline mp_real one() { mp_real one = 1.0; return one; }
1102  static inline mp_real conjugate(mp_real x) { return x; }
1103  static inline mp_real real(mp_real x) { return(x); }
1104  static inline mp_real imag(mp_real x) { return zero(); }
1105  static inline bool isnaninf(mp_real x) { return false; } // ToDo: Change this?
1106  static inline void seedrandom(unsigned int s) {
1107  long int seedVal = static_cast<long int>(s);
1108  srand48(seedVal);
1109  }
1110  static inline mp_real random() { return mp_rand(); }
1111  static inline std::string name() { return "mp_real"; }
1112  static inline mp_real squareroot(mp_real x) { return sqrt(x); }
1113  static inline mp_real pow(mp_real x, mp_real y) { return pow(x,y); }
1114  static inline mp_real pi() { return 3.14159265358979323846; }
1115  // Todo: RAB: 2004/05/28: Add nan() and isnaninf() functions when needed!
1116 };
1117 
1118 
1119 #endif // HAVE_TEUCHOS_ARPREC
1120 
1121 
1122 #ifdef HAVE_TEUCHOS_COMPLEX
1123 
1124 
1125 // Partial specialization for std::complex numbers templated on real type T
1126 template<class T>
1127 struct ScalarTraits<
1128  std::complex<T>
1129 >
1130 {
1131  typedef std::complex<T> ComplexT;
1132  typedef std::complex<typename ScalarTraits<T>::halfPrecision> halfPrecision;
1133  typedef std::complex<typename ScalarTraits<T>::doublePrecision> doublePrecision;
1134  typedef typename ScalarTraits<T>::magnitudeType magnitudeType;
1135  typedef typename ScalarTraits<T>::coordinateType coordinateType;
1136  static const bool isComplex = true;
1137  static const bool isOrdinal = ScalarTraits<T>::isOrdinal;
1138  static const bool isComparable = false;
1139  static const bool hasMachineParameters = true;
1140  static inline magnitudeType eps() { return ScalarTraits<magnitudeType>::eps(); }
1141  static inline magnitudeType sfmin() { return ScalarTraits<magnitudeType>::sfmin(); }
1142  static inline magnitudeType base() { return ScalarTraits<magnitudeType>::base(); }
1143  static inline magnitudeType prec() { return ScalarTraits<magnitudeType>::prec(); }
1144  static inline magnitudeType t() { return ScalarTraits<magnitudeType>::t(); }
1145  static inline magnitudeType rnd() { return ScalarTraits<magnitudeType>::rnd(); }
1146  static inline magnitudeType emin() { return ScalarTraits<magnitudeType>::emin(); }
1147  static inline magnitudeType rmin() { return ScalarTraits<magnitudeType>::rmin(); }
1148  static inline magnitudeType emax() { return ScalarTraits<magnitudeType>::emax(); }
1149  static inline magnitudeType rmax() { return ScalarTraits<magnitudeType>::rmax(); }
1150  static magnitudeType magnitude(ComplexT a)
1151  {
1152 #ifdef TEUCHOS_DEBUG
1153  TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
1154  a, "Error, the input value to magnitude(...) a = " << a << " can not be NaN!" );
1155 #endif
1156  return std::abs(a);
1157  }
1158  static inline ComplexT zero() { return ComplexT(ScalarTraits<magnitudeType>::zero(),ScalarTraits<magnitudeType>::zero()); }
1159  static inline ComplexT one() { return ComplexT(ScalarTraits<magnitudeType>::one(),ScalarTraits<magnitudeType>::zero()); }
1160  static inline ComplexT conjugate(ComplexT a){ return ComplexT(a.real(),-a.imag()); }
1161  static inline magnitudeType real(ComplexT a) { return a.real(); }
1162  static inline magnitudeType imag(ComplexT a) { return a.imag(); }
1163  static inline ComplexT nan() { return ComplexT(ScalarTraits<magnitudeType>::nan(),ScalarTraits<magnitudeType>::nan()); }
1164  static inline bool isnaninf(ComplexT x) { return ScalarTraits<magnitudeType>::isnaninf(x.real()) || ScalarTraits<magnitudeType>::isnaninf(x.imag()); }
1165  static inline void seedrandom(unsigned int s) { ScalarTraits<magnitudeType>::seedrandom(s); }
1166  static inline ComplexT random()
1167  {
1168  const T rnd1 = ScalarTraits<magnitudeType>::random();
1169  const T rnd2 = ScalarTraits<magnitudeType>::random();
1170  return ComplexT(rnd1,rnd2);
1171  }
1172  static inline std::string name() { return std::string("std::complex<")+std::string(ScalarTraits<magnitudeType>::name())+std::string(">"); }
1173  // This will only return one of the square roots of x, the other can be obtained by taking its conjugate
1174  static inline ComplexT squareroot(ComplexT x)
1175  {
1176 #ifdef TEUCHOS_DEBUG
1177  TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
1178  x, "Error, the input value to squareroot(...) x = " << x << " can not be NaN!" );
1179 #endif
1180  typedef ScalarTraits<magnitudeType> STMT;
1181  const T r = x.real(), i = x.imag(), zero = STMT::zero(), two = 2.0;
1182  const T a = STMT::squareroot((r*r)+(i*i));
1183  const T nr = STMT::squareroot((a+r)/two);
1184  const T ni = ( i == zero ? zero : STMT::squareroot((a-r)/two) );
1185  return ComplexT(nr,ni);
1186  }
1187  // 2010/03/19: rabartl: Above, I had to add the check for i == zero
1188  // to avoid a returned NaN on the Intel 10.1 compiler. For some
1189  // reason, having these two squareroot calls in a row produce a NaN
1190  // in an optimized build with this compiler. Amazingly, when I put
1191  // in print statements (i.e. std::cout << ...) the NaN went away and
1192  // the second squareroot((a-r)/two) returned zero correctly. I
1193  // guess that calling the output routine flushed the registers or
1194  // something and restarted the squareroot rountine for this compiler
1195  // or something similar. Actually, due to roundoff, it is possible that a-r
1196  // might be slightly less than zero (i.e. -1e-16) and that would cause
1197  // a possbile NaN return. THe above if test is the right thing to do
1198  // I think and is very cheap.
1199  static inline ComplexT pow(ComplexT x, ComplexT y) { return pow(x,y); }
1200  static inline ComplexT pi() { return ScalarTraits<T>::pi(); }
1201 };
1202 
1203 #endif // HAVE_TEUCHOS_COMPLEX
1204 #endif // DOXYGEN_SHOULD_SKIP_THIS
1205 
1206 } // Teuchos namespace
1207 
1208 #endif // _TEUCHOS_SCALARTRAITS_HPP_
Teuchos::ScalarTraits::magnitudeType
T magnitudeType
Mandatory typedef for result of magnitude.
Definition: Teuchos_ScalarTraitsDecl.hpp:93
Teuchos::ScalarTraits::eps
static magnitudeType eps()
Returns relative machine precision.
Definition: Teuchos_ScalarTraitsDecl.hpp:112
Teuchos::ScalarTraits::squareroot
static T squareroot(T x)
Returns a number of magnitudeType that is the square root of this scalar type x.
Definition: Teuchos_ScalarTraitsDecl.hpp:154
Teuchos::ScalarTraits::isComparable
static const bool isComparable
Determines if scalar type supports relational operators such as &lt;, &gt;, &lt;=, &gt;=.
Definition: Teuchos_ScalarTraitsDecl.hpp:105
Teuchos::ScalarTraits::real
static magnitudeType real(T a)
Returns the real part of the scalar type a.
Definition: Teuchos_ScalarTraitsDecl.hpp:138
Teuchos_ConfigDefs.hpp
Teuchos header file which uses auto-configuration information to include necessary C++ headers...
Teuchos_ScalarTraitsDecl.hpp
Declaration and default implementation for basic traits for the scalar field type.
Teuchos::ScalarTraits::pow
static T pow(T x, T y)
Returns the result of raising one scalar x to the power y.
Definition: Teuchos_ScalarTraitsDecl.hpp:156
Teuchos::ScalarTraits::emax
static magnitudeType emax()
Returns the largest exponent before overflow.
Definition: Teuchos_ScalarTraitsDecl.hpp:128
Teuchos::ScalarTraits::base
static magnitudeType base()
Returns the base of the machine.
Definition: Teuchos_ScalarTraitsDecl.hpp:116
Teuchos::ScalarTraits::hasMachineParameters
static const bool hasMachineParameters
Determines if scalar type have machine-specific parameters (i.e. eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax() are supported).
Definition: Teuchos_ScalarTraitsDecl.hpp:110
Teuchos::ScalarTraits::name
static std::string name()
Returns the name of this scalar type.
Definition: Teuchos_ScalarTraitsDecl.hpp:152
Teuchos::ScalarTraits::rmax
static magnitudeType rmax()
Overflow theshold - (base^emax)*(1-eps)
Definition: Teuchos_ScalarTraitsDecl.hpp:130
Teuchos::ScalarTraits::sfmin
static magnitudeType sfmin()
Returns safe minimum (sfmin), such that 1/sfmin does not overflow.
Definition: Teuchos_ScalarTraitsDecl.hpp:114
Teuchos::ScalarTraits::pi
static T pi()
Returns the value of PI.
Definition: Teuchos_ScalarTraitsDecl.hpp:158
Teuchos::ScalarTraits::coordinateType
T coordinateType
Typedef for coordinates.
Definition: Teuchos_ScalarTraitsDecl.hpp:99
Teuchos::ScalarTraits::conjugate
static T conjugate(T a)
Returns the conjugate of the scalar type a.
Definition: Teuchos_ScalarTraitsDecl.hpp:142
Teuchos::ScalarTraits::isOrdinal
static const bool isOrdinal
Determines if scalar type is an ordinal type.
Definition: Teuchos_ScalarTraitsDecl.hpp:103
Teuchos::ScalarTraits::prec
static magnitudeType prec()
Returns eps*base.
Definition: Teuchos_ScalarTraitsDecl.hpp:118
Teuchos::ScalarTraits::t
static magnitudeType t()
Returns the number of (base) digits in the mantissa.
Definition: Teuchos_ScalarTraitsDecl.hpp:120
Teuchos::ScalarTraits::rmin
static magnitudeType rmin()
Returns the underflow threshold - base^(emin-1)
Definition: Teuchos_ScalarTraitsDecl.hpp:126
Teuchos::ScalarTraits::doublePrecision
T doublePrecision
Typedef for double precision.
Definition: Teuchos_ScalarTraitsDecl.hpp:97
Teuchos::ScalarTraits::seedrandom
static void seedrandom(unsigned int s)
Seed the random number generator returned by random().
Definition: Teuchos_ScalarTraitsDecl.hpp:148
Teuchos::ScalarTraits::imag
static magnitudeType imag(T a)
Returns the imaginary part of the scalar type a.
Definition: Teuchos_ScalarTraitsDecl.hpp:140
Teuchos::ScalarTraits::isnaninf
static bool isnaninf(const T &x)
Returns true if x is NaN or Inf.
Definition: Teuchos_ScalarTraitsDecl.hpp:146
Teuchos::ScalarTraits::magnitude
static magnitudeType magnitude(T a)
Returns the magnitudeType of the scalar type a.
Definition: Teuchos_ScalarTraitsDecl.hpp:132
Teuchos::ScalarTraits::nan
static T nan()
Returns a number that represents NaN.
Definition: Teuchos_ScalarTraitsDecl.hpp:144
Teuchos::ScalarTraits::zero
static T zero()
Returns representation of zero for this scalar type.
Definition: Teuchos_ScalarTraitsDecl.hpp:134
Teuchos::ScalarTraits::random
static T random()
Returns a random number (between -one() and +one()) of this scalar type.
Definition: Teuchos_ScalarTraitsDecl.hpp:150
Teuchos::ScalarTraits::halfPrecision
T halfPrecision
Typedef for half precision.
Definition: Teuchos_ScalarTraitsDecl.hpp:95
Teuchos::ScalarTraits::isComplex
static const bool isComplex
Determines if scalar type is std::complex.
Definition: Teuchos_ScalarTraitsDecl.hpp:101
Teuchos::ScalarTraits::emin
static magnitudeType emin()
Returns the minimum exponent before (gradual) underflow.
Definition: Teuchos_ScalarTraitsDecl.hpp:124
Teuchos::ScalarTraits::rnd
static magnitudeType rnd()
Returns 1.0 when rounding occurs in addition, 0.0 otherwise.
Definition: Teuchos_ScalarTraitsDecl.hpp:122
Teuchos::ScalarTraits::one
static T one()
Returns representation of one for this scalar type.
Definition: Teuchos_ScalarTraitsDecl.hpp:136
Generated by   doxygen 1.8.5