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