Kokkos_ArithTraits.hpp

Go to the documentation of this file.

/*
//@HEADER
// ************************************************************************
//
//   Kokkos: Manycore Performance-Portable Multidimensional Arrays
//              Copyright (2012) Sandia Corporation
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact  H. Carter Edwards (hcedwar@sandia.gov)
//
// ************************************************************************
//@HEADER
*/

#ifndef KOKKOS_ARITHTRAITS_HPP
#define KOKKOS_ARITHTRAITS_HPP


#include <Kokkos_complex.hpp>

#include <cfloat>
#include <climits>
#include <cmath>
#include <cstdlib> // strtof, strtod, strtold
#include <complex> // std::complex
#include <limits> // std::numeric_limits
#ifdef __CUDACC__
#include <math_constants.h>
#endif
//
// mfh 24 Dec 2013: Temporary measure for testing; will go away.
//
#ifndef KOKKOS_FORCEINLINE_FUNCTION
#  ifdef __CUDA_ARCH__
#    define KOKKOS_FORCEINLINE_FUNCTION inline __host__ __device__
#  else
#    define KOKKOS_FORCEINLINE_FUNCTION
#  endif // __CUDA_ARCH__
#endif // KOKKOS_FORCEINLINE_FUNCTION

namespace { // anonymous

template<class IntType>
KOKKOS_FORCEINLINE_FUNCTION IntType
intPowImpl (const IntType x, const IntType y)
{
  // Recursion (unrolled into while loop): pow(x, 2y) = (x^y)^2
  IntType prod = x;
  IntType y_cur = 1;
  // If y == 1, then prod stays x.
  while (y_cur < y) {
    prod = prod * prod;
    y_cur = y_cur << 1;
  }
  // abs(y - y_cur) < floor(log2(y)), so it won't hurt asymptotic run
  // time to finish the remainder in a linear iteration.
  if (y > y_cur) {
    const IntType left = y - y_cur;
    for (IntType k = 0; k < left; ++k) {
      prod = prod * x;
    }
  }
  else if (y < y_cur) {
    // There's probably a better way to do this in order to avoid the
    // (expensive) integer division, but I'm not motivated to think of
    // it at the moment.
    const IntType left = y_cur - y;
    for (IntType k = 0; k < left; ++k) {
      prod = prod / x;
    }
  }
  return prod;

  // y = 8:
  //
  // x,1   -> x^2,2
  // x^2,2 -> x^4,4
  // x^4,4 -> x^8,8
  //
  // y = 9:
  //
  // x,1   -> x^2,2
  // x^2,2 -> x^4,4
  // x^4,4 -> x^8,8
  //
  // y - y_cur is what's left over.  Just do it one at a time.
  //
  // y = 3:
  // x,1   -> x^2,2
  // x^2,2 -> x^4,4
}



template<class IntType>
KOKKOS_FORCEINLINE_FUNCTION IntType
intPowSigned (const IntType x, const IntType y)
{
  // It's not entirely clear what to return if x and y are both zero.
  // In the case of floating-point numbers, 0^0 is NaN.  Here, though,
  // I think it's safe to return 0.
  if (x == 0) {
    return 0;
  } else if (y == 0) {
    return 1;
  } else if (y < 0) {
    if (x == 1) {
      return 1;
    }
    else if (x == -1) {
      return (y % 2 == 0) ? 1 : -1;
    }
    else {
      return 0; // round the fraction to zero
    }
  } else {
    return intPowImpl<IntType> (x, y);
  }
}

template<class IntType>
KOKKOS_FORCEINLINE_FUNCTION IntType
intPowUnsigned (const IntType x, const IntType y)
{
  // It's not entirely clear what to return if x and y are both zero.
  // In the case of floating-point numbers, 0^0 is NaN.  Here, though,
  // I think it's safe to return 0.
  if (x == 0) {
    return 0;
  } else if (y == 0) {
    return 1;
  } else {
    return intPowImpl<IntType> (x, y);
  }
}

// It might make sense to use special sqrt() approximations for
// integer arguments, like those presented on the following web site:
//
// http://www.azillionmonkeys.com/qed/sqroot.html#implementations
//
// Note that some of the implementations on the above page break ANSI
// C(++) aliasing rules (by assigning to the results of
// reinterpret_cast-ing between int and float).  It's also just a
// performance optimization and not required for a reasonable
// implementation.

} // namespace (anonymous)

namespace Kokkos {
namespace Details {

template<class T>
class ArithTraits {
public:
  typedef T val_type;
  typedef T mag_type;

  static const bool is_specialized = false;
  static const bool is_signed = false;
  static const bool is_integer = false;
  static const bool is_exact = false;
  static const bool is_complex = false;

  static KOKKOS_FORCEINLINE_FUNCTION bool isInf (const T& x);

  static KOKKOS_FORCEINLINE_FUNCTION bool isNan (const T& x);

  static KOKKOS_FORCEINLINE_FUNCTION mag_type abs (const T& x);

  static KOKKOS_FORCEINLINE_FUNCTION T zero ();

  static KOKKOS_FORCEINLINE_FUNCTION T one ();

  static KOKKOS_FORCEINLINE_FUNCTION T min ();

  static KOKKOS_FORCEINLINE_FUNCTION T max ();

  static KOKKOS_FORCEINLINE_FUNCTION mag_type real (const T& x);

  static KOKKOS_FORCEINLINE_FUNCTION mag_type imag (const T&);

  static KOKKOS_FORCEINLINE_FUNCTION T conj (const T&);

  static KOKKOS_FORCEINLINE_FUNCTION T pow (const T& x, const T& y);

  static KOKKOS_FORCEINLINE_FUNCTION T sqrt (const T& x);

  static KOKKOS_FORCEINLINE_FUNCTION T log (const T& x);

  static KOKKOS_FORCEINLINE_FUNCTION T log10 (const T& x);

  static KOKKOS_FORCEINLINE_FUNCTION T nan ();

  static KOKKOS_FORCEINLINE_FUNCTION mag_type epsilon ();


  typedef T magnitudeType;

  typedef T halfPrecision;

  typedef T doublePrecision;

  static const bool isComplex = false;
  static const bool isOrdinal = false;
  static const bool isComparable = false;

  static const bool hasMachineParameters = false;

  static KOKKOS_FORCEINLINE_FUNCTION mag_type eps ();

  static KOKKOS_FORCEINLINE_FUNCTION mag_type sfmin ();

  static KOKKOS_FORCEINLINE_FUNCTION int base ();

  static KOKKOS_FORCEINLINE_FUNCTION mag_type prec ();

  static KOKKOS_FORCEINLINE_FUNCTION int t ();

  static KOKKOS_FORCEINLINE_FUNCTION mag_type rnd ();

  static KOKKOS_FORCEINLINE_FUNCTION int emin ();

  static KOKKOS_FORCEINLINE_FUNCTION mag_type rmin ();

  static KOKKOS_FORCEINLINE_FUNCTION int emax ();

  static KOKKOS_FORCEINLINE_FUNCTION mag_type rmax ();

  static KOKKOS_FORCEINLINE_FUNCTION magnitudeType magnitude (const T& x);

  static KOKKOS_FORCEINLINE_FUNCTION T conjugate (const T& x);

  static KOKKOS_FORCEINLINE_FUNCTION bool isnaninf (const T& x);

  static std::string name ();

  static KOKKOS_FORCEINLINE_FUNCTION T squareroot (const T& x);
};


template<>
class ArithTraits<float> {
public:
  typedef float val_type;
  typedef val_type mag_type;

  static const bool is_specialized = true;
  static const bool is_signed = true;
  static const bool is_integer = false;
  static const bool is_exact = false;
  static const bool is_complex = false;

  static KOKKOS_FORCEINLINE_FUNCTION bool isInf (const float x) {
#ifdef __APPLE__
    return std::isinf (x);
#else
    return isinf (x);
#endif // __APPLE__
  }
  static KOKKOS_FORCEINLINE_FUNCTION bool isNan (const float x) {
#ifdef __APPLE__
    return std::isnan (x);
#else
    return isnan (x);
#endif // __APPLE__
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type abs (const float x) {
    return ::fabs (x);
  }
  static KOKKOS_FORCEINLINE_FUNCTION float zero () {
    return 0.0;
  }
  static KOKKOS_FORCEINLINE_FUNCTION float one () {
    return 1.0;
  }
  static KOKKOS_FORCEINLINE_FUNCTION float min () {
    return FLT_MIN;
  }
  static KOKKOS_FORCEINLINE_FUNCTION float max () {
    return FLT_MAX;
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type real (const float x) {
    return x;
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type imag (const float) {
    return 0.0;
  }
  static KOKKOS_FORCEINLINE_FUNCTION float conj (const float x) {
    return x;
  }
  static KOKKOS_FORCEINLINE_FUNCTION float pow (const float x, const float y) {
    return ::pow (x, y);
  }
  static KOKKOS_FORCEINLINE_FUNCTION float sqrt (const float x) {
    return ::sqrt (x);
  }
  static KOKKOS_FORCEINLINE_FUNCTION float log (const float x) {
    return ::log (x);
  }
  static KOKKOS_FORCEINLINE_FUNCTION float log10 (const float x) {
    return ::log10 (x);
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type epsilon () {
    return FLT_EPSILON;
  }

  // Backwards compatibility with Teuchos::ScalarTraits.
  typedef mag_type magnitudeType;
  // C++ doesn't have a standard "half-float" type.
  typedef float halfPrecision;
  typedef double doublePrecision;

  static const bool isComplex = false;
  static const bool isOrdinal = false;
  static const bool isComparable = true;
  static const bool hasMachineParameters = true;
  static KOKKOS_FORCEINLINE_FUNCTION bool isnaninf (const float x) {
    return isNan (x) || isInf (x);
  }
  static KOKKOS_FORCEINLINE_FUNCTION magnitudeType magnitude (const float x) {
    return abs (x);
  }
  static KOKKOS_FORCEINLINE_FUNCTION float conjugate (const float x) {
    return conj (x);
  }
  static std::string name () {
    return "float";
  }
  static KOKKOS_FORCEINLINE_FUNCTION float squareroot (const float x) {
    return sqrt (x);
  }
  static KOKKOS_FORCEINLINE_FUNCTION float nan () {
#ifdef __CUDA_ARCH__
    return CUDART_NAN_F;
    //return nan (); //this returns 0???
#else
    // http://pubs.opengroup.org/onlinepubs/009696899/functions/nan.html
    return strtof ("NAN()", (char**) NULL);
#endif // __CUDA_ARCH__
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type eps () {
    return epsilon ();
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type sfmin () {
    return FLT_MIN; // ???
  }
  static KOKKOS_FORCEINLINE_FUNCTION int base () {
    return FLT_RADIX;
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type prec () {
    return eps () * static_cast<mag_type> (base ());
  }
  static KOKKOS_FORCEINLINE_FUNCTION int t () {
    return FLT_MANT_DIG;
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type rnd () {
    return 1.0;
  }
  static KOKKOS_FORCEINLINE_FUNCTION int emin () {
    return FLT_MIN_EXP;
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type rmin () {
    return FLT_MIN; // ??? // should be base^(emin-1)
  }
  static KOKKOS_FORCEINLINE_FUNCTION int emax () {
    return FLT_MAX_EXP;
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type rmax () {
    return FLT_MAX; // ??? // should be (base^emax)*(1-eps)
  }
};


// The C++ Standard Library (with C++03 at least) only allows
// std::complex<T> for T = float, double, or long double.
template<class RealFloatType>
class ArithTraits<std::complex<RealFloatType> > {
public:
  typedef ::Kokkos::complex<RealFloatType> val_type;
  typedef RealFloatType mag_type;

  static const bool is_specialized = true;
  static const bool is_signed = true;
  static const bool is_integer = false;
  static const bool is_exact = false;
  static const bool is_complex = true;

  static bool isInf (const std::complex<RealFloatType>& x) {
#if defined(__CUDACC__) || defined(_CRAYC)
    return isinf (real (x)) || isinf (imag (x));
#else
    return std::isinf (real (x)) || std::isinf (imag (x));
#endif // __CUDACC__
  }
  static bool isNan (const std::complex<RealFloatType>& x) {
#if defined(__CUDACC__) || defined(_CRAYC)
    return isnan (real (x)) || isnan (imag (x));
#else
    return std::isnan (real (x)) || std::isnan (imag (x));
#endif // __CUDACC__
  }
  static mag_type abs (const std::complex<RealFloatType>& x) {
    return std::abs (x);
  }
  static std::complex<RealFloatType> zero () {
    return std::complex<RealFloatType> (ArithTraits<mag_type>::zero (), ArithTraits<mag_type>::zero ());
  }
  static std::complex<RealFloatType> one () {
    return std::complex<RealFloatType> (ArithTraits<mag_type>::one (), ArithTraits<mag_type>::zero ());
  }
  static std::complex<RealFloatType> min () {
    return std::complex<RealFloatType> (ArithTraits<mag_type>::min (), ArithTraits<mag_type>::zero ());
  }
  static std::complex<RealFloatType> max () {
    return std::complex<RealFloatType> (ArithTraits<mag_type>::max (), ArithTraits<mag_type>::zero ());
  }
  static mag_type real (const std::complex<RealFloatType>& x) {
    return std::real (x);
  }
  static mag_type imag (const std::complex<RealFloatType>& x) {
    return std::imag (x);
  }
  static std::complex<RealFloatType> conj (const std::complex<RealFloatType>& x) {
    return std::conj (x);
  }
  static std::complex<RealFloatType>
  pow (const std::complex<RealFloatType>& x, const std::complex<RealFloatType>& y) {
    // Fix for some weird gcc 4.2.1 inaccuracy.
    if (y == one ()) {
      return x;
    } else if (y == one () + one ()) {
      return x * x;
    } else {
      return std::pow (x, y);
    }
  }
  static std::complex<RealFloatType> sqrt (const std::complex<RealFloatType>& x) {
    return std::sqrt (x);
  }
  static std::complex<RealFloatType> log (const std::complex<RealFloatType>& x) {
    return std::log (x);
  }
  static std::complex<RealFloatType> log10 (const std::complex<RealFloatType>& x) {
    return std::log10 (x);
  }
  static std::complex<RealFloatType> nan () {
    const mag_type mag_nan = ArithTraits<mag_type>::nan ();
    return std::complex<RealFloatType> (mag_nan, mag_nan);
  }
  static mag_type epsilon () {
    return ArithTraits<mag_type>::epsilon ();
  }

  // Backwards compatibility with Teuchos::ScalarTraits.
  typedef mag_type magnitudeType;
  typedef std::complex<typename ArithTraits<mag_type>::halfPrecision> halfPrecision;
  typedef std::complex<typename ArithTraits<mag_type>::doublePrecision> doublePrecision;

  static const bool isComplex = true;
  static const bool isOrdinal = false;
  static const bool isComparable = false;
  static const bool hasMachineParameters = true;
  static bool isnaninf (const std::complex<RealFloatType>& x) {
    return isNan (x) || isInf (x);
  }
  static mag_type magnitude (const std::complex<RealFloatType>& x) {
    return abs (x);
  }
  static std::complex<RealFloatType> conjugate (const std::complex<RealFloatType>& x) {
    return conj (x);
  }
  static std::string name () {
    return std::string ("std::complex<") + ArithTraits<mag_type>::name () + ">";
  }
  static std::complex<RealFloatType> squareroot (const std::complex<RealFloatType>& x) {
    return sqrt (x);
  }
  static mag_type eps () {
    return epsilon ();
  }
  static mag_type sfmin () {
    return ArithTraits<mag_type>::sfmin ();
  }
  static int base () {
    return ArithTraits<mag_type>::base ();
  }
  static mag_type prec () {
    return ArithTraits<mag_type>::prec ();
  }
  static int t () {
    return ArithTraits<mag_type>::t ();
  }
  static mag_type rnd () {
    return ArithTraits<mag_type>::one ();
  }
  static int emin () {
    return ArithTraits<mag_type>::emin ();
  }
  static mag_type rmin () {
    return ArithTraits<mag_type>::rmin ();
  }
  static int emax () {
    return ArithTraits<mag_type>::emax ();
  }
  static mag_type rmax () {
    return ArithTraits<mag_type>::rmax ();
  }
};


template<>
class ArithTraits<double> {
public:
  typedef double val_type;
  typedef val_type mag_type;

  static const bool is_specialized = true;
  static const bool is_signed = true;
  static const bool is_integer = false;
  static const bool is_exact = false;
  static const bool is_complex = false;

  static KOKKOS_FORCEINLINE_FUNCTION bool isInf (const val_type x) {
#ifdef __APPLE__
    return std::isinf (x);
#else
    return isinf (x);
#endif // __APPLE__
  }
  static KOKKOS_FORCEINLINE_FUNCTION bool isNan (const val_type x) {
#ifdef __APPLE__
    return std::isnan (x);
#else
    return isnan (x);
#endif // __APPLE__
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type abs (const val_type x) {
    return ::fabs (x);
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type zero () {
    return 0.0;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type one () {
    return 1.0;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type min () {
    return DBL_MIN;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type max () {
    return DBL_MAX;
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type real (const val_type x) {
    return x;
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type imag (const val_type) {
    return 0.0;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type conj (const val_type x) {
    return x;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type pow (const val_type x, const val_type y) {
    return ::pow (x, y);
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type sqrt (const val_type x) {
    return ::sqrt (x);
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type log (const val_type x) {
    return ::log (x);
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type log10 (const val_type x) {
    return ::log10 (x);
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type nan () {
#ifdef __CUDA_ARCH__
    return CUDART_NAN;
    //return nan (); // this returns 0 ???
#else
    // http://pubs.opengroup.org/onlinepubs/009696899/functions/nan.html
    return strtod ("NAN", (char**) NULL);
#endif // __CUDA_ARCH__
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type epsilon () {
    return DBL_EPSILON;
  }

  // Backwards compatibility with Teuchos::ScalarTraits.
  typedef mag_type magnitudeType;
  typedef float halfPrecision;
#ifdef __CUDA_ARCH__
  typedef double doublePrecision; // CUDA doesn't support long double, unfortunately
#else
  typedef long double doublePrecision;
#endif // __CUDA_ARCH__
  static const bool isComplex = false;
  static const bool isOrdinal = false;
  static const bool isComparable = true;
  static const bool hasMachineParameters = true;
  static bool isnaninf (const val_type& x) {
    return isNan (x) || isInf (x);
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type magnitude (const val_type x) {
    return abs (x);
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type conjugate (const val_type x) {
    return conj (x);
  }
  static std::string name () {
    return "double";
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type squareroot (const val_type x) {
    return sqrt (x);
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type eps () {
    return epsilon ();
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type sfmin () {
    return DBL_MIN; // ???
  }
  static KOKKOS_FORCEINLINE_FUNCTION int base () {
    return FLT_RADIX; // same for float as for double
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type prec () {
    return eps () * static_cast<mag_type> (base ());
  }
  static KOKKOS_FORCEINLINE_FUNCTION int t () {
    return DBL_MANT_DIG;
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type rnd () {
    return 1.0;
  }
  static KOKKOS_FORCEINLINE_FUNCTION int emin () {
    return DBL_MIN_EXP;
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type rmin () {
    return DBL_MIN; // ??? // should be base^(emin-1)
  }
  static KOKKOS_FORCEINLINE_FUNCTION int emax () {
    return DBL_MAX_EXP;
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type rmax () {
    return DBL_MAX; // ??? // should be (base^emax)*(1-eps)
  }
};


// CUDA does not support long double in device functions, so none of
// the class methods in this specialization are marked as device
// functions.
template<>
class ArithTraits<long double> {
public:
  typedef long double val_type;
  typedef long double mag_type;

  static const bool is_specialized = true;
  static const bool is_signed = true;
  static const bool is_integer = false;
  static const bool is_exact = false;
  static const bool is_complex = false;

  static bool isInf (const val_type& x) {
#if defined(__CUDACC__) || defined(_CRAYC)
    return isinf (x);
#else
    return std::isinf (x);
#endif // __CUDACC__
  }
  static bool isNan (const val_type& x) {
#if defined(__CUDACC__) || defined(_CRAYC)
    return isnan (x);
#else
    return std::isnan (x);
#endif // __CUDACC__
  }
  static mag_type abs (const val_type& x) {
    return ::fabs (x);
  }
  static val_type zero () {
    return 0.0;
  }
  static val_type one () {
    return 1.0;
  }
  static val_type min () {
    return LDBL_MIN;
  }
  static val_type max () {
    return LDBL_MAX;
  }
  static mag_type real (const val_type& x) {
    return x;
  }
  static mag_type imag (const val_type&) {
    return zero ();
  }
  static val_type conj (const val_type& x) {
    return x;
  }
  static val_type pow (const val_type& x, const val_type& y) {
    return ::pow (x, y);
  }
  static val_type sqrt (const val_type& x) {
    return ::sqrt (x);
  }
  static val_type log (const val_type& x) {
    return ::log (x);
  }
  static val_type log10 (const val_type& x) {
    return ::log10 (x);
  }
  static val_type nan () {
    return strtold ("NAN()", (char**) NULL);
  }
  static mag_type epsilon () {
    return LDBL_EPSILON;
  }

  // Backwards compatibility with Teuchos::ScalarTraits.
  typedef mag_type magnitudeType;
  typedef double halfPrecision;
  // It might be appropriate to use QD's qd_real here.
  // For now, long double is the most you get.
  typedef val_type doublePrecision;

  static const bool isComplex = false;
  static const bool isOrdinal = false;
  static const bool isComparable = true;
  static const bool hasMachineParameters = true;
  static bool isnaninf (const val_type& x) {
    return isNan (x) || isInf (x);
  }
  static mag_type magnitude (const val_type& x) {
    return abs (x);
  }
  static val_type conjugate (const val_type& x) {
    return conj (x);
  }
  static std::string name () {
    return "long double";
  }
  static val_type squareroot (const val_type& x) {
    return sqrt (x);
  }
  static mag_type eps () {
    return epsilon ();
  }
  static mag_type sfmin () {
    return LDBL_MIN; // ???
  }
  static int base () {
    return FLT_RADIX; // same for float as for double or long double
  }
  static mag_type prec () {
    return eps () * static_cast<mag_type> (base ());
  }
  static int t () {
    return LDBL_MANT_DIG;
  }
  static mag_type rnd () {
    return one ();
  }
  static int emin () {
    return LDBL_MIN_EXP;
  }
  static mag_type rmin () {
    return LDBL_MIN;
  }
  static int emax () {
    return LDBL_MAX_EXP;
  }
  static mag_type rmax () {
    return LDBL_MAX;
  }
};


template<>
class ArithTraits< ::Kokkos::complex<float> > {
public:
  typedef ::Kokkos::complex<float> val_type;
  typedef float mag_type;

  static const bool is_specialized = true;
  static const bool is_signed = true;
  static const bool is_integer = false;
  static const bool is_exact = false;
  static const bool is_complex = true;

  static KOKKOS_FORCEINLINE_FUNCTION bool isInf (const val_type x) {
    return ArithTraits<mag_type>::isInf (x.real ()) ||
      ArithTraits<mag_type>::isInf (x.imag ());
  }
  static KOKKOS_FORCEINLINE_FUNCTION bool isNan (const val_type x) {
    return ArithTraits<mag_type>::isNan (x.real ()) ||
      ArithTraits<mag_type>::isNan (x.imag ());
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type abs (const val_type x) {
    return sqrt (::Kokkos::real (x) * ::Kokkos::real (x) +
                 ::Kokkos::imag (x) * ::Kokkos::imag (x));
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type zero () {
    return val_type (ArithTraits<mag_type>::zero (), ArithTraits<mag_type>::zero ());
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type one () {
    return val_type (ArithTraits<mag_type>::one (), ArithTraits<mag_type>::zero ());
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type min () {
    return val_type (ArithTraits<mag_type>::min (), ArithTraits<mag_type>::min ()); // ???
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type max () {
    return val_type (ArithTraits<mag_type>::max (), ArithTraits<mag_type>::max ()); // ???
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type real (const val_type x) {
    return x.real ();
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type imag (const val_type x) {
    return x.imag ();
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type conj (const val_type x) {
    return ::Kokkos::conj (x);
  }
  // static KOKKOS_FORCEINLINE_FUNCTION val_type pow (const val_type x, const val_type y) {
  //   return ::pow (x, y);
  // }
  // static KOKKOS_FORCEINLINE_FUNCTION val_type sqrt (const val_type x) {
  //   return ::sqrt (x);
  // }
  // static KOKKOS_FORCEINLINE_FUNCTION val_type log (const val_type x) {
  //   return ::log (x);
  // }
  // static KOKKOS_FORCEINLINE_FUNCTION val_type log10 (const val_type x) {
  //   return ::log10 (x);
  // }
  static KOKKOS_FORCEINLINE_FUNCTION val_type nan () {
    // ???
    return val_type (ArithTraits<mag_type>::nan (), ArithTraits<mag_type>::nan ());
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type epsilon () {
    return ArithTraits<mag_type>::epsilon (); // ???
  }

  // Backwards compatibility with Teuchos::ScalarTraits.
  typedef mag_type magnitudeType;
  typedef ::Kokkos::complex<ArithTraits<mag_type>::halfPrecision> halfPrecision;
  typedef ::Kokkos::complex<ArithTraits<mag_type>::doublePrecision> doublePrecision;

  static const bool isComplex = true;
  static const bool isOrdinal = false;
  static const bool isComparable = false;
  static const bool hasMachineParameters = ArithTraits<mag_type>::hasMachineParameters;
  static bool isnaninf (const val_type& x) {
    return isNan (x) || isInf (x);
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type magnitude (const val_type x) {
    return abs (x);
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type conjugate (const val_type x) {
    return conj (x);
  }
  static std::string name () {
    return "Kokkos::complex<float>";
  }
  // static KOKKOS_FORCEINLINE_FUNCTION val_type squareroot (const val_type x) {
  //   return sqrt (x);
  // }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type eps () {
    return epsilon ();
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type sfmin () {
    return ArithTraits<mag_type>::sfmin (); // ???
  }
  static KOKKOS_FORCEINLINE_FUNCTION int base () {
    return ArithTraits<mag_type>::base ();
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type prec () {
    return ArithTraits<mag_type>::prec (); // ???
  }
  static KOKKOS_FORCEINLINE_FUNCTION int t () {
    return ArithTraits<mag_type>::t ();
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type rnd () {
    return ArithTraits<mag_type>::rnd ();
  }
  static KOKKOS_FORCEINLINE_FUNCTION int emin () {
    return ArithTraits<mag_type>::emin ();
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type rmin () {
    return ArithTraits<mag_type>::rmin ();
  }
  static KOKKOS_FORCEINLINE_FUNCTION int emax () {
    return ArithTraits<mag_type>::emax ();
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type rmax () {
    return ArithTraits<mag_type>::rmax ();
  }
};


template<>
class ArithTraits< ::Kokkos::complex<double> > {
public:
  typedef ::Kokkos::complex<double> val_type;
  typedef double mag_type;

  static const bool is_specialized = true;
  static const bool is_signed = true;
  static const bool is_integer = false;
  static const bool is_exact = false;
  static const bool is_complex = true;

  static KOKKOS_FORCEINLINE_FUNCTION bool isInf (const val_type x) {
    return ArithTraits<mag_type>::isInf (x.real ()) ||
      ArithTraits<mag_type>::isInf (x.imag ());
  }
  static KOKKOS_FORCEINLINE_FUNCTION bool isNan (const val_type x) {
    return ArithTraits<mag_type>::isNan (x.real ()) ||
      ArithTraits<mag_type>::isNan (x.imag ());
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type abs (const val_type x) {
    return ::Kokkos::abs (x);
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type zero () {
    return val_type (ArithTraits<mag_type>::zero (), ArithTraits<mag_type>::zero ());
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type one () {
    return val_type (ArithTraits<mag_type>::one (), ArithTraits<mag_type>::zero ());
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type min () {
    return val_type (ArithTraits<mag_type>::min (), ArithTraits<mag_type>::min ()); // ???
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type max () {
    return val_type (ArithTraits<mag_type>::max (), ArithTraits<mag_type>::max ()); // ???
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type real (const val_type x) {
    return x.real ();
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type imag (const val_type x) {
    return x.imag ();
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type conj (const val_type x) {
    return ::Kokkos::conj (x);
  }
  // static KOKKOS_FORCEINLINE_FUNCTION val_type pow (const val_type x, const val_type y) {
  //   return ::pow (x, y);
  // }
  // static KOKKOS_FORCEINLINE_FUNCTION val_type sqrt (const val_type x) {
  //   return ::sqrt (x);
  // }
  // static KOKKOS_FORCEINLINE_FUNCTION val_type log (const val_type x) {
  //   return ::log (x);
  // }
  // static KOKKOS_FORCEINLINE_FUNCTION val_type log10 (const val_type x) {
  //   return ::log10 (x);
  // }
  static KOKKOS_FORCEINLINE_FUNCTION val_type nan () {
    // ???
    return val_type (ArithTraits<mag_type>::nan (), ArithTraits<mag_type>::nan ());
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type epsilon () {
    return ArithTraits<mag_type>::epsilon (); // ???
  }

  // Backwards compatibility with Teuchos::ScalarTraits.
  typedef mag_type magnitudeType;
  typedef ::Kokkos::complex<ArithTraits<mag_type>::halfPrecision> halfPrecision;
  typedef ::Kokkos::complex<ArithTraits<mag_type>::doublePrecision> doublePrecision;

  static const bool isComplex = true;
  static const bool isOrdinal = false;
  static const bool isComparable = false;
  static const bool hasMachineParameters = ArithTraits<mag_type>::hasMachineParameters;
  static bool isnaninf (const val_type& x) {
    return isNan (x) || isInf (x);
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type magnitude (const val_type x) {
    return abs (x);
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type conjugate (const val_type x) {
    return conj (x);
  }
  static std::string name () {
    return "Kokkos::complex<double>";
  }
  // static KOKKOS_FORCEINLINE_FUNCTION val_type squareroot (const val_type x) {
  //   return sqrt (x);
  // }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type eps () {
    return epsilon ();
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type sfmin () {
    return ArithTraits<mag_type>::sfmin (); // ???
  }
  static KOKKOS_FORCEINLINE_FUNCTION int base () {
    return ArithTraits<mag_type>::base ();
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type prec () {
    return ArithTraits<mag_type>::prec (); // ???
  }
  static KOKKOS_FORCEINLINE_FUNCTION int t () {
    return ArithTraits<mag_type>::t ();
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type rnd () {
    return ArithTraits<mag_type>::rnd ();
  }
  static KOKKOS_FORCEINLINE_FUNCTION int emin () {
    return ArithTraits<mag_type>::emin ();
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type rmin () {
    return ArithTraits<mag_type>::rmin ();
  }
  static KOKKOS_FORCEINLINE_FUNCTION int emax () {
    return ArithTraits<mag_type>::emax ();
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type rmax () {
    return ArithTraits<mag_type>::rmax ();
  }
};


template<>
class ArithTraits<char> {
public:
  typedef char val_type;
  typedef val_type mag_type;

  static const bool is_specialized = true;
  // The C(++) standard does not require that char be signed.  In
  // fact, signed char, unsigned char, and char are distinct types.
  // We can use std::numeric_limits here because it's a const bool,
  // not a class method.
  static const bool is_signed = std::numeric_limits<char>::is_signed;
  static const bool is_integer = true;
  static const bool is_exact = true;
  static const bool is_complex = false;

  static KOKKOS_FORCEINLINE_FUNCTION bool isInf (const val_type x) {
    return false;
  }
  static KOKKOS_FORCEINLINE_FUNCTION bool isNan (const val_type x) {
    return false;
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type abs (const val_type x) {
    // This may trigger a compiler warning if char is unsigned.  On
    // all platforms I have encountered, char is signed, but the C(++)
    // standard does not require this.
    return x >= 0 ? x : -x;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type zero () {
    return 0;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type one () {
    return 1;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type min () {
    return CHAR_MIN;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type max () {
    return CHAR_MAX;
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type real (const val_type x) {
    return x;
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type imag (const val_type x) {
    return 0;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type conj (const val_type x) {
    return x;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type
  pow (const val_type x, const val_type y) {
    if (is_signed) {
      return intPowSigned<val_type> (x, y);
    } else {
      return intPowUnsigned<val_type> (x, y);
    }
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type sqrt (const val_type x) {
    // C++11 defines std::sqrt for integer arguments.  However, we
    // currently can't assume C++11.
    //
    // This cast will result in no loss of accuracy, though it might
    // be more expensive than it should, if we were clever about using
    // bit operations.
    //
    // We take the absolute value first to avoid negative arguments.
    // Negative real arguments to sqrt(float) return (float) NaN, but
    // built-in integer types do not have an equivalent to NaN.
    // Casting NaN to an integer type will thus result in some integer
    // value which appears valid, but is not.  We cannot raise an
    // exception in device functions.  Thus, we prefer to take the
    // absolute value of x first, to avoid issues.  Another
    // possibility would be to test for a NaN output and convert it to
    // some reasonable value (like 0), though this might be more
    // expensive than the absolute value interpreted using the ternary
    // operator.
    return static_cast<val_type> ( ::sqrt (static_cast<float> (abs (x))));
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type log (const val_type x) {
    return static_cast<val_type> ( ::log (static_cast<float> (abs (x))));
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type log10 (const val_type x) {
    return static_cast<val_type> ( ::log10 (static_cast<float> (abs (x))));
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type epsilon () {
    return zero ();
  }

  // Backwards compatibility with Teuchos::ScalarTraits.
  typedef mag_type magnitudeType;
  typedef val_type halfPrecision;
  typedef val_type doublePrecision;

  static const bool isComplex = false;
  static const bool isOrdinal = true;
  static const bool isComparable = true;
  static const bool hasMachineParameters = false;
  static KOKKOS_FORCEINLINE_FUNCTION magnitudeType magnitude (const val_type x) {
    return abs (x);
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type conjugate (const val_type x) {
    return conj (x);
  }
  static KOKKOS_FORCEINLINE_FUNCTION bool isnaninf (const val_type) {
    return false;
  }
  static std::string name () {
    return "char";
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type squareroot (const val_type x) {
    return sqrt (x);
  }
};


template<>
class ArithTraits<signed char> {
public:
  typedef signed char val_type;
  typedef val_type mag_type;

  static const bool is_specialized = true;
  static const bool is_signed = true;
  static const bool is_integer = true;
  static const bool is_exact = true;
  static const bool is_complex = false;

  static KOKKOS_FORCEINLINE_FUNCTION bool isInf (const val_type x) {
    return false;
  }
  static KOKKOS_FORCEINLINE_FUNCTION bool isNan (const val_type x) {
    return false;
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type abs (const val_type x) {
    return x >= 0 ? x : -x;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type zero () {
    return 0;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type one () {
    return 1;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type min () {
    return SCHAR_MIN;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type max () {
    return SCHAR_MAX;
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type real (const val_type x) {
    return x;
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type imag (const val_type x) {
    return 0;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type conj (const val_type x) {
    return x;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type
  pow (const val_type x, const val_type y) {
    return intPowSigned<val_type> (x, y);
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type sqrt (const val_type x) {
    return static_cast<val_type> ( ::sqrt (static_cast<float> (abs (x))));
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type log (const val_type x) {
    return static_cast<val_type> ( ::log (static_cast<float> (abs (x))));
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type log10 (const val_type x) {
    return static_cast<val_type> ( ::log10 (static_cast<float> (abs (x))));
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type epsilon () {
    return zero ();
  }

  // Backwards compatibility with Teuchos::ScalarTraits.
  typedef mag_type magnitudeType;
  typedef val_type halfPrecision;
  typedef val_type doublePrecision;

  static const bool isComplex = false;
  static const bool isOrdinal = true;
  static const bool isComparable = true;
  static const bool hasMachineParameters = false;
  static KOKKOS_FORCEINLINE_FUNCTION magnitudeType magnitude (const val_type x) {
    return abs (x);
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type conjugate (const val_type x) {
    return conj (x);
  }
  static KOKKOS_FORCEINLINE_FUNCTION bool isnaninf (const val_type) {
    return false;
  }
  static std::string name () {
    return "signed char";
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type squareroot (const val_type x) {
    return sqrt (x);
  }
};


template<>
class ArithTraits<unsigned char> {
public:
  typedef unsigned char val_type;
  typedef val_type mag_type;

  static const bool is_specialized = true;
  static const bool is_signed = false;
  static const bool is_integer = true;
  static const bool is_exact = true;
  static const bool is_complex = false;

  static KOKKOS_FORCEINLINE_FUNCTION bool isInf (const val_type x) {
    return false;
  }
  static KOKKOS_FORCEINLINE_FUNCTION bool isNan (const val_type x) {
    return false;
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type abs (const val_type x) {
    return x; // it's unsigned, so it's positive
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type zero () {
    return 0;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type one () {
    return 1;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type min () {
    return 0;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type max () {
    return UCHAR_MAX;
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type real (const val_type x) {
    return x;
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type imag (const val_type x) {
    return 0;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type conj (const val_type x) {
    return x;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type
  pow (const val_type x, const val_type y) {
    return intPowUnsigned<val_type> (x, y);
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type sqrt (const val_type x) {
    // This will result in no loss of accuracy, though it might be
    // more expensive than it should, if we were clever about using
    // bit operations.
    return static_cast<val_type> ( ::sqrt (static_cast<float> (x)));
  }

  static KOKKOS_FORCEINLINE_FUNCTION val_type log (const val_type x) {
    return static_cast<val_type> ( ::log (static_cast<float> (x)));
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type log10 (const val_type x) {
    return static_cast<val_type> ( ::log10 (static_cast<float> (x)));
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type epsilon () {
    return zero ();
  }

  // Backwards compatibility with Teuchos::ScalarTraits.
  typedef mag_type magnitudeType;
  typedef val_type halfPrecision;
  typedef val_type doublePrecision;

  static const bool isComplex = false;
  static const bool isOrdinal = true;
  static const bool isComparable = true;
  static const bool hasMachineParameters = false;
  static KOKKOS_FORCEINLINE_FUNCTION magnitudeType magnitude (const val_type x) {
    return abs (x);
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type conjugate (const val_type x) {
    return conj (x);
  }
  static KOKKOS_FORCEINLINE_FUNCTION bool isnaninf (const val_type) {
    return false;
  }
  static std::string name () {
    return "unsigned char";
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type squareroot (const val_type x) {
    return sqrt (x);
  }
};


template<>
class ArithTraits<short> {
public:
  typedef short val_type;
  typedef val_type mag_type;

  static const bool is_specialized = true;
  static const bool is_signed = true;
  static const bool is_integer = true;
  static const bool is_exact = true;
  static const bool is_complex = false;

  static KOKKOS_FORCEINLINE_FUNCTION bool isInf (const val_type x) {
    return false;
  }
  static KOKKOS_FORCEINLINE_FUNCTION bool isNan (const val_type x) {
    return false;
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type abs (const val_type x) {
    // std::abs appears to work with CUDA 5.5 at least, but I'll use
    // the ternary expression for maximum generality.  Note that this
    // expression does not necessarily obey the rules for fabs() with
    // NaN input, so it should not be used for floating-point types.
    // It's perfectly fine for signed integer types, though.
    return x >= 0 ? x : -x;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type zero () {
    return 0;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type one () {
    return 1;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type min () {
    // Macros like this work with CUDA, but
    // std::numeric_limits<val_type>::min() does not, because it is
    // not marked as a __device__ function.
    return SHRT_MIN;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type max () {
    return SHRT_MAX;
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type real (const val_type x) {
    return x;
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type imag (const val_type x) {
    return 0;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type conj (const val_type x) {
    return x;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type pow (const val_type x, const val_type y) {
    return intPowSigned<val_type> (x, y);
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type sqrt (const val_type x) {
    // This will result in no loss of accuracy, though it might be
    // more expensive than it should, if we were clever about using
    // bit operations.
    return static_cast<val_type> ( ::sqrt (static_cast<float> (abs (x))));
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type log (const val_type x) {
    return static_cast<val_type> ( ::log (static_cast<float> (abs (x))));
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type log10 (const val_type x) {
    return static_cast<val_type> ( ::log10 (static_cast<float> (abs (x))));
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type epsilon () {
    return zero ();
  }

  // Backwards compatibility with Teuchos::ScalarTraits.
  typedef mag_type magnitudeType;
  typedef val_type halfPrecision;
  typedef val_type doublePrecision;

  static const bool isComplex = false;
  static const bool isOrdinal = true;
  static const bool isComparable = true;
  static const bool hasMachineParameters = false;
  static KOKKOS_FORCEINLINE_FUNCTION magnitudeType magnitude (const val_type x) {
    return abs (x);
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type conjugate (const val_type x) {
    return conj (x);
  }
  static KOKKOS_FORCEINLINE_FUNCTION bool isnaninf (const val_type) {
    return false;
  }
  static std::string name () {
    return "short";
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type squareroot (const val_type x) {
    return sqrt (x);
  }
};


template<>
class ArithTraits<unsigned short> {
public:
  typedef unsigned short val_type;
  typedef val_type mag_type;

  static const bool is_specialized = true;
  static const bool is_signed = false;
  static const bool is_integer = true;
  static const bool is_exact = true;
  static const bool is_complex = false;

  static KOKKOS_FORCEINLINE_FUNCTION bool isInf (const val_type x) {
    return false;
  }
  static KOKKOS_FORCEINLINE_FUNCTION bool isNan (const val_type x) {
    return false;
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type abs (const val_type x) {
    return x; // it's unsigned, so it's positive
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type zero () {
    return 0;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type one () {
    return 1;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type min () {
    return 0;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type max () {
    return USHRT_MAX;
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type real (const val_type x) {
    return x;
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type imag (const val_type x) {
    return 0;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type conj (const val_type x) {
    return x;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type
  pow (const val_type x, const val_type y) {
    return intPowUnsigned<val_type> (x, y);
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type sqrt (const val_type x) {
    // This will result in no loss of accuracy, though it might be
    // more expensive than it should, if we were clever about using
    // bit operations.
    return static_cast<val_type> ( ::sqrt (static_cast<float> (x)));
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type log (const val_type x) {
    return static_cast<val_type> ( ::log (static_cast<float> (x)));
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type log10 (const val_type x) {
    return static_cast<val_type> ( ::log10 (static_cast<float> (x)));
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type epsilon () {
    return zero ();
  }

  // Backwards compatibility with Teuchos::ScalarTraits.
  typedef mag_type magnitudeType;
  typedef val_type halfPrecision;
  typedef val_type doublePrecision;

  static const bool isComplex = false;
  static const bool isOrdinal = true;
  static const bool isComparable = true;
  static const bool hasMachineParameters = false;
  static KOKKOS_FORCEINLINE_FUNCTION magnitudeType magnitude (const val_type x) {
    return abs (x);
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type conjugate (const val_type x) {
    return conj (x);
  }
  static KOKKOS_FORCEINLINE_FUNCTION bool isnaninf (const val_type) {
    return false;
  }
  static std::string name () {
    return "unsigned short";
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type squareroot (const val_type x) {
    return sqrt (x);
  }
};


template<>
class ArithTraits<int> {
public:
  typedef int val_type;
  typedef val_type mag_type;

  static const bool is_specialized = true;
  static const bool is_signed = true;
  static const bool is_integer = true;
  static const bool is_exact = true;
  static const bool is_complex = false;

  static KOKKOS_FORCEINLINE_FUNCTION bool isInf (const val_type x) {
    return false;
  }
  static KOKKOS_FORCEINLINE_FUNCTION bool isNan (const val_type x) {
    return false;
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type abs (const val_type x) {
    // std::abs appears to work with CUDA 5.5 at least, but I'll use
    // the ternary expression for maximum generality.  Note that this
    // expression does not necessarily obey the rules for fabs() with
    // NaN input, so it should not be used for floating-point types.
    // It's perfectly fine for signed integer types, though.
    return x >= 0 ? x : -x;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type zero () {
    return 0;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type one () {
    return 1;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type min () {
    // Macros like INT_MIN work with CUDA, but
    // std::numeric_limits<val_type>::min() does not, because it is
    // not marked as a __device__ function.
    return INT_MIN;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type max () {
    return INT_MAX;
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type real (const val_type x) {
    return x;
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type imag (const val_type x) {
    return 0;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type conj (const val_type x) {
    return x;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type
  pow (const val_type x, const val_type y) {
    return intPowSigned<val_type> (x, y);
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type sqrt (const val_type x) {
    // This will result in no loss of accuracy, though it might be
    // more expensive than it should, if we were clever about using
    // bit operations.
    return static_cast<val_type> ( ::sqrt (static_cast<double> (abs (x))));
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type log (const val_type x) {
    return static_cast<val_type> ( ::log (static_cast<double> (abs (x))));
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type log10 (const val_type x) {
    return static_cast<val_type> ( ::log10 (static_cast<double> (abs (x))));
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type epsilon () {
    return zero ();
  }

  // Backwards compatibility with Teuchos::ScalarTraits.
  typedef mag_type magnitudeType;
  typedef val_type halfPrecision;
  typedef val_type doublePrecision;

  static const bool isComplex = false;
  static const bool isOrdinal = true;
  static const bool isComparable = true;
  static const bool hasMachineParameters = false;
  static KOKKOS_FORCEINLINE_FUNCTION magnitudeType magnitude (const val_type x) {
    return abs (x);
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type conjugate (const val_type x) {
    return conj (x);
  }
  static KOKKOS_FORCEINLINE_FUNCTION bool isnaninf (const val_type) {
    return false;
  }
  static std::string name () {
    return "int";
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type squareroot (const val_type x) {
    return sqrt (x);
  }
};


template<>
class ArithTraits<unsigned int> {
public:
  typedef unsigned int val_type;
  typedef val_type mag_type;

  static const bool is_specialized = true;
  static const bool is_signed = false;
  static const bool is_integer = true;
  static const bool is_exact = true;
  static const bool is_complex = false;

  static KOKKOS_FORCEINLINE_FUNCTION bool isInf (const val_type x) {
    return false;
  }
  static KOKKOS_FORCEINLINE_FUNCTION bool isNan (const val_type x) {
    return false;
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type abs (const val_type x) {
    return x; // it's unsigned, so it's positive
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type zero () {
    return 0;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type one () {
    return 1;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type min () {
    return 0;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type max () {
    return UINT_MAX;
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type real (const val_type x) {
    return x;
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type imag (const val_type x) {
    return 0;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type conj (const val_type x) {
    return x;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type
  pow (const val_type x, const val_type y) {
    return intPowUnsigned<val_type> (x, y);
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type sqrt (const val_type x) {
    // This will result in no loss of accuracy, though it might be
    // more expensive than it should, if we were clever about using
    // bit operations.
    return static_cast<val_type> ( ::sqrt (static_cast<double> (x)));
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type log (const val_type x) {
    return static_cast<val_type> ( ::log (static_cast<double> (x)));
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type log10 (const val_type x) {
    return static_cast<val_type> ( ::log10 (static_cast<double> (x)));
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type epsilon () {
    return zero ();
  }

  // Backwards compatibility with Teuchos::ScalarTraits.
  typedef mag_type magnitudeType;
  typedef val_type halfPrecision;
  typedef val_type doublePrecision;

  static const bool isComplex = false;
  static const bool isOrdinal = true;
  static const bool isComparable = true;
  static const bool hasMachineParameters = false;
  static KOKKOS_FORCEINLINE_FUNCTION magnitudeType magnitude (const val_type x) {
    return abs (x);
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type conjugate (const val_type x) {
    return conj (x);
  }
  static KOKKOS_FORCEINLINE_FUNCTION bool isnaninf (const val_type) {
    return false;
  }
  static std::string name () {
    return "unsigned int";
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type squareroot (const val_type x) {
    return sqrt (x);
  }
};


template<>
class ArithTraits<long> {
public:
  typedef long val_type;
  typedef val_type mag_type;

  static const bool is_specialized = true;
  static const bool is_signed = true;
  static const bool is_integer = true;
  static const bool is_exact = true;
  static const bool is_complex = false;

  static KOKKOS_FORCEINLINE_FUNCTION bool isInf (const val_type x) {
    return false;
  }
  static KOKKOS_FORCEINLINE_FUNCTION bool isNan (const val_type x) {
    return false;
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type abs (const val_type x) {
    return x >= 0 ? x : -x;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type zero () {
    return 0;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type one () {
    return 1;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type min () {
    return LONG_MIN;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type max () {
    return LONG_MAX;
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type real (const val_type x) {
    return x;
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type imag (const val_type x) {
    return 0;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type conj (const val_type x) {
    return x;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type
  pow (const val_type x, const val_type y) {
    return intPowSigned<val_type> (x, y);
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type sqrt (const val_type x) {
#ifdef __CUDA_ARCH__
    return static_cast<val_type> ( ::sqrt (static_cast<double> (abs (x))));
#else
    return static_cast<val_type> ( ::sqrt (static_cast<long double> (abs (x))));
#endif // __CUDA_ARCH__
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type log (const val_type x) {
    return static_cast<val_type> ( ::log (static_cast<double> (abs (x))));
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type log10 (const val_type x) {
    return static_cast<val_type> ( ::log10 (static_cast<double> (abs (x))));
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type epsilon () {
    return zero ();
  }

  // Backwards compatibility with Teuchos::ScalarTraits.
  typedef mag_type magnitudeType;
  typedef val_type halfPrecision;
  typedef val_type doublePrecision;

  static const bool isComplex = false;
  static const bool isOrdinal = true;
  static const bool isComparable = true;
  static const bool hasMachineParameters = false;
  static KOKKOS_FORCEINLINE_FUNCTION magnitudeType magnitude (const val_type x) {
    return abs (x);
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type conjugate (const val_type x) {
    return conj (x);
  }
  static KOKKOS_FORCEINLINE_FUNCTION bool isnaninf (const val_type) {
    return false;
  }
  static std::string name () {
    return "long";
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type squareroot (const val_type x) {
    return sqrt (x);
  }
};


template<>
class ArithTraits<unsigned long> {
public:
  typedef unsigned long val_type;
  typedef val_type mag_type;

  static const bool is_specialized = true;
  static const bool is_signed = false;
  static const bool is_integer = true;
  static const bool is_exact = true;
  static const bool is_complex = false;

  static KOKKOS_FORCEINLINE_FUNCTION bool isInf (const val_type x) {
    return false;
  }
  static KOKKOS_FORCEINLINE_FUNCTION bool isNan (const val_type x) {
    return false;
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type abs (const val_type x) {
    return x;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type zero () {
    return 0;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type one () {
    return 1;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type min () {
    return 0;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type max () {
    return ULONG_MAX;
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type real (const val_type x) {
    return x;
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type imag (const val_type x) {
    return 0;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type conj (const val_type x) {
    return x;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type pow (const val_type x, const val_type y) {
    return intPowUnsigned<val_type> (x, y);
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type sqrt (const val_type x) {
#ifdef __CUDA_ARCH__
    return static_cast<val_type> ( ::sqrt (static_cast<double> (x)));
#else
    return static_cast<val_type> ( ::sqrt (static_cast<long double> (x)));
#endif // __CUDA_ARCH__
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type log (const val_type x) {
    return static_cast<long> ( ::log (static_cast<double> (x)));
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type log10 (const val_type x) {
    return static_cast<long> ( ::log10 (static_cast<double> (x)));
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type epsilon () {
    return zero ();
  }

  // Backwards compatibility with Teuchos::ScalarTraits.
  typedef mag_type magnitudeType;
  typedef val_type halfPrecision;
  typedef val_type doublePrecision;

  static const bool isComplex = false;
  static const bool isOrdinal = true;
  static const bool isComparable = true;
  static const bool hasMachineParameters = false;
  static KOKKOS_FORCEINLINE_FUNCTION magnitudeType magnitude (const val_type x) {
    return abs (x);
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type conjugate (const val_type x) {
    return conj (x);
  }
  static KOKKOS_FORCEINLINE_FUNCTION bool isnaninf (const val_type) {
    return false;
  }
  static std::string name () {
    return "unsigned long";
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type squareroot (const val_type x) {
    return sqrt (x);
  }
};


template<>
class ArithTraits<long long> {
public:
  typedef long long val_type;
  typedef val_type mag_type;

  static const bool is_specialized = true;
  static const bool is_signed = true;
  static const bool is_integer = true;
  static const bool is_exact = true;
  static const bool is_complex = false;

  static KOKKOS_FORCEINLINE_FUNCTION bool isInf (const val_type x) {
    return false;
  }
  static KOKKOS_FORCEINLINE_FUNCTION bool isNan (const val_type x) {
    return false;
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type abs (const val_type x) {
    return x >= 0 ? x : -x;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type zero () {
    return 0;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type one () {
    return 1;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type min () {
    return LLONG_MIN;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type max () {
    return LLONG_MAX;
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type real (const val_type x) {
    return x;
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type imag (const val_type x) {
    return 0;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type conj (const val_type x) {
    return x;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type
  pow (const val_type x, const val_type y) {
    return intPowSigned<val_type> (x, y);
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type sqrt (const val_type x) {
#ifdef __CUDA_ARCH__
    // Casting from a 64-bit integer type to double does result in a
    // loss of accuracy.  However, it gives us a good first
    // approximation.  For very large numbers, we may lose some
    // significand bits, but will always get within a factor of two
    // (assuming correct rounding) of the exact double-precision
    // number.  We could then binary search between half the result
    // and twice the result (assuming the latter is <= INT64_MAX,
    // which it has to be, so we don't have to check) to ensure
    // correctness.  It actually should suffice to check numbers
    // within 1 of the result.
    return static_cast<val_type> ( ::sqrt (static_cast<double> (abs (x))));
#else
    // IEEE 754 promises that long double has at least 64 significand
    // bits, so we can use it to represent any signed or unsigned
    // 64-bit integer type exactly.  However, CUDA does not implement
    // long double for device functions.
    return static_cast<val_type> ( ::sqrt (static_cast<long double> (abs (x))));
#endif // __CUDA_ARCH__
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type log (const val_type x) {
    return static_cast<val_type> ( ::log (static_cast<double> (abs (x))));
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type log10 (const val_type x) {
    return static_cast<val_type> ( ::log10 (static_cast<double> (abs (x))));
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type epsilon () {
    return zero ();
  }

  // Backwards compatibility with Teuchos::ScalarTraits.
  typedef mag_type magnitudeType;
  typedef val_type halfPrecision;
  typedef val_type doublePrecision;

  static const bool isComplex = false;
  static const bool isOrdinal = true;
  static const bool isComparable = true;
  static const bool hasMachineParameters = false;
  static KOKKOS_FORCEINLINE_FUNCTION magnitudeType magnitude (const val_type x) {
    return abs (x);
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type conjugate (const val_type x) {
    return conj (x);
  }
  static KOKKOS_FORCEINLINE_FUNCTION bool isnaninf (const val_type) {
    return false;
  }
  static std::string name () {
    return "long long";
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type squareroot (const val_type x) {
    return sqrt (x);
  }
};


template<>
class ArithTraits<unsigned long long> {
public:
  typedef unsigned long long val_type;
  typedef val_type mag_type;

  static const bool is_specialized = true;
  static const bool is_signed = false;
  static const bool is_integer = true;
  static const bool is_exact = true;
  static const bool is_complex = false;

  static KOKKOS_FORCEINLINE_FUNCTION bool isInf (const val_type x) {
    return false;
  }
  static KOKKOS_FORCEINLINE_FUNCTION bool isNan (const val_type x) {
    return false;
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type abs (const val_type x) {
    return x; // unsigned integers are always nonnegative
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type zero () {
    return 0;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type one () {
    return 1;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type min () {
    return 0;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type max () {
    return ULLONG_MAX ;
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type real (const val_type x) {
    return x;
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type imag (const val_type x) {
    return 0;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type conj (const val_type x) {
    return x;
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type
  pow (const val_type x, const val_type y) {
    return intPowUnsigned<val_type> (x, y);
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type sqrt (const val_type x) {
#ifdef __CUDA_ARCH__
    return static_cast<val_type> ( ::sqrt (static_cast<double> (x)));
#else
    return static_cast<val_type> ( ::sqrt (static_cast<long double> (x)));
#endif // __CUDA_ARCH__
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type log (const val_type x) {
    return static_cast<val_type> ( ::log (static_cast<double> (x)));
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type log10 (const val_type x) {
    return static_cast<val_type> ( ::log10 (static_cast<double> (x)));
  }
  static KOKKOS_FORCEINLINE_FUNCTION mag_type epsilon () {
    return zero ();
  }

  // Backwards compatibility with Teuchos::ScalarTraits.
  typedef mag_type magnitudeType;
  typedef val_type halfPrecision;
  typedef val_type doublePrecision;

  static const bool isComplex = false;
  static const bool isOrdinal = true;
  static const bool isComparable = true;
  static const bool hasMachineParameters = false;
  static KOKKOS_FORCEINLINE_FUNCTION magnitudeType magnitude (const val_type x) {
    return abs (x);
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type conjugate (const val_type x) {
    return conj (x);
  }
  static KOKKOS_FORCEINLINE_FUNCTION bool isnaninf (const val_type) {
    return false;
  }
  static std::string name () {
    return "unsigned long long";
  }
  static KOKKOS_FORCEINLINE_FUNCTION val_type squareroot (const val_type x) {
    return sqrt (x);
  }
};

// dd_real and qd_real are floating-point types provided by the QD
// library of David Bailey (LBNL):
//
// http://crd-legacy.lbl.gov/~dhbailey/mpdist/
//
// dd_real uses two doubles (128 bits), and qd_real uses four doubles
// (256 bits).
//
// Kokkos does <i>not</i> currently support these types in device
// functions.  It should be possible to use Kokkos' support for
// aggregate types to implement device function support for dd_real
// and qd_real, but we have not done this yet (as of 07 Jan 2014).
// Hence, the class methods of the ArithTraits specializations for
// dd_real and qd_real are not marked as device functions.
#ifdef HAVE_KOKKOS_QD
template<>
struct ScalarTraits<dd_real>
{
  typedef dd_real val_type;
  typedef dd_real mag_type;

  static const bool is_specialized = true;
  static const bool is_signed = true;
  static const bool is_integer = false;
  static const bool is_exact = false;
  static const bool is_complex = false;

  static inline bool isInf (const val_type& x) {
    return isinf (x);
  }
  static inline bool isNan (const val_type& x) {
    return isnan (x);
  }
  static inline mag_type abs (const val_type& x) {
    return ::abs (x);
  }
  static inline val_type zero () {
    return val_type (0.0);
  }
  static inline val_type one () {
    return val_type (1.0);
  }
  static inline val_type min () {
    return std::numeric_limits<val_type>::min ();
  }
  static inline val_type max () {
    return std::numeric_limits<val_type>::max ();
  }
  static inline mag_type real (const val_type& x) {
    return x;
  }
  static inline mag_type imag (const val_type&) {
    return zero ();
  }
  static inline val_type conj (const val_type& x) {
    return x;
  }
  static inline val_type pow (const val_type& x, const val_type& y) {
    return ::pow(x,y);
  }
  static inline val_type sqrt (const val_type& x) {
      return ::sqrt (x);
  }
  static inline val_type log (const val_type& x) {
    // dd_real puts its transcendental functions in the global namespace.
    return ::log (x);
  }
  static inline val_type log10 (const val_type& x) {
    return ::log10 (x);
  }
  static inline val_type nan () {
    return val_type::_nan;
  }
  static val_type epsilon () {
    return std::numeric_limits<val_type>::epsilon ();
  }

  typedef dd_real magnitudeType;
  typedef double halfPrecision;
  typedef qd_real doublePrecision;

  static const bool isComplex = false;
  static const bool isOrdinal = false;
  static const bool isComparable = true;
  static const bool hasMachineParameters = true;

  static mag_type eps () {
    return epsilon ();
  }
  static mag_type sfmin () {
    return min ();
  }
  static int base ()  {
    return std::numeric_limits<val_type>::radix;
  }
  static mag_type prec () {
    return eps () * base ();
  }
  static int t () {
    return std::numeric_limits<val_type>::digits;
  }
  static mag_type rnd () {
    return std::numeric_limits<val_type>::round_style == std::round_to_nearest ?
      one () :
      zero ();
  }
  static int emin () {
    return std::numeric_limits<val_type>::min_exponent;
  }
  static mag_type rmin () {
    return std::numeric_limits<val_type>::min ();
  }
  static int emax () {
    return std::numeric_limits<val_type>::max_exponent;
  }
  static mag_type rmax () {
    return std::numeric_limits<val_type>::max ();
  }
  static mag_type magnitude (const val_type& x) {
    return ::abs (x);
  }
  static val_type conjugate (const val_type& x) {
    return conj (x);
  }
  static bool isnaninf (const val_type& x) {
    return isNan (x) || isInf (x);
  }
  static std::string name () { return "dd_real"; }
  static val_type squareroot (const val_type& x) {
    return ::sqrt (x);
  }
};


template<>
struct ScalarTraits<qd_real>
{
  typedef qd_real val_type;
  typedef qd_real mag_type;

  static const bool is_specialized = true;
  static const bool is_signed = true;
  static const bool is_integer = false;
  static const bool is_exact = false;
  static const bool is_complex = false;

  static inline bool isInf (const val_type& x) {
    return isinf (x);
  }
  static inline bool isNan (const val_type& x) {
    return isnan (x);
  }
  static inline mag_type abs (const val_type& x) {
    return ::abs (x);
  }
  static inline val_type zero () {
    return val_type (0.0);
  }
  static inline val_type one () {
    return val_type (1.0);
  }
  static inline val_type min () {
    return std::numeric_limits<val_type>::min ();
  }
  static inline val_type max () {
    return std::numeric_limits<val_type>::max ();
  }
  static inline mag_type real (const val_type& x) {
    return x;
  }
  static inline mag_type imag (const val_type&) {
    return zero ();
  }
  static inline val_type conj (const val_type& x) {
    return x;
  }
  static inline val_type pow (const val_type& x, const val_type& y) {
    return ::pow (x, y);
  }
  static inline val_type sqrt (const val_type& x) {
      return ::sqrt (x);
  }
  static inline val_type log (const val_type& x) {
    // val_type puts its transcendental functions in the global namespace.
    return ::log (x);
  }
  static inline val_type log10 (const val_type& x) {
    return ::log10 (x);
  }
  static inline val_type nan () {
    return val_type::_nan;
  }
  static inline val_type epsilon () {
    return std::numeric_limits<val_type>::epsilon ();
  }

  typedef qd_real magnitudeType;
  typedef dd_real halfPrecision;
  // The QD library does not have an "oct-double real" class.  One
  // could use an arbitrary-precision library like MPFR or ARPREC,
  // with the precision set appropriately, to get an
  // extended-precision type for qd_real.
  typedef qd_real doublePrecision;

  static const bool isComplex = false;
  static const bool isOrdinal = false;
  static const bool isComparable = true;
  static const bool hasMachineParameters = true;

  static mag_type eps () {
    return epsilon ();
  }
  static mag_type sfmin () {
    return min ();
  }
  static int base ()  {
    return std::numeric_limits<val_type>::radix;
  }
  static mag_type prec () {
    return eps () * base ();
  }
  static int t () {
    return std::numeric_limits<val_type>::digits;
  }
  static mag_type rnd () {
    return std::numeric_limits<val_type>::round_style == std::round_to_nearest ?
      one () :
      zero ();
  }
  static int emin () {
    return std::numeric_limits<val_type>::min_exponent;
  }
  static mag_type rmin () {
    return std::numeric_limits<val_type>::min ();
  }
  static int emax () {
    return std::numeric_limits<val_type>::max_exponent;
  }
  static mag_type rmax () {
    return std::numeric_limits<val_type>::max ();
  }
  static mag_type magnitude (const val_type& x) {
    return ::abs (x);
  }
  static val_type conjugate (const val_type& x) {
    return conj (x);
  }
  static bool isnaninf (const val_type& x) {
    return isNan (x) || isInf (x);
  }
  static std::string name () { return "qd_real"; }
  static val_type squareroot (const val_type& x) {
    return ::sqrt (x);
  }
};
#endif // HAVE_KOKKOS_QD

} // namespace Details
} // namespace Kokkos

#endif // KOKKOS_ARITHTRAITS_HPP