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