[Trilinos-Users] Tpetra and boost::interval

[Sensei]

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