quda/lib/blas_quda.cu

Go to the documentation of this file.

#include <stdlib.h>
#include <stdio.h>

#include <quda_internal.h>
#include <blas_quda.h>
#include <color_spinor_field.h>

#include <cuComplex.h>

#define REDUCE_MAX_BLOCKS 65536

#define REDUCE_DOUBLE 64
#define REDUCE_KAHAN 32

#if (__CUDA_ARCH__ >= 130)
#define REDUCE_TYPE REDUCE_DOUBLE
#define QudaSumFloat double
#define QudaSumComplex cuDoubleComplex
#define QudaSumFloat3 double3
#else
#define REDUCE_TYPE REDUCE_KAHAN
#define QudaSumFloat float
#define QudaSumComplex cuComplex
#define QudaSumFloat3 float3
#endif

// Required for the reduction kernels
#ifdef __DEVICE_EMULATION__
#define EMUSYNC __syncthreads()
#else
#define EMUSYNC
#endif

// These are used for reduction kernels
static QudaSumFloat *d_reduceFloat=0;
static QudaSumComplex *d_reduceComplex=0;
static QudaSumFloat3 *d_reduceFloat3=0;

static QudaSumFloat *h_reduceFloat=0;
static QudaSumComplex *h_reduceComplex=0;
static QudaSumFloat3 *h_reduceFloat3=0;

unsigned long long blas_quda_flops;
unsigned long long blas_quda_bytes;

static dim3 blasBlock;
static dim3 blasGrid;

// generated by blas_test
#include <blas_param.h>

double2 operator+(const double2& x, const double2 &y) {
  return make_double2(x.x + y.x, x.y + y.y);
}

double3 operator+(const double3& x, const double3 &y) {
  double3 z;
  z.x = x.x + y.x; z.y = x.y + y.y; z.z = x.z + y.z;
  return z;
}

__device__ float2 operator*(const float a, const float2 x) {
  float2 y;
  y.x = a*x.x;
  y.y = a*x.y;
  return y;
}

template <typename Float2>
__device__ Float2 operator+(const Float2 x, const Float2 y) {
  Float2 z;
  z.x = x.x + y.x;
  z.y = x.y + y.y;
  return z;
}

template <typename Float2>
__device__ Float2 operator+=(Float2 &x, const Float2 y) {
  x.x += y.x;
  x.y += y.y;
  return x;
}

template <typename Float2>
__device__ Float2 operator-=(Float2 &x, const Float2 y) {
  x.x -= y.x;
  x.y -= y.y;
  return x;
}

template <typename Float, typename Float2>
__device__ Float2 operator*=(Float2 &x, const Float a) {
  x.x *= a;
  x.y *= a;
  return x;
}

template <typename Float>
__device__ float4 operator*=(float4 &a, const Float &b) {
  a.x *= b;
  a.y *= b;
  a.z *= b;
  a.w *= b;
  return a;
}



void zeroCuda(cudaColorSpinorField &a) { a.zero(); }

void initBlas(void)
{  
  if (!d_reduceFloat) {
    if (cudaMalloc((void**) &d_reduceFloat, REDUCE_MAX_BLOCKS*sizeof(QudaSumFloat)) == cudaErrorMemoryAllocation) {
      errorQuda("Error allocating device reduction array");
    }
  }

  if (!d_reduceComplex) {
    if (cudaMalloc((void**) &d_reduceComplex, REDUCE_MAX_BLOCKS*sizeof(QudaSumComplex)) == cudaErrorMemoryAllocation) {
      errorQuda("Error allocating device reduction array");
    }
  }
  
  if (!d_reduceFloat3) {
    if (cudaMalloc((void**) &d_reduceFloat3, REDUCE_MAX_BLOCKS*sizeof(QudaSumFloat3)) == cudaErrorMemoryAllocation) {
      errorQuda("Error allocating device reduction array");
    }
  }

  if (!h_reduceFloat) {
    if (cudaMallocHost((void**) &h_reduceFloat, REDUCE_MAX_BLOCKS*sizeof(QudaSumFloat)) == cudaErrorMemoryAllocation) {
      errorQuda("Error allocating host reduction array");
    }
  }

  if (!h_reduceComplex) {
    if (cudaMallocHost((void**) &h_reduceComplex, REDUCE_MAX_BLOCKS*sizeof(QudaSumComplex)) == cudaErrorMemoryAllocation) {
      errorQuda("Error allocating host reduction array");
    }
  }
  
  if (!h_reduceFloat3) {
    if (cudaMallocHost((void**) &h_reduceFloat3, REDUCE_MAX_BLOCKS*sizeof(QudaSumFloat3)) == cudaErrorMemoryAllocation) {
      errorQuda("Error allocating host reduction array");
    }
  }
}

void endBlas(void)
{
  if (d_reduceFloat) cudaFree(d_reduceFloat);
  if (d_reduceComplex) cudaFree(d_reduceComplex);
  if (d_reduceFloat3) cudaFree(d_reduceFloat3);
  if (h_reduceFloat) cudaFreeHost(h_reduceFloat);
  if (h_reduceComplex) cudaFreeHost(h_reduceComplex);
  if (h_reduceFloat3) cudaFreeHost(h_reduceFloat3);
}

// blasTuning = 1 turns off error checking
static int blasTuning = 0;

void setBlasTuning(int tuning)
{
  blasTuning = tuning;
}

void setBlasParam(int kernel, int prec, int threads, int blocks)
{
  blas_threads[kernel][prec] = threads;
  blas_blocks[kernel][prec] = blocks;
}

void setBlock(int kernel, int length, QudaPrecision precision)
{
  int prec;
  switch(precision) {
  case QUDA_HALF_PRECISION:
    prec = 0;
    break;
  case QUDA_SINGLE_PRECISION:
    prec = 1;
    break;
  case QUDA_DOUBLE_PRECISION:
    prec = 2;
    break;
  }

  int blocks = min(blas_blocks[kernel][prec], max(length/blas_threads[kernel][prec], 1));
  blasBlock.x = blas_threads[kernel][prec];
  blasBlock.y = 1;
  blasBlock.z = 1;

  blasGrid.x = blocks;
  blasGrid.y = 1;
  blasGrid.z = 1;
}

#if (__CUDA_ARCH__ >= 130)
static __inline__ __device__ double2 fetch_double2(texture<int4, 1> t, int i)
{
  int4 v = tex1Dfetch(t,i);
  return make_double2(__hiloint2double(v.y, v.x), __hiloint2double(v.w, v.z));
}
#else
static __inline__ __device__ double2 fetch_double2(texture<int4, 1> t, int i)
{
  // do nothing
  return make_double2(0.0, 0.0);
}
#endif

float2 __device__ read_Float2(float2 *x, int i) {
  return make_float2(x[i].x, x[i].y);
}

double2 __device__ read_Float2(double2 *x, int i) {
  return make_double2(x[i].x, x[i].y);
}

#define READ_DOUBLE2_TEXTURE(x, i) \
  fetch_double2(x##TexDouble2, i)

#define READ_FLOAT2_TEXTURE(x, i) \
  tex1Dfetch(x##TexSingle2, i)

float2 __device__ make_Float2(float2 x) {
  return make_float2(x.x, x.y);
}

double2 __device__ make_Float2(double2 x) {
  return make_double2(x.x, x.y);
}

#define RECONSTRUCT_HALF_SPINOR(a, texHalf, texNorm, length)            \
  float a##c = tex1Dfetch(texNorm, i);                                  \
  float4 a##0 = tex1Dfetch(texHalf, i + 0*length);                      \
  float4 a##1 = tex1Dfetch(texHalf, i + 1*length);                      \
  float4 a##2 = tex1Dfetch(texHalf, i + 2*length);                      \
  float4 a##3 = tex1Dfetch(texHalf, i + 3*length);                      \
  float4 a##4 = tex1Dfetch(texHalf, i + 4*length);                      \
  float4 a##5 = tex1Dfetch(texHalf, i + 5*length);                      \
  a##0 *= a##c;                                                         \
  a##1 *= a##c;                                                         \
  a##2 *= a##c;                                                         \
  a##3 *= a##c;                                                         \
  a##4 *= a##c;                                                         \
  a##5 *= a##c;

#define RECONSTRUCT_HALF_SPINOR_ST(a, texHalf, texNorm, length)         \
  float a##c = tex1Dfetch(texNorm, i);                                  \
  float2 a##0 = tex1Dfetch(texHalf, i + 0*length);                      \
  float2 a##1 = tex1Dfetch(texHalf, i + 1*length);                      \
  float2 a##2 = tex1Dfetch(texHalf, i + 2*length);                      \
  (a##0) *= a##c;                                                               \
  (a##1) *= a##c;                                                               \
  (a##2) *= a##c;


// Some musings on how to clean up the blas code using Boost
/*#define BOOST_RECONSTRUCT_HALF_SPINOR(z, j, a, texHalf, length)       \
  float4 a##k tex1Dfetch(texHalf, i + j*length);        \
  a##k *= a##c;

#define RECONSTRUCT_HALF_SPINOR(a, texHalf, texNorm, length)            \
  BOOST_PP_REPEAT(6, BOOST_RECONSTRUCT_HALF_SPINOR, a, texHalf, length) \
*/

#define READ_HALF_SPINOR_TEX(a, tex, texNorm, length)                   \
  float a##c = tex1Dfetch(texNorm, i);                                  \
  float4 a##0 = tex1Dfetch(tex, i + 0*length);                          \
  float4 a##1 = tex1Dfetch(tex, i + 1*length);                          \
  float4 a##2 = tex1Dfetch(tex, i + 2*length);                          \
  float4 a##3 = tex1Dfetch(tex, i + 3*length);                          \
  float4 a##4 = tex1Dfetch(tex, i + 4*length);                          \
  float4 a##5 = tex1Dfetch(tex, i + 5*length);                          \

#define READ_HALF_SPINOR(a, tex, length)                                \
  float4 a##0 = tex1Dfetch(tex, i + 0*length);                          \
  float4 a##1 = tex1Dfetch(tex, i + 1*length);                          \
  float4 a##2 = tex1Dfetch(tex, i + 2*length);                          \
  float4 a##3 = tex1Dfetch(tex, i + 3*length);                          \
  float4 a##4 = tex1Dfetch(tex, i + 4*length);                          \
  float4 a##5 = tex1Dfetch(tex, i + 5*length);                          \
  float a##c = a##N[i];

#define READ_HALF_SPINOR_ST(a, tex, length)                             \
  float2 a##0 = tex1Dfetch(tex, i + 0*length);                          \
  float2 a##1 = tex1Dfetch(tex, i + 1*length);                          \
  float2 a##2 = tex1Dfetch(tex, i + 2*length);                          \
  float a##c = a##N[i];

#define SHORT_LENGTH 65536
#define SCALE_FLOAT ((SHORT_LENGTH-1) * 0.5)
#define SHIFT_FLOAT (-1.f / (SHORT_LENGTH-1))

__device__ short float2short(float c, float a) {
  //return (short)(a*MAX_SHORT);
  short rtn = (short)((a+SHIFT_FLOAT)*SCALE_FLOAT*c);
  return rtn;
}

__device__ float short2float(short a) {
  return (float)a/SCALE_FLOAT - SHIFT_FLOAT;
}

__device__ short4 float42short4(float c, float4 a) {
  return make_short4(float2short(c, a.x), float2short(c, a.y), float2short(c, a.z), float2short(c, a.w));
}

#define FAST_ABS_MAX(a, b) fmaxf(fabsf(a), fabsf(b));
#define FAST_MAX(a, b) fmaxf(a, b);

__device__ float fast_abs_max(float4 a) {
  float c0 = FAST_ABS_MAX(a.x, a.y);
  float c1 = FAST_ABS_MAX(a.z, a.w);
  return FAST_MAX(c0, c1);
}

#define CONSTRUCT_HALF_SPINOR_FROM_SINGLE(h, n, a, length) {            \
    float c0 = fast_abs_max(a##0);                                      \
    float c1 = fast_abs_max(a##1);                                      \
    c0 = FAST_MAX(c0, c1);                                              \
    float c2 = fast_abs_max(a##2);                                      \
    float c3 = fast_abs_max(a##3);                                      \
    c1 = FAST_MAX(c2, c3);                                              \
    c0 = FAST_MAX(c0, c1);                                              \
    c2 = fast_abs_max(a##4);                                            \
    c3 = fast_abs_max(a##5);                                            \
    c1 = FAST_MAX(c2, c3);                                              \
    c0 = FAST_MAX(c0, c1);                                              \
    n[i] = c0;                                                          \
    float C = __fdividef(MAX_SHORT, c0);                                \
    h[i+0*length] = make_short4((short)(C*(float)(a##0).x), (short)(C*(float)(a##0).y), \
                                (short)(C*(float)(a##0).z), (short)(C*(float)(a##0).w)); \
    h[i+1*length] = make_short4((short)(C*(float)(a##1).x), (short)(C*(float)(a##1).y), \
                                (short)(C*(float)(a##1).z), (short)(C*(float)(a##1).w)); \
    h[i+2*length] = make_short4((short)(C*(float)(a##2).x), (short)(C*(float)(a##2).y), \
                                (short)(C*(float)(a##2).z), (short)(C*(float)(a##2).w)); \
    h[i+3*length] = make_short4((short)(C*(float)(a##3).x), (short)(C*(float)(a##3).y), \
                                (short)(C*(float)(a##3).z), (short)(C*(float)(a##3).w)); \
    h[i+4*length] = make_short4((short)(C*(float)(a##4).x), (short)(C*(float)(a##4).y), \
                                (short)(C*(float)(a##4).z), (short)(C*(float)(a##4).w)); \
    h[i+5*length] = make_short4((short)(C*(float)(a##5).x), (short)(C*(float)(a##5).y), \
                                (short)(C*(float)(a##5).z), (short)(C*(float)(a##5).w));}

#define CONSTRUCT_HALF_SPINOR_FROM_DOUBLE(h, n, a, length)              \
  {float c0 = fmaxf(fabsf((a##0).x), fabsf((a##0).y));                  \
  float c1 = fmaxf(fabsf((a##1).x), fabsf((a##1).y));                   \
  float c2 = fmaxf(fabsf((a##2).x), fabsf((a##2).y));                   \
  float c3 = fmaxf(fabsf((a##3).x), fabsf((a##3).y));                   \
  float c4 = fmaxf(fabsf((a##4).x), fabsf((a##4).y));                   \
  float c5 = fmaxf(fabsf((a##5).x), fabsf((a##5).y));                   \
  float c6 = fmaxf(fabsf((a##6).x), fabsf((a##6).y));                   \
  float c7 = fmaxf(fabsf((a##7).x), fabsf((a##7).y));                   \
  float c8 = fmaxf(fabsf((a##8).x), fabsf((a##8).y));                   \
  float c9 = fmaxf(fabsf((a##9).x), fabsf((a##9).y));                   \
  float c10 = fmaxf(fabsf((a##10).x), fabsf((a##10).y));                \
  float c11 = fmaxf(fabsf((a##11).x), fabsf((a##11).y));                \
  c0 = fmaxf(c0, c1); c1 = fmaxf(c2, c3);  c2 = fmaxf(c4, c5); c3 = fmaxf(c6, c7); \
  c4 = fmaxf(c8, c9); c5 = fmaxf(c10, c11); c0 = fmaxf(c0, c1); c1 = fmaxf(c2, c3); \
  c2 = fmaxf(c4, c5); c0 = fmaxf(c0, c1); c0 = fmaxf(c0, c2);           \
  n[i] = c0;                                                            \
  float C = __fdividef(MAX_SHORT, c0);                                  \
  h[i+0*length] = make_short4((short)(C*(float)(a##0).x), (short)(C*(float)(a##0).y), \
                              (short)(C*(float)(a##1).x), (short)(C*(float)(a##1).y)); \
  h[i+1*length] = make_short4((short)(C*(float)(a##2).x), (short)(C*(float)(a##2).y), \
                              (short)(C*(float)(a##3).x), (short)(C*(float)(a##3).y)); \
  h[i+2*length] = make_short4((short)(C*(float)(a##4).x), (short)(C*(float)(a##4).y), \
                              (short)(C*(float)(a##5).x), (short)(C*(float)(a##5).y)); \
  h[i+3*length] = make_short4((short)(C*(float)(a##6).x), (short)(C*(float)(a##6).y), \
                              (short)(C*(float)(a##7).x), (short)(C*(float)(a##7).y)); \
  h[i+4*length] = make_short4((short)(C*(float)(a##8).x), (short)(C*(float)(a##8).y), \
                              (short)(C*(float)(a##9).x), (short)(C*(float)(a##9).y)); \
  h[i+5*length] = make_short4((short)(C*(float)(a##10).x), (short)(C*(float)(a##10).y), \
                              (short)(C*(float)(a##11).x), (short)(C*(float)(a##11).y));}

#define CONSTRUCT_HALF_SPINOR_FROM_SINGLE_ST(h, n, a, length)           \
    {float c0 = fmaxf(fabsf((a##0).x), fabsf((a##0).y));                \
        float c1 = fmaxf(fabsf((a##1).x), fabsf((a##1).y));             \
        float c2 = fmaxf(fabsf((a##2).x), fabsf((a##2).y));             \
        c0 = fmaxf(c0, c1); c0 = fmaxf(c0, c2);                         \
        n[i] = c0;                                                      \
        float C = __fdividef(MAX_SHORT, c0);                            \
        h[i+0*length] = make_short2((short)(C*(float)(a##0).x), (short)(C*(float)(a##0).y)); \
        h[i+1*length] = make_short2((short)(C*(float)(a##1).x), (short)(C*(float)(a##1).y)); \
        h[i+2*length] = make_short2((short)(C*(float)(a##2).x), (short)(C*(float)(a##2).y));}

#define CONSTRUCT_HALF_SPINOR_FROM_DOUBLE_ST(h, n, a, length)           \
    {float c0 = fmaxf(fabsf((a##0).x), fabsf((a##0).y));                \
        float c1 = fmaxf(fabsf((a##1).x), fabsf((a##1).y));             \
        float c2 = fmaxf(fabsf((a##2).x), fabsf((a##2).y));             \
        c0 = fmaxf(c0, c1); c0 = fmaxf(c0, c2);                         \
        n[i] = c0;                                                      \
        float C = __fdividef(MAX_SHORT, c0);                            \
        h[i+0*length] = make_short2((short)(C*(float)(a##0).x), (short)(C*(float)(a##0).y)); \
        h[i+1*length] = make_short2((short)(C*(float)(a##1).x), (short)(C*(float)(a##1).y)); \
        h[i+2*length] = make_short2((short)(C*(float)(a##2).x), (short)(C*(float)(a##2).y));}


#define SUM_FLOAT4(sum, a)                      \
  float sum = a.x + a.y + a.z + a.w;

#define SUM_FLOAT2(sum, a)                      \
  float sum = a.x + a.y;

#if (__CUDA_ARCH__ < 200) 
#define REAL_DOT_FLOAT4(dot, a, b) \
  float dot = a.x*b.x + a.y*b.y + a.z*b.z + a.w*b.w;
#else
#define REAL_DOT_FLOAT4(dot, a, b)              \
  float dot = fmaf(a.x, b.x, 0.0f);             \
  dot = fmaf(a.y, b.y, dot);                    \
  dot = fmaf(a.z, b.z, dot);                    \
  dot = fmaf(a.w, b.w, dot)
#endif

#define REAL_DOT_FLOAT2(dot, a, b) \
  float dot = a.x*b.x + a.y*b.y;

#if (__CUDA_ARCH__ < 200) 
#define IMAG_DOT_FLOAT4(dot, a, b)                      \
  float dot = a.x*b.y - a.y*b.x + a.z*b.w - a.w*b.z;
#else
#define IMAG_DOT_FLOAT4(dot, a, b)              \
  float dot = fmaf(a.x, b.y, 0.0f);             \
  dot = fmaf(-a.y, b.x, dot);                   \
  dot = fmaf(a.z, b.w, dot);                    \
  dot = fmaf(-a.w, b.z, dot)
#endif

#define IMAG_DOT_FLOAT2(dot, a, b)                      \
  float dot = a.x*b.y - a.y*b.x;

#define AX_FLOAT4(a, X)                         \
  X.x *= a; X.y *= a; X.z *= a; X.w *= a;

#define AX_FLOAT2(a, X)                         \
  X.x *= a; X.y *= a;

#define XPY_FLOAT4(X, Y)                     \
  Y.x += X.x; Y.y += X.y; Y.z += X.z; Y.w += X.w;

#define XPY_FLOAT2(X, Y)                     \
  Y.x += X.x; Y.y += X.y;

#define XMY_FLOAT4(X, Y)                     \
  Y.x = X.x - Y.x; Y.y = X.y - Y.y; Y.z = X.z - Y.z; Y.w = X.w - Y.w;

#define XMY_FLOAT2(X, Y)                     \
  Y.x = X.x - Y.x; Y.y = X.y - Y.y;

#define MXPY_FLOAT4(X, Y)                    \
  Y.x -= X.x; Y.y -= X.y; Y.z -= X.z; Y.w -= X.w;

#define MXPY_FLOAT2(X, Y)                    \
  Y.x -= X.x; Y.y -= X.y; 

#if (__CUDA_ARCH__ < 200) 
#define AXPY_FLOAT4(a, X, Y)                 \
  Y.x += a*X.x; Y.y += a*X.y;                \
  Y.z += a*X.z; Y.w += a*X.w;
#else
#define AXPY_FLOAT4(a, X, Y)                            \
  Y.x = fmaf(a, X.x, Y.x); Y.y = fmaf(a, X.y, Y.y);     \
  Y.z = fmaf(a, X.z, Y.z); Y.w = fmaf(a, X.w, Y.w);
#endif

#define AXPY_FLOAT2(a, X, Y)                 \
  Y.x += a*X.x; Y.y += a*X.y;                

#define AXPBY_FLOAT4(a, X, b, Y)                                \
  Y.x = b*Y.x; Y.x += a*X.x; Y.y = b*Y.y; Y.y += a*X.y;         \
  Y.z = b*Y.z; Y.z += a*X.z; Y.w = b*Y.w; Y.w += a*X.w; 

#define AXPBY_FLOAT2(a, X, b, Y)                        \
  Y.x = b*Y.x; Y.x += a*X.x; Y.y = b*Y.y; Y.y += a*X.y;         \

#if (__CUDA_ARCH__ < 200)
#define XPAY_FLOAT4(X, a, Y)                         \
  Y.x = X.x + a*Y.x; Y.y = X.y + a*Y.y;              \
  Y.z = X.z + a*Y.z; Y.w = X.w + a*Y.w;
#else
#define XPAY_FLOAT4(X, a, Y)                         \
  Y.x = fmaf(a, Y.x, X.x); Y.y = fmaf(a, Y.y, X.y);  \
  Y.z = fmaf(a, Y.z, X.z); Y.w = fmaf(a, Y.w, X.w);
#endif

#define XPAY_FLOAT2(X, a, Y)                         \
  Y.x = X.x + a*Y.x; Y.y = X.y + a*Y.y;              

#if (__CUDA_ARCH__ < 200)
#define CAXPY_FLOAT4(a, X, 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;
#else
#define CAXPY_FLOAT4(a, X, Y)                                   \
  Y.x = fmaf(a.x, X.x, Y.x); Y.x = fmaf(-a.y, X.y, Y.x);        \
  Y.y = fmaf(a.y, X.x, Y.y); Y.y = fmaf( a.x, X.y, Y.y);        \
  Y.z = fmaf(a.x, X.z, Y.z); Y.z = fmaf(-a.y, X.w, Y.z);        \
  Y.w = fmaf(a.y, X.z, Y.w); Y.w = fmaf( a.x, X.w, Y.w);
#endif // (__CUDA_ARCH__ < 200)

#if (__CUDA_ARCH__ < 200)
#define CAXPY_FLOAT2(a, X, 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;
#else
#define CAXPY_FLOAT2(a, X, Y)                                   \
  Y.x = fmaf(a.x, X.x, Y.x); Y.x = fmaf(-a.y, X.y, Y.x);        \
  Y.y = fmaf(a.y, X.x, Y.y); Y.y = fmaf( a.x, X.y, Y.y);
#endif // (__CUDA_ARCH__ < 200)

#define CMAXPY_FLOAT4(a, X, 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;

#define CAXPBY_FLOAT4(a, X, b, Y)                                       \
  { float2 y;                                                           \
  y.x = a.x*X.x; y.x -= a.y*X.y; y.x += b.x*Y.x; y.x -= b.y*Y.y;        \
  y.y = a.y*X.x; y.y += a.x*X.y; y.y += b.y*Y.x; y.y += b.x*Y.y;        \
  Y.x = y.x; Y.y = y.y;                                                 \
  y.x = a.x*X.z; y.x -= a.y*X.w; y.x += b.x*Y.z; y.x -= b.y*Y.w;        \
  y.y = a.y*X.z; y.y += a.x*X.w; y.y += b.y*Y.z; y.y += b.x*Y.w;        \
  Y.z = y.x; Y.w = y.y;}

#define CAXPBY_FLOAT2(a, X, b, Y)                                       \
  { float2 y;                                                           \
  y.x = a.x*X.x; y.x -= a.y*X.y; y.x += b.x*Y.x; y.x -= b.y*Y.y;        \
  y.y = a.y*X.x; y.y += a.x*X.y; y.y += b.y*Y.x; y.y += b.x*Y.y;        \
  Y.x = y.x; Y.y = y.y;}


#define CXPAYPBZ_FLOAT4(X, a, Y, b, Z)                                 \
  {float2 z;                                                           \
  z.x = X.x + a.x*Y.x; z.x -= a.y*Y.y; z.x += b.x*Z.x; z.x -= b.y*Z.y; \
  z.y = X.y + a.y*Y.x; z.y += a.x*Y.y; z.y += b.y*Z.x; z.y += b.x*Z.y; \
  Z.x = z.x; Z.y = z.y;                                                \
  z.x = X.z + a.x*Y.z; z.x -= a.y*Y.w; z.x += b.x*Z.z; z.x -= b.y*Z.w; \
  z.y = X.w + a.y*Y.z; z.y += a.x*Y.w; z.y += b.y*Z.z; z.y += b.x*Z.w; \
  Z.z = z.x; Z.w = z.y;}

#define CXPAYPBZ_FLOAT2(X, a, Y, b, Z)                                 \
  {float2 z;                                                           \
  z.x = X.x + a.x*Y.x; z.x -= a.y*Y.y; z.x += b.x*Z.x; z.x -= b.y*Z.y; \
  z.y = X.y + a.y*Y.x; z.y += a.x*Y.y; z.y += b.y*Z.x; z.y += b.x*Z.y; \
  Z.x = z.x; Z.y = z.y;}

#if (__CUDA_ARCH__ < 200)
#define CAXPBYPZ_FLOAT4(a, X, b, Y, Z)            \
  Z.x += a.x*X.x - a.y*X.y + b.x*Y.x - b.y*Y.y;   \
  Z.y += a.y*X.x + a.x*X.y + b.y*Y.x + b.x*Y.y;   \
  Z.z += a.x*X.z - a.y*X.w + b.x*Y.z - b.y*Y.w;   \
  Z.w += a.y*X.z + a.x*X.w + b.y*Y.z + b.x*Y.w;
#else
#define CAXPBYPZ_FLOAT4(a, X, b, Y, Z)                          \
  Z.x = fmaf(a.x, X.x, Z.x); Z.x = fmaf(-a.y, X.y, Z.x); Z.x = fmaf(b.x, Y.x, Z.x); Z.x = fmaf(-b.y, Y.y, Z.x); \
  Z.y = fmaf(a.y, X.x, Z.y); Z.y = fmaf( a.x, X.y, Z.y); Z.y = fmaf(b.y, Y.x, Z.y); Z.y = fmaf( b.x, Y.y, Z.y); \
  Z.z = fmaf(a.x, X.z, Z.z); Z.z = fmaf(-a.y, X.w, Z.z); Z.z = fmaf(b.x, Y.z, Z.z); Z.z = fmaf(-b.y, Y.w, Z.z); \
  Z.w = fmaf(a.y, X.z, Z.w); Z.w = fmaf( a.x, X.w, Z.w); Z.w = fmaf(b.y, Y.z, Z.w); Z.w = fmaf( b.x, Y.w, Z.w);
#endif // (__CUDA_ARCH__ < 200)

#if (__CUDA_ARCH__ < 200)
#define CAXPBYPZ_FLOAT2(a, X, b, Y, Z)            \
  Z.x += a.x*X.x - a.y*X.y + b.x*Y.x - b.y*Y.y;   \
  Z.y += a.y*X.x + a.x*X.y + b.y*Y.x + b.x*Y.y;
#else
#define CAXPBYPZ_FLOAT2(a, X, b, Y, Z)                          \
  Z.x = fmaf(a.x, X.x, Z.x); Z.x = fmaf(-a.y, X.y, Z.x); Z.x = fmaf(b.x, Y.x, Z.x); Z.x = fmaf(-b.y, Y.y, Z.x); \
  Z.y = fmaf(a.y, X.x, Z.y); Z.y = fmaf( a.x, X.y, Z.y); Z.y = fmaf(b.y, Y.x, Z.y); Z.y = fmaf( b.x, Y.y, Z.y);
#endif // (__CUDA_ARCH__ < 200)

// Double precision input spinor field
texture<int4, 1> xTexDouble2;
texture<int4, 1> yTexDouble2;
texture<int4, 1> zTexDouble2;
texture<int4, 1> wTexDouble2;
texture<int4, 1> uTexDouble2;

// Single precision input spinor field
texture<float2, 1> xTexSingle2;
texture<float2, 1> yTexSingle2;

texture<float4, 1> xTexSingle4;

// Half precision input spinor field
texture<short4, 1, cudaReadModeNormalizedFloat> texHalf1;
texture<short2, 1, cudaReadModeNormalizedFloat> texHalfSt1;
texture<float, 1, cudaReadModeElementType> texNorm1;

// Half precision input spinor field
texture<short4, 1, cudaReadModeNormalizedFloat> texHalf2;
texture<short2, 1, cudaReadModeNormalizedFloat> texHalfSt2;
texture<float, 1, cudaReadModeElementType> texNorm2;

// Half precision input spinor field
texture<short4, 1, cudaReadModeNormalizedFloat> texHalf3;
texture<short2, 1, cudaReadModeNormalizedFloat> texHalfSt3;
texture<float, 1, cudaReadModeElementType> texNorm3;

// Half precision input spinor field
texture<short4, 1, cudaReadModeNormalizedFloat> texHalf4;
texture<short2, 1, cudaReadModeNormalizedFloat> texHalfSt4;
texture<float, 1, cudaReadModeElementType> texNorm4;

// Half precision input spinor field
texture<short4, 1, cudaReadModeNormalizedFloat> texHalf5;
texture<short2, 1, cudaReadModeNormalizedFloat> texHalfSt5;
texture<float, 1, cudaReadModeElementType> texNorm5;

#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()); \
  }

// For kernels with precision conversion built in
#define checkSpinorLength(a, b)                                         \
  {                                                                     \
    if (a.Length() != b.Length()) {                                     \
      errorQuda("engths do not match: %d %d", a.Length(), b.Length());  \
    }                                                                   

__global__ void convertDSKernel(double2 *dst, float4 *src, int length) {
  unsigned int i = blockIdx.x*(blockDim.x) + threadIdx.x;
  unsigned int gridSize = gridDim.x*blockDim.x;
  while (i < length) {
    for (int k=0; k<6; k++) {
      dst[2*k*length+i].x = src[k*length+i].x;
      dst[2*k*length+i].y = src[k*length+i].y;
      dst[(2*k+1)*length+i].x = src[k*length+i].z;
      dst[(2*k+1)*length+i].y = src[k*length+i].w;
    }
    i += gridSize;
  }   
}

__global__ void convertDSKernel(double2 *dst, float2 *src, int length) {
  unsigned int i = blockIdx.x*(blockDim.x) + threadIdx.x;
  unsigned int gridSize = gridDim.x*blockDim.x;
  while (i < length) {
      for (int k=0; k<3; k++) {
          dst[k*length+i].x = src[k*length+i].x;
          dst[k*length+i].y = src[k*length+i].y;
      }
      i += gridSize;
  }   
}

__global__ void convertSDKernel(float4 *dst, double2 *src, int length) {
  unsigned int i = blockIdx.x*(blockDim.x) + threadIdx.x;
  unsigned int gridSize = gridDim.x*blockDim.x;
  while (i < length) {
    for (int k=0; k<6; k++) {
      dst[k*length+i].x = src[2*k*length+i].x;
      dst[k*length+i].y = src[2*k*length+i].y;
      dst[k*length+i].z = src[(2*k+1)*length+i].x;
      dst[k*length+i].w = src[(2*k+1)*length+i].y;
    }
    i += gridSize;
  }   
}

__global__ void convertSDKernel(float2 *dst, double2 *src, int length) {
    unsigned int i = blockIdx.x*(blockDim.x) + threadIdx.x;
  unsigned int gridSize = gridDim.x*blockDim.x;
  while (i < length) {
    for (int k=0; k<3; k++) {
        dst[k*length+i].x = src[k*length+i].x;
        dst[k*length+i].y = src[k*length+i].y;
    }
    i += gridSize;
  }   
}

__global__ void convertHSKernel(short4 *h, float *norm, int length, int real_length) {

  int i = blockIdx.x*(blockDim.x) + threadIdx.x;
  unsigned int gridSize = gridDim.x*blockDim.x;

  while(i < real_length) {
    float4 F0 = tex1Dfetch(xTexSingle4, i + 0*length);
    float4 F1 = tex1Dfetch(xTexSingle4, i + 1*length);
    float4 F2 = tex1Dfetch(xTexSingle4, i + 2*length);
    float4 F3 = tex1Dfetch(xTexSingle4, i + 3*length);
    float4 F4 = tex1Dfetch(xTexSingle4, i + 4*length);
    float4 F5 = tex1Dfetch(xTexSingle4, i + 5*length);
    CONSTRUCT_HALF_SPINOR_FROM_SINGLE(h, norm, F, length);
    i += gridSize;
  }

}

__global__ void convertHSKernel(short2 *h, float *norm, int length, int real_length) {

  int i = blockIdx.x*(blockDim.x) + threadIdx.x;
  unsigned int gridSize = gridDim.x*blockDim.x;

  while(i < real_length) {
      float2 F0 = tex1Dfetch(xTexSingle2, i + 0*length);
      float2 F1 = tex1Dfetch(xTexSingle2, i + 1*length);
      float2 F2 = tex1Dfetch(xTexSingle2, i + 2*length);
      CONSTRUCT_HALF_SPINOR_FROM_SINGLE_ST(h, norm, F, length);
      i += gridSize;
  }

}


__global__ void convertSHKernel(float4 *res, int length, int real_length) {

  int i = blockIdx.x*(blockDim.x) + threadIdx.x;
  unsigned int gridSize = gridDim.x*blockDim.x;

  while (i<real_length) {
    RECONSTRUCT_HALF_SPINOR(I, texHalf1, texNorm1, length);
    res[0*length+i] = I0;
    res[1*length+i] = I1;
    res[2*length+i] = I2;
    res[3*length+i] = I3;
    res[4*length+i] = I4;
    res[5*length+i] = I5;
    i += gridSize;
  }
}

__global__ void convertSHKernel(float2 *res, int length, int real_length) {

    int i = blockIdx.x*(blockDim.x) + threadIdx.x;
    unsigned int gridSize = gridDim.x*blockDim.x;
    
    while (i<real_length) {
        RECONSTRUCT_HALF_SPINOR_ST(I, texHalfSt1, texNorm1, length);
        res[0*length+i] = I0;
        res[1*length+i] = I1;
        res[2*length+i] = I2;
        i += gridSize;
    }
}

__global__ void convertHDKernel(short4 *h, float *norm, int length, int real_length) {

  int i = blockIdx.x*(blockDim.x) + threadIdx.x; 
  unsigned int gridSize = gridDim.x*blockDim.x;

  while(i < real_length) {
    double2 F0 = fetch_double2(xTexDouble2, i+0*length);
    double2 F1 = fetch_double2(xTexDouble2, i+1*length);
    double2 F2 = fetch_double2(xTexDouble2, i+2*length);
    double2 F3 = fetch_double2(xTexDouble2, i+3*length);
    double2 F4 = fetch_double2(xTexDouble2, i+4*length);
    double2 F5 = fetch_double2(xTexDouble2, i+5*length);
    double2 F6 = fetch_double2(xTexDouble2, i+6*length);
    double2 F7 = fetch_double2(xTexDouble2, i+7*length);
    double2 F8 = fetch_double2(xTexDouble2, i+8*length);
    double2 F9 = fetch_double2(xTexDouble2, i+9*length);
    double2 F10 = fetch_double2(xTexDouble2, i+10*length);
    double2 F11 = fetch_double2(xTexDouble2, i+11*length);
    CONSTRUCT_HALF_SPINOR_FROM_DOUBLE(h, norm, F, length);
    i += gridSize;
  }
}

__global__ void convertHDKernel(short2 *h, float *norm, int length, int real_length) {

  int i = blockIdx.x*(blockDim.x) + threadIdx.x; 
  unsigned int gridSize = gridDim.x*blockDim.x;

  while(i < real_length) {
      double2 F0 = fetch_double2(xTexDouble2, i+0*length);
      double2 F1 = fetch_double2(xTexDouble2, i+1*length);
      double2 F2 = fetch_double2(xTexDouble2, i+2*length);
      CONSTRUCT_HALF_SPINOR_FROM_DOUBLE_ST(h, norm, F, length);
      i += gridSize;
  }
}

__global__ void convertDHKernel(double2 *res, int length, int real_length) {

  int i = blockIdx.x*(blockDim.x) + threadIdx.x; 
  unsigned int gridSize = gridDim.x*blockDim.x;

  while(i < real_length) {
    RECONSTRUCT_HALF_SPINOR(I, texHalf1, texNorm1, length);
    res[0*length+i] = make_double2(I0.x, I0.y);
    res[1*length+i] = make_double2(I0.z, I0.w);
    res[2*length+i] = make_double2(I1.x, I1.y);
    res[3*length+i] = make_double2(I1.z, I1.w);
    res[4*length+i] = make_double2(I2.x, I2.y);
    res[5*length+i] = make_double2(I2.z, I2.w);
    res[6*length+i] = make_double2(I3.x, I3.y);
    res[7*length+i] = make_double2(I3.z, I3.w);
    res[8*length+i] = make_double2(I4.x, I4.y);
    res[9*length+i] = make_double2(I4.z, I4.w);
    res[10*length+i] = make_double2(I5.x, I5.y);
    res[11*length+i] = make_double2(I5.z, I5.w);
    i += gridSize;
  }

}

__global__ void convertDHKernelSt(double2 *res, int length, int real_length) {

  int i = blockIdx.x*(blockDim.x) + threadIdx.x; 
  unsigned int gridSize = gridDim.x*blockDim.x;
  
  while(i < real_length) {
      RECONSTRUCT_HALF_SPINOR_ST(I, texHalfSt1, texNorm1, length);
      res[0*length+i] = make_double2(I0.x, I0.y);
      res[1*length+i] = make_double2(I1.x, I1.y);
      res[2*length+i] = make_double2(I2.x, I2.y);
      i += gridSize;
  }

}



void copyCuda(cudaColorSpinorField &dst, const cudaColorSpinorField &src) {

  if (src.nSpin != 1 && src.nSpin != 4){
    errorQuda("nSpin(%d) not supported in function %s, line %d\n", src.nSpin, __FUNCTION__, __LINE__);  
  }
  if ((dst.precision == QUDA_HALF_PRECISION || src.precision == QUDA_HALF_PRECISION) &&
      (dst.siteSubset == QUDA_FULL_SITE_SUBSET || src.siteSubset == QUDA_FULL_SITE_SUBSET)) {
    copyCuda(dst.Even(), src.Even());
    copyCuda(dst.Odd(), src.Odd());
    return;
  }

  // For a given dst precision, there are two non-trivial possibilities for the
  // src precision.  The higher one corresponds to kernel index 0 (in the table
  // of block and grid dimensions), while the lower one corresponds to index 1.
  int id;
  if (src.Precision() == QUDA_DOUBLE_PRECISION ||
      dst.Precision() == QUDA_DOUBLE_PRECISION && src.Precision() == QUDA_SINGLE_PRECISION) {
    id = 0;
  } else {
    id = 1;
  }
  setBlock(id, dst.stride, dst.Precision());

  blas_quda_bytes += src.real_length*((int)src.precision + (int)dst.precision);

  if (dst.precision == QUDA_DOUBLE_PRECISION && src.precision == QUDA_SINGLE_PRECISION) {
      if (src.nSpin == 4){
          convertDSKernel<<<blasGrid, blasBlock>>>((double2*)dst.v, (float4*)src.v, src.stride);
      }else{ //src.nSpin == 1
          convertDSKernel<<<blasGrid, blasBlock>>>((double2*)dst.v, (float2*)src.v, src.stride);          
      }

  } else if (dst.precision == QUDA_SINGLE_PRECISION && src.precision == QUDA_DOUBLE_PRECISION) {
      if (src.nSpin == 4){
          convertSDKernel<<<blasGrid, blasBlock>>>((float4*)dst.v, (double2*)src.v, src.stride);
      }else{ //src.nSpin ==1
          convertSDKernel<<<blasGrid, blasBlock>>>((float2*)dst.v, (double2*)src.v, src.stride);
      }
  } else if (dst.precision == QUDA_SINGLE_PRECISION && src.precision == QUDA_HALF_PRECISION) {
      blas_quda_bytes += src.volume*sizeof(float);
      int spinor_bytes = src.length*sizeof(short);
      if (src.nSpin == 4){      
        cudaBindTexture(0, texHalf1, src.v, spinor_bytes); 
        cudaBindTexture(0, texNorm1, src.norm, spinor_bytes/12);
        convertSHKernel<<<blasGrid, blasBlock>>>((float4*)dst.v, src.stride, src.volume);
        cudaUnbindTexture(texHalf1); 
        cudaUnbindTexture(texNorm1);
      }else{ //nSpin== 1;
        cudaBindTexture(0, texHalfSt1, src.v, spinor_bytes); 
        cudaBindTexture(0, texNorm1, src.norm, spinor_bytes/3);
        convertSHKernel<<<blasGrid, blasBlock>>>((float2*)dst.v, src.stride, src.volume);
        cudaUnbindTexture(texHalfSt1); 
        cudaUnbindTexture(texNorm1);
      }
  } else if (dst.precision == QUDA_HALF_PRECISION && src.precision == QUDA_SINGLE_PRECISION) {
    blas_quda_bytes += dst.volume*sizeof(float);
    int spinor_bytes = src.length*sizeof(float);
    if (src.nSpin == 4){
        cudaBindTexture(0, xTexSingle4, src.v, spinor_bytes); 
        convertHSKernel<<<blasGrid, blasBlock>>>((short4*)dst.v, (float*)dst.norm, src.stride, src.volume);
        cudaUnbindTexture(xTexSingle4); 
    }else{ //nSpinr == 1
        cudaBindTexture(0, xTexSingle2, src.v, spinor_bytes); 
        convertHSKernel<<<blasGrid, blasBlock>>>((short2*)dst.v, (float*)dst.norm, src.stride, src.volume);     
        cudaUnbindTexture(xTexSingle2); 
    }
  } else if (dst.precision == QUDA_DOUBLE_PRECISION && src.precision == QUDA_HALF_PRECISION) {
    blas_quda_bytes += src.volume*sizeof(float);
    int spinor_bytes = src.length*sizeof(short);
    if (src.nSpin == 4){
        cudaBindTexture(0, texHalf1, src.v, spinor_bytes); 
        cudaBindTexture(0, texNorm1, src.norm, spinor_bytes/12);
        convertDHKernel<<<blasGrid, blasBlock>>>((double2*)dst.v, src.stride, src.volume);
        cudaUnbindTexture(texHalf1); 
        cudaUnbindTexture(texNorm1);
    }else{//nSpinr == 1
        cudaBindTexture(0, texHalfSt1, src.v, spinor_bytes); 
        cudaBindTexture(0, texNorm1, src.norm, spinor_bytes/3);
        convertDHKernelSt<<<blasGrid, blasBlock>>>((double2*)dst.v, src.stride, src.volume);
        cudaUnbindTexture(texHalfSt1); 
        cudaUnbindTexture(texNorm1);
    }
  } else if (dst.precision == QUDA_HALF_PRECISION && src.precision == QUDA_DOUBLE_PRECISION) {
    blas_quda_bytes += dst.volume*sizeof(float);
    int spinor_bytes = src.length*sizeof(double);
    cudaBindTexture(0, xTexDouble2, src.v, spinor_bytes); 
    if (src.nSpin == 4){
      convertHDKernel<<<blasGrid, blasBlock>>>((short4*)dst.v, (float*)dst.norm, src.stride, src.volume);
    }else{ //nSpinr == 1
      convertHDKernel<<<blasGrid, blasBlock>>>((short2*)dst.v, (float*)dst.norm, src.stride, src.volume);
    }
    cudaUnbindTexture(xTexDouble2); 
  } else {
    cudaMemcpy(dst.v, src.v, dst.bytes, cudaMemcpyDeviceToDevice);
    if (dst.precision == QUDA_HALF_PRECISION) {
      cudaMemcpy(dst.norm, src.norm, dst.bytes/(dst.nColor*dst.nSpin), cudaMemcpyDeviceToDevice);
      blas_quda_bytes += 2*dst.real_length*sizeof(float);
    }
  }
  
  cudaThreadSynchronize();
  if (!blasTuning) checkCudaError();

}


template <typename Float, typename Float2>
__global__ void axpbyKernel(Float a, Float2 *x, Float b, Float2 *y, int length) {
  unsigned int i = blockIdx.x*(blockDim.x) + threadIdx.x;
  unsigned int gridSize = gridDim.x*blockDim.x;
  while (i < length) {
    y[i] = a*x[i] + b*y[i];
    i += gridSize;
  } 
}

__global__ void axpbyHKernel(float a, float b, short4 *yH, float *yN, int stride, int length) {
  unsigned int i = blockIdx.x*(blockDim.x) + threadIdx.x;
  unsigned int gridSize = gridDim.x*blockDim.x;
  while (i < length) {
    RECONSTRUCT_HALF_SPINOR(x, texHalf1, texNorm1, stride);
    RECONSTRUCT_HALF_SPINOR(y, texHalf2, texNorm2, stride);
    AXPBY_FLOAT4(a, x0, b, y0);
    AXPBY_FLOAT4(a, x1, b, y1);
    AXPBY_FLOAT4(a, x2, b, y2);
    AXPBY_FLOAT4(a, x3, b, y3);
    AXPBY_FLOAT4(a, x4, b, y4);
    AXPBY_FLOAT4(a, x5, b, y5);
    CONSTRUCT_HALF_SPINOR_FROM_SINGLE(yH, yN, y, stride);
    i += gridSize;
  } 
  
}
__global__ void axpbyHKernel(float a, float b, short2 *yH, float *yN, int stride, int length) {
  unsigned int i = blockIdx.x*(blockDim.x) + threadIdx.x;
  unsigned int gridSize = gridDim.x*blockDim.x;
  while (i < length) {
    RECONSTRUCT_HALF_SPINOR_ST(x, texHalfSt1, texNorm1, stride);
    RECONSTRUCT_HALF_SPINOR_ST(y, texHalfSt2, texNorm2, stride);
    AXPBY_FLOAT2(a, x0, b, y0);
    AXPBY_FLOAT2(a, x1, b, y1);
    AXPBY_FLOAT2(a, x2, b, y2);
    CONSTRUCT_HALF_SPINOR_FROM_SINGLE_ST(yH, yN, y, stride);
    i += gridSize;
  } 
  
}


// performs the operation y[i] = a*x[i] + b*y[i]
void axpbyCuda(const double &a, cudaColorSpinorField &x, const double &b, cudaColorSpinorField &y) {
  setBlock(2, x.length, x.precision);
  checkSpinor(x, y);
  if (x.precision == QUDA_DOUBLE_PRECISION) {
    axpbyKernel<<<blasGrid, blasBlock>>>(a, (double*)x.v, b, (double*)y.v, x.length);
  } else if (x.precision == QUDA_SINGLE_PRECISION) {
    axpbyKernel<<<blasGrid, blasBlock>>>((float)a, (float2*)x.v, (float)b, (float2*)y.v, x.length/2);
  } else {
    if (x.siteSubset == QUDA_FULL_SITE_SUBSET) {
      axpbyCuda(a, x.Even(), b, y.Even());
      axpbyCuda(a, x.Odd(), b, y.Odd());
      return;
    }
    int spinor_bytes = x.length*sizeof(short);
    if (x.nSpin == 4){ //wilson
        cudaBindTexture(0, texHalf1, x.v, spinor_bytes); 
        cudaBindTexture(0, texNorm1, x.norm, spinor_bytes/12);    
        cudaBindTexture(0, texHalf2, y.v, spinor_bytes); 
        cudaBindTexture(0, texNorm2, y.norm, spinor_bytes/12);    
        axpbyHKernel<<<blasGrid, blasBlock>>>((float)a, (float)b, (short4*)y.v, 
                                              (float*)y.norm, y.stride, y.volume);
    }else if (x.nSpin == 1) {//staggered
        cudaBindTexture(0, texHalfSt1, x.v, spinor_bytes); 
        cudaBindTexture(0, texNorm1, x.norm, spinor_bytes/3);    
        cudaBindTexture(0, texHalfSt2, y.v, spinor_bytes); 
        cudaBindTexture(0, texNorm2, y.norm, spinor_bytes/3);    
        axpbyHKernel<<<blasGrid, blasBlock>>>((float)a, (float)b, (short2*)y.v, 
                                              (float*)y.norm, y.stride, y.volume);
    }else{
      errorQuda("ERROR: nSpin=%d is not supported\n", x.nSpin);
    }
    blas_quda_bytes += 3*x.volume*sizeof(float);
  }
  blas_quda_bytes += 3*x.real_length*x.precision;
  blas_quda_flops += 3*x.real_length;

  if (!blasTuning) checkCudaError();
}

template <typename Float>
__global__ void xpyKernel(Float *x, Float *y, int len) {
  unsigned int i = blockIdx.x*(blockDim.x) + threadIdx.x;
  unsigned int gridSize = gridDim.x*blockDim.x;
  while (i < len) {
    y[i] += x[i];
    i += gridSize;
  } 
}

__global__ void xpyHKernel(short4 *yH, float *yN, int stride, int length) {
  unsigned int i = blockIdx.x*(blockDim.x) + threadIdx.x;
  unsigned int gridSize = gridDim.x*blockDim.x;
  while (i < length) {
    RECONSTRUCT_HALF_SPINOR(x, texHalf1, texNorm1, stride);
    RECONSTRUCT_HALF_SPINOR(y, texHalf2, texNorm2, stride);
    XPY_FLOAT4(x0, y0);
    XPY_FLOAT4(x1, y1);
    XPY_FLOAT4(x2, y2);
    XPY_FLOAT4(x3, y3);
    XPY_FLOAT4(x4, y4);
    XPY_FLOAT4(x5, y5);
    CONSTRUCT_HALF_SPINOR_FROM_SINGLE(yH, yN, y, stride);
    i += gridSize;
  } 
  
}

__global__ void xpyHKernel(short2 *yH, float *yN, int stride, int length) {
  unsigned int i = blockIdx.x*(blockDim.x) + threadIdx.x;
  unsigned int gridSize = gridDim.x*blockDim.x;
  while (i < length) {
    RECONSTRUCT_HALF_SPINOR_ST(x, texHalfSt1, texNorm1, stride);
    RECONSTRUCT_HALF_SPINOR_ST(y, texHalfSt2, texNorm2, stride);
    XPY_FLOAT2(x0, y0);
    XPY_FLOAT2(x1, y1);
    XPY_FLOAT2(x2, y2);
    CONSTRUCT_HALF_SPINOR_FROM_SINGLE_ST(yH, yN, y, stride);
    i += gridSize;
  } 
  
}

// performs the operation y[i] = x[i] + y[i]
void xpyCuda(cudaColorSpinorField &x, cudaColorSpinorField &y) {
  checkSpinor(x,y);
  setBlock(3, x.length, x.precision);
  if (x.precision == QUDA_DOUBLE_PRECISION) {
    xpyKernel<<<blasGrid, blasBlock>>>((double*)x.v, (double*)y.v, x.length);
  } else if (x.precision == QUDA_SINGLE_PRECISION) {
    xpyKernel<<<blasGrid, blasBlock>>>((float2*)x.v, (float2*)y.v, x.length/2);
  } else {
    if (x.siteSubset == QUDA_FULL_SITE_SUBSET) {
      xpyCuda(x.Even(), y.Even());
      xpyCuda(x.Odd(), y.Odd());
      return;
    }
    int spinor_bytes = x.length*sizeof(short);
    if (x.nSpin == 4){ //wilson
      cudaBindTexture(0, texHalf1, x.v, spinor_bytes); 
      cudaBindTexture(0, texNorm1, x.norm, spinor_bytes/12);    
      cudaBindTexture(0, texHalf2, y.v, spinor_bytes); 
      cudaBindTexture(0, texNorm2, y.norm, spinor_bytes/12);    
      xpyHKernel<<<blasGrid, blasBlock>>>((short4*)y.v, (float*)y.norm, y.stride, y.volume);
    }else if (x.nSpin == 1){ //staggered
      cudaBindTexture(0, texHalfSt1, x.v, spinor_bytes); 
      cudaBindTexture(0, texNorm1, x.norm, spinor_bytes/3);    
      cudaBindTexture(0, texHalfSt2, y.v, spinor_bytes); 
      cudaBindTexture(0, texNorm2, y.norm, spinor_bytes/3);    
      xpyHKernel<<<blasGrid, blasBlock>>>((short2*)y.v, (float*)y.norm, y.stride, y.volume);
    }else{
      errorQuda("ERROR: nSpin=%d is not supported\n", x.nSpin);
    }
    blas_quda_bytes += 3*x.volume*sizeof(float);
  }
  blas_quda_bytes += 3*x.real_length*x.precision;
  blas_quda_flops += x.real_length;

  if (!blasTuning) checkCudaError();
}

template <typename Float, typename Float2>
__global__ void axpyKernel(Float a, Float2 *x, Float2 *y, int len) {
  unsigned int i = blockIdx.x*(blockDim.x) + threadIdx.x;
  unsigned int gridSize = gridDim.x*blockDim.x;
  while (i < len) {
    y[i] += a*x[i];
    i += gridSize;
  } 
}

__global__ void axpyHKernel(float a, short4 *yH, float *yN, int stride, int length) {
  unsigned int i = blockIdx.x*(blockDim.x) + threadIdx.x;
  unsigned int gridSize = gridDim.x*blockDim.x;
  while (i < length) {
    RECONSTRUCT_HALF_SPINOR(x, texHalf1, texNorm1, stride);
    RECONSTRUCT_HALF_SPINOR(y, texHalf2, texNorm2, stride);
    AXPY_FLOAT4(a, x0, y0);
    AXPY_FLOAT4(a, x1, y1);
    AXPY_FLOAT4(a, x2, y2);
    AXPY_FLOAT4(a, x3, y3);
    AXPY_FLOAT4(a, x4, y4);
    AXPY_FLOAT4(a, x5, y5);
    CONSTRUCT_HALF_SPINOR_FROM_SINGLE(yH, yN, y, stride);
    i += gridSize;
  } 
  
}


__global__ void axpyHKernel(float a, short2 *yH, float *yN, int stride, int length) {
  unsigned int i = blockIdx.x*(blockDim.x) + threadIdx.x;
  unsigned int gridSize = gridDim.x*blockDim.x;
  while (i < length) {
    RECONSTRUCT_HALF_SPINOR_ST(x, texHalfSt1, texNorm1, stride);
    RECONSTRUCT_HALF_SPINOR_ST(y, texHalfSt2, texNorm2, stride);
    AXPY_FLOAT2(a, x0, y0);
    AXPY_FLOAT2(a, x1, y1);
    AXPY_FLOAT2(a, x2, y2);    
    CONSTRUCT_HALF_SPINOR_FROM_SINGLE_ST(yH, yN, y, stride);
    i += gridSize;
  } 
  
}

// performs the operation y[i] = a*x[i] + y[i]
void axpyCuda(const double &a, cudaColorSpinorField &x, cudaColorSpinorField &y) {
  checkSpinor(x,y);
  setBlock(4, x.length, x.precision);
  if (x.precision == QUDA_DOUBLE_PRECISION) {
    axpyKernel<<<blasGrid, blasBlock>>>(a, (double*)x.v, (double*)y.v, x.length);
  } else if (x.precision == QUDA_SINGLE_PRECISION) {
    axpyKernel<<<blasGrid, blasBlock>>>((float)a, (float2*)x.v, (float2*)y.v, x.length/2);
  } else {
    if (x.siteSubset == QUDA_FULL_SITE_SUBSET) {
      axpyCuda(a, x.Even(), y.Even());
      axpyCuda(a, x.Odd(), y.Odd());
      return;
    }
    int spinor_bytes = x.length*sizeof(short);
    if (x.nSpin == 4){ //wilson
      cudaBindTexture(0, texHalf1, x.v, spinor_bytes); 
      cudaBindTexture(0, texNorm1, x.norm, spinor_bytes/12);    
      cudaBindTexture(0, texHalf2, y.v, spinor_bytes); 
      cudaBindTexture(0, texNorm2, y.norm, spinor_bytes/12);    
      axpyHKernel<<<blasGrid, blasBlock>>>((float)a, (short4*)y.v, (float*)y.norm, y.stride, y.volume);
    }else if (x.nSpin == 1){ //staggered
      cudaBindTexture(0, texHalfSt1, x.v, spinor_bytes); 
      cudaBindTexture(0, texNorm1, x.norm, spinor_bytes/3);    
      cudaBindTexture(0, texHalfSt2, y.v, spinor_bytes); 
      cudaBindTexture(0, texNorm2, y.norm, spinor_bytes/3);    
      axpyHKernel<<<blasGrid, blasBlock>>>((float)a, (short2*)y.v, (float*)y.norm, y.stride, y.volume);
    }else{
      errorQuda("ERROR: nSpin=%d is not supported\n", x.nSpin);
    }
    blas_quda_bytes += 3*x.volume*sizeof(float);
  }
  blas_quda_bytes += 3*x.real_length*x.precision;
  blas_quda_flops += 2*x.real_length;

  if (!blasTuning) checkCudaError();
}

template <typename Float, typename Float2>
__global__ void xpayKernel(const Float2 *x, Float a, Float2 *y, int len) {
  unsigned int i = blockIdx.x*(blockDim.x) + threadIdx.x;
  unsigned int gridSize = gridDim.x*blockDim.x;
  while (i < len) {
    y[i] = x[i] + a*y[i];
    i += gridSize;
  } 
}

__global__ void xpayHKernel(float a, short4 *yH, float *yN, int stride, int length) {
  unsigned int i = blockIdx.x*(blockDim.x) + threadIdx.x;
  unsigned int gridSize = gridDim.x*blockDim.x;
  while (i < length) {
    RECONSTRUCT_HALF_SPINOR(x, texHalf1, texNorm1, stride);
    RECONSTRUCT_HALF_SPINOR(y, texHalf2, texNorm2, stride);
    XPAY_FLOAT4(x0, a, y0);
    XPAY_FLOAT4(x1, a, y1);
    XPAY_FLOAT4(x2, a, y2);
    XPAY_FLOAT4(x3, a, y3);
    XPAY_FLOAT4(x4, a, y4);
    XPAY_FLOAT4(x5, a, y5);
    CONSTRUCT_HALF_SPINOR_FROM_SINGLE(yH, yN, y, stride);
    i += gridSize;
  }   
}

__global__ void xpayHKernel(float a, short2 *yH, float *yN, int stride, int length) {
  unsigned int i = blockIdx.x*(blockDim.x) + threadIdx.x;
  unsigned int gridSize = gridDim.x*blockDim.x;
  while (i < length) {
    RECONSTRUCT_HALF_SPINOR_ST(x, texHalfSt1, texNorm1, stride);
    RECONSTRUCT_HALF_SPINOR_ST(y, texHalfSt2, texNorm2, stride);
    XPAY_FLOAT2(x0, a, y0);
    XPAY_FLOAT2(x1, a, y1);
    XPAY_FLOAT2(x2, a, y2);
    CONSTRUCT_HALF_SPINOR_FROM_SINGLE_ST(yH, yN, y, stride);
    i += gridSize;
  }   
}


// performs the operation y[i] = x[i] + a*y[i]
void xpayCuda(const cudaColorSpinorField &x, const double &a, cudaColorSpinorField &y) {
  checkSpinor(x,y);
  setBlock(5, x.length, x.precision);
  if (x.precision == QUDA_DOUBLE_PRECISION) {
    xpayKernel<<<blasGrid, blasBlock>>>((double*)x.v, a, (double*)y.v, x.length);
  } else if (x.precision == QUDA_SINGLE_PRECISION) {
    xpayKernel<<<blasGrid, blasBlock>>>((float2*)x.v, (float)a, (float2*)y.v, x.length/2);
  } else {
    if (x.siteSubset == QUDA_FULL_SITE_SUBSET) {
      xpayCuda(x.Even(), a, y.Even());
      xpayCuda(x.Odd(), a, y.Odd());
      return;
    }
    int spinor_bytes = x.length*sizeof(short);
    if (x.nSpin == 4){ //wilson
      cudaBindTexture(0, texHalf1, x.v, spinor_bytes); 
      cudaBindTexture(0, texNorm1, x.norm, spinor_bytes/12);    
      cudaBindTexture(0, texHalf2, y.v, spinor_bytes); 
      cudaBindTexture(0, texNorm2, y.norm, spinor_bytes/12);    
      xpayHKernel<<<blasGrid, blasBlock>>>((float)a, (short4*)y.v, (float*)y.norm, y.stride, y.volume);
    }else if (x.nSpin ==1){ //staggered
      cudaBindTexture(0, texHalfSt1, x.v, spinor_bytes); 
      cudaBindTexture(0, texNorm1, x.norm, spinor_bytes/3);    
      cudaBindTexture(0, texHalfSt2, y.v, spinor_bytes); 
      cudaBindTexture(0, texNorm2, y.norm, spinor_bytes/3);    
      xpayHKernel<<<blasGrid, blasBlock>>>((float)a, (short2*)y.v, (float*)y.norm, y.stride, y.volume);
    }else{
      errorQuda("ERROR: nSpin=%d is not supported\n", x.nSpin);      
    }
    blas_quda_bytes += 3*x.volume*sizeof(float);
  }
  blas_quda_bytes += 3*x.real_length*x.precision;
  blas_quda_flops += 2*x.real_length;

  if (!blasTuning) checkCudaError();
}

template <typename Float>
__global__ void mxpyKernel(Float *x, Float *y, int len) {
  unsigned int i = blockIdx.x*(blockDim.x) + threadIdx.x;
  unsigned int gridSize = gridDim.x*blockDim.x;
  while (i < len) {
    y[i] -= x[i];
    i += gridSize;
  } 
}

__global__ void mxpyHKernel(short4 *yH, float *yN, int stride, int length) {
  unsigned int i = blockIdx.x*(blockDim.x) + threadIdx.x;
  unsigned int gridSize = gridDim.x*blockDim.x;
  while (i < length) {
    RECONSTRUCT_HALF_SPINOR(x, texHalf1, texNorm1, stride);
    RECONSTRUCT_HALF_SPINOR(y, texHalf2, texNorm2, stride);
    MXPY_FLOAT4(x0, y0);
    MXPY_FLOAT4(x1, y1);
    MXPY_FLOAT4(x2, y2);
    MXPY_FLOAT4(x3, y3);
    MXPY_FLOAT4(x4, y4);
    MXPY_FLOAT4(x5, y5);
    CONSTRUCT_HALF_SPINOR_FROM_SINGLE(yH, yN, y, stride);
    i += gridSize;
  } 
  
}

__global__ void mxpyHKernel(short2 *yH, float *yN, int stride, int length) {
  unsigned int i = blockIdx.x*(blockDim.x) + threadIdx.x;
  unsigned int gridSize = gridDim.x*blockDim.x;
  while (i < length) {
    RECONSTRUCT_HALF_SPINOR_ST(x, texHalfSt1, texNorm1, stride);
    RECONSTRUCT_HALF_SPINOR_ST(y, texHalfSt2, texNorm2, stride);
    MXPY_FLOAT2(x0, y0);
    MXPY_FLOAT2(x1, y1);
    MXPY_FLOAT2(x2, y2);
    CONSTRUCT_HALF_SPINOR_FROM_SINGLE_ST(yH, yN, y, stride);
    i += gridSize;
  } 
  
}


// performs the operation y[i] -= x[i] (minus x plus y)
void mxpyCuda(cudaColorSpinorField &x, cudaColorSpinorField &y) {
  checkSpinor(x,y);
  setBlock(6, x.length, x.precision);
  if (x.precision == QUDA_DOUBLE_PRECISION) {
    mxpyKernel<<<blasGrid, blasBlock>>>((double*)x.v, (double*)y.v, x.length);
  } else if (x.precision == QUDA_SINGLE_PRECISION) {
    mxpyKernel<<<blasGrid, blasBlock>>>((float2*)x.v, (float2*)y.v, x.length/2);
  } else {
    if (x.siteSubset == QUDA_FULL_SITE_SUBSET) {
      mxpyCuda(x.Even(), y.Even());
      mxpyCuda(x.Odd(), y.Odd());
      return;
    }
    int spinor_bytes = x.length*sizeof(short);
    if (x.nSpin == 4){ //wilson
      cudaBindTexture(0, texHalf1, x.v, spinor_bytes); 
      cudaBindTexture(0, texNorm1, x.norm, spinor_bytes/12);    
      cudaBindTexture(0, texHalf2, y.v, spinor_bytes); 
      cudaBindTexture(0, texNorm2, y.norm, spinor_bytes/12);    
      mxpyHKernel<<<blasGrid, blasBlock>>>((short4*)y.v, (float*)y.norm, y.stride, y.volume);
    }else if (x.nSpin == 1) { //staggered
      cudaBindTexture(0, texHalfSt1, x.v, spinor_bytes); 
      cudaBindTexture(0, texNorm1, x.norm, spinor_bytes/3);    
      cudaBindTexture(0, texHalfSt2, y.v, spinor_bytes); 
      cudaBindTexture(0, texNorm2, y.norm, spinor_bytes/3);    
      mxpyHKernel<<<blasGrid, blasBlock>>>((short2*)y.v, (float*)y.norm, y.stride, y.volume);
    }else{
      errorQuda("ERROR: nSpin=%d is not supported\n", x.nSpin);      
    }
    blas_quda_bytes += 3*x.volume*sizeof(float);
  }
  blas_quda_bytes += 3*x.real_length*x.precision;
  blas_quda_flops += x.real_length;

  if (!blasTuning) checkCudaError();
}

template <typename Float, typename Float2>
__global__ void axKernel(Float a, Float2 *x, int len) {
  unsigned int i = blockIdx.x*(blockDim.x) + threadIdx.x;
  unsigned int gridSize = gridDim.x*blockDim.x;
  while (i < len) {
    x[i] *= a;
    i += gridSize;
  } 
}

__global__ void axHKernel(float a, short4 *xH, float *xN, int stride, int length) {
  unsigned int i = blockIdx.x*(blockDim.x) + threadIdx.x;
  unsigned int gridSize = gridDim.x*blockDim.x;
  while (i < length) {
    RECONSTRUCT_HALF_SPINOR(x, texHalf1, texNorm1, stride);
    AX_FLOAT4(a, x0); AX_FLOAT4(a, x1); AX_FLOAT4(a, x2);
    AX_FLOAT4(a, x3); AX_FLOAT4(a, x4); AX_FLOAT4(a, x5);
    CONSTRUCT_HALF_SPINOR_FROM_SINGLE(xH, xN, x, stride);
    i += gridSize;
  } 
  
}

__global__ void axHKernel(float a, short2 *xH, float *xN, int stride, int length) {
  unsigned int i = blockIdx.x*(blockDim.x) + threadIdx.x;
  unsigned int gridSize = gridDim.x*blockDim.x;
  while (i < length) {
    RECONSTRUCT_HALF_SPINOR_ST(x, texHalfSt1, texNorm1, stride);
    AX_FLOAT2(a, x0); AX_FLOAT2(a, x1); AX_FLOAT2(a, x2);
    CONSTRUCT_HALF_SPINOR_FROM_SINGLE_ST(xH, xN, x, stride);
    i += gridSize;
  } 
  
}

// performs the operation x[i] = a*x[i]
void axCuda(const double &a, cudaColorSpinorField &x) {
  setBlock(7, x.length, x.precision);
  if (x.precision == QUDA_DOUBLE_PRECISION) {
    axKernel<<<blasGrid, blasBlock>>>(a, (double*)x.v, x.length);
  } else if (x.precision == QUDA_SINGLE_PRECISION) {
    axKernel<<<blasGrid, blasBlock>>>((float)a, (float2*)x.v, x.length/2);
  } else {
    if (x.siteSubset == QUDA_FULL_SITE_SUBSET) {
      axCuda(a, x.Even());
      axCuda(a, x.Odd());
      return;
    }
    int spinor_bytes = x.length*sizeof(short);
    if (x.nSpin == 4){ //wilson
      cudaBindTexture(0, texHalf1, x.v, spinor_bytes); 
      cudaBindTexture(0, texNorm1, x.norm, spinor_bytes/12);    
      axHKernel<<<blasGrid, blasBlock>>>((float)a, (short4*)x.v, (float*)x.norm, x.stride, x.volume);
    }else if (x.nSpin ==1){ //staggered
      cudaBindTexture(0, texHalfSt1, x.v, spinor_bytes); 
      cudaBindTexture(0, texNorm1, x.norm, spinor_bytes/3);    
      axHKernel<<<blasGrid, blasBlock>>>((float)a, (short2*)x.v, (float*)x.norm, x.stride, x.volume);
    }else{
      errorQuda("ERROR: nSpin=%d is not supported\n", x.nSpin);           
    }
    blas_quda_bytes += 2*x.volume*sizeof(float);
  }
  blas_quda_bytes += 2*x.real_length*x.precision;
  blas_quda_flops += x.real_length;

  if (!blasTuning) checkCudaError();
}

template <typename Float2>
__global__ void caxpyDKernel(Float2 a, Float2 *x, Float2 *y, int len) {
  
  unsigned int i = blockIdx.x*(blockDim.x) + threadIdx.x;
  unsigned int gridSize = gridDim.x*blockDim.x;
  while (i < len) {
    Float2 Z = READ_DOUBLE2_TEXTURE(x, i);
    y[i].x += a.x*Z.x - a.y*Z.y;
    y[i].y += a.y*Z.x + a.x*Z.y;
    i += gridSize;
  } 
  
}

template <typename Float2>
__global__ void caxpySKernel(Float2 a, Float2 *x, Float2 *y, int len) {
  
  unsigned int i = blockIdx.x*(blockDim.x) + threadIdx.x;
  unsigned int gridSize = gridDim.x*blockDim.x;
  while (i < len) {
    Float2 Z = read_Float2(x, i);
    y[i].x += a.x*Z.x - a.y*Z.y;
    y[i].y += a.y*Z.x + a.x*Z.y;
    i += gridSize;
  } 
  
}

__global__ void caxpyHKernel(float2 a, short4 *yH, float *yN, int stride, int length) {
  
  unsigned int i = blockIdx.x*(blockDim.x) + threadIdx.x;
  unsigned int gridSize = gridDim.x*blockDim.x;
  while (i < length) {
    RECONSTRUCT_HALF_SPINOR(x, texHalf1, texNorm1, stride);
    RECONSTRUCT_HALF_SPINOR(y, texHalf2, texNorm2, stride);
    CAXPY_FLOAT4(a, x0, y0);
    CAXPY_FLOAT4(a, x1, y1);
    CAXPY_FLOAT4(a, x2, y2);
    CAXPY_FLOAT4(a, x3, y3);
    CAXPY_FLOAT4(a, x4, y4);
    CAXPY_FLOAT4(a, x5, y5);
    CONSTRUCT_HALF_SPINOR_FROM_SINGLE(yH, yN, y, stride);
    i += gridSize;
  } 
  
}

__global__ void caxpyHKernel(float2 a, short2 *yH, float *yN, int stride, int length) {
  
  unsigned int i = blockIdx.x*(blockDim.x) + threadIdx.x;
  unsigned int gridSize = gridDim.x*blockDim.x;
  while (i < length) {
    RECONSTRUCT_HALF_SPINOR_ST(x, texHalfSt1, texNorm1, stride);
    RECONSTRUCT_HALF_SPINOR_ST(y, texHalfSt2, texNorm2, stride);
    CAXPY_FLOAT2(a, x0, y0);
    CAXPY_FLOAT2(a, x1, y1);
    CAXPY_FLOAT2(a, x2, y2);
    CONSTRUCT_HALF_SPINOR_FROM_SINGLE_ST(yH, yN, y, stride);
    i += gridSize;
  } 
  
}

// performs the operation y[i] += a*x[i]
void caxpyCuda(const Complex &a, cudaColorSpinorField &x, cudaColorSpinorField &y) {
  checkSpinor(x,y);
  int length = x.length/2;
  setBlock(8, length, x.precision);
  blas_quda_bytes += 3*x.real_length*x.precision;
  blas_quda_flops += 4*x.real_length;
  if (x.precision == QUDA_DOUBLE_PRECISION) {
    int spinor_bytes = x.length*sizeof(double);
    cudaBindTexture(0, xTexDouble2, x.v, spinor_bytes); 
    cudaBindTexture(0, yTexDouble2, y.v, spinor_bytes); 
    double2 a2 = make_double2(real(a), imag(a));
    caxpyDKernel<<<blasGrid, blasBlock>>>(a2, (double2*)x.v, (double2*)y.v, length);
  } else if (x.precision == QUDA_SINGLE_PRECISION) {
    float2 a2 = make_float2(real(a), imag(a));
    caxpySKernel<<<blasGrid, blasBlock>>>(a2, (float2*)x.v, (float2*)y.v, length);
  } else {
    if (x.siteSubset == QUDA_FULL_SITE_SUBSET) {
      caxpyCuda(a, x.Even(), y.Even());
      caxpyCuda(a, x.Odd(), y.Odd());
      return;
    }
    int spinor_bytes = x.length*sizeof(short);
    if (x.nSpin == 4){ //wilson
      cudaBindTexture(0, texHalf1, x.v, spinor_bytes); 
      cudaBindTexture(0, texNorm1, x.norm, spinor_bytes/12);    
      cudaBindTexture(0, texHalf2, y.v, spinor_bytes); 
      cudaBindTexture(0, texNorm2, y.norm, spinor_bytes/12);    
      float2 a2 = make_float2(real(a), imag(a));
      caxpyHKernel<<<blasGrid, blasBlock>>>(a2, (short4*)y.v, (float*)y.norm, y.stride, y.volume);
    } else if (x.nSpin == 1){ //staggered
      cudaBindTexture(0, texHalfSt1, x.v, spinor_bytes); 
      cudaBindTexture(0, texNorm1, x.norm, spinor_bytes/3);    
      cudaBindTexture(0, texHalfSt2, y.v, spinor_bytes); 
      cudaBindTexture(0, texNorm2, y.norm, spinor_bytes/3);    
      float2 a2 = make_float2(real(a), imag(a));
      caxpyHKernel<<<blasGrid, blasBlock>>>(a2, (short2*)y.v, (float*)y.norm, y.stride, y.volume);
    }else{
      errorQuda("ERROR: nSpin=%d is not supported\n", x.nSpin);     
    }
    blas_quda_bytes += 3*x.volume*sizeof(float);
  }

  if (!blasTuning) checkCudaError();
}

template <typename Float2>
__global__ void caxpbyDKernel(Float2 a, Float2 *x, Float2 b, Float2 *y, int len) {
  
  unsigned int i = blockIdx.x*(blockDim.x) + threadIdx.x;
  unsigned int gridSize = gridDim.x*blockDim.x;
  while (i < len) {
    Float2 Z1 = READ_DOUBLE2_TEXTURE(x, i);
    Float2 Z2 = READ_DOUBLE2_TEXTURE(y, i);
    y[i].x = a.x*Z1.x + b.x*Z2.x - a.y*Z1.y - b.y*Z2.y;
    y[i].y = a.y*Z1.x + b.y*Z2.x + a.x*Z1.y + b.x*Z2.y;
    i += gridSize;
  } 
}

template <typename Float2>
__global__ void caxpbySKernel(Float2 a, Float2 *x, Float2 b, Float2 *y, int len) {
  
  unsigned int i = blockIdx.x*(blockDim.x) + threadIdx.x;
  unsigned int gridSize = gridDim.x*blockDim.x;
  while (i < len) {
    Float2 Z1 = read_Float2(x, i);
    Float2 Z2 = read_Float2(y, i);
    y[i].x = a.x*Z1.x + b.x*Z2.x - a.y*Z1.y - b.y*Z2.y;
    y[i].y = a.y*Z1.x + b.y*Z2.x + a.x*Z1.y + b.x*Z2.y;
    i += gridSize;
  } 
}

__global__ void caxpbyHKernel(float2 a, float2 b, short4 *yH, float *yN, int stride, int length) {
  
  unsigned int i = blockIdx.x*(blockDim.x) + threadIdx.x;
  unsigned int gridSize = gridDim.x*blockDim.x;
  while (i < length) {
    RECONSTRUCT_HALF_SPINOR(x, texHalf1, texNorm1, stride);
    RECONSTRUCT_HALF_SPINOR(y, texHalf2, texNorm2, stride);
    CAXPBY_FLOAT4(a, x0, b, y0);
    CAXPBY_FLOAT4(a, x1, b, y1);
    CAXPBY_FLOAT4(a, x2, b, y2);
    CAXPBY_FLOAT4(a, x3, b, y3);
    CAXPBY_FLOAT4(a, x4, b, y4);
    CAXPBY_FLOAT4(a, x5, b, y5);
    CONSTRUCT_HALF_SPINOR_FROM_SINGLE(yH, yN, y, stride);
    i += gridSize;
  }   
}

__global__ void caxpbyHKernel(float2 a, float2 b, short2 *yH, float *yN, int stride, int length) {
  
  unsigned int i = blockIdx.x*(blockDim.x) + threadIdx.x;
  unsigned int gridSize = gridDim.x*blockDim.x;
  while (i < length) {
    RECONSTRUCT_HALF_SPINOR_ST(x, texHalfSt1, texNorm1, stride);
    RECONSTRUCT_HALF_SPINOR_ST(y, texHalfSt2, texNorm2, stride);
    CAXPBY_FLOAT2(a, x0, b, y0);
    CAXPBY_FLOAT2(a, x1, b, y1);
    CAXPBY_FLOAT2(a, x2, b, y2);
    CONSTRUCT_HALF_SPINOR_FROM_SINGLE_ST(yH, yN, y, stride);
    i += gridSize;
  }   
}


// performs the operation y[i] = c*x[i] + b*y[i]
void caxpbyCuda(const Complex &a, cudaColorSpinorField &x, const Complex &b, cudaColorSpinorField &y) {
  checkSpinor(x,y);
  int length = x.length/2;
  setBlock(9, length, x.precision);
  blas_quda_bytes += 3*x.real_length*x.precision;
  blas_quda_flops += 7*x.real_length;
  if (x.precision == QUDA_DOUBLE_PRECISION) {
    int spinor_bytes = x.length*sizeof(double);
    cudaBindTexture(0, xTexDouble2, x.v, spinor_bytes); 
    cudaBindTexture(0, yTexDouble2, y.v, spinor_bytes); 
    double2 a2 = make_double2(real(a), imag(a));
    double2 b2 = make_double2(real(b), imag(b));
    caxpbyDKernel<<<blasGrid, blasBlock>>>(a2, (double2*)x.v, b2, (double2*)y.v, length);
  } else if (x.precision == QUDA_SINGLE_PRECISION) {
    float2 a2 = make_float2(real(a), imag(a));
    float2 b2 = make_float2(real(b), imag(b));
    caxpbySKernel<<<blasGrid, blasBlock>>>(a2, (float2*)x.v, b2, (float2*)y.v, length);
  } else {
    if (x.siteSubset == QUDA_FULL_SITE_SUBSET) {
      caxpbyCuda(a, x.Even(), b, y.Even());
      caxpbyCuda(a, x.Odd(), b, y.Odd());
      return;
    }
    int spinor_bytes = x.length*sizeof(short);
    if (x.nSpin == 4){ //wilson
      cudaBindTexture(0, texHalf1, x.v, spinor_bytes); 
      cudaBindTexture(0, texNorm1, x.norm, spinor_bytes/12);    
      cudaBindTexture(0, texHalf2, y.v, spinor_bytes); 
      cudaBindTexture(0, texNorm2, y.norm, spinor_bytes/12);    
      float2 a2 = make_float2(real(a), imag(a));
      float2 b2 = make_float2(real(b), imag(b));
      caxpbyHKernel<<<blasGrid, blasBlock>>>(a2, b2, (short4*)y.v, (float*)y.norm, y.stride, y.volume);
    }else if (x.nSpin == 1){ //staggered
      cudaBindTexture(0, texHalfSt1, x.v, spinor_bytes); 
      cudaBindTexture(0, texNorm1, x.norm, spinor_bytes/3);    
      cudaBindTexture(0, texHalfSt2, y.v, spinor_bytes); 
      cudaBindTexture(0, texNorm2, y.norm, spinor_bytes/3);    
      float2 a2 = make_float2(real(a), imag(a));
      float2 b2 = make_float2(real(b), imag(b));
      caxpbyHKernel<<<blasGrid, blasBlock>>>(a2, b2, (short2*)y.v, (float*)y.norm, y.stride, y.volume);
    }else{
      errorQuda("ERROR: nSpin=%d is not supported\n", x.nSpin);     
    }
    blas_quda_bytes += 3*x.volume*sizeof(float);
  }

  if (!blasTuning) checkCudaError();
}

template <typename Float2>
__global__ void cxpaypbzDKernel(Float2 *x, Float2 a, Float2 *y, Float2 b, Float2 *z, int len) {
  
  unsigned int i = blockIdx.x*(blockDim.x) + threadIdx.x;
  unsigned int gridSize = gridDim.x*blockDim.x;
  while (i < len) {
    Float2 T1 = READ_DOUBLE2_TEXTURE(x, i);
    Float2 T2 = READ_DOUBLE2_TEXTURE(y, i);
    Float2 T3 = read_Float2(z, i);

    T1.x += a.x*T2.x - a.y*T2.y;
    T1.y += a.y*T2.x + a.x*T2.y;
    T1.x += b.x*T3.x - b.y*T3.y;
    T1.y += b.y*T3.x + b.x*T3.y;
    
    z[i] = make_Float2(T1);
    i += gridSize;
  } 
  
}

template <typename Float2>
__global__ void cxpaypbzSKernel(Float2 *x, Float2 a, Float2 *y, Float2 b, Float2 *z, int len) {
  
  unsigned int i = blockIdx.x*(blockDim.x) + threadIdx.x;
  unsigned int gridSize = gridDim.x*blockDim.x;
  while (i < len) {
    Float2 T1 = read_Float2(x, i);
    Float2 T2 = read_Float2(y, i);
    Float2 T3 = read_Float2(z, i);

    T1.x += a.x*T2.x - a.y*T2.y;
    T1.y += a.y*T2.x + a.x*T2.y;
    T1.x += b.x*T3.x - b.y*T3.y;
    T1.y += b.y*T3.x + b.x*T3.y;
    
    z[i] = make_Float2(T1);
    i += gridSize;
  } 
  
}

__global__ void cxpaypbzHKernel(float2 a, float2 b, short4 *zH, float *zN, int stride, int length) {
  
  unsigned int i = blockIdx.x*(blockDim.x) + threadIdx.x;
  unsigned int gridSize = gridDim.x*blockDim.x;
  while (i < length) {
    RECONSTRUCT_HALF_SPINOR(x, texHalf1, texNorm1, stride);
    RECONSTRUCT_HALF_SPINOR(y, texHalf2, texNorm2, stride);
    RECONSTRUCT_HALF_SPINOR(z, texHalf3, texNorm3, stride);
    CXPAYPBZ_FLOAT4(x0, a, y0, b, z0);
    CXPAYPBZ_FLOAT4(x1, a, y1, b, z1);
    CXPAYPBZ_FLOAT4(x2, a, y2, b, z2);
    CXPAYPBZ_FLOAT4(x3, a, y3, b, z3);
    CXPAYPBZ_FLOAT4(x4, a, y4, b, z4);
    CXPAYPBZ_FLOAT4(x5, a, y5, b, z5);
    CONSTRUCT_HALF_SPINOR_FROM_SINGLE(zH, zN, z, stride);
    i += gridSize;
  }   
}


__global__ void cxpaypbzHKernel(float2 a, float2 b, short2 *zH, float *zN, int stride, int length) {
  
  unsigned int i = blockIdx.x*(blockDim.x) + threadIdx.x;
  unsigned int gridSize = gridDim.x*blockDim.x;
  while (i < length) {
    RECONSTRUCT_HALF_SPINOR_ST(x, texHalfSt1, texNorm1, stride);
    RECONSTRUCT_HALF_SPINOR_ST(y, texHalfSt2, texNorm2, stride);
    RECONSTRUCT_HALF_SPINOR_ST(z, texHalfSt3, texNorm3, stride);
    CXPAYPBZ_FLOAT2(x0, a, y0, b, z0);
    CXPAYPBZ_FLOAT2(x1, a, y1, b, z1);
    CXPAYPBZ_FLOAT2(x2, a, y2, b, z2);
    CONSTRUCT_HALF_SPINOR_FROM_SINGLE_ST(zH, zN, z, stride);
    i += gridSize;
  }   
}


// performs the operation z[i] = x[i] + a*y[i] + b*z[i]
void cxpaypbzCuda(cudaColorSpinorField &x, const Complex &a, cudaColorSpinorField &y, 
                  const Complex &b, cudaColorSpinorField &z) {
  checkSpinor(x,y);
  checkSpinor(x,z);
  int length = x.length/2;
  setBlock(10, length, x.precision);
  blas_quda_bytes += 4*x.real_length*x.precision;
  blas_quda_flops += 8*x.real_length;
  if (x.precision == QUDA_DOUBLE_PRECISION) {
    int spinor_bytes = x.length*sizeof(double);
    cudaBindTexture(0, xTexDouble2, x.v, spinor_bytes); 
    cudaBindTexture(0, yTexDouble2, y.v, spinor_bytes); 
    double2 a2 = make_double2(real(a), imag(a));
    double2 b2 = make_double2(real(b), imag(b));
    cxpaypbzDKernel<<<blasGrid, blasBlock>>>((double2*)x.v, a2, (double2*)y.v, b2, (double2*)z.v, length);
  } else if (x.precision == QUDA_SINGLE_PRECISION) {
    float2 a2 = make_float2(real(a), imag(a));
    float2 b2 = make_float2(real(b), imag(b));
    cxpaypbzSKernel<<<blasGrid, blasBlock>>>((float2*)x.v, a2, (float2*)y.v, b2, (float2*)z.v, length);
  } else {
    if (x.siteSubset == QUDA_FULL_SITE_SUBSET) {
      cxpaypbzCuda(x.Even(), a, y.Even(), b, z.Even());
      cxpaypbzCuda(x.Odd(), a, y.Odd(), b, z.Odd());
      return;
    }
    int spinor_bytes = x.length*sizeof(short);
    blas_quda_bytes += 4*x.volume*sizeof(float);
    if (x.nSpin ==4 ){//wilson 
      cudaBindTexture(0, texHalf1, x.v, spinor_bytes); 
      cudaBindTexture(0, texNorm1, x.norm, spinor_bytes/12);    
      cudaBindTexture(0, texHalf2, y.v, spinor_bytes); 
      cudaBindTexture(0, texNorm2, y.norm, spinor_bytes/12);    
      cudaBindTexture(0, texHalf3, z.v, spinor_bytes); 
      cudaBindTexture(0, texNorm3, z.norm, spinor_bytes/12);    
      float2 a2 = make_float2(real(a), imag(a));
      float2 b2 = make_float2(real(b), imag(b));
      cxpaypbzHKernel<<<blasGrid, blasBlock>>>(a2, b2, (short4*)z.v, (float*)z.norm, z.stride, z.volume);
    } else if (x.nSpin ==1 ){//staggered 
      cudaBindTexture(0, texHalfSt1, x.v, spinor_bytes); 
      cudaBindTexture(0, texNorm1, x.norm, spinor_bytes/3);    
      cudaBindTexture(0, texHalfSt2, y.v, spinor_bytes); 
      cudaBindTexture(0, texNorm2, y.norm, spinor_bytes/3);    
      cudaBindTexture(0, texHalfSt3, z.v, spinor_bytes); 
      cudaBindTexture(0, texNorm3, z.norm, spinor_bytes/3);    
      float2 a2 = make_float2(real(a), imag(a));
      float2 b2 = make_float2(real(b), imag(b));
      cxpaypbzHKernel<<<blasGrid, blasBlock>>>(a2, b2, (short2*)z.v, (float*)z.norm, z.stride, z.volume);
    }else{
      errorQuda("ERROR: nSpin=%d is not supported\n", x.nSpin);           
    }
  }

  if (!blasTuning) checkCudaError();
}


template <typename Float, typename Float2>
__global__ void axpyBzpcxDKernel(Float a, Float2 *x, Float2 *y, Float b, Float2 *z, Float c, int len)
{
  unsigned int i = blockIdx.x*(blockDim.x) + threadIdx.x;
  unsigned int gridSize = gridDim.x*blockDim.x;
  while (i < len) {
    Float2 x_i = READ_DOUBLE2_TEXTURE(x, i);
    Float2 z_i = READ_DOUBLE2_TEXTURE(z, i);

    y[i].x += a*x_i.x;
    y[i].y += a*x_i.y;
    
    x[i].x = b*z_i.x + c*x_i.x;
    x[i].y = b*z_i.y + c*x_i.y;
    
    i += gridSize;
  }
}


template <typename Float, typename Float2>
__global__ void axpyBzpcxSKernel(Float a, Float2 *x, Float2 *y, Float b, Float2 *z, Float c, int len)
{
  unsigned int i = blockIdx.x*(blockDim.x) + threadIdx.x;
  unsigned int gridSize = gridDim.x*blockDim.x;
  while (i < len) {
    Float2 x_i = read_Float2(x, i);
    Float2 z_i = read_Float2(z, i);

    y[i].x += a*x_i.x;
    y[i].y += a*x_i.y;
    
    x[i].x = b*z_i.x + c*x_i.x;
    x[i].y = b*z_i.y + c*x_i.y;
    
    i += gridSize;
  }
}

__global__ void axpyBzpcxHKernel(float a, float b, float c, short4 *xH, float *xN, short4 *yH, float *yN, int stride, int length)
{
  
  unsigned int i = blockIdx.x*(blockDim.x) + threadIdx.x;
  unsigned int gridSize = gridDim.x*blockDim.x;
  while (i < length) {
    RECONSTRUCT_HALF_SPINOR(x, texHalf1, texNorm1, stride);
    RECONSTRUCT_HALF_SPINOR(y, texHalf2, texNorm2, stride);
    RECONSTRUCT_HALF_SPINOR(z, texHalf3, texNorm3, stride);
    AXPY_FLOAT4(a, x0, y0);
    AXPBY_FLOAT4(b, z0, c, x0);
    AXPY_FLOAT4(a, x1, y1);
    AXPBY_FLOAT4(b, z1, c, x1);
    AXPY_FLOAT4(a, x2, y2);
    AXPBY_FLOAT4(b, z2, c, x2);
    AXPY_FLOAT4(a, x3, y3);
    AXPBY_FLOAT4(b, z3, c, x3);
    AXPY_FLOAT4(a, x4, y4);
    AXPBY_FLOAT4(b, z4, c, x4);
    AXPY_FLOAT4(a, x5, y5);
    AXPBY_FLOAT4(b, z5, c, x5);    
    CONSTRUCT_HALF_SPINOR_FROM_SINGLE(yH, yN, y, stride);
    CONSTRUCT_HALF_SPINOR_FROM_SINGLE(xH, xN, x, stride);
    i += gridSize;
  }
}


__global__ void axpyBzpcxHKernel(float a, float b, float c, short2 *xH, float *xN, short2 *yH, float *yN, int stride, int length)
{
  
  unsigned int i = blockIdx.x*(blockDim.x) + threadIdx.x;
  unsigned int gridSize = gridDim.x*blockDim.x;
  while (i < length) {
    RECONSTRUCT_HALF_SPINOR_ST(x, texHalfSt1, texNorm1, stride);
    RECONSTRUCT_HALF_SPINOR_ST(y, texHalfSt2, texNorm2, stride);
    RECONSTRUCT_HALF_SPINOR_ST(z, texHalfSt3, texNorm3, stride);
    AXPY_FLOAT2(a, x0, y0);
    AXPBY_FLOAT2(b, z0, c, x0);
    AXPY_FLOAT2(a, x1, y1);
    AXPBY_FLOAT2(b, z1, c, x1);
    AXPY_FLOAT2(a, x2, y2);
    AXPBY_FLOAT2(b, z2, c, x2);
    CONSTRUCT_HALF_SPINOR_FROM_SINGLE_ST(yH, yN, y, stride);
    CONSTRUCT_HALF_SPINOR_FROM_SINGLE_ST(xH, xN, x, stride);
    i += gridSize;
  }
}

// performs the operations: {y[i] = a*x[i] + y[i]; x[i] = b*z[i] + c*x[i]}
void axpyBzpcxCuda(const double &a, cudaColorSpinorField& x, cudaColorSpinorField& y, const double &b, 
                   cudaColorSpinorField& z, const double &c) 
{
  checkSpinor(x,y);
  checkSpinor(x,z);
  setBlock(11, x.length, x.precision);
  if (x.precision == QUDA_DOUBLE_PRECISION) {
    int spinor_bytes = x.length*sizeof(double);
    cudaBindTexture(0, xTexDouble2, x.v, spinor_bytes); 
    cudaBindTexture(0, zTexDouble2, z.v, spinor_bytes); 
    axpyBzpcxDKernel<<<blasGrid, blasBlock>>>(a, (double2*)x.v, (double2*)y.v, b, (double2*)z.v, c, x.length/2);
  } else if (x.precision == QUDA_SINGLE_PRECISION) {
    axpyBzpcxSKernel<<<blasGrid, blasBlock>>>((float)a, (float2*)x.v, (float2*)y.v, (float)b, (float2*)z.v, (float)c, x.length/2);
  } else {
    if (x.siteSubset == QUDA_FULL_SITE_SUBSET){
      axpyBzpcxCuda(a, x.Even(), y.Even(), b, z.Even(), c);
      axpyBzpcxCuda(a, x.Odd(), y.Odd(), b, z.Odd(), c);
      return ;
    }

    int spinor_bytes = x.length*sizeof(short);
    if (x.nSpin == 4){ //wilson
      cudaBindTexture(0, texHalf1, x.v, spinor_bytes);
      cudaBindTexture(0, texNorm1, x.norm, spinor_bytes/12);
      cudaBindTexture(0, texHalf2, y.v, spinor_bytes);
      cudaBindTexture(0, texNorm2, y.norm, spinor_bytes/12);
      cudaBindTexture(0, texHalf3, z.v, spinor_bytes);
      cudaBindTexture(0, texNorm3, z.norm, spinor_bytes/12);
      axpyBzpcxHKernel<<<blasGrid, blasBlock>>>((float)a, (float)b, (float)c, (short4*)x.v, (float*)x.norm,
                                              (short4*)y.v, (float*)y.norm, z.stride, z.volume);
    }else if (x.nSpin == 1){ //staggered
      cudaBindTexture(0, texHalfSt1, x.v, spinor_bytes);
      cudaBindTexture(0, texNorm1, x.norm, spinor_bytes/3);
      cudaBindTexture(0, texHalfSt2, y.v, spinor_bytes);
      cudaBindTexture(0, texNorm2, y.norm, spinor_bytes/3);
      cudaBindTexture(0, texHalfSt3, z.v, spinor_bytes);
      cudaBindTexture(0, texNorm3, z.norm, spinor_bytes/3);
      axpyBzpcxHKernel<<<blasGrid, blasBlock>>>((float)a, (float)b, (float)c, (short2*)x.v, (float*)x.norm,
                                              (short2*)y.v, (float*)y.norm, z.stride, z.volume);
    }else{
      errorQuda("ERROR: nSpin=%d is not supported\n", x.nSpin);    
    }
    blas_quda_bytes += 5*x.volume*sizeof(float);
  }
  blas_quda_bytes += 5*x.real_length*x.precision;
  blas_quda_flops += 10*x.real_length;

  if (!blasTuning) checkCudaError();
}



template <typename Float, typename Float2>
__global__ void axpyZpbxDKernel(Float a, Float2 *x, Float2 *y, Float2 *z, Float b, int len) {
  unsigned int i = blockIdx.x*(blockDim.x) + threadIdx.x;
  unsigned int gridSize = gridDim.x*blockDim.x;
  while (i < len) {
    Float2 x_i = READ_DOUBLE2_TEXTURE(x, i);
    Float2 z_i = READ_DOUBLE2_TEXTURE(z, i);
    y[i].x += a*x_i.x;
    y[i].y += a*x_i.y;
    x[i].x = z_i.x + b*x_i.x;
    x[i].y = z_i.y + b*x_i.y;
    i += gridSize;
  }
}

template <typename Float, typename Float2>
__global__ void axpyZpbxSKernel(Float a, Float2 *x, Float2 *y, Float2 *z, Float b, int len) {
  unsigned int i = blockIdx.x*(blockDim.x) + threadIdx.x;
  unsigned int gridSize = gridDim.x*blockDim.x;
  while (i < len) {
    Float2 x_i = read_Float2(x, i);
    Float2 z_i = read_Float2(z, i);
    y[i].x += a*x_i.x;
    y[i].y += a*x_i.y;
    x[i].x = z_i.x + b*x_i.x;
    x[i].y = z_i.y + b*x_i.y;
    i += gridSize;
  }
}

__global__ void axpyZpbxHKernel(float a, float b, short4 *xH, float *xN, short4 *yH, float *yN, int stride, int length) {
  
  unsigned int i = blockIdx.x*(blockDim.x) + threadIdx.x;
  unsigned int gridSize = gridDim.x*blockDim.x;
  while (i < length) {
    RECONSTRUCT_HALF_SPINOR(x, texHalf1, texNorm1, stride);
    RECONSTRUCT_HALF_SPINOR(y, texHalf2, texNorm2, stride);
    AXPY_FLOAT4(a, x0, y0);
    AXPY_FLOAT4(a, x1, y1);
    AXPY_FLOAT4(a, x2, y2);
    AXPY_FLOAT4(a, x3, y3);
    AXPY_FLOAT4(a, x4, y4);
    AXPY_FLOAT4(a, x5, y5);
    CONSTRUCT_HALF_SPINOR_FROM_SINGLE(yH, yN, y, stride);
    RECONSTRUCT_HALF_SPINOR(z, texHalf3, texNorm3, stride);
    XPAY_FLOAT4(z0, b, x0);
    XPAY_FLOAT4(z1, b, x1);
    XPAY_FLOAT4(z2, b, x2);
    XPAY_FLOAT4(z3, b, x3);
    XPAY_FLOAT4(z4, b, x4);
    XPAY_FLOAT4(z5, b, x5);
    CONSTRUCT_HALF_SPINOR_FROM_SINGLE(xH, xN, x, stride);
    i += gridSize;
  }   
}

__global__ void axpyZpbxHKernel(float a, float b, short2 *xH, float *xN, short2 *yH, float *yN, int stride, int length) {
  
  unsigned int i = blockIdx.x*(blockDim.x) + threadIdx.x;
  unsigned int gridSize = gridDim.x*blockDim.x;
  while (i < length) {
    RECONSTRUCT_HALF_SPINOR_ST(x, texHalfSt1, texNorm1, stride);
    RECONSTRUCT_HALF_SPINOR_ST(y, texHalfSt2, texNorm2, stride);
    RECONSTRUCT_HALF_SPINOR_ST(z, texHalfSt3, texNorm3, stride);
    AXPY_FLOAT2(a, x0, y0);
    XPAY_FLOAT2(z0, b, x0);
    AXPY_FLOAT2(a, x1, y1);
    XPAY_FLOAT2(z1, b, x1);
    AXPY_FLOAT2(a, x2, y2);
    XPAY_FLOAT2(z2, b, x2);
    CONSTRUCT_HALF_SPINOR_FROM_SINGLE_ST(yH, yN, y, stride);
    CONSTRUCT_HALF_SPINOR_FROM_SINGLE_ST(xH, xN, x, stride);
    i += gridSize;
  }   
}


// performs the operations: {y[i] = a*x[i] + y[i]; x[i] = z[i] + b*x[i]}
void axpyZpbxCuda(const double &a, cudaColorSpinorField &x, cudaColorSpinorField &y, 
                  cudaColorSpinorField &z, const double &b) {
  checkSpinor(x,y);
  checkSpinor(x,z);
  setBlock(12, x.length, x.precision);
  if (x.precision == QUDA_DOUBLE_PRECISION) {
    int spinor_bytes = x.length*sizeof(double);
    cudaBindTexture(0, xTexDouble2, x.v, spinor_bytes); 
    cudaBindTexture(0, zTexDouble2, z.v, spinor_bytes); 
    axpyZpbxDKernel<<<blasGrid, blasBlock>>>
      (a, (double2*)x.v, (double2*)y.v, (double2*)z.v, b, x.length/2);
  } else if (x.precision == QUDA_SINGLE_PRECISION) {
    axpyZpbxSKernel<<<blasGrid, blasBlock>>>
      ((float)a, (float2*)x.v, (float2*)y.v, (float2*)z.v, (float)b, x.length/2);
  } else {
    if (x.siteSubset == QUDA_FULL_SITE_SUBSET) {
      axpyZpbxCuda(a, x.Even(), y.Even(), z.Even(), b);
      axpyZpbxCuda(a, x.Odd(), y.Odd(), z.Odd(), b);
      return;
    }
    int spinor_bytes = x.length*sizeof(short);
    if (x.nSpin ==4){ //wilson 
      cudaBindTexture(0, texHalf1, x.v, spinor_bytes); 
      cudaBindTexture(0, texNorm1, x.norm, spinor_bytes/12);    
      cudaBindTexture(0, texHalf2, y.v, spinor_bytes); 
      cudaBindTexture(0, texNorm2, y.norm, spinor_bytes/12);    
      cudaBindTexture(0, texHalf3, z.v, spinor_bytes); 
      cudaBindTexture(0, texNorm3, z.norm, spinor_bytes/12);    
      axpyZpbxHKernel<<<blasGrid, blasBlock>>>((float)a, (float)b, (short4*)x.v, (float*)x.norm,
                                               (short4*)y.v, (float*)y.norm, z.stride, z.volume);
    }else if (x.nSpin == 1){ //staggered
      cudaBindTexture(0, texHalfSt1, x.v, spinor_bytes); 
      cudaBindTexture(0, texNorm1, x.norm, spinor_bytes/3);    
      cudaBindTexture(0, texHalfSt2, y.v, spinor_bytes); 
      cudaBindTexture(0, texNorm2, y.norm, spinor_bytes/3);    
      cudaBindTexture(0, texHalfSt3, z.v, spinor_bytes); 
      cudaBindTexture(0, texNorm3, z.norm, spinor_bytes/3);    
      axpyZpbxHKernel<<<blasGrid, blasBlock>>>((float)a, (float)b, (short2*)x.v, (float*)x.norm,
                                               (short2*)y.v, (float*)y.norm, z.stride, z.volume);
    }else{
      errorQuda("ERROR: nSpin=%d is not supported\n", x.nSpin);           
    }
    blas_quda_bytes += 5*x.volume*sizeof(float);
  }
  blas_quda_bytes += 5*x.real_length*x.precision;
  blas_quda_flops += 8*x.real_length;

  if (!blasTuning) checkCudaError();
}

template <typename Float2>
__global__ void caxpbypzYmbwDKernel(Float2 a, Float2 *x, Float2 b, Float2 *y, Float2 *z, Float2 *w, int len) {

  unsigned int i = blockIdx.x*(blockDim.x) + threadIdx.x;
  unsigned int gridSize = gridDim.x*blockDim.x;
  while (i < len) {
    Float2 X = READ_DOUBLE2_TEXTURE(x, i);
    Float2 Z = read_Float2(z, i);

    Z.x += a.x*X.x - a.y*X.y;
    Z.y += a.y*X.x + a.x*X.y;

    Float2 Y = READ_DOUBLE2_TEXTURE(y, i);
    Z.x += b.x*Y.x - b.y*Y.y;
    Z.y += b.y*Y.x + b.x*Y.y;
    z[i] = make_Float2(Z);

    Float2 W = read_Float2(w, i);

    Y.x -= b.x*W.x - b.y*W.y;
    Y.y -= b.y*W.x + b.x*W.y;   
    
    y[i] = make_Float2(Y);
    i += gridSize;
  } 
}

template <typename Float2>
__global__ void caxpbypzYmbwSKernel(Float2 a, Float2 *x, Float2 b, Float2 *y, Float2 *z, Float2 *w, int len) {

  unsigned int i = blockIdx.x*(blockDim.x) + threadIdx.x;
  unsigned int gridSize = gridDim.x*blockDim.x;
  while (i < len) {
    Float2 X = read_Float2(x, i);
    Float2 Z = read_Float2(z, i);

    Z.x += a.x*X.x - a.y*X.y;
    Z.y += a.y*X.x + a.x*X.y;

    Float2 Y = read_Float2(y, i);
    Z.x += b.x*Y.x - b.y*Y.y;
    Z.y += b.y*Y.x + b.x*Y.y;
    z[i] = make_Float2(Z);

    Float2 W = read_Float2(w, i);

    Y.x -= b.x*W.x - b.y*W.y;
    Y.y -= b.y*W.x + b.x*W.y;   
    
    y[i] = make_Float2(Y);
    i += gridSize;
  } 
}

__global__ void caxpbypzYmbwHKernel(float2 a, float2 b, float *xN, short4 *yH, float *yN, 
                                    short4 *zH, float *zN, float *wN, int stride, int length) {
  
  unsigned int i = blockIdx.x*(blockDim.x) + threadIdx.x;
  unsigned int gridSize = gridDim.x*blockDim.x;
  while (i < length) {
    RECONSTRUCT_HALF_SPINOR(x, texHalf1, texNorm1, stride);
    RECONSTRUCT_HALF_SPINOR(y, texHalf2, texNorm2, stride);
    RECONSTRUCT_HALF_SPINOR(z, texHalf3, texNorm3, stride);
    CAXPBYPZ_FLOAT4(a, x0, b, y0, z0);
    CAXPBYPZ_FLOAT4(a, x1, b, y1, z1);
    CAXPBYPZ_FLOAT4(a, x2, b, y2, z2);
    CAXPBYPZ_FLOAT4(a, x3, b, y3, z3);
    CAXPBYPZ_FLOAT4(a, x4, b, y4, z4);
    CAXPBYPZ_FLOAT4(a, x5, b, y5, z5);
    CONSTRUCT_HALF_SPINOR_FROM_SINGLE(zH, zN, z, stride);
    READ_HALF_SPINOR(w, texHalf4, stride);
    float2 b2 = -wc*b; 
    CAXPY_FLOAT4(b2, w0, y0);
    CAXPY_FLOAT4(b2, w1, y1);
    CAXPY_FLOAT4(b2, w2, y2);
    CAXPY_FLOAT4(b2, w3, y3);
    CAXPY_FLOAT4(b2, w4, y4);
    CAXPY_FLOAT4(b2, w5, y5);
    CONSTRUCT_HALF_SPINOR_FROM_SINGLE(yH, yN, y, stride);
    i += gridSize;
  }   
}

__global__ void caxpbypzYmbwHKernel(float2 a, float2 b, float *xN, short2 *yH, float *yN, 
                                    short2 *zH, float *zN, float *wN, int stride, int length) {
  
  unsigned int i = blockIdx.x*(blockDim.x) + threadIdx.x;
  unsigned int gridSize = gridDim.x*blockDim.x;
  while (i < length) {
    RECONSTRUCT_HALF_SPINOR_ST(x, texHalfSt1, texNorm1, stride);
    RECONSTRUCT_HALF_SPINOR_ST(y, texHalfSt2, texNorm2, stride);
    RECONSTRUCT_HALF_SPINOR_ST(z, texHalfSt3, texNorm3, stride);
    CAXPBYPZ_FLOAT2(a, x0, b, y0, z0);
    CAXPBYPZ_FLOAT2(a, x1, b, y1, z1);
    CAXPBYPZ_FLOAT2(a, x2, b, y2, z2);
    CONSTRUCT_HALF_SPINOR_FROM_SINGLE_ST(zH, zN, z, stride);
    READ_HALF_SPINOR_ST(w, texHalfSt4, stride);
    float2 b2 = -wc*b; 
    CAXPY_FLOAT2(b2, w0, y0);
    CAXPY_FLOAT2(b2, w1, y1);
    CAXPY_FLOAT2(b2, w2, y2);
    CONSTRUCT_HALF_SPINOR_FROM_SINGLE_ST(yH, yN, y, stride);
    i += gridSize;
  }   
}

// performs the operation z[i] = a*x[i] + b*y[i] + z[i] and y[i] -= b*w[i]
void caxpbypzYmbwCuda(const Complex &a, cudaColorSpinorField &x, const Complex &b, cudaColorSpinorField &y,
                      cudaColorSpinorField &z, cudaColorSpinorField &w) {
  checkSpinor(x,y);
  checkSpinor(x,z);
  checkSpinor(x,w);
  int length = x.length/2;
  setBlock(13, length, x.precision);
  if (x.precision == QUDA_DOUBLE_PRECISION) {
    int spinor_bytes = x.length*sizeof(double);
    cudaBindTexture(0, xTexDouble2, x.v, spinor_bytes); 
    cudaBindTexture(0, yTexDouble2, y.v, spinor_bytes); 
    cudaBindTexture(0, zTexDouble2, z.v, spinor_bytes); 
    double2 a2 = make_double2(real(a), imag(a));
    double2 b2 = make_double2(real(b), imag(b));
    caxpbypzYmbwDKernel<<<blasGrid, blasBlock>>>(a2, (double2*)x.v, b2, (double2*)y.v, 
                                          (double2*)z.v, (double2*)w.v, length); 
  } else if (x.precision == QUDA_SINGLE_PRECISION) {
    float2 a2 = make_float2(real(a), imag(a));
    float2 b2 = make_float2(real(b), imag(b));
    caxpbypzYmbwSKernel<<<blasGrid, blasBlock>>>(a2, (float2*)x.v, b2, (float2*)y.v, 
                                                 (float2*)z.v, (float2*)w.v, length); 
  } else {
    if (x.siteSubset == QUDA_FULL_SITE_SUBSET) {
      caxpbypzYmbwCuda(a, x.Even(), b, y.Even(), z.Even(), w.Even());
      caxpbypzYmbwCuda(a, x.Odd(), b, y.Odd(), z.Odd(), w.Odd());
      return;
    }
    int spinor_bytes = x.length*sizeof(short);
    blas_quda_bytes += 6*x.volume*sizeof(float);
    float2 a2 = make_float2(real(a), imag(a));
    float2 b2 = make_float2(real(b), imag(b));
    if (x.nSpin == 4){ //wilson
      cudaBindTexture(0, texHalf1, x.v, spinor_bytes); 
      cudaBindTexture(0, texNorm1, x.norm, spinor_bytes/12);    
      cudaBindTexture(0, texHalf2, y.v, spinor_bytes); 
      cudaBindTexture(0, texNorm2, y.norm, spinor_bytes/12);    
      cudaBindTexture(0, texHalf3, z.v, spinor_bytes); 
      cudaBindTexture(0, texNorm3, z.norm, spinor_bytes/12);    
      cudaBindTexture(0, texHalf4, w.v, spinor_bytes); 
      cudaBindTexture(0, texNorm4, w.norm, spinor_bytes/12); 
      caxpbypzYmbwHKernel<<<blasGrid, blasBlock>>>(a2, b2, (float*)x.norm, (short4*)y.v, (float*)y.norm,
                                                   (short4*)z.v, (float*)z.norm, (float*)w.norm,
                                                   z.stride, z.volume);
    } else if (x.nSpin == 1){ //staggered
      cudaBindTexture(0, texHalfSt1, x.v, spinor_bytes); 
      cudaBindTexture(0, texNorm1, x.norm, spinor_bytes/3);    
      cudaBindTexture(0, texHalfSt2, y.v, spinor_bytes); 
      cudaBindTexture(0, texNorm2, y.norm, spinor_bytes/3);    
      cudaBindTexture(0, texHalfSt3, z.v, spinor_bytes); 
      cudaBindTexture(0, texNorm3, z.norm, spinor_bytes/3);    
      cudaBindTexture(0, texHalfSt4, w.v, spinor_bytes);
      cudaBindTexture(0, texNorm4, w.norm, spinor_bytes/3); 
      caxpbypzYmbwHKernel<<<blasGrid, blasBlock>>>(a2, b2, (float*)x.norm, (short2*)y.v, (float*)y.norm,
                                                   (short2*)z.v, (float*)z.norm, (float*)w.norm,
                                                   z.stride, z.volume);
    }else{
     errorQuda("ERROR: nSpin=%d is not supported\n", x.nSpin);                 
    }
  }
  blas_quda_bytes += 6*x.real_length*x.precision;
  blas_quda_flops += 12*x.real_length;

  if (!blasTuning) checkCudaError();
}


// Computes c = a + b in "double single" precision.
__device__ void dsadd(volatile QudaSumFloat &c0, volatile QudaSumFloat &c1, const volatile QudaSumFloat &a0, 
                      const volatile QudaSumFloat &a1, const float b0, const float b1) {
  // Compute dsa + dsb using Knuth's trick.
  QudaSumFloat t1 = a0 + b0;
  QudaSumFloat e = t1 - a0;
  QudaSumFloat t2 = ((b0 - e) + (a0 - (t1 - e))) + a1 + b1;
  // The result is t1 + t2, after normalization.
  c0 = e = t1 + t2;
  c1 = t2 - (e - t1);
}

// Computes c = a + b in "double single" precision (complex version)
__device__ void zcadd(volatile QudaSumComplex &c0, volatile QudaSumComplex &c1, const volatile QudaSumComplex &a0, 
                      const volatile QudaSumComplex &a1, const volatile QudaSumComplex &b0, const volatile QudaSumComplex &b1) {
  // Compute dsa + dsb using Knuth's trick.
  QudaSumFloat t1 = a0.x + b0.x;
  QudaSumFloat e = t1 - a0.x;
  QudaSumFloat t2 = ((b0.x - e) + (a0.x - (t1 - e))) + a1.x + b1.x;
  // The result is t1 + t2, after normalization.
  c0.x = e = t1 + t2;
  c1.x = t2 - (e - t1);
  
  // Compute dsa + dsb using Knuth's trick.
  t1 = a0.y + b0.y;
  e = t1 - a0.y;
  t2 = ((b0.y - e) + (a0.y - (t1 - e))) + a1.y + b1.y;
  // The result is t1 + t2, after normalization.
  c0.y = e = t1 + t2;
  c1.y = t2 - (e - t1);
}

// Computes c = a + b in "double single" precision (float3 version)
__device__ void dsadd3(volatile QudaSumFloat3 &c0, volatile QudaSumFloat3 &c1, const volatile QudaSumFloat3 &a0, 
                       const volatile QudaSumFloat3 &a1, const volatile QudaSumFloat3 &b0, const volatile QudaSumFloat3 &b1) {
  // Compute dsa + dsb using Knuth's trick.
  QudaSumFloat t1 = a0.x + b0.x;
  QudaSumFloat e = t1 - a0.x;
  QudaSumFloat t2 = ((b0.x - e) + (a0.x - (t1 - e))) + a1.x + b1.x;
  // The result is t1 + t2, after normalization.
  c0.x = e = t1 + t2;
  c1.x = t2 - (e - t1);
  
  // Compute dsa + dsb using Knuth's trick.
  t1 = a0.y + b0.y;
  e = t1 - a0.y;
  t2 = ((b0.y - e) + (a0.y - (t1 - e))) + a1.y + b1.y;
  // The result is t1 + t2, after normalization.
  c0.y = e = t1 + t2;
  c1.y = t2 - (e - t1);
  
  // Compute dsa + dsb using Knuth's trick.
  t1 = a0.z + b0.z;
  e = t1 - a0.z;
  t2 = ((b0.z - e) + (a0.z - (t1 - e))) + a1.z + b1.z;
  // The result is t1 + t2, after normalization.
  c0.z = e = t1 + t2;
  c1.z = t2 - (e - t1);
}

//
// double sumCuda(float *a, int n) {}
//
template <unsigned int reduce_threads, typename Float>
#define REDUCE_FUNC_NAME(suffix) sumD##suffix
#define REDUCE_TYPES Float *a
#define REDUCE_PARAMS a
#define REDUCE_AUXILIARY(i)
#define REDUCE_OPERATION(i) a[i]
#include "reduce_core.h"
#undef REDUCE_FUNC_NAME
#undef REDUCE_TYPES
#undef REDUCE_PARAMS
#undef REDUCE_AUXILIARY
#undef REDUCE_OPERATION

template <unsigned int reduce_threads, typename Float>
#define REDUCE_FUNC_NAME(suffix) sumS##suffix
#define REDUCE_TYPES Float *a
#define REDUCE_PARAMS a
#define REDUCE_AUXILIARY(i)
#define REDUCE_OPERATION(i) a[i].x + a[i].y
#include "reduce_core.h"
#undef REDUCE_FUNC_NAME
#undef REDUCE_TYPES
#undef REDUCE_PARAMS
#undef REDUCE_AUXILIARY
#undef REDUCE_OPERATION

template <unsigned int reduce_threads, typename Float>
#define REDUCE_FUNC_NAME(suffix) sumH##suffix
#define REDUCE_TYPES Float *aN, int stride
#define REDUCE_PARAMS aN, stride
#define REDUCE_AUXILIARY(i)                                             \
  READ_HALF_SPINOR(a, texHalf1, stride);                                \
  SUM_FLOAT4(s0, a0);                                                   \
  SUM_FLOAT4(s1, a1);                                                   \
  SUM_FLOAT4(s2, a2);                                                   \
  SUM_FLOAT4(s3, a3);                                                   \
  SUM_FLOAT4(s4, a4);                                                   \
  SUM_FLOAT4(s5, a5);                                                   \
  s0 += s1; s2 += s3; s4 += s5; s0 += s2; s0 += s4;
#define REDUCE_OPERATION(i) (ac*s0)
#include "reduce_core.h"
#undef REDUCE_FUNC_NAME
#undef REDUCE_TYPES
#undef REDUCE_PARAMS
#undef REDUCE_AUXILIARY
#undef REDUCE_OPERATION

template <unsigned int reduce_threads, typename Float>
#define REDUCE_FUNC_NAME(suffix) sumHSt##suffix
#define REDUCE_TYPES Float *aN, int stride
#define REDUCE_PARAMS aN, stride
#define REDUCE_AUXILIARY(i)                                             \
  READ_HALF_SPINOR_ST(a, texHalfSt1, stride);                           \
  SUM_FLOAT2(s0, a0);                                                   \
  SUM_FLOAT2(s1, a1);                                                   \
  SUM_FLOAT2(s2, a2);                                                   \
  s0 += s1; s0 += s2;
#define REDUCE_OPERATION(i) (ac*s0)
#include "reduce_core.h"
#undef REDUCE_FUNC_NAME
#undef REDUCE_TYPES
#undef REDUCE_PARAMS
#undef REDUCE_AUXILIARY
#undef REDUCE_OPERATION

double sumCuda(cudaColorSpinorField &a) {
  const int id = 14;
  blas_quda_flops += a.real_length;
  blas_quda_bytes += a.real_length*a.precision;
  if (a.precision == QUDA_DOUBLE_PRECISION) {
    return sumDCuda((double*)a.v, a.length, id, a.precision);
  } else if (a.precision == QUDA_SINGLE_PRECISION) {
    return sumSCuda((float2*)a.v, a.length/2, id, a.precision);
  } else {
    if (a.siteSubset == QUDA_FULL_SITE_SUBSET) return sumCuda(a.Even()) + sumCuda(a.Odd());
    int spinor_bytes = a.length*sizeof(short);
    blas_quda_bytes += a.volume*sizeof(float);
    if (a.nSpin ==4) { // wilson
      cudaBindTexture(0, texHalf1, a.v, spinor_bytes); 
      cudaBindTexture(0, texNorm1, a.norm, spinor_bytes/12);    
      return sumHCuda((float*)a.norm, a.stride, a.volume, id, a.precision);
    } else if (a.nSpin ==1) { // staggered
      cudaBindTexture(0, texHalfSt1, a.v, spinor_bytes); 
      cudaBindTexture(0, texNorm1, a.norm, spinor_bytes/3);    
      return sumHStCuda((float*)a.norm, a.stride, a.volume, id, a.precision);
    } else {
      errorQuda("ERROR: nSpin=%d is not supported\n", a.nSpin);           
      return 0;
    }
  }
}

//
// double normCuda(float *a, int n) {}
//
template <unsigned int reduce_threads, typename Float>
#define REDUCE_FUNC_NAME(suffix) normD##suffix
#define REDUCE_TYPES Float *a
#define REDUCE_PARAMS a
#define REDUCE_AUXILIARY(i)
#define REDUCE_OPERATION(i) (a[i]*a[i])
#include "reduce_core.h"
#undef REDUCE_FUNC_NAME
#undef REDUCE_TYPES
#undef REDUCE_PARAMS
#undef REDUCE_AUXILIARY
#undef REDUCE_OPERATION

template <unsigned int reduce_threads, typename Float>
#define REDUCE_FUNC_NAME(suffix) normS##suffix
#define REDUCE_TYPES Float *a
#define REDUCE_PARAMS a
#define REDUCE_AUXILIARY(i)
#define REDUCE_OPERATION(i) (a[i].x*a[i].x + a[i].y*a[i].y)
#include "reduce_core.h"
#undef REDUCE_FUNC_NAME
#undef REDUCE_TYPES
#undef REDUCE_PARAMS
#undef REDUCE_AUXILIARY
#undef REDUCE_OPERATION

//
// double normHCuda(char *, int n) {}
//
template <unsigned int reduce_threads, typename Float>
#define REDUCE_FUNC_NAME(suffix) normH##suffix
#define REDUCE_TYPES Float *aN, int stride // dummy type
#define REDUCE_PARAMS aN, stride
#define REDUCE_AUXILIARY(i)                                             \
  READ_HALF_SPINOR(a, texHalf1, stride);                                \
  REAL_DOT_FLOAT4(norm0, a0, a0);                                       \
  REAL_DOT_FLOAT4(norm1, a1, a1);                                       \
  REAL_DOT_FLOAT4(norm2, a2, a2);                                       \
  REAL_DOT_FLOAT4(norm3, a3, a3);                                       \
  REAL_DOT_FLOAT4(norm4, a4, a4);                                       \
  REAL_DOT_FLOAT4(norm5, a5, a5);                                       \
  norm0 += norm1; norm2 += norm3; norm4 += norm5; norm0 += norm2, norm0 += norm4;
#define REDUCE_OPERATION(i) (ac*ac*norm0)
#include "reduce_core.h"
#undef REDUCE_FUNC_NAME
#undef REDUCE_TYPES
#undef REDUCE_PARAMS
#undef REDUCE_AUXILIARY
#undef REDUCE_OPERATION

template <unsigned int reduce_threads, typename Float>
#define REDUCE_FUNC_NAME(suffix) normHSt##suffix
#define REDUCE_TYPES Float *aN, int stride // dummy type
#define REDUCE_PARAMS aN, stride
#define REDUCE_AUXILIARY(i)                                             \
    READ_HALF_SPINOR_ST(a, texHalfSt1, stride);                         \
    REAL_DOT_FLOAT2(norm0, a0, a0);                                     \
    REAL_DOT_FLOAT2(norm1, a1, a1);                                     \
    REAL_DOT_FLOAT2(norm2, a2, a2);                                     \
    norm0 += norm1; norm0 += norm2; 
#define REDUCE_OPERATION(i) (ac*ac*norm0)
#include "reduce_core.h"
#undef REDUCE_FUNC_NAME
#undef REDUCE_TYPES
#undef REDUCE_PARAMS
#undef REDUCE_AUXILIARY
#undef REDUCE_OPERATION

double normCuda(const cudaColorSpinorField &a) {
  const int id = 15;
  blas_quda_flops += 2*a.real_length;
  blas_quda_bytes += a.real_length*a.precision;
  if (a.precision == QUDA_DOUBLE_PRECISION) {
    return normDCuda((double*)a.v, a.length, id, a.precision);
  } else if (a.precision == QUDA_SINGLE_PRECISION) {
    return normSCuda((float2*)a.v, a.length/2, id, a.precision);
  } else {
    if (a.siteSubset == QUDA_FULL_SITE_SUBSET) return normCuda(a.Even()) + normCuda(a.Odd());
    int spinor_bytes = a.length*sizeof(short);
    int half_norm_ratio = (a.nColor*a.nSpin*2*sizeof(short))/sizeof(float);
    blas_quda_bytes += (a.real_length*a.precision) / (a.nColor * a.nSpin);
    cudaBindTexture(0, texNorm1, a.norm, spinor_bytes/half_norm_ratio);    
    if (a.nSpin == 4){ //wilson
        cudaBindTexture(0, texHalf1, a.v, spinor_bytes); 
        return normHCuda((float*)a.norm, a.stride, a.volume, id, a.precision);
    }else if (a.nSpin == 1) { //staggered
        cudaBindTexture(0, texHalfSt1, a.v, spinor_bytes);
        return normHStCuda((float*)a.norm, a.stride, a.volume, id, a.precision);        
    }else{
      errorQuda("%s: nSpin(%d) is not supported\n", __FUNCTION__, a.nSpin);
      return 0;
    }
  }

}

//
// double reDotProductFCuda(float *a, float *b, int n) {}
//
template <unsigned int reduce_threads, typename Float>
#define REDUCE_FUNC_NAME(suffix) reDotProductD##suffix
#define REDUCE_TYPES Float *a, Float *b
#define REDUCE_PARAMS a, b
#define REDUCE_AUXILIARY(i)
#define REDUCE_OPERATION(i) (a[i]*b[i])
#include "reduce_core.h"
#undef REDUCE_FUNC_NAME
#undef REDUCE_TYPES
#undef REDUCE_PARAMS
#undef REDUCE_AUXILIARY
#undef REDUCE_OPERATION

template <unsigned int reduce_threads, typename Float>
#define REDUCE_FUNC_NAME(suffix) reDotProductS##suffix
#define REDUCE_TYPES Float *a, Float *b
#define REDUCE_PARAMS a, b
#define REDUCE_AUXILIARY(i)
#define REDUCE_OPERATION(i) (a[i].x*b[i].x + a[i].y*b[i].y)
#include "reduce_core.h"
#undef REDUCE_FUNC_NAME
#undef REDUCE_TYPES
#undef REDUCE_PARAMS
#undef REDUCE_AUXILIARY
#undef REDUCE_OPERATION

//
// double reDotProductHCuda(float *a, float *b, int n) {}
//
template <unsigned int reduce_threads, typename Float>
#define REDUCE_FUNC_NAME(suffix) reDotProductH##suffix
#define REDUCE_TYPES Float *aN, Float *bN, int stride
#define REDUCE_PARAMS aN, bN, stride
#define REDUCE_AUXILIARY(i)                                             \
  READ_HALF_SPINOR(a, texHalf1, stride);                                \
  READ_HALF_SPINOR(b, texHalf2, stride);                                \
  REAL_DOT_FLOAT4(rdot0, a0, b0);                                       \
  REAL_DOT_FLOAT4(rdot1, a1, b1);                                       \
  REAL_DOT_FLOAT4(rdot2, a2, b2);                                       \
  REAL_DOT_FLOAT4(rdot3, a3, b3);                                       \
  REAL_DOT_FLOAT4(rdot4, a4, b4);                                       \
  REAL_DOT_FLOAT4(rdot5, a5, b5);                                       \
  rdot0 += rdot1; rdot2 += rdot3; rdot4 += rdot5; rdot0 += rdot2; rdot0 += rdot4;
#define REDUCE_OPERATION(i) (ac*bc*rdot0)
#include "reduce_core.h"
#undef REDUCE_FUNC_NAME
#undef REDUCE_TYPES
#undef REDUCE_PARAMS
#undef REDUCE_AUXILIARY
#undef REDUCE_OPERATION

template <unsigned int reduce_threads, typename Float>
#define REDUCE_FUNC_NAME(suffix) reDotProductHSt##suffix
#define REDUCE_TYPES Float *aN, Float *bN, int stride
#define REDUCE_PARAMS aN, bN, stride
#define REDUCE_AUXILIARY(i)                                             \
  READ_HALF_SPINOR_ST(a, texHalfSt1, stride);                           \
  READ_HALF_SPINOR_ST(b, texHalfSt2, stride);                           \
  REAL_DOT_FLOAT2(rdot0, a0, b0);                                       \
  REAL_DOT_FLOAT2(rdot1, a1, b1);                                       \
  REAL_DOT_FLOAT2(rdot2, a2, b2);                                       \
  rdot0 += rdot1; rdot0 += rdot2; 
#define REDUCE_OPERATION(i) (ac*bc*rdot0)
#include "reduce_core.h"
#undef REDUCE_FUNC_NAME
#undef REDUCE_TYPES
#undef REDUCE_PARAMS
#undef REDUCE_AUXILIARY
#undef REDUCE_OPERATION  

double reDotProductCuda(cudaColorSpinorField &a, cudaColorSpinorField &b) {
  const int id = 16;
  blas_quda_flops += 2*a.real_length;
  checkSpinor(a, b);
  blas_quda_bytes += 2*a.real_length*a.precision;
  if (a.precision == QUDA_DOUBLE_PRECISION) {
    return reDotProductDCuda((double*)a.v, (double*)b.v, a.length, id, a.precision);
  } else if (a.precision == QUDA_SINGLE_PRECISION) {
    return reDotProductSCuda((float2*)a.v, (float2*)b.v, a.length/2, id, a.precision);
  } else {
    if (a.siteSubset == QUDA_FULL_SITE_SUBSET) {
      return reDotProductCuda(a.Even(), b.Even()) + reDotProductCuda(a.Odd(), b.Odd());
    }
    blas_quda_bytes += 2*a.volume*sizeof(float);
    int spinor_bytes = a.length*sizeof(short);
    if (a.nSpin == 4){ //wilson
      cudaBindTexture(0, texHalf1, a.v, spinor_bytes); 
      cudaBindTexture(0, texNorm1, a.norm, spinor_bytes/12);    
      cudaBindTexture(0, texHalf2, b.v, spinor_bytes); 
      cudaBindTexture(0, texNorm2, b.norm, spinor_bytes/12);    
      return reDotProductHCuda((float*)a.norm, (float*)b.norm, a.stride, a.volume, id, a.precision);
    }else if (a.nSpin == 1){ //staggered
      cudaBindTexture(0, texHalfSt1, a.v, spinor_bytes); 
      cudaBindTexture(0, texNorm1, a.norm, spinor_bytes/3);    
      cudaBindTexture(0, texHalfSt2, b.v, spinor_bytes); 
      cudaBindTexture(0, texNorm2, b.norm, spinor_bytes/3);    
      return reDotProductHStCuda((float*)a.norm, (float*)b.norm, a.stride, a.volume, id, a.precision);
    }else{
      errorQuda("%s: nSpin(%d) is not supported\n", __FUNCTION__, a.nSpin);      
      return 0;
    }
  }

}

//
// double axpyNormCuda(float a, float *x, float *y, n){}
//
// First performs the operation y[i] = a*x[i] + y[i]
// Second returns the norm of y
//

template <unsigned int reduce_threads, typename Float>
#define REDUCE_FUNC_NAME(suffix) axpyNormF##suffix
#define REDUCE_TYPES Float a, Float *x, Float *y
#define REDUCE_PARAMS a, x, y
#define REDUCE_AUXILIARY(i) y[i] = a*x[i] + y[i]
#define REDUCE_OPERATION(i) (y[i]*y[i])
#include "reduce_core.h"
#undef REDUCE_FUNC_NAME
#undef REDUCE_TYPES
#undef REDUCE_PARAMS
#undef REDUCE_AUXILIARY
#undef REDUCE_OPERATION

template <unsigned int reduce_threads, typename Float>
#define REDUCE_FUNC_NAME(suffix) axpyNormH##suffix
#define REDUCE_TYPES Float a, short4 *yH, float *yN, int stride
#define REDUCE_PARAMS a, yH, yN, stride
#define REDUCE_AUXILIARY(i)                                             \
  RECONSTRUCT_HALF_SPINOR(x, texHalf1, texNorm1, stride);               \
  RECONSTRUCT_HALF_SPINOR(y, texHalf2, texNorm2, stride);               \
  AXPY_FLOAT4(a, x0, y0);                                               \
  REAL_DOT_FLOAT4(norm0, y0, y0);                                       \
  AXPY_FLOAT4(a, x1, y1);                                               \
  REAL_DOT_FLOAT4(norm1, y1, y1);                                       \
  AXPY_FLOAT4(a, x2, y2);                                               \
  REAL_DOT_FLOAT4(norm2, y2, y2);                                       \
  AXPY_FLOAT4(a, x3, y3);                                               \
  REAL_DOT_FLOAT4(norm3, y3, y3);                                       \
  AXPY_FLOAT4(a, x4, y4);                                               \
  REAL_DOT_FLOAT4(norm4, y4, y4);                                       \
  AXPY_FLOAT4(a, x5, y5);                                               \
  REAL_DOT_FLOAT4(norm5, y5, y5);                                       \
  norm0 += norm1; norm2 += norm3; norm4 += norm5; norm0 += norm2; norm0 += norm4; \
  CONSTRUCT_HALF_SPINOR_FROM_SINGLE(yH, yN, y, stride);
#define REDUCE_OPERATION(i) (norm0)
#include "reduce_core.h"
#undef REDUCE_FUNC_NAME
#undef REDUCE_TYPES
#undef REDUCE_PARAMS
#undef REDUCE_AUXILIARY
#undef REDUCE_OPERATION

template <unsigned int reduce_threads, typename Float>
#define REDUCE_FUNC_NAME(suffix) axpyNormH##suffix
#define REDUCE_TYPES Float a, short2 *yH, float *yN, int stride
#define REDUCE_PARAMS a, yH, yN, stride
#define REDUCE_AUXILIARY(i)                                             \
  RECONSTRUCT_HALF_SPINOR_ST(x, texHalfSt1, texNorm1, stride);          \
  RECONSTRUCT_HALF_SPINOR_ST(y, texHalfSt2, texNorm2, stride);          \
  AXPY_FLOAT2(a, x0, y0);                                               \
  REAL_DOT_FLOAT2(norm0, y0, y0);                                       \
  AXPY_FLOAT2(a, x1, y1);                                               \
  REAL_DOT_FLOAT2(norm1, y1, y1);                                       \
  AXPY_FLOAT2(a, x2, y2);                                               \
  REAL_DOT_FLOAT2(norm2, y2, y2);                                       \
  norm0 += norm1; norm0 += norm2;                                       \
  CONSTRUCT_HALF_SPINOR_FROM_SINGLE_ST(yH, yN, y, stride);
#define REDUCE_OPERATION(i) (norm0)
#include "reduce_core.h"
#undef REDUCE_FUNC_NAME
#undef REDUCE_TYPES
#undef REDUCE_PARAMS
#undef REDUCE_AUXILIARY
#undef REDUCE_OPERATION

double axpyNormCuda(const double &a, cudaColorSpinorField &x, cudaColorSpinorField &y) {
  const int id = 17;
  blas_quda_flops += 4*x.real_length;
  checkSpinor(x,y);
  blas_quda_bytes += 3*x.real_length*x.precision;
  if (x.precision == QUDA_DOUBLE_PRECISION) {
    return axpyNormFCuda(a, (double*)x.v, (double*)y.v, x.length, id, x.precision);
  } else if (x.precision == QUDA_SINGLE_PRECISION) {
    return axpyNormFCuda((float)a, (float*)x.v, (float*)y.v, x.length, id, x.precision);
  } else {
    if (x.siteSubset == QUDA_FULL_SITE_SUBSET) 
      return axpyNormCuda(a, x.Even(), y.Even()) + axpyNormCuda(a, x.Odd(), y.Odd());

    cudaBindTexture(0, texNorm1, x.norm, x.bytes/(x.nColor*x.nSpin));    
    cudaBindTexture(0, texNorm2, y.norm, x.bytes/(x.nColor*x.nSpin));    
    blas_quda_bytes += 3*x.volume*sizeof(float);
    
    if (x.nSpin == 4){ //wilson
      cudaBindTexture(0, texHalf1, x.v, x.bytes); 
      cudaBindTexture(0, texHalf2, y.v, x.bytes); 
      return axpyNormHCuda((float)a, (short4*)y.v, (float*)y.norm, x.stride, x.volume, id, x.precision);
    }else if (x.nSpin == 1){ //staggered
      cudaBindTexture(0, texHalfSt1, x.v, x.bytes); 
      cudaBindTexture(0, texHalfSt2, y.v, x.bytes); 
      return axpyNormHCuda((float)a, (short2*)y.v, (float*)y.norm, x.stride, x.volume, id, x.precision);
    }else{
      errorQuda("%s: nSpin(%d) is not supported\n", __FUNCTION__, x.nSpin);            
      return 0;
    }
  }

}

//
// double xmyNormCuda(float a, float *x, float *y, n){}
//
// First performs the operation y[i] = x[i] - y[i]
// Second returns the norm of y
//

template <unsigned int reduce_threads, typename Float>
#define REDUCE_FUNC_NAME(suffix) xmyNormF##suffix
#define REDUCE_TYPES Float *x, Float *y
#define REDUCE_PARAMS x, y
#define REDUCE_AUXILIARY(i) y[i] = x[i] - y[i]
#define REDUCE_OPERATION(i) (y[i]*y[i])
#include "reduce_core.h"
#undef REDUCE_FUNC_NAME
#undef REDUCE_TYPES
#undef REDUCE_PARAMS
#undef REDUCE_AUXILIARY
#undef REDUCE_OPERATION

template <unsigned int reduce_threads, typename Float>
#define REDUCE_FUNC_NAME(suffix) xmyNormH##suffix
#define REDUCE_TYPES Float *d1, Float *d2, short4 *yH, float *yN, int stride
#define REDUCE_PARAMS d1, d2, yH, yN, stride
#define REDUCE_AUXILIARY(i)                                             \
  RECONSTRUCT_HALF_SPINOR(x, texHalf1, texNorm1, stride);               \
  RECONSTRUCT_HALF_SPINOR(y, texHalf2, texNorm2, stride);               \
  XMY_FLOAT4(x0, y0);                                                   \
  REAL_DOT_FLOAT4(norm0, y0, y0);                                       \
  XMY_FLOAT4(x1, y1);                                                   \
  REAL_DOT_FLOAT4(norm1, y1, y1);                                       \
  XMY_FLOAT4(x2, y2);                                                   \
  REAL_DOT_FLOAT4(norm2, y2, y2);                                       \
  XMY_FLOAT4(x3, y3);                                                   \
  REAL_DOT_FLOAT4(norm3, y3, y3);                                       \
  XMY_FLOAT4(x4, y4);                                                   \
  REAL_DOT_FLOAT4(norm4, y4, y4);                                       \
  XMY_FLOAT4(x5, y5);                                                   \
  REAL_DOT_FLOAT4(norm5, y5, y5);                                       \
  norm0 += norm1; norm2 += norm3; norm4 += norm5; norm0 += norm2; norm0 += norm4; \
  CONSTRUCT_HALF_SPINOR_FROM_SINGLE(yH, yN, y, stride);
#define REDUCE_OPERATION(i) (norm0)
#include "reduce_core.h"
#undef REDUCE_FUNC_NAME
#undef REDUCE_TYPES
#undef REDUCE_PARAMS
#undef REDUCE_AUXILIARY
#undef REDUCE_OPERATION

template <unsigned int reduce_threads, typename Float>
#define REDUCE_FUNC_NAME(suffix) xmyNormH##suffix
#define REDUCE_TYPES Float *d1, Float *d2, short2 *yH, float *yN, int stride
#define REDUCE_PARAMS d1, d2, yH, yN, stride
#define REDUCE_AUXILIARY(i)                                             \
  RECONSTRUCT_HALF_SPINOR_ST(x, texHalfSt1, texNorm1, stride);          \
  RECONSTRUCT_HALF_SPINOR_ST(y, texHalfSt2, texNorm2, stride);          \
  XMY_FLOAT2(x0, y0);                                                   \
  REAL_DOT_FLOAT2(norm0, y0, y0);                                       \
  XMY_FLOAT2(x1, y1);                                                   \
  REAL_DOT_FLOAT2(norm1, y1, y1);                                       \
  XMY_FLOAT2(x2, y2);                                                   \
  REAL_DOT_FLOAT2(norm2, y2, y2);                                       \
  norm0 += norm1; norm0 += norm2;                                       \
  CONSTRUCT_HALF_SPINOR_FROM_SINGLE_ST(yH, yN, y, stride);
#define REDUCE_OPERATION(i) (norm0)
#include "reduce_core.h"
#undef REDUCE_FUNC_NAME
#undef REDUCE_TYPES
#undef REDUCE_PARAMS
#undef REDUCE_AUXILIARY
#undef REDUCE_OPERATION


double xmyNormCuda(cudaColorSpinorField &x, cudaColorSpinorField &y) {
  const int id = 18;
  blas_quda_flops += 3*x.real_length;
  checkSpinor(x,y);
  blas_quda_bytes += 3*x.real_length*x.precision;
  if (x.precision == QUDA_DOUBLE_PRECISION) {
    return xmyNormFCuda((double*)x.v, (double*)y.v, x.length, id, x.precision);
  } else if (x.precision == QUDA_SINGLE_PRECISION) {
    return xmyNormFCuda((float*)x.v, (float*)y.v, x.length, id, x.precision);
  } else { 
    if (x.siteSubset == QUDA_FULL_SITE_SUBSET) 
      return xmyNormCuda(x.Even(), y.Even()) + xmyNormCuda(x.Odd(), y.Odd());

    cudaBindTexture(0, texNorm1, x.norm, x.bytes/(x.nColor*x.nSpin));    
    cudaBindTexture(0, texNorm2, y.norm, x.bytes/(x.nColor*x.nSpin));    
    blas_quda_bytes += 3*x.volume*sizeof(float);
    
    if (x.nSpin ==4 ){ //wilsin
      cudaBindTexture(0, texHalf1, x.v, x.bytes); 
      cudaBindTexture(0, texHalf2, y.v, x.bytes); 
      return xmyNormHCuda((char*)0, (char*)0, (short4*)y.v, (float*)y.norm, y.stride, y.volume, id, x.precision);
    }else if (x.nSpin == 1){
      cudaBindTexture(0, texHalfSt1, x.v, x.bytes); 
      cudaBindTexture(0, texHalfSt2, y.v, x.bytes); 
      return xmyNormHCuda((char*)0, (char*)0, (short2*)y.v, (float*)y.norm, y.stride, y.volume, id, x.precision);
    }else{
      errorQuda("%s: nSpin(%d) is not supported\n", __FUNCTION__, x.nSpin);                  
    }
  }

}


//
// double2 cDotProductCuda(float2 *x, float2 *y, int n) {}
//
template <unsigned int reduce_threads, typename Float, typename Float2>
#define REDUCE_FUNC_NAME(suffix) cDotProductD##suffix
#define REDUCE_TYPES Float2 *x, Float2 *y, Float c
#define REDUCE_PARAMS x, y, c
#define REDUCE_REAL_AUXILIARY(i) Float2 a = READ_DOUBLE2_TEXTURE(x, i);
#define REDUCE_IMAG_AUXILIARY(i) Float2 b = READ_DOUBLE2_TEXTURE(y, i);
#define REDUCE_REAL_OPERATION(i) (a.x*b.x + a.y*b.y)
#define REDUCE_IMAG_OPERATION(i) (a.x*b.y - a.y*b.x)
#include "reduce_complex_core.h"
#undef REDUCE_FUNC_NAME
#undef REDUCE_TYPES
#undef REDUCE_PARAMS
#undef REDUCE_REAL_AUXILIARY
#undef REDUCE_IMAG_AUXILIARY
#undef REDUCE_REAL_OPERATION
#undef REDUCE_IMAG_OPERATION

template <unsigned int reduce_threads, typename Float, typename Float2>
#define REDUCE_FUNC_NAME(suffix) cDotProductS##suffix
#define REDUCE_TYPES Float2 *x, Float2 *y, Float c
#define REDUCE_PARAMS x, y, c
#define REDUCE_REAL_AUXILIARY(i) Float2 a = read_Float2(x, i);
#define REDUCE_IMAG_AUXILIARY(i) Float2 b = read_Float2(y, i);
#define REDUCE_REAL_OPERATION(i) (a.x*b.x + a.y*b.y)
#define REDUCE_IMAG_OPERATION(i) (a.x*b.y - a.y*b.x)
#include "reduce_complex_core.h"
#undef REDUCE_FUNC_NAME
#undef REDUCE_TYPES
#undef REDUCE_PARAMS
#undef REDUCE_REAL_AUXILIARY
#undef REDUCE_IMAG_AUXILIARY
#undef REDUCE_REAL_OPERATION
#undef REDUCE_IMAG_OPERATION

template <unsigned int reduce_threads, typename Float, typename Float2>
#define REDUCE_FUNC_NAME(suffix) cDotProductH##suffix
#define REDUCE_TYPES Float *aN, Float2 *bN, int stride
#define REDUCE_PARAMS aN, bN, stride
#define REDUCE_REAL_AUXILIARY(i)                                        \
  READ_HALF_SPINOR(a, texHalf1, stride);                                \
  READ_HALF_SPINOR(b, texHalf2, stride);                                \
  REAL_DOT_FLOAT4(rdot0, a0, b0);                                       \
  REAL_DOT_FLOAT4(rdot1, a1, b1);                                       \
  REAL_DOT_FLOAT4(rdot2, a2, b2);                                       \
  REAL_DOT_FLOAT4(rdot3, a3, b3);                                       \
  REAL_DOT_FLOAT4(rdot4, a4, b4);                                       \
  REAL_DOT_FLOAT4(rdot5, a5, b5);                                       \
  rdot0 += rdot1; rdot2 += rdot3; rdot4 += rdot5; rdot0 += rdot2; rdot0 += rdot4;
#define REDUCE_IMAG_AUXILIARY(i)                                        \
  IMAG_DOT_FLOAT4(idot0, a0, b0);                                       \
  IMAG_DOT_FLOAT4(idot1, a1, b1);                                       \
  IMAG_DOT_FLOAT4(idot2, a2, b2);                                       \
  IMAG_DOT_FLOAT4(idot3, a3, b3);                                       \
  IMAG_DOT_FLOAT4(idot4, a4, b4);                                       \
  IMAG_DOT_FLOAT4(idot5, a5, b5);                                       \
  idot0 += idot1; idot2 += idot3; idot4 += idot5; idot0 += idot2; idot0 += idot4;
#define REDUCE_REAL_OPERATION(i) (ac*bc*rdot0)
#define REDUCE_IMAG_OPERATION(i) (ac*bc*idot0)
#include "reduce_complex_core.h"
#undef REDUCE_FUNC_NAME
#undef REDUCE_TYPES
#undef REDUCE_PARAMS
#undef REDUCE_REAL_AUXILIARY
#undef REDUCE_IMAG_AUXILIARY
#undef REDUCE_REAL_OPERATION
#undef REDUCE_IMAG_OPERATION

template <unsigned int reduce_threads, typename Float, typename Float2>
#define REDUCE_FUNC_NAME(suffix) cDotProductHSt##suffix
#define REDUCE_TYPES Float *aN, Float2 *bN, int stride
#define REDUCE_PARAMS aN, bN, stride
#define REDUCE_REAL_AUXILIARY(i)                                        \
  READ_HALF_SPINOR_ST(a, texHalfSt1, stride);                           \
  READ_HALF_SPINOR_ST(b, texHalfSt2, stride);                           \
  REAL_DOT_FLOAT2(rdot0, a0, b0);                                       \
  REAL_DOT_FLOAT2(rdot1, a1, b1);                                       \
  REAL_DOT_FLOAT2(rdot2, a2, b2);                                       \
  rdot0 += rdot1; rdot0 += rdot2;
#define REDUCE_IMAG_AUXILIARY(i)                                        \
  IMAG_DOT_FLOAT2(idot0, a0, b0);                                       \
  IMAG_DOT_FLOAT2(idot1, a1, b1);                                       \
  IMAG_DOT_FLOAT2(idot2, a2, b2);                                       \
  idot0 += idot1; idot0 += idot2;
#define REDUCE_REAL_OPERATION(i) (ac*bc*rdot0)
#define REDUCE_IMAG_OPERATION(i) (ac*bc*idot0)
#include "reduce_complex_core.h"
#undef REDUCE_FUNC_NAME
#undef REDUCE_TYPES
#undef REDUCE_PARAMS
#undef REDUCE_REAL_AUXILIARY
#undef REDUCE_IMAG_AUXILIARY
#undef REDUCE_REAL_OPERATION
#undef REDUCE_IMAG_OPERATION

Complex cDotProductCuda(cudaColorSpinorField &x, cudaColorSpinorField &y) {
  const int id = 19;
  blas_quda_flops += 4*x.real_length;
  checkSpinor(x,y);
  int length = x.length/2;
  blas_quda_bytes += 2*x.real_length*x.precision;
  double2 dot;
  if (x.precision == QUDA_DOUBLE_PRECISION) {
    char c = 0;
    int spinor_bytes = x.length*sizeof(double);
    cudaBindTexture(0, xTexDouble2, x.v, spinor_bytes); 
    cudaBindTexture(0, yTexDouble2, y.v, spinor_bytes); 
    dot = cDotProductDCuda((double2*)x.v, (double2*)y.v, c, length, id, x.precision);
  } else if (x.precision == QUDA_SINGLE_PRECISION) {
    char c = 0;
    int spinor_bytes = x.length*sizeof(float);
    cudaBindTexture(0, xTexSingle2, x.v, spinor_bytes); 
    cudaBindTexture(0, yTexSingle2, y.v, spinor_bytes); 
    dot = cDotProductSCuda((float2*)x.v, (float2*)y.v, c, length, id, x.precision);
  } else {
    if (x.siteSubset == QUDA_FULL_SITE_SUBSET) 
      return cDotProductCuda(x.Even(), y.Even()) + cDotProductCuda(x.Odd(), y.Odd());
    int spinor_bytes = x.length*sizeof(short);
    blas_quda_bytes += 2*x.volume*sizeof(float);
    if (x.nSpin == 4){
      cudaBindTexture(0, texHalf1, x.v, spinor_bytes); 
      cudaBindTexture(0, texNorm1, x.norm, spinor_bytes/12);    
      cudaBindTexture(0, texHalf2, y.v, spinor_bytes); 
      cudaBindTexture(0, texNorm2, y.norm, spinor_bytes/12);    
      dot = cDotProductHCuda((float*)x.norm, (float*)y.norm, x.stride, x.volume, id, x.precision);
    } else if (x.nSpin == 1){
      cudaBindTexture(0, texHalfSt1, x.v, spinor_bytes); 
      cudaBindTexture(0, texNorm1, x.norm, spinor_bytes/3);    
      cudaBindTexture(0, texHalfSt2, y.v, spinor_bytes); 
      cudaBindTexture(0, texNorm2, y.norm, spinor_bytes/3);    
      dot = cDotProductHStCuda((float*)x.norm, (float*)y.norm, x.stride, x.volume, id, x.precision);
    }else{
      errorQuda("%s: nSpin(%d) is not supported\n", __FUNCTION__, x.nSpin);      
    } 
  }

  return Complex(dot.x, dot.y);
}

//
// double2 xpaycDotzyCuda(float2 *x, float a, float2 *y, float2 *z, int n) {}
//
// First performs the operation y = x + a*y
// Second returns complex dot product (z,y)
//

template <unsigned int reduce_threads, typename Float, typename Float2>
#define REDUCE_FUNC_NAME(suffix) xpaycDotzyD##suffix
#define REDUCE_TYPES Float2 *x, Float a, Float2 *y, Float2 *z
#define REDUCE_PARAMS x, a, y, z
#define REDUCE_REAL_AUXILIARY(i)                \
  Float2 X = READ_DOUBLE2_TEXTURE(x, i);        \
  Float2 Y = READ_DOUBLE2_TEXTURE(y, i);        \
  Float2 Z = READ_DOUBLE2_TEXTURE(z, i);
#define REDUCE_IMAG_AUXILIARY(i) y[i].x = X.x + a*Y.x; y[i].y = X.y + a*Y.y
#define REDUCE_REAL_OPERATION(i) (Z.x*y[i].x + Z.y*y[i].y)
#define REDUCE_IMAG_OPERATION(i) (Z.x*y[i].y - Z.y*y[i].x)
#include "reduce_complex_core.h"
#undef REDUCE_FUNC_NAME
#undef REDUCE_TYPES
#undef REDUCE_PARAMS
#undef REDUCE_REAL_AUXILIARY
#undef REDUCE_IMAG_AUXILIARY
#undef REDUCE_REAL_OPERATION
#undef REDUCE_IMAG_OPERATION

template <unsigned int reduce_threads, typename Float, typename Float2>
#define REDUCE_FUNC_NAME(suffix) xpaycDotzyS##suffix
#define REDUCE_TYPES Float2 *x, Float a, Float2 *y, Float2 *z
#define REDUCE_PARAMS x, a, y, z
#define REDUCE_REAL_AUXILIARY(i) y[i].x = x[i].x + a*y[i].x
#define REDUCE_IMAG_AUXILIARY(i) y[i].y = x[i].y + a*y[i].y
#define REDUCE_REAL_OPERATION(i) (z[i].x*y[i].x + z[i].y*y[i].y)
#define REDUCE_IMAG_OPERATION(i) (z[i].x*y[i].y - z[i].y*y[i].x)
#include "reduce_complex_core.h"
#undef REDUCE_FUNC_NAME
#undef REDUCE_TYPES
#undef REDUCE_PARAMS
#undef REDUCE_REAL_AUXILIARY
#undef REDUCE_IMAG_AUXILIARY
#undef REDUCE_REAL_OPERATION
#undef REDUCE_IMAG_OPERATION

template <unsigned int reduce_threads, typename Float, typename Float2>
#define REDUCE_FUNC_NAME(suffix) xpaycDotzyH##suffix
#define REDUCE_TYPES Float a, short4 *yH, Float2 *yN, int stride
#define REDUCE_PARAMS a, yH, yN, stride
#define REDUCE_REAL_AUXILIARY(i)                                        \
  RECONSTRUCT_HALF_SPINOR(x, texHalf1, texNorm1, stride);               \
  RECONSTRUCT_HALF_SPINOR(y, texHalf2, texNorm2, stride);               \
  RECONSTRUCT_HALF_SPINOR(z, texHalf3, texNorm3, stride);               \
  XPAY_FLOAT4(x0, a, y0);                                               \
  XPAY_FLOAT4(x1, a, y1);                                               \
  XPAY_FLOAT4(x2, a, y2);                                               \
  XPAY_FLOAT4(x3, a, y3);                                               \
  XPAY_FLOAT4(x4, a, y4);                                               \
  XPAY_FLOAT4(x5, a, y5);                                               \
  REAL_DOT_FLOAT4(rdot0, z0, y0);                                       \
  REAL_DOT_FLOAT4(rdot1, z1, y1);                                       \
  REAL_DOT_FLOAT4(rdot2, z2, y2);                                       \
  REAL_DOT_FLOAT4(rdot3, z3, y3);                                       \
  REAL_DOT_FLOAT4(rdot4, z4, y4);                                       \
  REAL_DOT_FLOAT4(rdot5, z5, y5);                                       \
  rdot0 += rdot1; rdot2 += rdot3; rdot4 += rdot5; rdot0 += rdot2; rdot0 += rdot4;
#define REDUCE_IMAG_AUXILIARY(i)                                        \
  IMAG_DOT_FLOAT4(idot0, z0, y0);                                       \
  IMAG_DOT_FLOAT4(idot1, z1, y1);                                       \
  IMAG_DOT_FLOAT4(idot2, z2, y2);                                       \
  IMAG_DOT_FLOAT4(idot3, z3, y3);                                       \
  IMAG_DOT_FLOAT4(idot4, z4, y4);                                       \
  IMAG_DOT_FLOAT4(idot5, z5, y5);                                       \
  idot0 += idot1; idot2 += idot3; idot4 += idot5; idot0 += idot2; idot0 += idot4; \
  CONSTRUCT_HALF_SPINOR_FROM_SINGLE(yH, yN, y, stride);
#define REDUCE_REAL_OPERATION(i) (rdot0)
#define REDUCE_IMAG_OPERATION(i) (idot0)
#include "reduce_complex_core.h"
#undef REDUCE_FUNC_NAME
#undef REDUCE_TYPES
#undef REDUCE_PARAMS
#undef REDUCE_REAL_AUXILIARY
#undef REDUCE_IMAG_AUXILIARY
#undef REDUCE_REAL_OPERATION
#undef REDUCE_IMAG_OPERATION

template <unsigned int reduce_threads, typename Float, typename Float2>
#define REDUCE_FUNC_NAME(suffix) xpaycDotzyH##suffix
#define REDUCE_TYPES Float a, short2 *yH, Float2 *yN, int stride
#define REDUCE_PARAMS a, yH, yN, stride
#define REDUCE_REAL_AUXILIARY(i)                                        \
  RECONSTRUCT_HALF_SPINOR_ST(x, texHalfSt1, texNorm1, stride);          \
  RECONSTRUCT_HALF_SPINOR_ST(y, texHalfSt2, texNorm2, stride);          \
  RECONSTRUCT_HALF_SPINOR_ST(z, texHalfSt3, texNorm3, stride);          \
  XPAY_FLOAT2(x0, a, y0);                                               \
  XPAY_FLOAT2(x1, a, y1);                                               \
  XPAY_FLOAT2(x2, a, y2);                                               \
  REAL_DOT_FLOAT2(rdot0, z0, y0);                                       \
  REAL_DOT_FLOAT2(rdot1, z1, y1);                                       \
  REAL_DOT_FLOAT2(rdot2, z2, y2);                                       \
  rdot0 += rdot1; rdot0 += rdot2;
#define REDUCE_IMAG_AUXILIARY(i)                                        \
  IMAG_DOT_FLOAT2(idot0, z0, y0);                                       \
  IMAG_DOT_FLOAT2(idot1, z1, y1);                                       \
  IMAG_DOT_FLOAT2(idot2, z2, y2);                                       \
  idot0 += idot1; idot0 += idot2;                                       \
  CONSTRUCT_HALF_SPINOR_FROM_SINGLE_ST(yH, yN, y, stride);
#define REDUCE_REAL_OPERATION(i) (rdot0)
#define REDUCE_IMAG_OPERATION(i) (idot0)
#include "reduce_complex_core.h"
#undef REDUCE_FUNC_NAME
#undef REDUCE_TYPES
#undef REDUCE_PARAMS
#undef REDUCE_REAL_AUXILIARY
#undef REDUCE_IMAG_AUXILIARY
#undef REDUCE_REAL_OPERATION
#undef REDUCE_IMAG_OPERATION

Complex xpaycDotzyCuda(cudaColorSpinorField &x, const double &a, cudaColorSpinorField &y, cudaColorSpinorField &z) {
  const int id = 20;
  blas_quda_flops += 6*x.real_length;
  checkSpinor(x,y);
  checkSpinor(x,z);
  int length = x.length/2;
  blas_quda_bytes += 4*x.real_length*x.precision;
  double2 dot;
  if (x.precision == QUDA_DOUBLE_PRECISION) {
    int spinor_bytes = x.length*sizeof(double);
    cudaBindTexture(0, xTexDouble2, x.v, spinor_bytes); 
    cudaBindTexture(0, yTexDouble2, y.v, spinor_bytes); 
    cudaBindTexture(0, zTexDouble2, z.v, spinor_bytes); 
    dot = xpaycDotzyDCuda((double2*)x.v, a, (double2*)y.v, (double2*)z.v, length, id, x.precision);
  } else if (x.precision == QUDA_SINGLE_PRECISION) {
    dot = xpaycDotzySCuda((float2*)x.v, (float)a, (float2*)y.v, (float2*)z.v, length, id, x.precision);
  } else {
    if (x.siteSubset == QUDA_FULL_SITE_SUBSET) 
      return xpaycDotzyCuda(x.Even(), a, y.Even(), z.Even()) + xpaycDotzyCuda(x.Odd(), a, y.Odd(), z.Odd());
    int spinor_bytes = x.length*sizeof(short);
    blas_quda_bytes += 4*x.volume*sizeof(float);
    if (x.nSpin ==4 ){//wilson
      cudaBindTexture(0, texHalf1, x.v, spinor_bytes); 
      cudaBindTexture(0, texNorm1, x.norm, spinor_bytes/12);    
      cudaBindTexture(0, texHalf2, y.v, spinor_bytes); 
      cudaBindTexture(0, texNorm2, y.norm, spinor_bytes/12);    
      cudaBindTexture(0, texHalf3, z.v, spinor_bytes); 
      cudaBindTexture(0, texNorm3, z.norm, spinor_bytes/12);    
      dot = xpaycDotzyHCuda((float)a, (short4*)y.v, (float*)y.norm, x.stride, x.volume, id, x.precision);
    } else if (x.nSpin ==1 ){//wilson
      cudaBindTexture(0, texHalfSt1, x.v, spinor_bytes); 
      cudaBindTexture(0, texNorm1, x.norm, spinor_bytes/3);    
      cudaBindTexture(0, texHalfSt2, y.v, spinor_bytes); 
      cudaBindTexture(0, texNorm2, y.norm, spinor_bytes/3);    
      cudaBindTexture(0, texHalfSt3, z.v, spinor_bytes); 
      cudaBindTexture(0, texNorm3, z.norm, spinor_bytes/3);    
      dot = xpaycDotzyHCuda((float)a, (short2*)y.v, (float*)y.norm, x.stride, x.volume, id, x.precision);
    }else{
      errorQuda("%s: nSpin(%d) is not supported\n", __FUNCTION__, x.nSpin);            
    }
  }
  return Complex(dot.x, dot.y);
}


//
// double3 cDotProductNormACuda(float2 *a, float2 *b, int n) {}
//
template <unsigned int reduce_threads, typename Float2>
#define REDUCE_FUNC_NAME(suffix) cDotProductNormAD##suffix
#define REDUCE_TYPES Float2 *x, Float2 *y
#define REDUCE_PARAMS x, y
#define REDUCE_X_AUXILIARY(i) Float2 a = READ_DOUBLE2_TEXTURE(x, i);
#define REDUCE_Y_AUXILIARY(i) Float2 b = READ_DOUBLE2_TEXTURE(y, i);
#define REDUCE_Z_AUXILIARY(i)
#define REDUCE_X_OPERATION(i) (a.x*b.x + a.y*b.y)
#define REDUCE_Y_OPERATION(i) (a.x*b.y - a.y*b.x)
#define REDUCE_Z_OPERATION(i) (a.x*a.x + a.y*a.y)
#include "reduce_triple_core.h"
#undef REDUCE_FUNC_NAME
#undef REDUCE_TYPES
#undef REDUCE_PARAMS
#undef REDUCE_X_AUXILIARY
#undef REDUCE_Y_AUXILIARY
#undef REDUCE_Z_AUXILIARY
#undef REDUCE_X_OPERATION
#undef REDUCE_Y_OPERATION
#undef REDUCE_Z_OPERATION

template <unsigned int reduce_threads, typename Float2>
#define REDUCE_FUNC_NAME(suffix) cDotProductNormAS##suffix
#define REDUCE_TYPES Float2 *a, Float2 *b
#define REDUCE_PARAMS a, b
#define REDUCE_X_AUXILIARY(i)
#define REDUCE_Y_AUXILIARY(i)
#define REDUCE_Z_AUXILIARY(i)
#define REDUCE_X_OPERATION(i) (a[i].x*b[i].x + a[i].y*b[i].y)
#define REDUCE_Y_OPERATION(i) (a[i].x*b[i].y - a[i].y*b[i].x)
#define REDUCE_Z_OPERATION(i) (a[i].x*a[i].x + a[i].y*a[i].y)
#include "reduce_triple_core.h"
#undef REDUCE_FUNC_NAME
#undef REDUCE_TYPES
#undef REDUCE_PARAMS
#undef REDUCE_X_AUXILIARY
#undef REDUCE_Y_AUXILIARY
#undef REDUCE_Z_AUXILIARY
#undef REDUCE_X_OPERATION
#undef REDUCE_Y_OPERATION
#undef REDUCE_Z_OPERATION

template <unsigned int reduce_threads, typename Float2>
#define REDUCE_FUNC_NAME(suffix) cDotProductNormAH##suffix
#define REDUCE_TYPES Float2 *xN, Float2 *yN, int stride
#define REDUCE_PARAMS xN, yN, stride
#define REDUCE_X_AUXILIARY(i)                                           \
  READ_HALF_SPINOR(x, texHalf1, stride);                                \
  READ_HALF_SPINOR(y, texHalf2, stride);                                \
  REAL_DOT_FLOAT4(norm0, x0, x0);                                       \
  REAL_DOT_FLOAT4(norm1, x1, x1);                                       \
  REAL_DOT_FLOAT4(norm2, x2, x2);                                       \
  REAL_DOT_FLOAT4(norm3, x3, x3);                                       \
  REAL_DOT_FLOAT4(norm4, x4, x4);                                       \
  REAL_DOT_FLOAT4(norm5, x5, x5);                                       \
  norm0 += norm1; norm2 += norm3; norm4 += norm5; norm0 += norm2, norm0 += norm4;
#define REDUCE_Y_AUXILIARY(i)                                           \
  REAL_DOT_FLOAT4(rdot0, x0, y0);                                       \
  REAL_DOT_FLOAT4(rdot1, x1, y1);                                       \
  REAL_DOT_FLOAT4(rdot2, x2, y2);                                       \
  REAL_DOT_FLOAT4(rdot3, x3, y3);                                       \
  REAL_DOT_FLOAT4(rdot4, x4, y4);                                       \
  REAL_DOT_FLOAT4(rdot5, x5, y5);                                       \
  rdot0 += rdot1; rdot2 += rdot3; rdot4 += rdot5; rdot0 += rdot2; rdot0 += rdot4;
#define REDUCE_Z_AUXILIARY(i)                                           \
  IMAG_DOT_FLOAT4(idot0, x0, y0);                                       \
  IMAG_DOT_FLOAT4(idot1, x1, y1);                                       \
  IMAG_DOT_FLOAT4(idot2, x2, y2);                                       \
  IMAG_DOT_FLOAT4(idot3, x3, y3);                                       \
  IMAG_DOT_FLOAT4(idot4, x4, y4);                                       \
  IMAG_DOT_FLOAT4(idot5, x5, y5);                                       \
  idot0 += idot1; idot2 += idot3; idot4 += idot5; idot0 += idot2; idot0 += idot4;  
#define REDUCE_X_OPERATION(i) (xc*yc*rdot0)
#define REDUCE_Y_OPERATION(i) (xc*yc*idot0)
#define REDUCE_Z_OPERATION(i) (xc*xc*norm0)
#include "reduce_triple_core.h"
#undef REDUCE_FUNC_NAME
#undef REDUCE_TYPES
#undef REDUCE_PARAMS
#undef REDUCE_X_AUXILIARY
#undef REDUCE_Y_AUXILIARY
#undef REDUCE_Z_AUXILIARY
#undef REDUCE_X_OPERATION
#undef REDUCE_Y_OPERATION
#undef REDUCE_Z_OPERATION

template <unsigned int reduce_threads, typename Float2>
#define REDUCE_FUNC_NAME(suffix) cDotProductNormAHSt##suffix
#define REDUCE_TYPES Float2 *xN, Float2 *yN, int stride
#define REDUCE_PARAMS xN, yN, stride
#define REDUCE_X_AUXILIARY(i)                                           \
  READ_HALF_SPINOR_ST(x, texHalfSt1, stride);                           \
  READ_HALF_SPINOR_ST(y, texHalfSt2, stride);                           \
  REAL_DOT_FLOAT2(norm0, x0, x0);                                       \
  REAL_DOT_FLOAT2(norm1, x1, x1);                                       \
  REAL_DOT_FLOAT2(norm2, x2, x2);                                       \
  norm0 += norm1; norm0 += norm2;
#define REDUCE_Y_AUXILIARY(i)                                           \
  REAL_DOT_FLOAT2(rdot0, x0, y0);                                       \
  REAL_DOT_FLOAT2(rdot1, x1, y1);                                       \
  REAL_DOT_FLOAT2(rdot2, x2, y2);                                       \
  rdot0 += rdot1; rdot0 += rdot2;
#define REDUCE_Z_AUXILIARY(i)                                           \
  IMAG_DOT_FLOAT2(idot0, x0, y0);                                       \
  IMAG_DOT_FLOAT2(idot1, x1, y1);                                       \
  IMAG_DOT_FLOAT2(idot2, x2, y2);                                       \
  idot0 += idot1; idot0 += idot2;
#define REDUCE_X_OPERATION(i) (xc*yc*rdot0)
#define REDUCE_Y_OPERATION(i) (xc*yc*idot0)
#define REDUCE_Z_OPERATION(i) (xc*xc*norm0)
#include "reduce_triple_core.h"
#undef REDUCE_FUNC_NAME
#undef REDUCE_TYPES
#undef REDUCE_PARAMS
#undef REDUCE_X_AUXILIARY
#undef REDUCE_Y_AUXILIARY
#undef REDUCE_Z_AUXILIARY
#undef REDUCE_X_OPERATION
#undef REDUCE_Y_OPERATION
#undef REDUCE_Z_OPERATION

double3 cDotProductNormACuda(cudaColorSpinorField &x, cudaColorSpinorField &y) {
  const int id = 21;
  blas_quda_flops += 6*x.real_length;
  checkSpinor(x,y);
  int length = x.length/2;
  blas_quda_bytes += 2*x.real_length*x.precision;
  if (x.precision == QUDA_DOUBLE_PRECISION) {
    int spinor_bytes = x.length*sizeof(double);
    cudaBindTexture(0, xTexDouble2, x.v, spinor_bytes); 
    cudaBindTexture(0, yTexDouble2, y.v, spinor_bytes); 
    return cDotProductNormADCuda((double2*)x.v, (double2*)y.v, length, id, x.precision);
  } else if (x.precision == QUDA_SINGLE_PRECISION) {
    return cDotProductNormASCuda((float2*)x.v, (float2*)y.v, length, id, x.precision);
  } else {
    if (x.siteSubset == QUDA_FULL_SITE_SUBSET) 
      return cDotProductNormACuda(x.Even(), y.Even()) + cDotProductNormACuda(x.Odd(), y.Odd());
    int spinor_bytes = x.length*sizeof(short);
    blas_quda_bytes += 2*x.volume*sizeof(float);
    if (x.nSpin == 4){ //wilson
      cudaBindTexture(0, texHalf1, x.v, spinor_bytes); 
      cudaBindTexture(0, texNorm1, x.norm, spinor_bytes/12);    
      cudaBindTexture(0, texHalf2, y.v, spinor_bytes); 
      cudaBindTexture(0, texNorm2, y.norm, spinor_bytes/12);    
      return cDotProductNormAHCuda((float*)x.norm, (float*)y.norm, x.stride, x.volume, id, x.precision);
    } else if (x.nSpin == 1){ //staggered
      cudaBindTexture(0, texHalfSt1, x.v, spinor_bytes); 
      cudaBindTexture(0, texNorm1, x.norm, spinor_bytes/3);    
      cudaBindTexture(0, texHalfSt2, y.v, spinor_bytes); 
      cudaBindTexture(0, texNorm2, y.norm, spinor_bytes/3);    
      return cDotProductNormAHStCuda((float*)x.norm, (float*)y.norm, x.stride, x.volume, id, x.precision);
    }else{
      errorQuda("%s: nSpin(%d) is not supported\n", __FUNCTION__, x.nSpin);            
    }
  }

}

//
// double3 cDotProductNormBCuda(float2 *a, float2 *b, int n) {}
//
template <unsigned int reduce_threads, typename Float2>
#define REDUCE_FUNC_NAME(suffix) cDotProductNormBD##suffix
#define REDUCE_TYPES Float2 *x, Float2 *y
#define REDUCE_PARAMS x, y
#define REDUCE_X_AUXILIARY(i) Float2 a = READ_DOUBLE2_TEXTURE(x, i);
#define REDUCE_Y_AUXILIARY(i) Float2 b = READ_DOUBLE2_TEXTURE(y, i);
#define REDUCE_Z_AUXILIARY(i)
#define REDUCE_X_OPERATION(i) (a.x*b.x + a.y*b.y)
#define REDUCE_Y_OPERATION(i) (a.x*b.y - a.y*b.x)
#define REDUCE_Z_OPERATION(i) (b.x*b.x + b.y*b.y)
#include "reduce_triple_core.h"
#undef REDUCE_FUNC_NAME
#undef REDUCE_TYPES
#undef REDUCE_PARAMS
#undef REDUCE_X_AUXILIARY
#undef REDUCE_Y_AUXILIARY
#undef REDUCE_Z_AUXILIARY
#undef REDUCE_X_OPERATION
#undef REDUCE_Y_OPERATION
#undef REDUCE_Z_OPERATION

template <unsigned int reduce_threads, typename Float2>
#define REDUCE_FUNC_NAME(suffix) cDotProductNormBS##suffix
#define REDUCE_TYPES Float2 *a, Float2 *b
#define REDUCE_PARAMS a, b
#define REDUCE_X_AUXILIARY(i)
#define REDUCE_Y_AUXILIARY(i)
#define REDUCE_Z_AUXILIARY(i)
#define REDUCE_X_OPERATION(i) (a[i].x*b[i].x + a[i].y*b[i].y)
#define REDUCE_Y_OPERATION(i) (a[i].x*b[i].y - a[i].y*b[i].x)
#define REDUCE_Z_OPERATION(i) (b[i].x*b[i].x + b[i].y*b[i].y)
#include "reduce_triple_core.h"
#undef REDUCE_FUNC_NAME
#undef REDUCE_TYPES
#undef REDUCE_PARAMS
#undef REDUCE_X_AUXILIARY
#undef REDUCE_Y_AUXILIARY
#undef REDUCE_Z_AUXILIARY
#undef REDUCE_X_OPERATION
#undef REDUCE_Y_OPERATION
#undef REDUCE_Z_OPERATION

template <unsigned int reduce_threads, typename Float2>
#define REDUCE_FUNC_NAME(suffix) cDotProductNormBH##suffix
#define REDUCE_TYPES Float2 *xN, Float2 *yN, int stride
#define REDUCE_PARAMS xN, yN, stride
#define REDUCE_X_AUXILIARY(i)                                           \
  READ_HALF_SPINOR(x, texHalf1, stride);                                \
  READ_HALF_SPINOR(y, texHalf2, stride);                                \
  REAL_DOT_FLOAT4(norm0, y0, y0);                                       \
  REAL_DOT_FLOAT4(norm1, y1, y1);                                       \
  REAL_DOT_FLOAT4(norm2, y2, y2);                                       \
  REAL_DOT_FLOAT4(norm3, y3, y3);                                       \
  REAL_DOT_FLOAT4(norm4, y4, y4);                                       \
  REAL_DOT_FLOAT4(norm5, y5, y5);                                       \
  norm0 += norm1; norm2 += norm3; norm4 += norm5; norm0 += norm2, norm0 += norm4;
#define REDUCE_Y_AUXILIARY(i)                                           \
  REAL_DOT_FLOAT4(rdot0, x0, y0);                                       \
  REAL_DOT_FLOAT4(rdot1, x1, y1);                                       \
  REAL_DOT_FLOAT4(rdot2, x2, y2);                                       \
  REAL_DOT_FLOAT4(rdot3, x3, y3);                                       \
  REAL_DOT_FLOAT4(rdot4, x4, y4);                                       \
  REAL_DOT_FLOAT4(rdot5, x5, y5);                                       \
  rdot0 += rdot1; rdot2 += rdot3; rdot4 += rdot5; rdot0 += rdot2; rdot0 += rdot4;
#define REDUCE_Z_AUXILIARY(i)                                           \
  IMAG_DOT_FLOAT4(idot0, x0, y0);                                       \
  IMAG_DOT_FLOAT4(idot1, x1, y1);                                       \
  IMAG_DOT_FLOAT4(idot2, x2, y2);                                       \
  IMAG_DOT_FLOAT4(idot3, x3, y3);                                       \
  IMAG_DOT_FLOAT4(idot4, x4, y4);                                       \
  IMAG_DOT_FLOAT4(idot5, x5, y5);                                       \
  idot0 += idot1; idot2 += idot3; idot4 += idot5; idot0 += idot2; idot0 += idot4;
#define REDUCE_X_OPERATION(i) (xc*yc*rdot0)
#define REDUCE_Y_OPERATION(i) (xc*yc*idot0)
#define REDUCE_Z_OPERATION(i) (yc*yc*norm0)
#include "reduce_triple_core.h"
#undef REDUCE_FUNC_NAME
#undef REDUCE_TYPES
#undef REDUCE_PARAMS
#undef REDUCE_X_AUXILIARY
#undef REDUCE_Y_AUXILIARY
#undef REDUCE_Z_AUXILIARY
#undef REDUCE_X_OPERATION
#undef REDUCE_Y_OPERATION
#undef REDUCE_Z_OPERATION

template <unsigned int reduce_threads, typename Float2>
#define REDUCE_FUNC_NAME(suffix) cDotProductNormBHSt##suffix
#define REDUCE_TYPES Float2 *xN, Float2 *yN, int stride
#define REDUCE_PARAMS xN, yN, stride
#define REDUCE_X_AUXILIARY(i)                                           \
  READ_HALF_SPINOR_ST(x, texHalfSt1, stride);                           \
  READ_HALF_SPINOR_ST(y, texHalfSt2, stride);                           \
  REAL_DOT_FLOAT2(norm0, y0, y0);                                       \
  REAL_DOT_FLOAT2(norm1, y1, y1);                                       \
  REAL_DOT_FLOAT2(norm2, y2, y2);                                       \
  norm0 += norm1; norm0 += norm2;
#define REDUCE_Y_AUXILIARY(i)                                           \
  REAL_DOT_FLOAT2(rdot0, x0, y0);                                       \
  REAL_DOT_FLOAT2(rdot1, x1, y1);                                       \
  REAL_DOT_FLOAT2(rdot2, x2, y2);                                       \
  rdot0 += rdot1; rdot0 += rdot2;
#define REDUCE_Z_AUXILIARY(i)                                           \
  IMAG_DOT_FLOAT2(idot0, x0, y0);                                       \
  IMAG_DOT_FLOAT2(idot1, x1, y1);                                       \
  IMAG_DOT_FLOAT2(idot2, x2, y2);                                       \
  idot0 += idot1; idot0 += idot2;
#define REDUCE_X_OPERATION(i) (xc*yc*rdot0)
#define REDUCE_Y_OPERATION(i) (xc*yc*idot0)
#define REDUCE_Z_OPERATION(i) (yc*yc*norm0)
#include "reduce_triple_core.h"
#undef REDUCE_FUNC_NAME
#undef REDUCE_TYPES
#undef REDUCE_PARAMS
#undef REDUCE_X_AUXILIARY
#undef REDUCE_Y_AUXILIARY
#undef REDUCE_Z_AUXILIARY
#undef REDUCE_X_OPERATION
#undef REDUCE_Y_OPERATION
#undef REDUCE_Z_OPERATION

double3 cDotProductNormBCuda(cudaColorSpinorField &x, cudaColorSpinorField &y) {
  const int id = 22;
  blas_quda_flops += 6*x.real_length;
  checkSpinor(x,y);
  int length = x.length/2;
  blas_quda_bytes += 2*x.real_length*x.precision;
  if (x.precision == QUDA_DOUBLE_PRECISION) {
    int spinor_bytes = x.length*sizeof(double);
    cudaBindTexture(0, xTexDouble2, x.v, spinor_bytes); 
    cudaBindTexture(0, yTexDouble2, y.v, spinor_bytes); 
    return cDotProductNormBDCuda((double2*)x.v, (double2*)y.v, length, id, x.precision);
  } else if (x.precision == QUDA_SINGLE_PRECISION) {
    return cDotProductNormBSCuda((float2*)x.v, (float2*)y.v, length, id, x.precision);
  } else {
    if (x.siteSubset == QUDA_FULL_SITE_SUBSET) 
      return cDotProductNormBCuda(x.Even(), y.Even()) + cDotProductNormBCuda(x.Odd(), y.Odd());
    int spinor_bytes = x.length*sizeof(short);
    blas_quda_bytes += 2*x.volume*sizeof(float);
    if (x.nSpin == 4){ //wilson
      cudaBindTexture(0, texHalf1, x.v, spinor_bytes); 
      cudaBindTexture(0, texNorm1, x.norm, spinor_bytes/12);    
      cudaBindTexture(0, texHalf2, y.v, spinor_bytes); 
      cudaBindTexture(0, texNorm2, y.norm, spinor_bytes/12);    
      return cDotProductNormBHCuda((float*)x.norm, (float*)y.norm, x.stride, x.volume, id, x.precision);
    } else if (x.nSpin == 1){ //staggered
      cudaBindTexture(0, texHalfSt1, x.v, spinor_bytes); 
      cudaBindTexture(0, texNorm1, x.norm, spinor_bytes/3);    
      cudaBindTexture(0, texHalfSt2, y.v, spinor_bytes); 
      cudaBindTexture(0, texNorm2, y.norm, spinor_bytes/3);    
      return cDotProductNormBHStCuda((float*)x.norm, (float*)y.norm, x.stride, x.volume, id, x.precision);
    }else{
      errorQuda("%s: nSpin(%d) is not supported\n", __FUNCTION__, x.nSpin);            
    }
  }

}


//
// double3 caxpbypzYmbwcDotProductWYNormYCuda(float2 a, float2 *x, float2 b, float2 *y, 
// float2 *z, float2 *w, float2 *u, int len)
//
template <unsigned int reduce_threads, typename Float2>
#define REDUCE_FUNC_NAME(suffix) caxpbypzYmbwcDotProductWYNormYD##suffix
#define REDUCE_TYPES Float2 a, Float2 *x, Float2 b, Float2 *y, Float2 *z, Float2 *w, Float2 *u
#define REDUCE_PARAMS a, x, b, y, z, w, u
#define REDUCE_X_AUXILIARY(i)                           \
  Float2 X = READ_DOUBLE2_TEXTURE(x, i);                \
  Float2 Y = READ_DOUBLE2_TEXTURE(y, i);                \
  Float2 W = READ_DOUBLE2_TEXTURE(w, i);                                
#define REDUCE_Y_AUXILIARY(i)                   \
  Float2 Z = read_Float2(z, i);                 \
  Z.x += a.x*X.x - a.y*X.y;                     \
  Z.y += a.y*X.x + a.x*X.y;                     \
  Z.x += b.x*Y.x - b.y*Y.y;                     \
  Z.y += b.y*Y.x + b.x*Y.y;                     \
  Y.x -= b.x*W.x - b.y*W.y;                     \
  Y.y -= b.y*W.x + b.x*W.y;     
#define REDUCE_Z_AUXILIARY(i)         \
  z[i] = make_Float2(Z);              \
  y[i] = make_Float2(Y);              
#define REDUCE_X_OPERATION(i) (u[i].x*y[i].x + u[i].y*y[i].y)
#define REDUCE_Y_OPERATION(i) (u[i].x*y[i].y - u[i].y*y[i].x)
#define REDUCE_Z_OPERATION(i) (y[i].x*y[i].x + y[i].y*y[i].y)
#include "reduce_triple_core.h"
#undef REDUCE_FUNC_NAME
#undef REDUCE_TYPES
#undef REDUCE_PARAMS
#undef REDUCE_X_AUXILIARY
#undef REDUCE_Y_AUXILIARY
#undef REDUCE_Z_AUXILIARY
#undef REDUCE_X_OPERATION
#undef REDUCE_Y_OPERATION
#undef REDUCE_Z_OPERATION

template <unsigned int reduce_threads, typename Float2>
#define REDUCE_FUNC_NAME(suffix) caxpbypzYmbwcDotProductWYNormYS##suffix
#define REDUCE_TYPES Float2 a, Float2 *x, Float2 b, Float2 *y, Float2 *z, Float2 *w, Float2 *u
#define REDUCE_PARAMS a, x, b, y, z, w, u
#define REDUCE_X_AUXILIARY(i)                           \
  Float2 X = read_Float2(x, i);         \
  Float2 Y = read_Float2(y, i);         \
  Float2 W = read_Float2(w, i);         
#define REDUCE_Y_AUXILIARY(i)                   \
  Float2 Z = read_Float2(z, i); \
  Z.x += a.x*X.x - a.y*X.y;                     \
  Z.y += a.y*X.x + a.x*X.y;                     \
  Z.x += b.x*Y.x - b.y*Y.y;                     \
  Z.y += b.y*Y.x + b.x*Y.y;                     \
  Y.x -= b.x*W.x - b.y*W.y;                     \
  Y.y -= b.y*W.x + b.x*W.y;     
#define REDUCE_Z_AUXILIARY(i)         \
  z[i] = make_Float2(Z);              \
  y[i] = make_Float2(Y);              
#define REDUCE_X_OPERATION(i) (u[i].x*y[i].x + u[i].y*y[i].y)
#define REDUCE_Y_OPERATION(i) (u[i].x*y[i].y - u[i].y*y[i].x)
#define REDUCE_Z_OPERATION(i) (y[i].x*y[i].x + y[i].y*y[i].y)
#include "reduce_triple_core.h"
#undef REDUCE_FUNC_NAME
#undef REDUCE_TYPES
#undef REDUCE_PARAMS
#undef REDUCE_X_AUXILIARY
#undef REDUCE_Y_AUXILIARY
#undef REDUCE_Z_AUXILIARY
#undef REDUCE_X_OPERATION
#undef REDUCE_Y_OPERATION
#undef REDUCE_Z_OPERATION

//
// double3 caxpbypzYmbwcDotProductWYNormYCuda(float2 a, float2 *x, float2 b, float2 *y, 
// float2 *z, float2 *w, float2 *u, int len)
//
template <unsigned int reduce_threads, typename Float2>
#define REDUCE_FUNC_NAME(suffix) caxpbypzYmbwcDotProductWYNormYH##suffix
#define REDUCE_TYPES Float2 a, Float2 b, short4 *yH, float *yN, short4 *zH, float *zN, float *wN, float *uN, int stride
#define REDUCE_PARAMS a, b, yH, yN, zH, zN, wN, uN, stride
#define REDUCE_X_AUXILIARY(i)                                           \
  RECONSTRUCT_HALF_SPINOR(x, texHalf1, texNorm1, stride);               \
  RECONSTRUCT_HALF_SPINOR(y, texHalf2, texNorm2, stride);               \
  RECONSTRUCT_HALF_SPINOR(z, texHalf3, texNorm3, stride);               \
  CAXPBYPZ_FLOAT4(a, x0, b, y0, z0);                                    \
  CAXPBYPZ_FLOAT4(a, x1, b, y1, z1);                                    \
  CAXPBYPZ_FLOAT4(a, x2, b, y2, z2);                                    \
  CAXPBYPZ_FLOAT4(a, x3, b, y3, z3);                                    \
  CAXPBYPZ_FLOAT4(a, x4, b, y4, z4);                                    \
  CAXPBYPZ_FLOAT4(a, x5, b, y5, z5);                                    \
  CONSTRUCT_HALF_SPINOR_FROM_SINGLE(zH, zN, z, stride);                 \
  READ_HALF_SPINOR(w, texHalf4, stride);                                \
  float2 bwc = -wc*b;                                                   \
  CAXPY_FLOAT4(bwc, w0, y0);                                            \
  CAXPY_FLOAT4(bwc, w1, y1);                                            \
  CAXPY_FLOAT4(bwc, w2, y2);                                            \
  CAXPY_FLOAT4(bwc, w3, y3);                                            \
  CAXPY_FLOAT4(bwc, w4, y4);                                            \
  CAXPY_FLOAT4(bwc, w5, y5);                                            \
  REAL_DOT_FLOAT4(norm0, y0, y0);                                       \
  REAL_DOT_FLOAT4(norm1, y1, y1);                                       \
  REAL_DOT_FLOAT4(norm2, y2, y2);                                       \
  REAL_DOT_FLOAT4(norm3, y3, y3);                                       \
  REAL_DOT_FLOAT4(norm4, y4, y4);                                       \
  REAL_DOT_FLOAT4(norm5, y5, y5);                                       \
  CONSTRUCT_HALF_SPINOR_FROM_SINGLE(yH, yN, y, stride);                 
#define REDUCE_Y_AUXILIARY(i)                                           \
  READ_HALF_SPINOR(u, texHalf5, stride);                                \
  REAL_DOT_FLOAT4(rdot0, u0, y0);                                       \
  REAL_DOT_FLOAT4(rdot1, u1, y1);                                       \
  REAL_DOT_FLOAT4(rdot2, u2, y2);                                       \
  REAL_DOT_FLOAT4(rdot3, u3, y3);                                       \
  REAL_DOT_FLOAT4(rdot4, u4, y4);                                       \
  REAL_DOT_FLOAT4(rdot5, u5, y5);                                       \
  IMAG_DOT_FLOAT4(idot0, u0, y0);                                       \
  IMAG_DOT_FLOAT4(idot1, u1, y1);                                       \
  IMAG_DOT_FLOAT4(idot2, u2, y2);                                       \
  IMAG_DOT_FLOAT4(idot3, u3, y3);                                       \
  IMAG_DOT_FLOAT4(idot4, u4, y4);                                       \
  IMAG_DOT_FLOAT4(idot5, u5, y5);                                       
#define REDUCE_Z_AUXILIARY(i)                                           \
  norm0 += norm1; norm2 += norm3; norm4 += norm5; norm0 += norm2, norm0 += norm4; \
  rdot0 += rdot1; rdot2 += rdot3; rdot4 += rdot5; rdot0 += rdot2; rdot0 += rdot4; \
  idot0 += idot1; idot2 += idot3; idot4 += idot5; idot0 += idot2; idot0 += idot4;

#define REDUCE_X_OPERATION(i) (uc*rdot0)
#define REDUCE_Y_OPERATION(i) (uc*idot0)
#define REDUCE_Z_OPERATION(i) (norm0)

#include "reduce_triple_core.h"
#undef REDUCE_FUNC_NAME
#undef REDUCE_TYPES
#undef REDUCE_PARAMS
#undef REDUCE_X_AUXILIARY
#undef REDUCE_Y_AUXILIARY
#undef REDUCE_Z_AUXILIARY
#undef REDUCE_X_OPERATION
#undef REDUCE_Y_OPERATION
#undef REDUCE_Z_OPERATION


template <unsigned int reduce_threads, typename Float2>
#define REDUCE_FUNC_NAME(suffix) caxpbypzYmbwcDotProductWYNormYH##suffix
#define REDUCE_TYPES Float2 a, Float2 b, short2 *yH, float *yN, short2 *zH, float *zN, float *wN, float *uN, int stride
#define REDUCE_PARAMS a, b, yH, yN, zH, zN, wN, uN, stride
#define REDUCE_X_AUXILIARY(i)                                           \
  RECONSTRUCT_HALF_SPINOR_ST(x, texHalfSt1, texNorm1, stride);          \
  RECONSTRUCT_HALF_SPINOR_ST(y, texHalfSt2, texNorm2, stride);          \
  RECONSTRUCT_HALF_SPINOR_ST(z, texHalfSt3, texNorm3, stride);          \
  CAXPBYPZ_FLOAT2(a, x0, b, y0, z0);                                    \
  CAXPBYPZ_FLOAT2(a, x1, b, y1, z1);                                    \
  CAXPBYPZ_FLOAT2(a, x2, b, y2, z2);                                    \
  CONSTRUCT_HALF_SPINOR_FROM_SINGLE_ST(zH, zN, z, stride);                      \
  READ_HALF_SPINOR_ST(w, texHalfSt4, stride);                           \
  float2 bwc = -wc*b;                                                   \
  CAXPY_FLOAT2(bwc, w0, y0);                                            \
  CAXPY_FLOAT2(bwc, w1, y1);                                            \
  CAXPY_FLOAT2(bwc, w2, y2);                                            \
  REAL_DOT_FLOAT2(norm0, y0, y0);                                       \
  REAL_DOT_FLOAT2(norm1, y1, y1);                                       \
  REAL_DOT_FLOAT2(norm2, y2, y2);                                       \
  CONSTRUCT_HALF_SPINOR_FROM_SINGLE_ST(yH, yN, y, stride);                      
#define REDUCE_Y_AUXILIARY(i)                                           \
  READ_HALF_SPINOR_ST(u, texHalfSt5, stride);                           \
  REAL_DOT_FLOAT2(rdot0, u0, y0);                                       \
  REAL_DOT_FLOAT2(rdot1, u1, y1);                                       \
  REAL_DOT_FLOAT2(rdot2, u2, y2);                                       \
  IMAG_DOT_FLOAT2(idot0, u0, y0);                                       \
  IMAG_DOT_FLOAT2(idot1, u1, y1);                                       \
  IMAG_DOT_FLOAT2(idot2, u2, y2);                                       
#define REDUCE_Z_AUXILIARY(i)                                           \
  norm0 += norm1; norm0 += norm2; \
  rdot0 += rdot1; rdot0 += rdot2; \
  idot0 += idot1; idot0 += idot2; 

#define REDUCE_X_OPERATION(i) (uc*rdot0)
#define REDUCE_Y_OPERATION(i) (uc*idot0)
#define REDUCE_Z_OPERATION(i) (norm0)

#include "reduce_triple_core.h"
#undef REDUCE_FUNC_NAME
#undef REDUCE_TYPES
#undef REDUCE_PARAMS
#undef REDUCE_X_AUXILIARY
#undef REDUCE_Y_AUXILIARY
#undef REDUCE_Z_AUXILIARY
#undef REDUCE_X_OPERATION
#undef REDUCE_Y_OPERATION
#undef REDUCE_Z_OPERATION

// This convoluted kernel does the following: z += a*x + b*y, y -= b*w, norm = (y,y), dot = (u, y)
double3 caxpbypzYmbwcDotProductWYNormYCuda(const Complex &a, cudaColorSpinorField &x, const Complex &b, cudaColorSpinorField &y,
                                           cudaColorSpinorField &z, cudaColorSpinorField &w, cudaColorSpinorField &u) {
  const int id = 23;
  blas_quda_flops += 18*x.real_length;
  checkSpinor(x,y);
  checkSpinor(x,z);
  checkSpinor(x,w);
  checkSpinor(x,u);
  int length = x.length/2;
  blas_quda_bytes += 7*x.real_length*x.precision;
  if (x.precision == QUDA_DOUBLE_PRECISION) {
    int spinor_bytes = x.length*sizeof(double);
    cudaBindTexture(0, xTexDouble2, x.v, spinor_bytes); 
    cudaBindTexture(0, yTexDouble2, y.v, spinor_bytes); 
    cudaBindTexture(0, zTexDouble2, z.v, spinor_bytes); 
    cudaBindTexture(0, wTexDouble2, w.v, spinor_bytes); 
    cudaBindTexture(0, uTexDouble2, u.v, spinor_bytes); 
    double2 a2 = make_double2(real(a), imag(a));
    double2 b2 = make_double2(real(b), imag(b));
    return caxpbypzYmbwcDotProductWYNormYDCuda(a2, (double2*)x.v, b2, (double2*)y.v, (double2*)z.v, 
                                               (double2*)w.v, (double2*)u.v, length, id, x.precision);
  } else if (x.precision == QUDA_SINGLE_PRECISION) {
    float2 a2 = make_float2(real(a), imag(a));
    float2 b2 = make_float2(real(b), imag(b));
    return caxpbypzYmbwcDotProductWYNormYSCuda(a2, (float2*)x.v, b2, (float2*)y.v, (float2*)z.v,
                                               (float2*)w.v, (float2*)u.v, length, id, x.precision);
  } else {
    // fused nSpin=4 kernel is slow on Fermi
    // N.B. this introduces an extra half truncation so will affect convergence
    if (!blasTuning && (__CUDA_ARCH__ >= 200) && x.nSpin == 4) {
      caxpbypzYmbwCuda(a, x, b, y, z, w);
      return cDotProductNormBCuda(u, y);
    }
      
    if (x.siteSubset == QUDA_FULL_SITE_SUBSET) 
      return caxpbypzYmbwcDotProductWYNormYCuda(a, x.Even(), b, y.Even(), z.Even(), w.Even(), u.Even()) + 
        caxpbypzYmbwcDotProductWYNormYCuda(a, x.Odd(), b, y.Odd(), z.Odd(), w.Odd(), u.Odd());
    int spinor_bytes = x.length*sizeof(short);
    blas_quda_bytes += 7*x.volume*sizeof(float);
    if (x.nSpin == 4) { // wilson
      cudaBindTexture(0, texHalf1, x.v, spinor_bytes); 
      cudaBindTexture(0, texNorm1, x.norm, spinor_bytes/12);    
      cudaBindTexture(0, texHalf2, y.v, spinor_bytes); 
      cudaBindTexture(0, texNorm2, y.norm, spinor_bytes/12);    
      cudaBindTexture(0, texHalf3, z.v, spinor_bytes); 
      cudaBindTexture(0, texNorm3, z.norm, spinor_bytes/12);    
      cudaBindTexture(0, texHalf4, w.v, spinor_bytes); 
      cudaBindTexture(0, texNorm4, w.norm, spinor_bytes/12);    
      cudaBindTexture(0, texHalf5, u.v, spinor_bytes); 
      cudaBindTexture(0, texNorm5, u.norm, spinor_bytes/12);    
      float2 a2 = make_float2(real(a), imag(a));
      float2 b2 = make_float2(real(b), imag(b));
      return caxpbypzYmbwcDotProductWYNormYHCuda(a2, b2, (short4*)y.v, (float*)y.norm, 
                                                 (short4*)z.v, (float*)z.norm, (float*)w.norm, (float*)u.norm, 
                                                 y.stride, y.volume, id, x.precision);
    } else if (x.nSpin == 1){ // staggered
      cudaBindTexture(0, texHalfSt1, x.v, spinor_bytes); 
      cudaBindTexture(0, texNorm1, x.norm, spinor_bytes/3);    
      cudaBindTexture(0, texHalfSt2, y.v, spinor_bytes); 
      cudaBindTexture(0, texNorm2, y.norm, spinor_bytes/3);    
      cudaBindTexture(0, texHalfSt3, z.v, spinor_bytes); 
      cudaBindTexture(0, texNorm3, z.norm, spinor_bytes/3);    
      cudaBindTexture(0, texHalfSt4, w.v, spinor_bytes); 
      cudaBindTexture(0, texNorm4, w.norm, spinor_bytes/3);    
      cudaBindTexture(0, texHalfSt5, u.v, spinor_bytes); 
      cudaBindTexture(0, texNorm5, u.norm, spinor_bytes/3);    
      float2 a2 = make_float2(real(a), imag(a));
      float2 b2 = make_float2(real(b), imag(b));
      return caxpbypzYmbwcDotProductWYNormYHCuda(a2, b2, (short2*)y.v, (float*)y.norm, 
                                                 (short2*)z.v, (float*)z.norm, (float*)w.norm, (float*)u.norm, 
                                                 y.stride, y.volume, id, x.precision);
    } else {
      errorQuda("%s: nSpin(%d) is not supported\n", __FUNCTION__, x.nSpin);            
    }
  }

}