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