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 #ifdef HAVE_TEUCHOS_LONG_LONG_INT
409 template<>
410 struct ScalarTraits<long long int>
411 {
412  typedef long long int magnitudeType;
413  typedef long long int halfPrecision;
414  typedef long long int doublePrecision;
415  static const bool isComplex = false;
416  static const bool isOrdinal = true;
417  static const bool isComparable = true;
418  static const bool hasMachineParameters = false;
419  // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
420  static inline magnitudeType magnitude(long long int a) { return static_cast<long long int>(std::fabs(static_cast<double>(a))); }
421  static inline long long int zero() { return 0; }
422  static inline long long int one() { return 1; }
423  static inline long long int conjugate(long long int x) { return x; }
424  static inline long long int real(long long int x) { return x; }
425  static inline long long int imag(long long int) { return 0; }
426  static inline bool isnaninf(long long int) { return false; }
427  static inline void seedrandom(unsigned int s) {
428  std::srand(s);
429 #ifdef __APPLE__
430  // throw away first random number to address bug 3655
431  // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
432  random();
433 #endif
434  }
435  //static inline int random() { return (-1 + 2*rand()); } // RAB: This version should be used to be consistent with others
436  static inline long long int random() { return std::rand(); } // RAB: This version should be used for an unsigned int, not int
437  static inline std::string name() { return "long long int"; }
438  static inline long long int squareroot(long long int x) { return (long long int) std::sqrt((double) x); }
439  static inline long long int pow(long long int x, long long int y) { return (long long int) std::pow((double)x,(double)y); }
440  // Note: Depending on the number of bits in long long int, the cast
441  // from long long int to double may not be exact.
442  static inline long long int log(long long int x) { return static_cast<long long int> (std::log (static_cast<double> (x))); }
443  static inline long long int log10(long long int x) { return static_cast<long long int> (std::log10 (static_cast<double> (x))); }
444 };
445 
446 template<>
447 struct ScalarTraits<unsigned long long int>
448 {
449  typedef unsigned long long int magnitudeType;
450  typedef unsigned long long int halfPrecision;
451  typedef unsigned long long int doublePrecision;
452  static const bool isComplex = false;
453  static const bool isOrdinal = true;
454  static const bool isComparable = true;
455  static const bool hasMachineParameters = false;
456  // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
457  static inline magnitudeType magnitude(unsigned long long int a) { return static_cast<unsigned long long int>(std::fabs(static_cast<double>(a))); }
458  static inline unsigned long long int zero() { return 0; }
459  static inline unsigned long long int one() { return 1; }
460  static inline unsigned long long int conjugate(unsigned long long int x) { return x; }
461  static inline unsigned long long int real(unsigned long long int x) { return x; }
462  static inline unsigned long long int imag(unsigned long long int) { return 0; }
463  static inline bool isnaninf(unsigned long long int) { return false; }
464  static inline void seedrandom(unsigned int s) {
465  std::srand(s);
466 #ifdef __APPLE__
467  // throw away first random number to address bug 3655
468  // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
469  random();
470 #endif
471  }
472  //static inline int random() { return (-1 + 2*rand()); } // RAB: This version should be used to be consistent with others
473  static inline unsigned long long int random() { return std::rand(); } // RAB: This version should be used for an unsigned int, not int
474  static inline std::string name() { return "unsigned long long int"; }
475  static inline unsigned long long int squareroot(unsigned long long int x) { return (unsigned long long int) std::sqrt((double) x); }
476  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); }
477  // Note: Depending on the number of bits in unsigned long long int,
478  // the cast from unsigned long long int to double may not be exact.
479  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))); }
480  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))); }
481 };
482 #endif // HAVE_TEUCHOS_LONG_LONG_INT
483 
484 
485 #ifdef HAVE_TEUCHOS___INT64
486 
487 template<>
488 struct ScalarTraits<__int64>
489 {
490  typedef __int64 magnitudeType;
491  typedef __int64 halfPrecision;
492  typedef __int64 doublePrecision;
493  static const bool isComplex = false;
494  static const bool isOrdinal = true;
495  static const bool isComparable = true;
496  static const bool hasMachineParameters = false;
497  // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
498  static inline magnitudeType magnitude(__int64 a) { return static_cast<__int64>(std::fabs(static_cast<double>(a))); }
499  static inline __int64 zero() { return 0; }
500  static inline __int64 one() { return 1; }
501  static inline __int64 conjugate(__int64 x) { return x; }
502  static inline __int64 real(__int64 x) { return x; }
503  static inline __int64 imag(__int64) { return 0; }
504  static inline void seedrandom(unsigned int s) {
505  std::srand(s);
506 #ifdef __APPLE__
507  // throw away first random number to address bug 3655
508  // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
509  random();
510 #endif
511  }
512  //static inline int random() { return (-1 + 2*rand()); } // RAB: This version should be used to be consistent with others
513  static inline __int64 random() { return std::rand(); } // RAB: This version should be used for an unsigned int, not int
514  static inline std::string name() { return "__int64"; }
515  static inline __int64 squareroot(__int64 x) { return (__int64) std::sqrt((double) x); }
516  static inline __int64 pow(__int64 x, __int64 y) { return (__int64) std::pow((double)x,(double)y); }
517  // Note: Depending on the number of bits in __int64, the cast
518  // from __int64 to double may not be exact.
519  static inline __int64 log(__int64 x) { return static_cast<__int64> (std::log (static_cast<double> (x))); }
520  static inline __int64 log10(__int64 x) { return static_cast<__int64> (std::log10 (static_cast<double> (x))); }
521 };
522 
523 template<>
524 struct ScalarTraits<unsigned __int64>
525 {
526  typedef unsigned __int64 magnitudeType;
527  typedef unsigned __int64 halfPrecision;
528  typedef unsigned __int64 doublePrecision;
529  static const bool isComplex = false;
530  static const bool isOrdinal = true;
531  static const bool isComparable = true;
532  static const bool hasMachineParameters = false;
533  // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
534  static inline magnitudeType magnitude(unsigned __int64 a) { return static_cast<unsigned __int64>(std::fabs(static_cast<double>(a))); }
535  static inline unsigned __int64 zero() { return 0; }
536  static inline unsigned __int64 one() { return 1; }
537  static inline unsigned __int64 conjugate(unsigned __int64 x) { return x; }
538  static inline unsigned __int64 real(unsigned __int64 x) { return x; }
539  static inline unsigned __int64 imag(unsigned __int64) { return 0; }
540  static inline void seedrandom(unsigned int s) {
541  std::srand(s);
542 #ifdef __APPLE__
543  // throw away first random number to address bug 3655
544  // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
545  random();
546 #endif
547  }
548  //static inline int random() { return (-1 + 2*rand()); } // RAB: This version should be used to be consistent with others
549  static inline unsigned __int64 random() { return std::rand(); } // RAB: This version should be used for an unsigned int, not int
550  static inline std::string name() { return "unsigned __int64"; }
551  static inline unsigned __int64 squareroot(unsigned __int64 x) { return (unsigned __int64) std::sqrt((double) x); }
552  static inline unsigned __int64 pow(unsigned __int64 x, unsigned __int64 y) { return (unsigned __int64) std::pow((double)x,(double)y); }
553  // Note: Depending on the number of bits in unsigned __int64,
554  // the cast from unsigned __int64 to double may not be exact.
555  static inline unsigned __int64 log(unsigned __int64 x) { return static_cast<unsigned __int64> (std::log (static_cast<double> (x))); }
556  static inline unsigned __int64 log10(unsigned __int64 x) { return static_cast<unsigned __int64> (std::log10 (static_cast<double> (x))); }
557 };
558 
559 #endif // HAVE_TEUCHOS___INT64
560 
561 
562 #ifndef __sun
563 extern TEUCHOSCORE_LIB_DLL_EXPORT const float flt_nan;
564 #endif
565 
566 
567 template<>
568 struct ScalarTraits<float>
569 {
570  typedef float magnitudeType;
571  typedef float halfPrecision; // should become IEEE754-2008 binary16 or fp16 later, perhaps specified at configure according to architectural support
572  typedef double doublePrecision;
573  static const bool isComplex = false;
574  static const bool isOrdinal = false;
575  static const bool isComparable = true;
576  static const bool hasMachineParameters = true;
577  static inline float eps() {
578  return std::numeric_limits<float>::epsilon();
579  }
580  static inline float sfmin() {
581  return std::numeric_limits<float>::min();
582  }
583  static inline float base() {
584  return static_cast<float>(std::numeric_limits<float>::radix);
585  }
586  static inline float prec() {
587  return eps()*base();
588  }
589  static inline float t() {
590  return static_cast<float>(std::numeric_limits<float>::digits);
591  }
592  static inline float rnd() {
593  return ( std::numeric_limits<float>::round_style == std::round_to_nearest ? one() : zero() );
594  }
595  static inline float emin() {
596  return static_cast<float>(std::numeric_limits<float>::min_exponent);
597  }
598  static inline float rmin() {
599  return std::numeric_limits<float>::min();
600  }
601  static inline float emax() {
602  return static_cast<float>(std::numeric_limits<float>::max_exponent);
603  }
604  static inline float rmax() {
605  return std::numeric_limits<float>::max();
606  }
607  static inline magnitudeType magnitude(float a)
608  {
609 #ifdef TEUCHOS_DEBUG
610  TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
611  a, "Error, the input value to magnitude(...) a = " << a << " can not be NaN!" );
612 #endif
613  return std::fabs(a);
614  }
615  static inline float zero() { return(0.0f); }
616  static inline float one() { return(1.0f); }
617  static inline float conjugate(float x) { return(x); }
618  static inline float real(float x) { return x; }
619  static inline float imag(float) { return zero(); }
620  static inline float nan() {
621 #ifdef __sun
622  return 0.0f/std::sin(0.0f);
623 #else
624  return flt_nan;
625 #endif
626  }
627  static inline bool isnaninf(float x) {
628  return generic_real_isnaninf<float>(x);
629  }
630  static inline void seedrandom(unsigned int s) {
631  std::srand(s);
632 #ifdef __APPLE__
633  // throw away first random number to address bug 3655
634  // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
635  random();
636 #endif
637  }
638  static inline float random() { float rnd = (float) std::rand() / RAND_MAX; return (-1.0f + 2.0f * rnd); }
639  static inline std::string name() { return "float"; }
640  static inline float squareroot(float x)
641  {
642 #ifdef TEUCHOS_DEBUG
643  TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
644  x, "Error, the input value to squareroot(...) x = " << x << " can not be NaN!" );
645 #endif
646  errno = 0;
647  const float rtn = std::sqrt(x);
648  if (errno)
649  return nan();
650  return rtn;
651  }
652  static inline float pow(float x, float y) { return std::pow(x,y); }
653  static inline float pi() { return 3.14159265358979323846f; }
654  static inline float log(float x) { return std::log(x); }
655  static inline float log10(float x) { return std::log10(x); }
656 };
657 
658 
659 #ifndef __sun
660 extern TEUCHOSCORE_LIB_DLL_EXPORT const double dbl_nan;
661 #endif
662 
663 
664 template<>
665 struct ScalarTraits<double>
666 {
667  typedef double magnitudeType;
668  typedef float halfPrecision;
669  /* there are different options as to how to double "double"
670  - QD's DD (if available)
671  - ARPREC
672  - GNU MP
673  - a true hardware quad
674 
675  in the shortterm, this should be specified at configure time. I have inserted a configure-time option (--enable-teuchos-double-to-dd)
676  which uses QD's DD when available. This must be used alongside --enable-teuchos-qd.
677  */
678 #if defined(HAVE_TEUCHOS_DOUBLE_TO_QD)
679  typedef dd_real doublePrecision;
680 #elif defined(HAVE_TEUCHOS_DOUBLE_TO_ARPREC)
681  typedef mp_real doublePrecision;
682 #else
683  typedef double doublePrecision; // don't double "double" in this case
684 #endif
685  static const bool isComplex = false;
686  static const bool isOrdinal = false;
687  static const bool isComparable = true;
688  static const bool hasMachineParameters = true;
689  static inline double eps() {
690  return std::numeric_limits<double>::epsilon();
691  }
692  static inline double sfmin() {
693  return std::numeric_limits<double>::min();
694  }
695  static inline double base() {
696  return std::numeric_limits<double>::radix;
697  }
698  static inline double prec() {
699  return eps()*base();
700  }
701  static inline double t() {
702  return std::numeric_limits<double>::digits;
703  }
704  static inline double rnd() {
705  return ( std::numeric_limits<double>::round_style == std::round_to_nearest ? double(1.0) : double(0.0) );
706  }
707  static inline double emin() {
708  return std::numeric_limits<double>::min_exponent;
709  }
710  static inline double rmin() {
711  return std::numeric_limits<double>::min();
712  }
713  static inline double emax() {
714  return std::numeric_limits<double>::max_exponent;
715  }
716  static inline double rmax() {
717  return std::numeric_limits<double>::max();
718  }
719  static inline magnitudeType magnitude(double a)
720  {
721 #ifdef TEUCHOS_DEBUG
722  TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
723  a, "Error, the input value to magnitude(...) a = " << a << " can not be NaN!" );
724 #endif
725  return std::fabs(a);
726  }
727  static inline double zero() { return 0.0; }
728  static inline double one() { return 1.0; }
729  static inline double conjugate(double x) { return(x); }
730  static inline double real(double x) { return(x); }
731  static inline double imag(double) { return(0); }
732  static inline double nan() {
733 #ifdef __sun
734  return 0.0/std::sin(0.0);
735 #else
736  return dbl_nan;
737 #endif
738  }
739  static inline bool isnaninf(double x) {
740  return generic_real_isnaninf<double>(x);
741  }
742  static inline void seedrandom(unsigned int s) {
743  std::srand(s);
744 #ifdef __APPLE__
745  // throw away first random number to address bug 3655
746  // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
747  random();
748 #endif
749  }
750  static inline double random() { double rnd = (double) std::rand() / RAND_MAX; return (double)(-1.0 + 2.0 * rnd); }
751  static inline std::string name() { return "double"; }
752  static inline double squareroot(double x)
753  {
754 #ifdef TEUCHOS_DEBUG
755  TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
756  x, "Error, the input value to squareroot(...) x = " << x << " can not be NaN!" );
757 #endif
758  errno = 0;
759  const double rtn = std::sqrt(x);
760  if (errno)
761  return nan();
762  return rtn;
763  }
764  static inline double pow(double x, double y) { return std::pow(x,y); }
765  static inline double pi() { return 3.14159265358979323846; }
766  static inline double log(double x) { return std::log(x); }
767  static inline double log10(double x) { return std::log10(x); }
768 };
769 
770 
771 #ifdef HAVE_TEUCHOSCORE_QUADMATH
772 
773 template<>
774 struct ScalarTraits<__float128> {
775  typedef __float128 magnitudeType;
776  // Unfortunately, we can't rely on a standard __float256 type. It
777  // might make sense to use qd_real, but mixing quadmath and QD might
778  // cause unforeseen issues.
779  typedef __float128 doublePrecision;
780  typedef double halfPrecision;
781 
782  static const bool isComplex = false;
783  static const bool isOrdinal = false;
784  static const bool isComparable = true;
785  static const bool hasMachineParameters = true;
786 
787  static __float128 eps () {
788  return FLT128_EPSILON;
789  }
790  static __float128 sfmin () {
791  return FLT128_MIN; // ???
792  }
793  static __float128 base () {
794  return 2;
795  }
796  static __float128 prec () {
797  return eps () * base ();
798  }
799  static __float128 t () {
800  return FLT128_MANT_DIG;
801  }
802  static __float128 rnd () {
803  return 1.0;
804  }
805  static __float128 emin () {
806  return FLT128_MIN_EXP;
807  }
808  static __float128 rmin () {
809  return FLT128_MIN; // ??? // should be base^(emin-1)
810  }
811  static __float128 emax () {
812  return FLT128_MAX_EXP;
813  }
814  static __float128 rmax () {
815  return FLT128_MAX; // ??? // should be (base^emax)*(1-eps)
816  }
817  static magnitudeType magnitude (const __float128& x) {
818  return fabsq (x);
819  }
820  static __float128 zero () {
821  return 0.0;
822  }
823  static __float128 one () {
824  return 1.0;
825  }
826  static __float128 conjugate (const __float128& x) {
827  return x;
828  }
829  static __float128 real (const __float128& x) {
830  return x;
831  }
832  static __float128 imag (const __float128& /* x */) {
833  return 0.0;
834  }
835  static __float128 nan () {
836  return strtoflt128 ("NAN()", NULL); // ???
837  }
838  static bool isnaninf (const __float128& x) {
839  return isinfq (x) || isnanq (x);
840  }
841  static inline void seedrandom (unsigned int s) {
842  std::srand (s);
843 #ifdef __APPLE__
844  // throw away first random number to address bug 3655
845  // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
846  random ();
847 #endif
848  }
849  static __float128 random () {
850  // Half the smallest normalized double, is the scaling factor of
851  // the lower-order term in the double-double representation.
852  const __float128 scalingFactor =
853  static_cast<__float128> (std::numeric_limits<double>::min ()) /
854  static_cast<__float128> (2.0);
855  const __float128 higherOrderTerm =
856  static_cast<__float128> (ScalarTraits<double>::random ());
857  const __float128 lowerOrderTerm =
858  static_cast<__float128> (ScalarTraits<double>::random ()) *
859  scalingFactor;
860  return higherOrderTerm + lowerOrderTerm;
861  }
862  static std::string name () {
863  return "__float128";
864  }
865  static __float128 squareroot (const __float128& x) {
866  return sqrtq (x);
867  }
868  static __float128 pow (const __float128& x, const __float128& y) {
869  return powq (x, y);
870  }
871  static __float128 pi() { return 3.14159265358979323846; }
872  static __float128 log (const __float128& x) {
873  return logq (x);
874  }
875  static __float128 log10 (const __float128& x) {
876  return log10q (x);
877  }
878 };
879 #endif // HAVE_TEUCHOSCORE_QUADMATH
880 
881 
882 
883 #ifdef HAVE_TEUCHOS_QD
884 
885 bool operator&&(const dd_real &a, const dd_real &b);
886 bool operator&&(const qd_real &a, const qd_real &b);
887 
888 template<>
889 struct ScalarTraits<dd_real>
890 {
891  typedef dd_real magnitudeType;
892  typedef double halfPrecision;
893  typedef qd_real doublePrecision;
894  static const bool isComplex = false;
895  static const bool isOrdinal = false;
896  static const bool isComparable = true;
897  static const bool hasMachineParameters = true;
898  static inline dd_real eps() { return std::numeric_limits<dd_real>::epsilon(); }
899  static inline dd_real sfmin() { return std::numeric_limits<dd_real>::min(); }
900  static inline dd_real base() { return std::numeric_limits<dd_real>::radix; }
901  static inline dd_real prec() { return eps()*base(); }
902  static inline dd_real t() { return std::numeric_limits<dd_real>::digits; }
903  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) ); }
904  static inline dd_real emin() { return std::numeric_limits<dd_real>::min_exponent; }
905  static inline dd_real rmin() { return std::numeric_limits<dd_real>::min(); }
906  static inline dd_real emax() { return std::numeric_limits<dd_real>::max_exponent; }
907  static inline dd_real rmax() { return std::numeric_limits<dd_real>::max(); }
908  static inline magnitudeType magnitude(dd_real a)
909  {
910 #ifdef TEUCHOS_DEBUG
911  TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
912  a, "Error, the input value to magnitude(...) a = " << a << " can not be NaN!" );
913 #endif
914  return ::abs(a);
915  }
916  static inline dd_real zero() { return dd_real(0.0); }
917  static inline dd_real one() { return dd_real(1.0); }
918  static inline dd_real conjugate(dd_real x) { return(x); }
919  static inline dd_real real(dd_real x) { return x ; }
920  static inline dd_real imag(dd_real) { return zero(); }
921  static inline dd_real nan() { return dd_real::_nan; }
922  static inline bool isnaninf(dd_real x) { return isnan(x) || isinf(x); }
923  static inline void seedrandom(unsigned int s) {
924  // ddrand() uses std::rand(), so the std::srand() is our seed
925  std::srand(s);
926 #ifdef __APPLE__
927  // throw away first random number to address bug 3655
928  // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
929  random();
930 #endif
931  }
932  static inline dd_real random() { return ddrand(); }
933  static inline std::string name() { return "dd_real"; }
934  static inline dd_real squareroot(dd_real x)
935  {
936 #ifdef TEUCHOS_DEBUG
937  TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
938  x, "Error, the input value to squareroot(...) x = " << x << " can not be NaN!" );
939 #endif
940  return ::sqrt(x);
941  }
942  static inline dd_real pow(dd_real x, dd_real y) { return ::pow(x,y); }
943  static inline dd_real pi() { return 3.14159265358979323846; }
944  // dd_real puts its transcendental functions in the global namespace.
945  static inline dd_real log(dd_real x) { return ::log(x); }
946  static inline dd_real log10(dd_real x) { return ::log10(x); }
947 };
948 
949 
950 template<>
951 struct ScalarTraits<qd_real>
952 {
953  typedef qd_real magnitudeType;
954  typedef dd_real halfPrecision;
955  typedef qd_real doublePrecision;
956  static const bool isComplex = false;
957  static const bool isOrdinal = false;
958  static const bool isComparable = true;
959  static const bool hasMachineParameters = true;
960  static inline qd_real eps() { return std::numeric_limits<qd_real>::epsilon(); }
961  static inline qd_real sfmin() { return std::numeric_limits<qd_real>::min(); }
962  static inline qd_real base() { return std::numeric_limits<qd_real>::radix; }
963  static inline qd_real prec() { return eps()*base(); }
964  static inline qd_real t() { return std::numeric_limits<qd_real>::digits; }
965  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) ); }
966  static inline qd_real emin() { return std::numeric_limits<qd_real>::min_exponent; }
967  static inline qd_real rmin() { return std::numeric_limits<qd_real>::min(); }
968  static inline qd_real emax() { return std::numeric_limits<qd_real>::max_exponent; }
969  static inline qd_real rmax() { return std::numeric_limits<qd_real>::max(); }
970  static inline magnitudeType magnitude(qd_real a)
971  {
972 #ifdef TEUCHOS_DEBUG
973  TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
974  a, "Error, the input value to magnitude(...) a = " << a << " can not be NaN!" );
975 #endif
976  return ::abs(a);
977  }
978  static inline qd_real zero() { return qd_real(0.0); }
979  static inline qd_real one() { return qd_real(1.0); }
980  static inline qd_real conjugate(qd_real x) { return(x); }
981  static inline qd_real real(qd_real x) { return x ; }
982  static inline qd_real imag(qd_real) { return zero(); }
983  static inline qd_real nan() { return qd_real::_nan; }
984  static inline bool isnaninf(qd_real x) { return isnan(x) || isinf(x); }
985  static inline void seedrandom(unsigned int s) {
986  // qdrand() uses std::rand(), so the std::srand() is our seed
987  std::srand(s);
988 #ifdef __APPLE__
989  // throw away first random number to address bug 3655
990  // http://software.sandia.gov/bugzilla/show_bug.cgi?id=3655
991  random();
992 #endif
993  }
994  static inline qd_real random() { return qdrand(); }
995  static inline std::string name() { return "qd_real"; }
996  static inline qd_real squareroot(qd_real x)
997  {
998 #ifdef TEUCHOS_DEBUG
999  TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
1000  x, "Error, the input value to squareroot(...) x = " << x << " can not be NaN!" );
1001 #endif
1002  return ::sqrt(x);
1003  }
1004  static inline qd_real pow(qd_real x, qd_real y) { return ::pow(x,y); }
1005  static inline qd_real pi() { return 3.14159265358979323846; }
1006  // qd_real puts its transcendental functions in the global namespace.
1007  static inline qd_real log(qd_real x) { return ::log(x); }
1008  static inline qd_real log10(qd_real x) { return ::log10(x); }
1009 };
1010 
1011 
1012 #endif // HAVE_TEUCHOS_QD
1013 
1014 
1015 #ifdef HAVE_TEUCHOS_GNU_MP
1016 
1017 
1018 extern gmp_randclass gmp_rng;
1019 
1020 
1021 /* Regarding halfPrecision, doublePrecision and mpf_class:
1022  Because the precision of an mpf_class float is not determined by the data type,
1023  there is no way to fill the typedefs for this object.
1024 
1025  Instead, we could create new data classes (e.g., Teuchos::MPF128, Teuchos::MPF256) for
1026  commonly used levels of precision, and fill out ScalarTraits for these. This would allow us
1027  to typedef the promotions and demotions in the appropriate way. These classes would serve to
1028  wrap an mpf_class object, calling the constructor for the appropriate precision, exposing the
1029  arithmetic routines but hiding the precision-altering routines.
1030 
1031  Alternatively (perhaps, preferably), would could create a single class templated on the precision (e.g., Teuchos::MPF<N>).
1032  Then we have a single (partially-specialized) implementation of ScalarTraits. This class, as above, must expose all of the
1033  operations expected of a scalar type; however, most of these can be trivially stolen from the gmpcxx.h header file
1034 
1035  CGB/RAB, 01/05/2009
1036 */
1037 template<>
1038 struct ScalarTraits<mpf_class>
1039 {
1040  typedef mpf_class magnitudeType;
1041  typedef mpf_class halfPrecision;
1042  typedef mpf_class doublePrecision;
1043  static const bool isComplex = false;
1044  static const bool hasMachineParameters = false;
1045  // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
1046  static magnitudeType magnitude(mpf_class a) { return std::abs(a); }
1047  static inline mpf_class zero() { mpf_class zero = 0.0; return zero; }
1048  static inline mpf_class one() { mpf_class one = 1.0; return one; }
1049  static inline mpf_class conjugate(mpf_class x) { return x; }
1050  static inline mpf_class real(mpf_class x) { return(x); }
1051  static inline mpf_class imag(mpf_class x) { return(0); }
1052  static inline bool isnaninf(mpf_class x) { return false; } // mpf_class currently can't handle nan or inf!
1053  static inline void seedrandom(unsigned int s) {
1054  unsigned long int seedVal = static_cast<unsigned long int>(s);
1055  gmp_rng.seed( seedVal );
1056  }
1057  static inline mpf_class random() {
1058  return gmp_rng.get_f();
1059  }
1060  static inline std::string name() { return "mpf_class"; }
1061  static inline mpf_class squareroot(mpf_class x) { return std::sqrt(x); }
1062  static inline mpf_class pow(mpf_class x, mpf_class y) { return pow(x,y); }
1063  // Todo: RAB: 2004/05/28: Add nan() and isnaninf() functions when needed!
1064 };
1065 
1066 #endif // HAVE_TEUCHOS_GNU_MP
1067 
1068 #ifdef HAVE_TEUCHOS_ARPREC
1069 
1070 /* See discussion above for mpf_class, regarding halfPrecision and doublePrecision. Something similar will need to be done
1071  for ARPREC. */
1072 template<>
1073 struct ScalarTraits<mp_real>
1074 {
1075  typedef mp_real magnitudeType;
1076  typedef double halfPrecision;
1077  typedef mp_real doublePrecision;
1078  static const bool isComplex = false;
1079  static const bool isComparable = true;
1080  static const bool isOrdinal = false;
1081  static const bool hasMachineParameters = false;
1082  // Not defined: eps(), sfmin(), base(), prec(), t(), rnd(), emin(), rmin(), emax(), rmax()
1083  static magnitudeType magnitude(mp_real a) { return abs(a); }
1084  static inline mp_real zero() { mp_real zero = 0.0; return zero; }
1085  static inline mp_real one() { mp_real one = 1.0; return one; }
1086  static inline mp_real conjugate(mp_real x) { return x; }
1087  static inline mp_real real(mp_real x) { return(x); }
1088  static inline mp_real imag(mp_real x) { return zero(); }
1089  static inline bool isnaninf(mp_real x) { return false; } // ToDo: Change this?
1090  static inline void seedrandom(unsigned int s) {
1091  long int seedVal = static_cast<long int>(s);
1092  srand48(seedVal);
1093  }
1094  static inline mp_real random() { return mp_rand(); }
1095  static inline std::string name() { return "mp_real"; }
1096  static inline mp_real squareroot(mp_real x) { return sqrt(x); }
1097  static inline mp_real pow(mp_real x, mp_real y) { return pow(x,y); }
1098  static inline mp_real pi() { return 3.14159265358979323846; }
1099  // Todo: RAB: 2004/05/28: Add nan() and isnaninf() functions when needed!
1100 };
1101 
1102 
1103 #endif // HAVE_TEUCHOS_ARPREC
1104 
1105 
1106 #ifdef HAVE_TEUCHOS_COMPLEX
1107 
1108 
1109 // Partial specialization for std::complex numbers templated on real type T
1110 template<class T>
1111 struct ScalarTraits<
1112  std::complex<T>
1113 >
1114 {
1115  typedef std::complex<T> ComplexT;
1116  typedef std::complex<typename ScalarTraits<T>::halfPrecision> halfPrecision;
1117  typedef std::complex<typename ScalarTraits<T>::doublePrecision> doublePrecision;
1118  typedef typename ScalarTraits<T>::magnitudeType magnitudeType;
1119  static const bool isComplex = true;
1120  static const bool isOrdinal = ScalarTraits<T>::isOrdinal;
1121  static const bool isComparable = false;
1122  static const bool hasMachineParameters = true;
1123  static inline magnitudeType eps() { return ScalarTraits<magnitudeType>::eps(); }
1124  static inline magnitudeType sfmin() { return ScalarTraits<magnitudeType>::sfmin(); }
1125  static inline magnitudeType base() { return ScalarTraits<magnitudeType>::base(); }
1126  static inline magnitudeType prec() { return ScalarTraits<magnitudeType>::prec(); }
1127  static inline magnitudeType t() { return ScalarTraits<magnitudeType>::t(); }
1128  static inline magnitudeType rnd() { return ScalarTraits<magnitudeType>::rnd(); }
1129  static inline magnitudeType emin() { return ScalarTraits<magnitudeType>::emin(); }
1130  static inline magnitudeType rmin() { return ScalarTraits<magnitudeType>::rmin(); }
1131  static inline magnitudeType emax() { return ScalarTraits<magnitudeType>::emax(); }
1132  static inline magnitudeType rmax() { return ScalarTraits<magnitudeType>::rmax(); }
1133  static magnitudeType magnitude(ComplexT a)
1134  {
1135 #ifdef TEUCHOS_DEBUG
1136  TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
1137  a, "Error, the input value to magnitude(...) a = " << a << " can not be NaN!" );
1138 #endif
1139  return std::abs(a);
1140  }
1141  static inline ComplexT zero() { return ComplexT(ScalarTraits<magnitudeType>::zero(),ScalarTraits<magnitudeType>::zero()); }
1142  static inline ComplexT one() { return ComplexT(ScalarTraits<magnitudeType>::one(),ScalarTraits<magnitudeType>::zero()); }
1143  static inline ComplexT conjugate(ComplexT a){ return ComplexT(a.real(),-a.imag()); }
1144  static inline magnitudeType real(ComplexT a) { return a.real(); }
1145  static inline magnitudeType imag(ComplexT a) { return a.imag(); }
1146  static inline ComplexT nan() { return ComplexT(ScalarTraits<magnitudeType>::nan(),ScalarTraits<magnitudeType>::nan()); }
1147  static inline bool isnaninf(ComplexT x) { return ScalarTraits<magnitudeType>::isnaninf(x.real()) || ScalarTraits<magnitudeType>::isnaninf(x.imag()); }
1148  static inline void seedrandom(unsigned int s) { ScalarTraits<magnitudeType>::seedrandom(s); }
1149  static inline ComplexT random()
1150  {
1151  const T rnd1 = ScalarTraits<magnitudeType>::random();
1152  const T rnd2 = ScalarTraits<magnitudeType>::random();
1153  return ComplexT(rnd1,rnd2);
1154  }
1155  static inline std::string name() { return std::string("std::complex<")+std::string(ScalarTraits<magnitudeType>::name())+std::string(">"); }
1156  // This will only return one of the square roots of x, the other can be obtained by taking its conjugate
1157  static inline ComplexT squareroot(ComplexT x)
1158  {
1159 #ifdef TEUCHOS_DEBUG
1160  TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
1161  x, "Error, the input value to squareroot(...) x = " << x << " can not be NaN!" );
1162 #endif
1163  typedef ScalarTraits<magnitudeType> STMT;
1164  const T r = x.real(), i = x.imag(), zero = STMT::zero(), two = 2.0;
1165  const T a = STMT::squareroot((r*r)+(i*i));
1166  const T nr = STMT::squareroot((a+r)/two);
1167  const T ni = ( i == zero ? zero : STMT::squareroot((a-r)/two) );
1168  return ComplexT(nr,ni);
1169  }
1170  // 2010/03/19: rabartl: Above, I had to add the check for i == zero
1171  // to avoid a returned NaN on the Intel 10.1 compiler. For some
1172  // reason, having these two squareroot calls in a row produce a NaN
1173  // in an optimized build with this compiler. Amazingly, when I put
1174  // in print statements (i.e. std::cout << ...) the NaN went away and
1175  // the second squareroot((a-r)/two) returned zero correctly. I
1176  // guess that calling the output routine flushed the registers or
1177  // something and restarted the squareroot rountine for this compiler
1178  // or something similar. Actually, due to roundoff, it is possible that a-r
1179  // might be slightly less than zero (i.e. -1e-16) and that would cause
1180  // a possbile NaN return. THe above if test is the right thing to do
1181  // I think and is very cheap.
1182  static inline ComplexT pow(ComplexT x, ComplexT y) { return pow(x,y); }
1183  static inline ComplexT pi() { return ScalarTraits<T>::pi(); }
1184 };
1185 
1186 #endif // HAVE_TEUCHOS_COMPLEX
1187 #endif // DOXYGEN_SHOULD_SKIP_THIS
1188 
1189 } // Teuchos namespace
1190 
1191 #endif // _TEUCHOS_SCALARTRAITS_HPP_
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
The Teuchos namespace contains all of the classes, structs and enums used by Teuchos, as well as a number of utility routines.
Generated by   doxygen 1.8.13