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