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