Teuchos - Trilinos Tools Package  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Teuchos_ScalarTraits.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Teuchos: Common Tools Package
5 // Copyright (2004) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
42 // NOTE: Before adding specializations of ScalarTraits, make sure that they do
43 // not duplicate specializations already present in PyTrilinos (see
44 // packages/PyTrilinos/src/Teuchos_Traits.i)
45 
46 // NOTE: halfPrecision and doublePrecision are not currently implemented for
47 // ARPREC, GMP or the ordinal types (e.g., int, char)
48 
49 #ifndef _TEUCHOS_SCALARTRAITS_HPP_
50 #define _TEUCHOS_SCALARTRAITS_HPP_
51 
56 #include "Teuchos_ConfigDefs.hpp"
57 
58 #ifdef HAVE_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 
113 #include "Teuchos_ScalarTraitsDecl.hpp"
114 
115 
116 namespace Teuchos {
117 
118 
119 #ifndef DOXYGEN_SHOULD_SKIP_THIS
120 
121 
122 TEUCHOSCORE_LIB_DLL_EXPORT
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() / 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_TEUCHOSCORE_QUADMATH
787 
788 template<>
789 struct ScalarTraits<__float128> {
790  typedef __float128 magnitudeType;
791  // Unfortunately, we can't rely on a standard __float256 type. It
792  // might make sense to use qd_real, but mixing quadmath and QD might
793  // cause unforeseen issues.
794  typedef __float128 doublePrecision;
795  typedef double halfPrecision;
796  typedef __float128 coordinateType;
797 
798  static const bool isComplex = false;
799  static const bool isOrdinal = false;
800  static const bool isComparable = true;
801  static const bool hasMachineParameters = true;
802 
803  static __float128 eps () {
804  return FLT128_EPSILON;
805  }
806  static __float128 sfmin () {
807  return FLT128_MIN; // ???
808  }
809  static __float128 base () {
810  return 2;
811  }
812  static __float128 prec () {
813  return eps () * base ();
814  }
815  static __float128 t () {
816  return FLT128_MANT_DIG;
817  }
818  static __float128 rnd () {
819  return 1.0;
820  }
821  static __float128 emin () {
822  return FLT128_MIN_EXP;
823  }
824  static __float128 rmin () {
825  return FLT128_MIN; // ??? // should be base^(emin-1)
826  }
827  static __float128 emax () {
828  return FLT128_MAX_EXP;
829  }
830  static __float128 rmax () {
831  return FLT128_MAX; // ??? // should be (base^emax)*(1-eps)
832  }
833  static magnitudeType magnitude (const __float128& x) {
834  return fabsq (x);
835  }
836  static __float128 zero () {
837  return 0.0;
838  }
839  static __float128 one () {
840  return 1.0;
841  }
842  static __float128 conjugate (const __float128& x) {
843  return x;
844  }
845  static __float128 real (const __float128& x) {
846  return x;
847  }
848  static __float128 imag (const __float128& /* x */) {
849  return 0.0;
850  }
851  static __float128 nan () {
852  return strtoflt128 ("NAN()", NULL); // ???
853  }
854  static bool isnaninf (const __float128& x) {
855  return isinfq (x) || isnanq (x);
856  }
857  static inline void seedrandom (unsigned int s) {
858  std::srand (s);
859 #ifdef __APPLE__
860  // throw away first random number to address bug 3655
861  // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
862  random ();
863 #endif
864  }
865  static __float128 random () {
866  // Half the smallest normalized double, is the scaling factor of
867  // the lower-order term in the double-double representation.
868  const __float128 scalingFactor =
869  static_cast<__float128> (std::numeric_limits<double>::min ()) /
870  static_cast<__float128> (2.0);
871  const __float128 higherOrderTerm =
872  static_cast<__float128> (ScalarTraits<double>::random ());
873  const __float128 lowerOrderTerm =
874  static_cast<__float128> (ScalarTraits<double>::random ()) *
875  scalingFactor;
876  return higherOrderTerm + lowerOrderTerm;
877  }
878  static std::string name () {
879  return "__float128";
880  }
881  static __float128 squareroot (const __float128& x) {
882  return sqrtq (x);
883  }
884  static __float128 pow (const __float128& x, const __float128& y) {
885  return powq (x, y);
886  }
887  static __float128 pi() { return 3.14159265358979323846; }
888  static __float128 log (const __float128& x) {
889  return logq (x);
890  }
891  static __float128 log10 (const __float128& x) {
892  return log10q (x);
893  }
894 };
895 #endif // HAVE_TEUCHOSCORE_QUADMATH
896 
897 
898 
899 #ifdef HAVE_TEUCHOS_QD
900 
901 bool operator&&(const dd_real &a, const dd_real &b);
902 bool operator&&(const qd_real &a, const qd_real &b);
903 
904 template<>
905 struct ScalarTraits<dd_real>
906 {
907  typedef dd_real magnitudeType;
908  typedef double halfPrecision;
909  typedef qd_real doublePrecision;
910  typedef dd_real coordinateType;
911  static const bool isComplex = false;
912  static const bool isOrdinal = false;
913  static const bool isComparable = true;
914  static const bool hasMachineParameters = true;
915  static inline dd_real eps() { return std::numeric_limits<dd_real>::epsilon(); }
916  static inline dd_real sfmin() { return std::numeric_limits<dd_real>::min(); }
917  static inline dd_real base() { return std::numeric_limits<dd_real>::radix; }
918  static inline dd_real prec() { return eps()*base(); }
919  static inline dd_real t() { return std::numeric_limits<dd_real>::digits; }
920  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) ); }
921  static inline dd_real emin() { return std::numeric_limits<dd_real>::min_exponent; }
922  static inline dd_real rmin() { return std::numeric_limits<dd_real>::min(); }
923  static inline dd_real emax() { return std::numeric_limits<dd_real>::max_exponent; }
924  static inline dd_real rmax() { return std::numeric_limits<dd_real>::max(); }
925  static inline magnitudeType magnitude(dd_real a)
926  {
927 #ifdef TEUCHOS_DEBUG
928  TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
929  a, "Error, the input value to magnitude(...) a = " << a << " can not be NaN!" );
930 #endif
931  return ::abs(a);
932  }
933  static inline dd_real zero() { return dd_real(0.0); }
934  static inline dd_real one() { return dd_real(1.0); }
935  static inline dd_real conjugate(dd_real x) { return(x); }
936  static inline dd_real real(dd_real x) { return x ; }
937  static inline dd_real imag(dd_real) { return zero(); }
938  static inline dd_real nan() { return dd_real::_nan; }
939  static inline bool isnaninf(dd_real x) { return isnan(x) || isinf(x); }
940  static inline void seedrandom(unsigned int s) {
941  // ddrand() uses std::rand(), so the std::srand() is our seed
942  std::srand(s);
943 #ifdef __APPLE__
944  // throw away first random number to address bug 3655
945  // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
946  random();
947 #endif
948  }
949  static inline dd_real random() { return ddrand(); }
950  static inline std::string name() { return "dd_real"; }
951  static inline dd_real squareroot(dd_real x)
952  {
953 #ifdef TEUCHOS_DEBUG
954  TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
955  x, "Error, the input value to squareroot(...) x = " << x << " can not be NaN!" );
956 #endif
957  return ::sqrt(x);
958  }
959  static inline dd_real pow(dd_real x, dd_real y) { return ::pow(x,y); }
960  static inline dd_real pi() { return 3.14159265358979323846; }
961  // dd_real puts its transcendental functions in the global namespace.
962  static inline dd_real log(dd_real x) { return ::log(x); }
963  static inline dd_real log10(dd_real x) { return ::log10(x); }
964 };
965 
966 
967 template<>
968 struct ScalarTraits<qd_real>
969 {
970  typedef qd_real magnitudeType;
971  typedef dd_real halfPrecision;
972  typedef qd_real doublePrecision;
973  typedef qd_real coordinateType;
974  static const bool isComplex = false;
975  static const bool isOrdinal = false;
976  static const bool isComparable = true;
977  static const bool hasMachineParameters = true;
978  static inline qd_real eps() { return std::numeric_limits<qd_real>::epsilon(); }
979  static inline qd_real sfmin() { return std::numeric_limits<qd_real>::min(); }
980  static inline qd_real base() { return std::numeric_limits<qd_real>::radix; }
981  static inline qd_real prec() { return eps()*base(); }
982  static inline qd_real t() { return std::numeric_limits<qd_real>::digits; }
983  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) ); }
984  static inline qd_real emin() { return std::numeric_limits<qd_real>::min_exponent; }
985  static inline qd_real rmin() { return std::numeric_limits<qd_real>::min(); }
986  static inline qd_real emax() { return std::numeric_limits<qd_real>::max_exponent; }
987  static inline qd_real rmax() { return std::numeric_limits<qd_real>::max(); }
988  static inline magnitudeType magnitude(qd_real a)
989  {
990 #ifdef TEUCHOS_DEBUG
991  TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
992  a, "Error, the input value to magnitude(...) a = " << a << " can not be NaN!" );
993 #endif
994  return ::abs(a);
995  }
996  static inline qd_real zero() { return qd_real(0.0); }
997  static inline qd_real one() { return qd_real(1.0); }
998  static inline qd_real conjugate(qd_real x) { return(x); }
999  static inline qd_real real(qd_real x) { return x ; }
1000  static inline qd_real imag(qd_real) { return zero(); }
1001  static inline qd_real nan() { return qd_real::_nan; }
1002  static inline bool isnaninf(qd_real x) { return isnan(x) || isinf(x); }
1003  static inline void seedrandom(unsigned int s) {
1004  // qdrand() uses std::rand(), so the std::srand() is our seed
1005  std::srand(s);
1006 #ifdef __APPLE__
1007  // throw away first random number to address bug 3655
1008  // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
1009  random();
1010 #endif
1011  }
1012  static inline qd_real random() { return qdrand(); }
1013  static inline std::string name() { return "qd_real"; }
1014  static inline qd_real squareroot(qd_real x)
1015  {
1016 #ifdef TEUCHOS_DEBUG
1017  TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
1018  x, "Error, the input value to squareroot(...) x = " << x << " can not be NaN!" );
1019 #endif
1020  return ::sqrt(x);
1021  }
1022  static inline qd_real pow(qd_real x, qd_real y) { return ::pow(x,y); }
1023  static inline qd_real pi() { return 3.14159265358979323846; }
1024  // qd_real puts its transcendental functions in the global namespace.
1025  static inline qd_real log(qd_real x) { return ::log(x); }
1026  static inline qd_real log10(qd_real x) { return ::log10(x); }
1027 };
1028 
1029 
1030 #endif // HAVE_TEUCHOS_QD
1031 
1032 
1033 #ifdef HAVE_TEUCHOS_GNU_MP
1034 
1035 
1036 extern gmp_randclass gmp_rng;
1037 
1038 
1039 /* Regarding halfPrecision, doublePrecision and mpf_class:
1040  Because the precision of an mpf_class float is not determined by the data type,
1041  there is no way to fill the typedefs for this object.
1042 
1043  Instead, we could create new data classes (e.g., Teuchos::MPF128, Teuchos::MPF256) for
1044  commonly used levels of precision, and fill out ScalarTraits for these. This would allow us
1045  to typedef the promotions and demotions in the appropriate way. These classes would serve to
1046  wrap an mpf_class object, calling the constructor for the appropriate precision, exposing the
1047  arithmetic routines but hiding the precision-altering routines.
1048 
1049  Alternatively (perhaps, preferably), would could create a single class templated on the precision (e.g., Teuchos::MPF<N>).
1050  Then we have a single (partially-specialized) implementation of ScalarTraits. This class, as above, must expose all of the
1051  operations expected of a scalar type; however, most of these can be trivially stolen from the gmpcxx.h header file
1052 
1053  CGB/RAB, 01/05/2009
1054 */
1055 template<>
1056 struct ScalarTraits<mpf_class>
1057 {
1058  typedef mpf_class magnitudeType;
1059  typedef mpf_class halfPrecision;
1060  typedef mpf_class doublePrecision;
1061  typedef mpf_class coordinateType;
1062  static const bool isComplex = false;
1063  static const bool hasMachineParameters = false;
1064  // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
1065  static magnitudeType magnitude(mpf_class a) { return std::abs(a); }
1066  static inline mpf_class zero() { mpf_class zero = 0.0; return zero; }
1067  static inline mpf_class one() { mpf_class one = 1.0; return one; }
1068  static inline mpf_class conjugate(mpf_class x) { return x; }
1069  static inline mpf_class real(mpf_class x) { return(x); }
1070  static inline mpf_class imag(mpf_class x) { return(0); }
1071  static inline bool isnaninf(mpf_class x) { return false; } // mpf_class currently can't handle nan or inf!
1072  static inline void seedrandom(unsigned int s) {
1073  unsigned long int seedVal = static_cast<unsigned long int>(s);
1074  gmp_rng.seed( seedVal );
1075  }
1076  static inline mpf_class random() {
1077  return gmp_rng.get_f();
1078  }
1079  static inline std::string name() { return "mpf_class"; }
1080  static inline mpf_class squareroot(mpf_class x) { return std::sqrt(x); }
1081  static inline mpf_class pow(mpf_class x, mpf_class y) { return pow(x,y); }
1082  // Todo: RAB: 2004/05/28: Add nan() and isnaninf() functions when needed!
1083 };
1084 
1085 #endif // HAVE_TEUCHOS_GNU_MP
1086 
1087 #ifdef HAVE_TEUCHOS_ARPREC
1088 
1089 /* See discussion above for mpf_class, regarding halfPrecision and doublePrecision. Something similar will need to be done
1090  for ARPREC. */
1091 template<>
1092 struct ScalarTraits<mp_real>
1093 {
1094  typedef mp_real magnitudeType;
1095  typedef double halfPrecision;
1096  typedef mp_real doublePrecision;
1097  typedef mp_real coordinateType;
1098  static const bool isComplex = false;
1099  static const bool isComparable = true;
1100  static const bool isOrdinal = false;
1101  static const bool hasMachineParameters = false;
1102  // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
1103  static magnitudeType magnitude(mp_real a) { return abs(a); }
1104  static inline mp_real zero() { mp_real zero = 0.0; return zero; }
1105  static inline mp_real one() { mp_real one = 1.0; return one; }
1106  static inline mp_real conjugate(mp_real x) { return x; }
1107  static inline mp_real real(mp_real x) { return(x); }
1108  static inline mp_real imag(mp_real x) { return zero(); }
1109  static inline bool isnaninf(mp_real x) { return false; } // ToDo: Change this?
1110  static inline void seedrandom(unsigned int s) {
1111  long int seedVal = static_cast<long int>(s);
1112  srand48(seedVal);
1113  }
1114  static inline mp_real random() { return mp_rand(); }
1115  static inline std::string name() { return "mp_real"; }
1116  static inline mp_real squareroot(mp_real x) { return sqrt(x); }
1117  static inline mp_real pow(mp_real x, mp_real y) { return pow(x,y); }
1118  static inline mp_real pi() { return 3.14159265358979323846; }
1119  // Todo: RAB: 2004/05/28: Add nan() and isnaninf() functions when needed!
1120 };
1121 
1122 
1123 #endif // HAVE_TEUCHOS_ARPREC
1124 
1125 
1126 #ifdef HAVE_TEUCHOS_COMPLEX
1127 
1128 
1129 // Partial specialization for std::complex numbers templated on real type T
1130 template<class T>
1131 struct ScalarTraits<
1132  std::complex<T>
1133 >
1134 {
1135  typedef std::complex<T> ComplexT;
1136  typedef std::complex<typename ScalarTraits<T>::halfPrecision> halfPrecision;
1137  typedef std::complex<typename ScalarTraits<T>::doublePrecision> doublePrecision;
1138  typedef typename ScalarTraits<T>::magnitudeType magnitudeType;
1139  typedef typename ScalarTraits<T>::coordinateType coordinateType;
1140  static const bool isComplex = true;
1141  static const bool isOrdinal = ScalarTraits<T>::isOrdinal;
1142  static const bool isComparable = false;
1143  static const bool hasMachineParameters = true;
1144  static inline magnitudeType eps() { return ScalarTraits<magnitudeType>::eps(); }
1145  static inline magnitudeType sfmin() { return ScalarTraits<magnitudeType>::sfmin(); }
1146  static inline magnitudeType base() { return ScalarTraits<magnitudeType>::base(); }
1147  static inline magnitudeType prec() { return ScalarTraits<magnitudeType>::prec(); }
1148  static inline magnitudeType t() { return ScalarTraits<magnitudeType>::t(); }
1149  static inline magnitudeType rnd() { return ScalarTraits<magnitudeType>::rnd(); }
1150  static inline magnitudeType emin() { return ScalarTraits<magnitudeType>::emin(); }
1151  static inline magnitudeType rmin() { return ScalarTraits<magnitudeType>::rmin(); }
1152  static inline magnitudeType emax() { return ScalarTraits<magnitudeType>::emax(); }
1153  static inline magnitudeType rmax() { return ScalarTraits<magnitudeType>::rmax(); }
1154  static magnitudeType magnitude(ComplexT a)
1155  {
1156 #ifdef TEUCHOS_DEBUG
1157  TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
1158  a, "Error, the input value to magnitude(...) a = " << a << " can not be NaN!" );
1159 #endif
1160  return std::abs(a);
1161  }
1162  static inline ComplexT zero() { return ComplexT(ScalarTraits<magnitudeType>::zero(),ScalarTraits<magnitudeType>::zero()); }
1163  static inline ComplexT one() { return ComplexT(ScalarTraits<magnitudeType>::one(),ScalarTraits<magnitudeType>::zero()); }
1164  static inline ComplexT conjugate(ComplexT a){ return ComplexT(a.real(),-a.imag()); }
1165  static inline magnitudeType real(ComplexT a) { return a.real(); }
1166  static inline magnitudeType imag(ComplexT a) { return a.imag(); }
1167  static inline ComplexT nan() { return ComplexT(ScalarTraits<magnitudeType>::nan(),ScalarTraits<magnitudeType>::nan()); }
1168  static inline bool isnaninf(ComplexT x) { return ScalarTraits<magnitudeType>::isnaninf(x.real()) || ScalarTraits<magnitudeType>::isnaninf(x.imag()); }
1169  static inline void seedrandom(unsigned int s) { ScalarTraits<magnitudeType>::seedrandom(s); }
1170  static inline ComplexT random()
1171  {
1172  const T rnd1 = ScalarTraits<magnitudeType>::random();
1173  const T rnd2 = ScalarTraits<magnitudeType>::random();
1174  return ComplexT(rnd1,rnd2);
1175  }
1176  static inline std::string name() { return std::string("std::complex<")+std::string(ScalarTraits<magnitudeType>::name())+std::string(">"); }
1177  // This will only return one of the square roots of x, the other can be obtained by taking its conjugate
1178  static inline ComplexT squareroot(ComplexT x)
1179  {
1180 #ifdef TEUCHOS_DEBUG
1181  TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
1182  x, "Error, the input value to squareroot(...) x = " << x << " can not be NaN!" );
1183 #endif
1184  typedef ScalarTraits<magnitudeType> STMT;
1185  const T r = x.real(), i = x.imag(), zero = STMT::zero(), two = 2.0;
1186  const T a = STMT::squareroot((r*r)+(i*i));
1187  const T nr = STMT::squareroot((a+r)/two);
1188  const T ni = ( i == zero ? zero : STMT::squareroot((a-r)/two) );
1189  return ComplexT(nr,ni);
1190  }
1191  // 2010/03/19: rabartl: Above, I had to add the check for i == zero
1192  // to avoid a returned NaN on the Intel 10.1 compiler. For some
1193  // reason, having these two squareroot calls in a row produce a NaN
1194  // in an optimized build with this compiler. Amazingly, when I put
1195  // in print statements (i.e. std::cout << ...) the NaN went away and
1196  // the second squareroot((a-r)/two) returned zero correctly. I
1197  // guess that calling the output routine flushed the registers or
1198  // something and restarted the squareroot rountine for this compiler
1199  // or something similar. Actually, due to roundoff, it is possible that a-r
1200  // might be slightly less than zero (i.e. -1e-16) and that would cause
1201  // a possbile NaN return. THe above if test is the right thing to do
1202  // I think and is very cheap.
1203  static inline ComplexT pow(ComplexT x, ComplexT y) { return pow(x,y); }
1204  static inline ComplexT pi() { return ScalarTraits<T>::pi(); }
1205 };
1206 #endif // HAVE_TEUCHOS_COMPLEX
1207 
1208 #ifdef HAVE_TEUCHOSCORE_KOKKOSCORE
1209 // Partial specialization for Kokkos::complex<T>
1210 template<class T>
1211 struct ScalarTraits<
1212  Kokkos::complex<T>
1213 >
1214 {
1215  typedef Kokkos::complex<T> ComplexT;
1216  typedef Kokkos::complex<typename ScalarTraits<T>::halfPrecision> halfPrecision;
1217  typedef Kokkos::complex<typename ScalarTraits<T>::doublePrecision> doublePrecision;
1218  typedef typename ScalarTraits<T>::magnitudeType magnitudeType;
1219  typedef typename ScalarTraits<T>::coordinateType coordinateType;
1220  static const bool isComplex = true;
1221  static const bool isOrdinal = ScalarTraits<T>::isOrdinal;
1222  static const bool isComparable = false;
1223  static const bool hasMachineParameters = true;
1224  static inline magnitudeType eps() { return ScalarTraits<magnitudeType>::eps(); }
1225  static inline magnitudeType sfmin() { return ScalarTraits<magnitudeType>::sfmin(); }
1226  static inline magnitudeType base() { return ScalarTraits<magnitudeType>::base(); }
1227  static inline magnitudeType prec() { return ScalarTraits<magnitudeType>::prec(); }
1228  static inline magnitudeType t() { return ScalarTraits<magnitudeType>::t(); }
1229  static inline magnitudeType rnd() { return ScalarTraits<magnitudeType>::rnd(); }
1230  static inline magnitudeType emin() { return ScalarTraits<magnitudeType>::emin(); }
1231  static inline magnitudeType rmin() { return ScalarTraits<magnitudeType>::rmin(); }
1232  static inline magnitudeType emax() { return ScalarTraits<magnitudeType>::emax(); }
1233  static inline magnitudeType rmax() { return ScalarTraits<magnitudeType>::rmax(); }
1234  static magnitudeType magnitude(ComplexT a)
1235  {
1236 #ifdef TEUCHOS_DEBUG
1237  TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
1238  a, "Error, the input value to magnitude(...) a = " << a << " can not be NaN!" );
1239 #endif
1240  return std::abs(std::complex<T>(a));
1241  }
1242  static inline ComplexT zero() { return ComplexT(ScalarTraits<magnitudeType>::zero(),ScalarTraits<magnitudeType>::zero()); }
1243  static inline ComplexT one() { return ComplexT(ScalarTraits<magnitudeType>::one(),ScalarTraits<magnitudeType>::zero()); }
1244  static inline ComplexT conjugate(ComplexT a){ return ComplexT(a.real(),-a.imag()); }
1245  static inline magnitudeType real(ComplexT a) { return a.real(); }
1246  static inline magnitudeType imag(ComplexT a) { return a.imag(); }
1247  static inline ComplexT nan() { return ComplexT(ScalarTraits<magnitudeType>::nan(),ScalarTraits<magnitudeType>::nan()); }
1248  static inline bool isnaninf(ComplexT x) { return ScalarTraits<magnitudeType>::isnaninf(x.real()) || ScalarTraits<magnitudeType>::isnaninf(x.imag()); }
1249  static inline void seedrandom(unsigned int s) { ScalarTraits<magnitudeType>::seedrandom(s); }
1250  static inline ComplexT random()
1251  {
1252  const T rnd1 = ScalarTraits<magnitudeType>::random();
1253  const T rnd2 = ScalarTraits<magnitudeType>::random();
1254  return ComplexT(rnd1,rnd2);
1255  }
1256  static inline std::string name() { return std::string("Kokkos::complex<")+std::string(ScalarTraits<magnitudeType>::name())+std::string(">"); }
1257  // This will only return one of the square roots of x, the other can be obtained by taking its conjugate
1258  static inline ComplexT squareroot(ComplexT x)
1259  {
1260 #ifdef TEUCHOS_DEBUG
1261  TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
1262  x, "Error, the input value to squareroot(...) x = " << x << " can not be NaN!" );
1263 #endif
1264  typedef ScalarTraits<magnitudeType> STMT;
1265  const T r = x.real(), i = x.imag(), zero = STMT::zero(), two = 2.0;
1266  const T a = STMT::squareroot((r*r)+(i*i));
1267  const T nr = STMT::squareroot((a+r)/two);
1268  const T ni = ( i == zero ? zero : STMT::squareroot((a-r)/two) );
1269  return ComplexT(nr,ni);
1270  }
1271  static inline ComplexT pow(ComplexT x, ComplexT y) { return pow(std::complex<T>(x), std::complex<T>(y)); }
1272  static inline ComplexT pi() { return ScalarTraits<T>::pi(); }
1273 };
1274 #endif // HAVE_TEUCHOSCORE_KOKKOSCORE
1275 
1276 #endif // DOXYGEN_SHOULD_SKIP_THIS
1277 
1278 } // Teuchos namespace
1279 
1280 #endif // _TEUCHOS_SCALARTRAITS_HPP_
Teuchos::ScalarTraits::magnitudeType
T magnitudeType
Mandatory typedef for result of magnitude.
Definition: Teuchos_ScalarTraitsDecl.hpp:93
Teuchos::ScalarTraits::eps
static magnitudeType eps()
Returns relative machine precision.
Definition: Teuchos_ScalarTraitsDecl.hpp:112
Teuchos::ScalarTraits::squareroot
static T squareroot(T x)
Returns a number of magnitudeType that is the square root of this scalar type x.
Definition: Teuchos_ScalarTraitsDecl.hpp:154
Teuchos::ScalarTraits::isComparable
static const bool isComparable
Determines if scalar type supports relational operators such as &lt;, &gt;, &lt;=, &gt;=.
Definition: Teuchos_ScalarTraitsDecl.hpp:105
Teuchos::ScalarTraits::real
static magnitudeType real(T a)
Returns the real part of the scalar type a.
Definition: Teuchos_ScalarTraitsDecl.hpp:138
Teuchos_ConfigDefs.hpp
Teuchos header file which uses auto-configuration information to include necessary C++ headers...
Teuchos_ScalarTraitsDecl.hpp
Declaration and default implementation for basic traits for the scalar field type.
Teuchos::ScalarTraits::pow
static T pow(T x, T y)
Returns the result of raising one scalar x to the power y.
Definition: Teuchos_ScalarTraitsDecl.hpp:156
Teuchos::ScalarTraits::emax
static magnitudeType emax()
Returns the largest exponent before overflow.
Definition: Teuchos_ScalarTraitsDecl.hpp:128
Teuchos::ScalarTraits::base
static magnitudeType base()
Returns the base of the machine.
Definition: Teuchos_ScalarTraitsDecl.hpp:116
Teuchos::ScalarTraits::hasMachineParameters
static const bool hasMachineParameters
Determines if scalar type have machine-specific parameters (i.e. eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax() are supported).
Definition: Teuchos_ScalarTraitsDecl.hpp:110
Teuchos::ScalarTraits::name
static std::string name()
Returns the name of this scalar type.
Definition: Teuchos_ScalarTraitsDecl.hpp:152
Teuchos::ScalarTraits::rmax
static magnitudeType rmax()
Overflow theshold - (base^emax)*(1-eps)
Definition: Teuchos_ScalarTraitsDecl.hpp:130
Teuchos::ScalarTraits::sfmin
static magnitudeType sfmin()
Returns safe minimum (sfmin), such that 1/sfmin does not overflow.
Definition: Teuchos_ScalarTraitsDecl.hpp:114
Teuchos::ScalarTraits::pi
static T pi()
Returns the value of PI.
Definition: Teuchos_ScalarTraitsDecl.hpp:158
Teuchos::ScalarTraits::coordinateType
T coordinateType
Typedef for coordinates.
Definition: Teuchos_ScalarTraitsDecl.hpp:99
Teuchos::ScalarTraits::conjugate
static T conjugate(T a)
Returns the conjugate of the scalar type a.
Definition: Teuchos_ScalarTraitsDecl.hpp:142
Teuchos::ScalarTraits::isOrdinal
static const bool isOrdinal
Determines if scalar type is an ordinal type.
Definition: Teuchos_ScalarTraitsDecl.hpp:103
Teuchos::ScalarTraits::prec
static magnitudeType prec()
Returns eps*base.
Definition: Teuchos_ScalarTraitsDecl.hpp:118
Teuchos::ScalarTraits::t
static magnitudeType t()
Returns the number of (base) digits in the mantissa.
Definition: Teuchos_ScalarTraitsDecl.hpp:120
Teuchos::ScalarTraits::rmin
static magnitudeType rmin()
Returns the underflow threshold - base^(emin-1)
Definition: Teuchos_ScalarTraitsDecl.hpp:126
Teuchos::ScalarTraits::doublePrecision
T doublePrecision
Typedef for double precision.
Definition: Teuchos_ScalarTraitsDecl.hpp:97
Teuchos::ScalarTraits::seedrandom
static void seedrandom(unsigned int s)
Seed the random number generator returned by random().
Definition: Teuchos_ScalarTraitsDecl.hpp:148
Teuchos::ScalarTraits::imag
static magnitudeType imag(T a)
Returns the imaginary part of the scalar type a.
Definition: Teuchos_ScalarTraitsDecl.hpp:140
Teuchos::ScalarTraits::isnaninf
static bool isnaninf(const T &x)
Returns true if x is NaN or Inf.
Definition: Teuchos_ScalarTraitsDecl.hpp:146
Teuchos::ScalarTraits::magnitude
static magnitudeType magnitude(T a)
Returns the magnitudeType of the scalar type a.
Definition: Teuchos_ScalarTraitsDecl.hpp:132
Teuchos::ScalarTraits::nan
static T nan()
Returns a number that represents NaN.
Definition: Teuchos_ScalarTraitsDecl.hpp:144
Teuchos::ScalarTraits::zero
static T zero()
Returns representation of zero for this scalar type.
Definition: Teuchos_ScalarTraitsDecl.hpp:134
Teuchos::ScalarTraits::random
static T random()
Returns a random number (between -one() and +one()) of this scalar type.
Definition: Teuchos_ScalarTraitsDecl.hpp:150
Teuchos::ScalarTraits::halfPrecision
T halfPrecision
Typedef for half precision.
Definition: Teuchos_ScalarTraitsDecl.hpp:95
Teuchos::ScalarTraits::isComplex
static const bool isComplex
Determines if scalar type is std::complex.
Definition: Teuchos_ScalarTraitsDecl.hpp:101
Teuchos::ScalarTraits::emin
static magnitudeType emin()
Returns the minimum exponent before (gradual) underflow.
Definition: Teuchos_ScalarTraitsDecl.hpp:124
Teuchos::ScalarTraits::rnd
static magnitudeType rnd()
Returns 1.0 when rounding occurs in addition, 0.0 otherwise.
Definition: Teuchos_ScalarTraitsDecl.hpp:122
Teuchos::ScalarTraits::one
static T one()
Returns representation of one for this scalar type.
Definition: Teuchos_ScalarTraitsDecl.hpp:136
Generated by   doxygen 1.8.5