# [Trilinos-Users] Tpetra and boost::interval

Sensei senseiwa at gmail.com
Wed Jul 22 09:23:19 EDT 2015

```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()
{
}

//! 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));
}
};
}

```