[Trilinos-Users] [EXTERNAL] Re: Tpetra and boost::interval
Phipps, Eric T
etphipp at sandia.gov
Wed Jul 22 13:58:06 EDT 2015
Well, you need to write down mathematically what it means to solve the
problem with an interval scalar type. For example, what does it mean for
an interval linear system solved by an iterative solver method to have
converged. That will tell you what the magnitude type should be (since
tolerances are magnitudes). Related to this is what comparison operators
mean, e.g., a < b when a and or b are intervals.
You can see examples of implementing scalar traits in the Stokhos package,
which defines several scalar types for use within uncertainty
quantification, e.g.,
Trilinos/packages/stokhos/src/sacado/kokkos/vector/Sacado_MP_VectorTraits.h
pp
Note, for Tpetra to work out-of-the-box with any scalar type, it needs to
be a POD type. It can be made to work with non-POD scalar types (there
are examples of this in Stokhos), but it requires further specialization
that requires expert knowledge of the internals of Tpetra and Kokkos.
-Eric
On 7/22/15, 7:23 AM, "Sensei" <senseiwa at gmail.com> wrote:
>Ok, I've set up a sort of specialization for Teuchos::ScalarTraits, but
>I am not really skilled in interval logic, nor I can see if some of the
>following methods should be double or interval.
>
>As I understand, all methods like min/max/base & Co., must return a
>double (in my case), some others a scalar type.
>
>However, magnitudeType isn't clear to me. The result of a magniture
>means a double, but what should I put there? The absolute value of an
>interval, i.e., another interval?
>
>
>
>Thanks for any hints you can give me!
>
>
>
>namespace Teuchos
>{
> template <>
> struct ScalarTraits<scalar>
> {
> //! Mandatory typedef for result of magnitude
> typedef double magnitudeType;
>
> //! Typedef for half precision
> typedef double halfPrecision;
>
> //! Typedef for double precision
> typedef double doublePrecision;
>
> //! Determines if scalar type is std::complex
> static const bool isComplex = false;
>
> //! Determines if scalar type is an ordinal type
> static const bool isOrdinal = false;
>
> //! Determines if scalar type supports relational operators
>such as <, >, <=, >=.
> static const bool isComparable = true;
>
> /** \brief Determines if scalar type have machine-specific
>parameters
> * (i.e. eps(), sfmin(), base(), prec(), t(), rnd(), emin(),
>rmin(), emax(),
> * rmax() are supported).
> */
> static const bool hasMachineParameters = false;
>
> //! Returns relative machine precision.
> static inline magnitudeType eps()
> {
> return std::numeric_limits<double>::epsilon();
> }
>
> //! Returns safe minimum (sfmin), such that 1/sfmin does not
>overflow.
> static inline magnitudeType sfmin()
> {
> return std::numeric_limits<double>::min();
> }
>
> //! Returns the base of the machine.
> static inline magnitudeType base()
> {
> return std::numeric_limits<double>::radix;
> }
>
> //! Returns \c eps*base.
> static inline magnitudeType prec()
> {
> return eps() * base();
> }
>
> //! Returns the number of (base) digits in the mantissa.
> static inline magnitudeType t()
> {
> return std::numeric_limits<double>::digits;
> }
>
> //! Returns 1.0 when rounding occurs in addition, 0.0 otherwise
> static inline magnitudeType rnd()
> {
> return (std::numeric_limits<double>::round_style ==
>std::round_to_nearest ? double(1.0) : double(0.0));
> }
>
> //! Returns the minimum exponent before (gradual) underflow.
> static inline magnitudeType emin()
> {
> return std::numeric_limits<double>::min_exponent;
> }
>
> //! Returns the underflow threshold - \c base^(emin-1)
> static inline magnitudeType rmin()
> {
> return std::numeric_limits<double>::min();
> }
>
> //! Returns the largest exponent before overflow.
> static inline magnitudeType emax()
> {
> return std::numeric_limits<double>::max_exponent;
> }
>
> //! Overflow theshold - \c (base^emax)*(1-eps)
> static inline magnitudeType rmax()
> {
> return std::numeric_limits<double>::max();
> }
>
> //! Returns the magnitudeType of the scalar type \c a.
> static inline magnitudeType magnitude(scalar a)
> {
> // ????
> return std::fabs(0.0);
> }
>
> //! Returns representation of zero for this scalar type.
> static inline scalar zero()
> {
> return scalar(0.0, 0.0);
> }
>
> //! Returns representation of one for this scalar type.
> static inline scalar one()
> {
> return scalar(1.0, 1.0);
> }
>
> //! Returns the real part of the scalar type \c a.
> static inline scalar real(scalar a)
> {
> return a;
> }
>
> //! Returns the imaginary part of the scalar type \c a.
> static inline scalar imag(scalar a)
> {
> return zero();
> }
>
> //! Returns the conjugate of the scalar type \c a.
> static inline scalar conjugate(scalar a)
> {
> return a;
> }
>
> //! Returns a number that represents NaN.
> static inline scalar nan()
> {
> // SHOULD THIS BE THE ARCHETYPE FOR REAL/IMAG/CONJUGATE ?
> return scalar(dbl_nan, dbl_nan);
> }
>
> //! Returns <tt>true</tt> if <tt>x</tt> is NaN or Inf.
> static inline bool isnaninf(const scalar& x)
> {
> return generic_real_isnaninf(x.lower()) ||
>generic_real_isnaninf(x.upper());
> }
>
> //! Seed the random number generator returned by
><tt>random()</tt>.
> static inline void seedrandom(unsigned int s)
> {
> // Copied from teuchos
> std::srand(s);
>#if __APPLE__
> random();
>#endif
> }
>
> //! Returns a random number (between -one() and +one()) of this
>scalar type.
> static inline scalar random()
> {
> // Copied from teuchos, now we just have two random numbers
> double d1, d2;
>
> d1 = -1.0 + 2.0 * ((double) std::rand() / RAND_MAX);
> d2 = -1.0 + 2.0 * ((double) std::rand() / RAND_MAX);
>
> return scalar(d1, d2);
> }
>
> //! Returns the name of this scalar type.
> static inline std::string name()
> {
> return "scalar aka 'boost::numeric::interval<double>'";
> }
>
> //! Returns a number of magnitudeType that is the square root
>of this scalar type \c x.
> static inline scalar squareroot(scalar x)
> {
> // Defined by boost
> return boost::numeric::sqrt(x);;
> }
>
> //! Returns the result of raising one scalar \c x to the power
>\c y.
> static inline scalar pow(scalar x, scalar y)
> {
> // ????
> return zero(); // ERROR: boost::numeric::exp(y *
>boost::numeric::log(x));
> }
> };
>}
>
>
>
>
>
>
>
>
>
>_______________________________________________
>Trilinos-Users mailing list
>Trilinos-Users at trilinos.org
>https://trilinos.org/mailman/listinfo/trilinos-users
More information about the Trilinos-Users
mailing list