reduce_quda.cu

Go to the documentation of this file.

#include <blas_quda.h>
#include <tune_quda.h>
#include <float_vector.h>

#if (__COMPUTE_CAPABILITY__ >= 130)
#define QudaSumFloat double
#define QudaSumFloat2 double2
#define QudaSumFloat3 double3
#else
#define QudaSumFloat doublesingle
#define QudaSumFloat2 doublesingle2
#define QudaSumFloat3 doublesingle3
#include <double_single.h>
#endif

#define REDUCE_MAX_BLOCKS 65536

#define checkSpinor(a, b)                                               \
  {                                                                     \
    if (a.Precision() != b.Precision())                                 \
      errorQuda("precisions do not match: %d %d", a.Precision(), b.Precision()); \
    if (a.Length() != b.Length())                                       \
      errorQuda("lengths do not match: %d %d", a.Length(), b.Length()); \
    if (a.Stride() != b.Stride())                                       \
      errorQuda("strides do not match: %d %d", a.Stride(), b.Stride()); \
  }

#define checkLength(a, b)                                               \
  {                                                                     \
    if (a.Length() != b.Length())                                       \
      errorQuda("lengths do not match: %d %d", a.Length(), b.Length()); \
    if (a.Stride() != b.Stride())                                       \
      errorQuda("strides do not match: %d %d", a.Stride(), b.Stride()); \
  }

static struct {
  const char *vol_str;
  const char *aux_str;
  char aux_tmp[quda::TuneKey::aux_n];
} blasStrings;

// These are used for reduction kernels
static QudaSumFloat *d_reduce=0;
static QudaSumFloat *h_reduce=0;
static QudaSumFloat *hd_reduce=0;
static cudaEvent_t reduceEnd;
    
namespace quda {

  cudaStream_t* getBlasStream();
    
  void initReduce()
  { 

    const int MaxReduce = 12;
    // reduction buffer size
    size_t bytes = MaxReduce*3*REDUCE_MAX_BLOCKS*sizeof(QudaSumFloat); // Factor of N for composite reductions


    if (!d_reduce) d_reduce = (QudaSumFloat *) device_malloc(bytes);
    
    // these arrays are actually oversized currently (only needs to be QudaSumFloat3)
    
    // if the device supports host-mapped memory then use a host-mapped array for the reduction
    if (!h_reduce) {
      // only use zero copy reductions when using 64-bit
#if (defined(_MSC_VER) && defined(_WIN64)) || defined(__LP64__)
      if(deviceProp.canMapHostMemory) {
        h_reduce = (QudaSumFloat *) mapped_malloc(bytes);       
        cudaHostGetDevicePointer(&hd_reduce, h_reduce, 0); // set the matching device pointer
      } else 
#endif 
      {
        h_reduce = (QudaSumFloat *) pinned_malloc(bytes);
        hd_reduce = d_reduce;
      }
      memset(h_reduce, 0, bytes); // added to ensure that valgrind doesn't report h_reduce is unitialised
    }
    
    cudaEventCreateWithFlags(&reduceEnd, cudaEventDisableTiming);
    
    checkCudaError();
  }

  void endReduce(void)
  {
    if (d_reduce) {
      device_free(d_reduce);
      d_reduce = 0;
    }
    if (h_reduce) {
      host_free(h_reduce);
      h_reduce = 0;
    }
    hd_reduce = 0;
    
    cudaEventDestroy(reduceEnd);
  }

  namespace reduce {

#include <texture.h>
#include <reduce_core.h>
#include <reduce_mixed_core.h>
    
  } // namespace reduce

  template <typename ReduceType, typename Float2, typename FloatN>
  struct ReduceFunctor {
    
    virtual __device__ void pre() { ; }
    
    virtual __device__ void operator()(ReduceType &sum, FloatN &x, FloatN &y, 
                                       FloatN &z, FloatN &w, FloatN &v) = 0;
    
    virtual __device__ void post(ReduceType &sum) { ; }
    
  };

  __device__ double norm2_(const double2 &a) { return a.x*a.x + a.y*a.y; }
  __device__ float norm2_(const float2 &a) { return a.x*a.x + a.y*a.y; }
  __device__ float norm2_(const float4 &a) { return a.x*a.x + a.y*a.y + a.z*a.z + a.w*a.w; }

  template <typename ReduceType, typename Float2, typename FloatN>
#if (__COMPUTE_CAPABILITY__ >= 200)
  struct Norm2 : public ReduceFunctor<ReduceType, Float2, FloatN> {
#else
  struct Norm2 {
#endif
    Norm2(const Float2 &a, const Float2 &b) { ; }
    __device__ void operator()(ReduceType &sum, FloatN &x, FloatN &y, FloatN &z,FloatN  &w, FloatN &v) { sum += norm2_(x); }
    static int streams() { return 1; } 
    static int flops() { return 2; } 
  };

  double normCuda(const cudaColorSpinorField &x) {
    cudaColorSpinorField &y = (cudaColorSpinorField&)x; // FIXME
    return reduce::reduceCuda<double,QudaSumFloat,QudaSumFloat,Norm2,0,0,0,0,0,false>
      (make_double2(0.0, 0.0), make_double2(0.0, 0.0), y, y, y, y, y);
  }

  __device__ double dot_(const double2 &a, const double2 &b) { return a.x*b.x + a.y*b.y; }
  __device__ float dot_(const float2 &a, const float2 &b) { return a.x*b.x + a.y*b.y; }
  __device__ float dot_(const float4 &a, const float4 &b) { return a.x*b.x + a.y*b.y + a.z*b.z + a.w*b.w; }

  template <typename ReduceType, typename Float2, typename FloatN>
#if (__COMPUTE_CAPABILITY__ >= 200)
  struct Dot : public ReduceFunctor<ReduceType, Float2, FloatN> {
#else
  struct Dot {
#endif
    Dot(const Float2 &a, const Float2 &b) { ; }
    __device__ void operator()(ReduceType &sum, FloatN &x, FloatN &y, FloatN &z, FloatN &w, FloatN &v) { sum += dot_(x,y); }
    static int streams() { return 2; } 
    static int flops() { return 2; } 
  };

  double reDotProductCuda(cudaColorSpinorField &x, cudaColorSpinorField &y) {
    return reduce::reduceCuda<double,QudaSumFloat,QudaSumFloat,Dot,0,0,0,0,0,false>
      (make_double2(0.0, 0.0), make_double2(0.0, 0.0), x, y, x, x, x);
  }

 
  void reDotProductCuda(double* result, std::vector<cudaColorSpinorField*>& x, std::vector<cudaColorSpinorField*>& y){
#ifndef SSTEP
    errorQuda("S-step code not built\n");
#else 
    switch(x.size()){
      case 1:
        reduce::multiReduceCuda<1,double,QudaSumFloat,QudaSumFloat,Dot,0,0,0,0,0,false>
        (result, make_double2(0.0, 0.0), make_double2(0.0, 0.0), x, y, x, x, x);
        break;
      case 2:
        reduce::multiReduceCuda<2,double,QudaSumFloat,QudaSumFloat,Dot,0,0,0,0,0,false>
        (result, make_double2(0.0, 0.0), make_double2(0.0, 0.0), x, y, x, x, x);
        break; 
      case 3:
        reduce::multiReduceCuda<3,double,QudaSumFloat,QudaSumFloat,Dot,0,0,0,0,0,false>
        (result, make_double2(0.0, 0.0), make_double2(0.0, 0.0), x, y, x, x, x);
        break;
      case 4:
        reduce::multiReduceCuda<4,double,QudaSumFloat,QudaSumFloat,Dot,0,0,0,0,0,false>
        (result, make_double2(0.0, 0.0), make_double2(0.0, 0.0), x, y, x, x, x);
        break;
      case 5:
        reduce::multiReduceCuda<5,double,QudaSumFloat,QudaSumFloat,Dot,0,0,0,0,0,false>
        (result, make_double2(0.0, 0.0), make_double2(0.0, 0.0), x, y, x, x, x);
        break;
      case 6:
        reduce::multiReduceCuda<6,double,QudaSumFloat,QudaSumFloat,Dot,0,0,0,0,0,false>
        (result, make_double2(0.0, 0.0), make_double2(0.0, 0.0), x, y, x, x, x);
        break;
      case 7:
        reduce::multiReduceCuda<7,double,QudaSumFloat,QudaSumFloat,Dot,0,0,0,0,0,false>
        (result, make_double2(0.0, 0.0), make_double2(0.0, 0.0), x, y, x, x, x);
        break;
      case 8:
        reduce::multiReduceCuda<8,double,QudaSumFloat,QudaSumFloat,Dot,0,0,0,0,0,false>
        (result, make_double2(0.0, 0.0), make_double2(0.0, 0.0), x, y, x, x, x);
        break;
      case 9:
        reduce::multiReduceCuda<9,double,QudaSumFloat,QudaSumFloat,Dot,0,0,0,0,0,false>
        (result, make_double2(0.0, 0.0), make_double2(0.0, 0.0), x, y, x, x, x);
        break;
      case 10:
        reduce::multiReduceCuda<10,double,QudaSumFloat,QudaSumFloat,Dot,0,0,0,0,0,false>
        (result, make_double2(0.0, 0.0), make_double2(0.0, 0.0), x, y, x, x, x);
        break;
      case 11:
        reduce::multiReduceCuda<11,double,QudaSumFloat,QudaSumFloat,Dot,0,0,0,0,0,false>
        (result, make_double2(0.0, 0.0), make_double2(0.0, 0.0), x, y, x, x, x);
        break;
      case 12:
        reduce::multiReduceCuda<12,double,QudaSumFloat,QudaSumFloat,Dot,0,0,0,0,0,false>
        (result, make_double2(0.0, 0.0), make_double2(0.0, 0.0), x, y, x, x, x);
        break;
      case 13:
        reduce::multiReduceCuda<13,double,QudaSumFloat,QudaSumFloat,Dot,0,0,0,0,0,false>
        (result, make_double2(0.0, 0.0), make_double2(0.0, 0.0), x, y, x, x, x);
        break;
      case 14:
        reduce::multiReduceCuda<14,double,QudaSumFloat,QudaSumFloat,Dot,0,0,0,0,0,false>
        (result, make_double2(0.0, 0.0), make_double2(0.0, 0.0), x, y, x, x, x);
        break;
      case 15:
        reduce::multiReduceCuda<15,double,QudaSumFloat,QudaSumFloat,Dot,0,0,0,0,0,false>
        (result, make_double2(0.0, 0.0), make_double2(0.0, 0.0), x, y, x, x, x);
        break;
      case 16:
        reduce::multiReduceCuda<16,double,QudaSumFloat,QudaSumFloat,Dot,0,0,0,0,0,false>
        (result, make_double2(0.0, 0.0), make_double2(0.0, 0.0), x, y, x, x, x);
        break;
      case 17:
        reduce::multiReduceCuda<17,double,QudaSumFloat,QudaSumFloat,Dot,0,0,0,0,0,false>
        (result, make_double2(0.0, 0.0), make_double2(0.0, 0.0), x, y, x, x, x);
        break;
      case 18:
        reduce::multiReduceCuda<18,double,QudaSumFloat,QudaSumFloat,Dot,0,0,0,0,0,false>
        (result, make_double2(0.0, 0.0), make_double2(0.0, 0.0), x, y, x, x, x);
        break;
      case 19:
        reduce::multiReduceCuda<19,double,QudaSumFloat,QudaSumFloat,Dot,0,0,0,0,0,false>
        (result, make_double2(0.0, 0.0), make_double2(0.0, 0.0), x, y, x, x, x);
        break;
      case 20:
        reduce::multiReduceCuda<20,double,QudaSumFloat,QudaSumFloat,Dot,0,0,0,0,0,false>
        (result, make_double2(0.0, 0.0), make_double2(0.0, 0.0), x, y, x, x, x);
        break;
      default:
        errorQuda("Unsupported vector size");
        break;
    }
#endif // SSTEP
  } 


    /* 
      returns the real component of the dot product of a and b 
      and the norm of a
    */
  __device__ double2 dotNormA_(const double2 &a, const double2 &b)
  { return make_double2(a.x*b.x + a.y*b.y, a.x*a.x + a.y*a.y); }
 
  __device__ double2 dotNormA_(const float2 &a, const float2 &b)
  { return make_double2(a.x*b.x + a.y*b.y, a.x*a.x + a.y*a.y); }

 
  __device__ double2 dotNormA_(const float4 &a, const float4 & b)
  { return make_double2(a.x*b.x + a.y*b.y + a.z*b.z + a.w*b.w, a.x*a.x + a.y*a.y + a.z*a.z +     a.w*a.w); }



  template <typename ReduceType, typename Float2, typename FloatN>
#if (__COMPUTE_CAPABILITY__ >= 200)
  struct DotNormA : public ReduceFunctor<ReduceType, Float2, FloatN> {
#else
  struct DotNormA {
#endif
    DotNormA(const Float2 &a, const Float2 &b){}
    __device__ void operator()(ReduceType &sum, FloatN &x, FloatN &y, FloatN &z,  FloatN &w, FloatN &v){sum += dotNormA_(x,y);}
    static int streams() { return 2; }
    static int flops() { return 4; }
  };

  double2 reDotProductNormACuda(cudaColorSpinorField &x,cudaColorSpinorField &y){
    return reduce::reduceCuda<double2,QudaSumFloat2,QudaSumFloat,DotNormA,0,0,0,0,0,false>
      (make_double2(0.0, 0.0), make_double2(0.0, 0.0), x, y, x, x, x);
  }


  template <typename ReduceType, typename Float2, typename FloatN>
#if (__COMPUTE_CAPABILITY__ >= 200)
  struct axpyNorm2 : public ReduceFunctor<ReduceType, Float2, FloatN> {
#else
  struct axpyNorm2 {
#endif
    Float2 a;
    axpyNorm2(const Float2 &a, const Float2 &b) : a(a) { ; }
    __device__ void operator()(ReduceType &sum, FloatN &x, FloatN &y, FloatN &z, FloatN &w, FloatN &v) { 
      y += a.x*x; sum += norm2_(y); }
    static int streams() { return 3; } 
    static int flops() { return 4; } 
  };

  double axpyNormCuda(const double &a, cudaColorSpinorField &x, cudaColorSpinorField &y) {
    return reduce::reduceCuda<double,QudaSumFloat,QudaSumFloat,axpyNorm2,0,1,0,0,0,false>
      (make_double2(a, 0.0), make_double2(0.0, 0.0), x, y, x, x, x);
  }

  template <typename ReduceType, typename Float2, typename FloatN>
#if (__COMPUTE_CAPABILITY__ >= 200)
  struct xmyNorm2 : public ReduceFunctor<ReduceType, Float2, FloatN> {
#else
  struct xmyNorm2 {
#endif
    xmyNorm2(const Float2 &a, const Float2 &b) { ; }
    __device__ void operator()(ReduceType &sum, FloatN &x, FloatN &y, FloatN &z, FloatN &w, FloatN &v) { 
      y = x - y; sum += norm2_(y); }
    static int streams() { return 3; } 
    static int flops() { return 3; } 
  };

  double xmyNormCuda(cudaColorSpinorField &x, cudaColorSpinorField &y) {
    return reduce::reduceCuda<double,QudaSumFloat,QudaSumFloat,xmyNorm2,0,1,0,0,0,false>
      (make_double2(0.0, 0.0), make_double2(0.0, 0.0), x, y, x, x, x);
  }


  __device__ void Caxpy_(const float2 &a, const float4 &x, float4 &y) {
    y.x += a.x*x.x; y.x -= a.y*x.y;
    y.y += a.y*x.x; y.y += a.x*x.y;
    y.z += a.x*x.z; y.z -= a.y*x.w;
    y.w += a.y*x.z; y.w += a.x*x.w;
  }

  __device__ void Caxpy_(const float2 &a, const float2 &x, float2 &y) {
    y.x += a.x*x.x; y.x -= a.y*x.y;
    y.y += a.y*x.x; y.y += a.x*x.y;
  }

  __device__ void Caxpy_(const double2 &a, const double2 &x, double2 &y) {
    y.x += a.x*x.x; y.x -= a.y*x.y;
    y.y += a.y*x.x; y.y += a.x*x.y;
  }

  template <typename ReduceType, typename Float2, typename FloatN>
#if (__COMPUTE_CAPABILITY__ >= 200)
  struct caxpyNorm2 : public ReduceFunctor<ReduceType, Float2, FloatN> {
#else
  struct caxpyNorm2 {
#endif
    Float2 a;
    caxpyNorm2(const Float2 &a, const Float2 &b) : a(a) { ; }
    __device__ void operator()(ReduceType &sum, FloatN &x, FloatN &y, FloatN &z, FloatN &w, FloatN &v) { 
      Caxpy_(a, x, y); sum += norm2_(y); }
    static int streams() { return 3; } 
    static int flops() { return 6; } 
  };

  double caxpyNormCuda(const Complex &a, cudaColorSpinorField &x, cudaColorSpinorField &y) {
    return reduce::reduceCuda<double,QudaSumFloat,QudaSumFloat,caxpyNorm2,0,1,0,0,0,false>
      (make_double2(REAL(a), IMAG(a)), make_double2(0.0, 0.0), x, y, x, x, x);
  }

  template <typename ReduceType, typename Float2, typename FloatN>
#if (__COMPUTE_CAPABILITY__ >= 200)
  struct caxpyxmaznormx : public ReduceFunctor<ReduceType, Float2, FloatN> {
#else
  struct caxpyxmaznormx {
#endif
    Float2 a;
    caxpyxmaznormx(const Float2 &a, const Float2 &b) : a(a) { ; }
    __device__ void operator()(ReduceType &sum, FloatN &x, FloatN &y, FloatN &z, FloatN &w, FloatN &v) { Caxpy_(a, x, y); x-= a.x*z; sum += norm2_(x); }
    static int streams() { return 5; } 
    static int flops() { return 10; } 
  };

  double caxpyXmazNormXCuda(const Complex &a, cudaColorSpinorField &x, 
                            cudaColorSpinorField &y, cudaColorSpinorField &z) {
    return reduce::reduceCuda<double,QudaSumFloat,QudaSumFloat,caxpyxmaznormx,1,1,0,0,0,false>
      (make_double2(REAL(a), IMAG(a)), make_double2(0.0, 0.0), x, y, z, x, x);
  }

    template <typename ReduceType, typename Float2, typename FloatN>
#if (__COMPUTE_CAPABILITY__ >= 200)
    struct cabxpyaxnorm : public ReduceFunctor<ReduceType, Float2, FloatN> {
#else
    struct cabxpyaxnorm {
#endif
    Float2 a;
    Float2 b;
    cabxpyaxnorm(const Float2 &a, const Float2 &b) : a(a), b(b) { ; }
      __device__ void operator()(ReduceType &sum, FloatN &x, FloatN &y, FloatN &z, FloatN &w, FloatN &v) { x *= a.x; Caxpy_(b, x, y); sum += norm2_(y); }
    static int streams() { return 4; } 
    static int flops() { return 10; } 
  };

  double cabxpyAxNormCuda(const double &a, const Complex &b, 
                          cudaColorSpinorField &x, cudaColorSpinorField &y) {
    return reduce::reduceCuda<double,QudaSumFloat,QudaSumFloat,cabxpyaxnorm,1,1,0,0,0,false>
      (make_double2(a, 0.0), make_double2(REAL(b), IMAG(b)), x, y, x, x, x);
  }

  __device__ double2 cdot_(const double2 &a, const double2 &b) 
  { return make_double2(a.x*b.x + a.y*b.y, a.x*b.y - a.y*b.x); }
  __device__ double2 cdot_(const float2 &a, const float2 &b) 
  { return make_double2(a.x*b.x + a.y*b.y, a.x*b.y - a.y*b.x); }
  __device__ double2 cdot_(const float4 &a, const float4 &b) 
  { return make_double2(a.x*b.x + a.y*b.y + a.z*b.z + a.w*b.w, a.x*b.y - a.y*b.x + a.z*b.w - a.w*b.z); }

  template <typename ReduceType, typename Float2, typename FloatN>
#if (__COMPUTE_CAPABILITY__ >= 200)
  struct Cdot : public ReduceFunctor<ReduceType, Float2, FloatN> {
#else
  struct Cdot {
#endif
    Cdot(const Float2 &a, const Float2 &b) { ; }
    __device__ void operator()(ReduceType &sum, FloatN &x, FloatN &y, FloatN &z, FloatN &w, FloatN &v) { sum += cdot_(x,y); }
    static int streams() { return 2; } 
    static int flops() { return 4; } 
  };

  Complex cDotProductCuda(cudaColorSpinorField &x, cudaColorSpinorField &y) {
    double2 cdot = reduce::reduceCuda<double2,QudaSumFloat2,QudaSumFloat,Cdot,0,0,0,0,0,false>
      (make_double2(0.0, 0.0), make_double2(0.0, 0.0), x, y, x, x, x);
    return Complex(cdot.x, cdot.y);
  }

  void cDotProductCuda(Complex* result, std::vector<cudaColorSpinorField*>& x, std::vector<cudaColorSpinorField*>& y){
#ifndef SSTEP
    errorQuda("S-step code not built\n");
#else
    double2* cdot = new double2[x.size()];

    switch(x.size()){
      case 1:
        reduce::multiReduceCuda<1,double2,QudaSumFloat2,QudaSumFloat,Cdot,0,0,0,0,0,false>
        (cdot, make_double2(0.0, 0.0), make_double2(0.0, 0.0), x, y, x, x, x);
        break;
      case 6:
        reduce::multiReduceCuda<6,double2,QudaSumFloat2,QudaSumFloat,Cdot,0,0,0,0,0,false>
        (cdot, make_double2(0.0, 0.0), make_double2(0.0, 0.0), x, y, x, x, x);
        break;
      case 10:
        reduce::multiReduceCuda<10,double2,QudaSumFloat2,QudaSumFloat,Cdot,0,0,0,0,0,false>
        (cdot, make_double2(0.0, 0.0), make_double2(0.0, 0.0), x, y, x, x, x);
        break;
      case 14:
        reduce::multiReduceCuda<14,double2,QudaSumFloat2,QudaSumFloat,Cdot,0,0,0,0,0,false>
        (cdot, make_double2(0.0, 0.0), make_double2(0.0, 0.0), x, y, x, x, x);
        break;
      case 18:
        reduce::multiReduceCuda<18,double2,QudaSumFloat2,QudaSumFloat,Cdot,0,0,0,0,0,false>
        (cdot, make_double2(0.0, 0.0), make_double2(0.0, 0.0), x, y, x, x, x);
        break;
      case 22:
        reduce::multiReduceCuda<22,double2,QudaSumFloat2,QudaSumFloat,Cdot,0,0,0,0,0,false>
        (cdot, make_double2(0.0, 0.0), make_double2(0.0, 0.0), x, y, x, x, x);
        break;
      default:
        errorQuda("Unsupported vector size\n");
        break;
    }

    for(int i=0; i<x.size(); ++i) result[i] = Complex(cdot[i].x,cdot[i].y);
    delete[] cdot;
#endif
  }

  template <typename ReduceType, typename Float2, typename FloatN>
#if (__COMPUTE_CAPABILITY__ >= 200)
  struct xpaycdotzy : public ReduceFunctor<ReduceType, Float2, FloatN> {
#else
  struct xpaycdotzy {
#endif
    Float2 a;
    xpaycdotzy(const Float2 &a, const Float2 &b) : a(a) { ; }
    __device__ void operator()(ReduceType &sum, FloatN &x, FloatN &y, FloatN &z, FloatN &w, FloatN &v) { y = x + a.x*y; sum += cdot_(z,y); }
    static int streams() { return 4; } 
    static int flops() { return 6; } 
  };

  Complex xpaycDotzyCuda(cudaColorSpinorField &x, const double &a, cudaColorSpinorField &y, cudaColorSpinorField &z) {
    double2 cdot = reduce::reduceCuda<double2,QudaSumFloat2,QudaSumFloat,xpaycdotzy,0,1,0,0,0,false>
      (make_double2(a, 0.0), make_double2(0.0, 0.0), x, y, z, x, x);
    return Complex(cdot.x, cdot.y);
  }

  template <typename ReduceType, typename Float2, typename FloatN>
#if (__COMPUTE_CAPABILITY__ >= 200)
  struct caxpydotzy : public ReduceFunctor<ReduceType, Float2, FloatN> {
#else
  struct caxpydotzy {
#endif
    Float2 a;
    caxpydotzy(const Float2 &a, const Float2 &b) : a(a) { ; }
    __device__ void operator()(ReduceType &sum, FloatN &x, FloatN &y, FloatN &z, FloatN &w, FloatN &v) { Caxpy_(a, x, y); sum += cdot_(z,y); }
    static int streams() { return 4; } 
    static int flops() { return 8; } 
  };

  Complex caxpyDotzyCuda(const Complex &a, cudaColorSpinorField &x, cudaColorSpinorField &y,
                         cudaColorSpinorField &z) {
    double2 cdot = reduce::reduceCuda<double2,QudaSumFloat2,QudaSumFloat,caxpydotzy,0,1,0,0,0,false>
      (make_double2(REAL(a), IMAG(a)), make_double2(0.0, 0.0), x, y, z, x, x);
    return Complex(cdot.x, cdot.y);
  }

  __device__ double3 cdotNormA_(const double2 &a, const double2 &b) 
  { return make_double3(a.x*b.x + a.y*b.y, a.x*b.y - a.y*b.x, a.x*a.x + a.y*a.y); }
  __device__ double3 cdotNormA_(const float2 &a, const float2 &b) 
  { return make_double3(a.x*b.x + a.y*b.y, a.x*b.y - a.y*b.x, a.x*a.x + a.y*a.y); }
  __device__ double3 cdotNormA_(const float4 &a, const float4 &b) 
  { return make_double3(a.x*b.x + a.y*b.y + a.z*b.z + a.w*b.w, 
                        a.x*b.y - a.y*b.x + a.z*b.w - a.w*b.z,
                        a.x*a.x + a.y*a.y + a.z*a.z + a.w*a.w); }

  template <typename ReduceType, typename Float2, typename FloatN>
#if (__COMPUTE_CAPABILITY__ >= 200)
  struct CdotNormA : public ReduceFunctor<ReduceType, Float2, FloatN> {
#else
  struct CdotNormA {
#endif
    CdotNormA(const Float2 &a, const Float2 &b) { ; }
    __device__ void operator()(ReduceType &sum, FloatN &x, FloatN &y, FloatN &z, FloatN &w, FloatN &v) { sum += cdotNormA_(x,y); }
    static int streams() { return 2; } 
    static int flops() { return 6; } 
  };

  double3 cDotProductNormACuda(cudaColorSpinorField &x, cudaColorSpinorField &y) {
    return reduce::reduceCuda<double3,QudaSumFloat3,QudaSumFloat,CdotNormA,0,0,0,0,0,false>
      (make_double2(0.0, 0.0), make_double2(0.0, 0.0), x, y, x, x, x);
  }

  __device__ double3 cdotNormB_(const double2 &a, const double2 &b) 
  { return make_double3(a.x*b.x + a.y*b.y, a.x*b.y - a.y*b.x, b.x*b.x + b.y*b.y); }
  __device__ double3 cdotNormB_(const float2 &a, const float2 &b) 
  { return make_double3(a.x*b.x + a.y*b.y, a.x*b.y - a.y*b.x, b.x*b.x + b.y*b.y); }
  __device__ double3 cdotNormB_(const float4 &a, const float4 &b) 
  { return make_double3(a.x*b.x + a.y*b.y + a.z*b.z + a.w*b.w, a.x*b.y - a.y*b.x + a.z*b.w - a.w*b.z,
                        b.x*b.x + b.y*b.y + b.z*b.z + b.w*b.w); }

  template <typename ReduceType, typename Float2, typename FloatN>
#if (__COMPUTE_CAPABILITY__ >= 200)
  struct CdotNormB : public ReduceFunctor<ReduceType, Float2, FloatN> {
#else
  struct CdotNormB {
#endif
    CdotNormB(const Float2 &a, const Float2 &b) { ; }
    __device__ void operator()(ReduceType &sum, FloatN &x, FloatN &y, FloatN &z, FloatN &w, FloatN &v) { sum += cdotNormB_(x,y); }
    static int streams() { return 2; } 
    static int flops() { return 6; } 
  };

  double3 cDotProductNormBCuda(cudaColorSpinorField &x, cudaColorSpinorField &y) {
    return reduce::reduceCuda<double3,QudaSumFloat3,QudaSumFloat,CdotNormB,0,0,0,0,0,false>
      (make_double2(0.0, 0.0), make_double2(0.0, 0.0), x, y, x, x, x);
  }

  template <typename ReduceType, typename Float2, typename FloatN>
#if (__COMPUTE_CAPABILITY__ >= 200)
  struct caxpbypzYmbwcDotProductUYNormY : public ReduceFunctor<ReduceType, Float2, FloatN> {
#else
  struct caxpbypzYmbwcDotProductUYNormY {
#endif
    Float2 a;
    Float2 b;
    caxpbypzYmbwcDotProductUYNormY(const Float2 &a, const Float2 &b) : a(a), b(b) { ; }
    __device__ void operator()(ReduceType &sum, FloatN &x, FloatN &y, FloatN &z, FloatN &w, FloatN &v) { Caxpy_(a, x, z); Caxpy_(b, y, z); Caxpy_(-b, w, y); sum += cdotNormB_(v,y); }
    static int streams() { return 7; } 
    static int flops() { return 18; } 
  };

  double3 caxpbypzYmbwcDotProductUYNormYCuda(const Complex &a, cudaColorSpinorField &x, 
                                             const Complex &b, cudaColorSpinorField &y,
                                             cudaColorSpinorField &z, cudaColorSpinorField &w,
                                             cudaColorSpinorField &u) {
    if (x.Precision() != z.Precision()) {
      return reduce::mixed::reduceCuda<double3,QudaSumFloat3,QudaSumFloat,caxpbypzYmbwcDotProductUYNormY,0,1,1,0,0,false>
      (make_double2(REAL(a), IMAG(a)), make_double2(REAL(b), IMAG(b)), x, y, z, w, u);

    } else {
      return reduce::reduceCuda<double3,QudaSumFloat3,QudaSumFloat,caxpbypzYmbwcDotProductUYNormY,0,1,1,0,0,false>
      (make_double2(REAL(a), IMAG(a)), make_double2(REAL(b), IMAG(b)), x, y, z, w, u);
    }
  }


  template <typename ReduceType, typename Float2, typename FloatN>
#if (__COMPUTE_CAPABILITY__ >= 200)
  struct axpyCGNorm2 : public ReduceFunctor<ReduceType, Float2, FloatN> {
#else
  struct axpyCGNorm2 {
#endif
    Float2 a;
    axpyCGNorm2(const Float2 &a, const Float2 &b) : a(a) { ; }
    __device__ void operator()(ReduceType &sum, FloatN &x, FloatN &y, FloatN &z, FloatN &w, FloatN &v) { 
      FloatN y_new = y + a.x*x;
      sum.x += norm2_(y_new); 
      sum.y += dot_(y_new, y_new-y);
      y = y_new;
    }
    static int streams() { return 3; } 
    static int flops() { return 6; } 
  };

  Complex axpyCGNormCuda(const double &a, cudaColorSpinorField &x, cudaColorSpinorField &y) {
    double2 cg_norm = reduce::reduceCuda<double2,QudaSumFloat2,QudaSumFloat,axpyCGNorm2,0,1,0,0,0,false>
      (make_double2(a, 0.0), make_double2(0.0, 0.0), x, y, x, x, x);
    return Complex(cg_norm.x, cg_norm.y);
  }

#if (__COMPUTE_CAPABILITY__ >= 200)

  template <typename ReduceType, typename Float2, typename FloatN>
  struct HeavyQuarkResidualNorm : public ReduceFunctor<ReduceType, Float2, FloatN> {
    Float2 a;
    Float2 b;
    ReduceType aux;
    HeavyQuarkResidualNorm(const Float2 &a, const Float2 &b) : a(a), b(b) { ; }
    
    __device__ void pre() { aux.x = 0; aux.y = 0; }
    
    __device__ void operator()(ReduceType &sum, FloatN &x, FloatN &y, FloatN &z, FloatN &w, FloatN &v) { aux.x += norm2_(x); aux.y += norm2_(y); }
    
    __device__ void post(ReduceType &sum) 
    { 
      sum.x += aux.x; sum.y += aux.y; sum.z += (aux.x > 0.0) ? (aux.y / aux.x) : 1.0; 
    }
    
    static int streams() { return 2; } 
    static int flops() { return 4; } 
  };
  
  double3 HeavyQuarkResidualNormCuda(cudaColorSpinorField &x, cudaColorSpinorField &r) {
    double3 rtn = reduce::reduceCuda<double3,QudaSumFloat3,QudaSumFloat,HeavyQuarkResidualNorm,0,0,0,0,0,true>
      (make_double2(0.0, 0.0), make_double2(0.0, 0.0), x, r, r, r, r);
#ifdef MULTI_GPU
    rtn.z /= (x.Volume()*comm_size());
#else
    rtn.z /= x.Volume();
#endif
    return rtn;
  }
  
  template <typename ReduceType, typename Float2, typename FloatN>
  struct xpyHeavyQuarkResidualNorm : public ReduceFunctor<ReduceType, Float2, FloatN> {
    Float2 a;
    Float2 b;
    ReduceType aux;
    xpyHeavyQuarkResidualNorm(const Float2 &a, const Float2 &b) : a(a), b(b) { ; }
    
    __device__ void pre() { aux.x = 0; aux.y = 0; }
    
    __device__ void operator()(ReduceType &sum, FloatN &x, FloatN &y, FloatN &z, FloatN &w, FloatN &v) 
    { aux.x += norm2_(x + y); aux.y += norm2_(z); }
    
    __device__ void post(ReduceType &sum) 
    { 
      sum.x += aux.x; sum.y += aux.y; sum.z += (aux.x > 0.0) ? (aux.y / aux.x) : 1.0; 
    }
    
    static int streams() { return 3; } 
    static int flops() { return 5; }
  };
  
  double3 xpyHeavyQuarkResidualNormCuda(cudaColorSpinorField &x, cudaColorSpinorField &y,
                                        cudaColorSpinorField &r) {
    double3 rtn = reduce::reduceCuda<double3,QudaSumFloat3,QudaSumFloat,xpyHeavyQuarkResidualNorm,0,0,0,0,0,true>
      (make_double2(0.0, 0.0), make_double2(0.0, 0.0), x, y, r, r, r);
#ifdef MULTI_GPU
    rtn.z /= (x.Volume()*comm_size());
#else
    rtn.z /= x.Volume();
#endif
    return rtn;
  }

#else
  
  double3 HeavyQuarkResidualNormCuda(cudaColorSpinorField &x, cudaColorSpinorField &r) {
    errorQuda("Not supported on pre-Fermi architectures");
    return make_double3(0.0,0.0,0.0);
  }

  double3 xpyHeavyQuarkResidualNormCuda(cudaColorSpinorField &x, cudaColorSpinorField &y,
                                        cudaColorSpinorField &r) {
    errorQuda("Not supported on pre-Fermi architectures");
    return make_double3(0.0,0.0,0.0);
  }

#endif
    
  template <typename ReduceType, typename Float2, typename FloatN>
#if (__COMPUTE_CAPABILITY__ >= 200)
  struct tripleCGReduction : public ReduceFunctor<ReduceType, Float2, FloatN> {
#else
  struct tripleCGReduction {
#endif
    tripleCGReduction(const Float2 &a, const Float2 &b) { ; }
    __device__ void operator()(ReduceType &sum, FloatN &x, FloatN &y, FloatN &z, FloatN &w, FloatN &v) 
    { sum.x += norm2_(x); sum.y += norm2_(y); sum.z += dot_(y,z); }
    static int streams() { return 3; } 
    static int flops() { return 6; } 
  };

  double3 tripleCGReductionCuda(cudaColorSpinorField &x, cudaColorSpinorField &y, cudaColorSpinorField &z) {
    return reduce::reduceCuda<double3,QudaSumFloat3,QudaSumFloat,tripleCGReduction,0,0,0,0,0,false>
      (make_double2(0.0, 0.0), make_double2(0.0, 0.0), x, y, z, x, x);
  }

} // namespace quda