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