clover_field_order.h

Go to the documentation of this file.

#ifndef _CLOVER_ORDER_H
#define _CLOVER_ORDER_H
 
#include <register_traits.h>
#include <convert.h>
#include <clover_field.h>
#include <complex_quda.h>
#include <quda_matrix.h>
#include <color_spinor.h>
#include <trove_helper.cuh>
#include <transform_reduce.h>
 
namespace quda {
 
  template <typename Float, typename T>
    struct clover_wrapper {
      T &field;
      const int x_cb;
      const int parity;
      const int chirality;
 
      __device__ __host__ inline clover_wrapper<Float,T>(T &field, int x_cb, int parity, int chirality)
        : field(field), x_cb(x_cb), parity(parity), chirality(chirality) { }
 
      template<typename C>
      __device__ __host__ inline void operator=(const C &a) {
        field.save(a.data, x_cb, parity, chirality);
      }
    };
 
  template <typename T, int N>
    template <typename S>
    __device__ __host__ inline void HMatrix<T,N>::operator=(const clover_wrapper<T,S> &a) {
    a.field.load(data, a.x_cb, a.parity, a.chirality);
  }
 
  template <typename T, int N>
    template <typename S>
    __device__ __host__ inline HMatrix<T,N>::HMatrix(const clover_wrapper<T,S> &a) {
    a.field.load(data, a.x_cb, a.parity, a.chirality);
  }
 
  namespace clover {
 
    template<typename ReduceType, typename Float> struct square_ {
      __host__ __device__ inline ReduceType operator()(const quda::complex<Float> &x)
      { return static_cast<ReduceType>(norm(x)); }
    };
 
    template<typename ReduceType, typename Float> struct abs_ {
      __host__ __device__ inline ReduceType operator()(const quda::complex<Float> &x)
      { return static_cast<ReduceType>(abs(x)); }
    };
 
    template<typename Float, int nColor, int nSpin, QudaCloverFieldOrder order> struct Accessor {
      mutable complex<Float> dummy;
      Accessor(const CloverField &A, bool inverse=false) {
        errorQuda("Not implemented for order %d", order);
      }
 
      __device__ __host__ inline complex<Float>& operator()(int parity, int x, int s_row, int s_col,
                                                            int c_row, int c_col) const {
        return dummy;
      }
 
      template <typename helper, typename reducer>
      __host__ double transform_reduce(QudaFieldLocation location, helper h, double i, reducer r) const
      {
        return 0.0;
      }
    };
 
    template<typename Float, int nColor, int nSpin>
      struct Accessor<Float,nColor,nSpin,QUDA_FLOAT2_CLOVER_ORDER> {
      Float *a;
      int stride;
      size_t offset_cb;
      static constexpr int N = nSpin * nColor / 2;
    Accessor(const CloverField &A, bool inverse=false)
      : a(static_cast<Float*>(const_cast<void*>(A.V(inverse)))), stride(A.Stride()),
        offset_cb(A.Bytes()/(2*sizeof(Float))) { }
 
      __device__ __host__ inline complex<Float> operator()(int parity, int x, int s_row, int s_col, int c_row, int c_col) const {
        // if not in the diagonal chiral block then return 0.0
        if (s_col / 2 != s_row / 2) { return complex<Float>(0.0); }
 
        const int chirality = s_col / 2;
 
        int row = s_row%2 * nColor + c_row;
        int col = s_col%2 * nColor + c_col;
        Float *a_ = a+parity*offset_cb+stride*chirality*N*N;
 
        if (row == col) {
          return 2*a_[ row*stride+x ];
        } else if (col < row) {
          // switch coordinates to count from bottom right instead of top left of matrix
          int k = N*(N-1)/2 - (N-col)*(N-col-1)/2 + row - col - 1;
          complex<Float> *off = reinterpret_cast<complex<Float>*>(a_ + N);
 
          return 2*off[k*stride + x];
        } else {
          // requesting upper triangular so return conjugate transpose
          // switch coordinates to count from bottom right instead of top left of matrix
          int k = N*(N-1)/2 - (N-row)*(N-row-1)/2 + col - row - 1;
          complex<Float> *off = reinterpret_cast<complex<Float>*>(a_ + N);
          return 2*conj(off[k*stride + x]);
        }
 
      }
 
      template <typename helper, typename reducer>
      __host__ double transform_reduce(QudaFieldLocation location, helper h, double init, reducer r) const
      {
        // just use offset_cb, since factor of two from parity is equivalent to complexity
        double result = ::quda::transform_reduce(location, reinterpret_cast<complex<Float> *>(a), offset_cb, h, init, r);
        return 2.0 * result; // factor of two is normalization
      }
    };
 
    template<int N>
      __device__ __host__ inline int indexFloatN(int k, int stride, int x) {
      int j = k / N;
      int i = k % N;
      return (j*stride+x)*N + i;
    };
 
    template<typename Float, int nColor, int nSpin>
      struct Accessor<Float,nColor,nSpin,QUDA_FLOAT4_CLOVER_ORDER> {
      Float *a;
      int stride;
      size_t offset_cb;
      static constexpr int N = nSpin * nColor / 2;
    Accessor(const CloverField &A, bool inverse=false)
      : a(static_cast<Float*>(const_cast<void*>(A.V(inverse)))), stride(A.Stride()),
        offset_cb(A.Bytes()/(2*sizeof(Float))) { }
 
      __device__ __host__ inline complex<Float> operator()(int parity, int x, int s_row, int s_col, int c_row, int c_col) const {
        // if not in the diagonal chiral block then return 0.0
        if (s_col / 2 != s_row / 2) { return complex<Float>(0.0); }
 
        const int chirality = s_col / 2;
 
        int row = s_row%2 * nColor + c_row;
        int col = s_col%2 * nColor + c_col;
        Float *a_ = a+parity*offset_cb+stride*chirality*N*N;
 
        if (row == col) {
          return 2*a_[ indexFloatN<QUDA_FLOAT4_CLOVER_ORDER>(row, stride, x) ];
        } else if (col < row) {
          // switch coordinates to count from bottom right instead of top left of matrix
          int k = N*(N-1)/2 - (N-col)*(N-col-1)/2 + row - col - 1;
          int idx = N + 2*k;
 
          return 2*complex<Float>(a_[ indexFloatN<QUDA_FLOAT4_CLOVER_ORDER>(idx+0,stride,x) ],
                                  a_[ indexFloatN<QUDA_FLOAT4_CLOVER_ORDER>(idx+1,stride,x) ]);
        } else {
          // requesting upper triangular so return conjugate transpose
          // switch coordinates to count from bottom right instead of top left of matrix
          int k = N*(N-1)/2 - (N-row)*(N-row-1)/2 + col - row - 1;
          int idx = N + 2*k;
 
          return 2*complex<Float>( a_[ indexFloatN<QUDA_FLOAT4_CLOVER_ORDER>(idx+0,stride,x) ],
                                  -a_[ indexFloatN<QUDA_FLOAT4_CLOVER_ORDER>(idx+1,stride,x) ]);
        }
 
      }
 
      template <typename helper, typename reducer>
      __host__ double transform_reduce(QudaFieldLocation location, helper h, double init, reducer r) const
      {
        // just use offset_cb, since factor of two from parity is equivalent to complexity
        double result = ::quda::transform_reduce(location, reinterpret_cast<complex<Float> *>(a), offset_cb, h, init, r);
        return 2.0 * result; // factor of two is normalization
      }
    };
 
    template<typename Float, int nColor, int nSpin> 
      struct Accessor<Float,nColor,nSpin,QUDA_PACKED_CLOVER_ORDER> { 
      Float *a[2];
      const int N = nSpin * nColor / 2;
      complex<Float> zero;
      Accessor(const CloverField &A, bool inverse=false) {
        // even
        a[0] = static_cast<Float*>(const_cast<void*>(A.V(inverse)));
        // odd
        a[1] = static_cast<Float*>(const_cast<void*>(A.V(inverse))) + A.Bytes()/(2*sizeof(Float));
        zero = complex<Float>(0.0,0.0);
      }
 
      __device__ __host__ inline complex<Float> operator()(int parity, int x, int s_row, int s_col, int c_row, int c_col) const {
        // if not in the diagonal chiral block then return 0.0
        if (s_col / 2 != s_row / 2) { return zero; }
 
        const int chirality = s_col / 2;
 
        unsigned int row = s_row%2 * nColor + c_row;
        unsigned int col = s_col%2 * nColor + c_col;
 
        if (row == col) {
          complex<Float> tmp = a[parity][(x*2 + chirality)*N*N + row];
          return tmp;
        } else if (col < row) {
          // switch coordinates to count from bottom right instead of top left of matrix
          int k = N*(N-1)/2 - (N-col)*(N-col-1)/2 + row - col - 1;
          int idx = (x*2 + chirality)*N*N + N + 2*k;
          return complex<Float>(a[parity][idx], a[parity][idx+1]);
        } else {
          // switch coordinates to count from bottom right instead of top left of matrix
          int k = N*(N-1)/2 - (N-row)*(N-row-1)/2 + col - row - 1;
          int idx = (x*2 + chirality)*N*N + N + 2*k;
          return complex<Float>(a[parity][idx], -a[parity][idx+1]);
        }
      }
 
      template <typename helper, typename reducer>
      __host__ double transform_reduce(QudaFieldLocation location, helper h, double init, reducer r) const
      {
        errorQuda("Not implemented");
        return 0.0;
      }
    };
 
    /*
      FIXME the below is the old optimization used for reading the
      clover field, making use of the symmetry to reduce the number of
      reads.
 
#define READ_CLOVER2_DOUBLE_STR(clover_, chi)                           \
    double2 C0, C1, C2, C3, C4, C5, C6, C7, C8, C9;                       \
    double2 C10, C11, C12, C13, C14, C15, C16, C17;                       \
    double2* clover = (double2*)clover_;                                  \
    load_streaming_double2(C0, &clover[sid + (18*chi+0)*param.cl_stride]); \
    load_streaming_double2(C1, &clover[sid + (18*chi+1)*param.cl_stride]); \
    double diag = 0.5*(C0.x + C1.y);                                      \
    double diag_inv = 1.0/diag;                                           \
    C2 = make_double2(diag*(2-C0.y*diag_inv), diag*(2-C1.x*diag_inv));    \
    load_streaming_double2(C3, &clover[sid + (18*chi+3)*param.cl_stride]);        \
    load_streaming_double2(C4, &clover[sid + (18*chi+4)*param.cl_stride]);        \
    load_streaming_double2(C5, &clover[sid + (18*chi+5)*param.cl_stride]);        \
    load_streaming_double2(C6, &clover[sid + (18*chi+6)*param.cl_stride]);        \
    load_streaming_double2(C7, &clover[sid + (18*chi+7)*param.cl_stride]);        \
    load_streaming_double2(C8, &clover[sid + (18*chi+8)*param.cl_stride]);        \
    load_streaming_double2(C9, &clover[sid + (18*chi+9)*param.cl_stride]);        \
    load_streaming_double2(C10, &clover[sid + (18*chi+10)*param.cl_stride]);      \
    load_streaming_double2(C11, &clover[sid + (18*chi+11)*param.cl_stride]);      \
    load_streaming_double2(C12, &clover[sid + (18*chi+12)*param.cl_stride]);      \
    load_streaming_double2(C13, &clover[sid + (18*chi+13)*param.cl_stride]);      \
    load_streaming_double2(C14, &clover[sid + (18*chi+14)*param.cl_stride]); \
    C15 = make_double2(-C3.x,-C3.y);                                      \
    C16 = make_double2(-C4.x,-C4.y);                                      \
    C17 = make_double2(-C8.x,-C8.y);                                      \
    */
 
    template <typename Float, int nColor, int nSpin, QudaCloverFieldOrder order>
      struct FieldOrder {
 
      static constexpr bool supports_ghost_zone = false;
 
    protected:
      CloverField &A;
      const int volumeCB;
      const Accessor<Float, nColor, nSpin, order> accessor;
      bool inverse;
      const QudaFieldLocation location;
 
    public:
      FieldOrder(CloverField &A, bool inverse=false)
      : A(A), volumeCB(A.VolumeCB()), accessor(A,inverse), inverse(inverse), location(A.Location())
        { }
        
        CloverField& Field() { return A; }
        
        __device__ __host__ inline const complex<Float> operator()(int parity, int x, int s_row,
                                                                   int s_col, int c_row, int c_col) const {
          return accessor(parity, x, s_row, s_col, c_row, c_col);
        }
        
        __device__ __host__ inline complex<Float> operator()(int dummy, int parity, int x, int s_row,
                                                             int s_col, int c_row, int c_col) const {
          return accessor(parity,x,s_row,s_col,c_row,c_col);
        }
 
        /*
        __device__ __host__ inline complex<Float>& operator()(int parity, int x, int s_row,
                                                             int s_col, int c_row, int c_col) {
          //errorQuda("Clover accessor not implemented as a lvalue");
          return accessor(parity, x, s_row, s_col, c_row, c_col);
          }
        */
        
        __device__ __host__ inline int Ncolor() const { return nColor; }
 
        __device__ __host__ inline int Volume() const { return 2*volumeCB; }
 
        __device__ __host__ inline int VolumeCB() const { return volumeCB; }
 
        size_t Bytes() const {
          constexpr int n = (nSpin * nColor) / 2;
          constexpr int chiral_block = n * n / 2;
          return static_cast<size_t>(volumeCB) * chiral_block * 2ll * 2ll * sizeof(Float); // 2 from complex, 2 from chirality
        }
 
        __host__ double norm1(int dim=-1, bool global=true) const {
          double nrm1 = accessor.transform_reduce(location, abs_<double, Float>(), 0.0, plus<double>());
          if (global) comm_allreduce(&nrm1);
          return nrm1;
        }
 
        __host__ double norm2(int dim=-1, bool global=true) const {
          double nrm2 = accessor.transform_reduce(location, square_<double, Float>(), 0.0, plus<double>());
          if (global) comm_allreduce(&nrm2);
          return nrm2;
        }
 
        __host__ double abs_max(int dim=-1, bool global=true) const {
          double absmax = accessor.transform_reduce(location, abs_<Float, Float>(), 0.0, maximum<Float>());
          if (global) comm_allreduce_max(&absmax);
          return absmax;
        }
 
        __host__ double abs_min(int dim=-1, bool global=true) const {
          double absmax = accessor.transform_reduce(location, abs_<Float, Float>(), std::numeric_limits<double>::max(),
                                                    minimum<Float>());
          if (global) comm_allreduce_min(&absmax);
          return absmax;
        }
      };
 
    template <typename Float, int length, int N, bool add_rho=false, bool huge_alloc=false>
    struct FloatNOrder {
      using Accessor = FloatNOrder<Float, length, N, add_rho, huge_alloc>;
      using real = typename mapper<Float>::type;
      typedef typename VectorType<Float, N>::type Vector;
      typedef typename AllocType<huge_alloc>::type AllocInt;
      typedef float norm_type;
      static const int M = length / (N * 2); // number of short vectors per chiral block
      static const int block = length / 2;   // chiral block size
      Float *clover;
      norm_type *norm;
      const AllocInt offset; // offset can be 32-bit or 64-bit
      const AllocInt norm_offset;
        const int volumeCB;
        const int stride;
 
        const bool twisted;
        const real mu2;
        const real rho;
 
        size_t bytes;
        size_t norm_bytes;
        void *backup_h; 
        void *backup_norm_h; 
 
        FloatNOrder(const CloverField &clover, bool is_inverse, Float *clover_ = 0, norm_type *norm_ = 0,
                    bool override = false) :
          offset(clover.Bytes() / (2 * sizeof(Float) * N)),
          norm_offset(clover.NormBytes() / (2 * sizeof(norm_type))),
          volumeCB(clover.VolumeCB()),
          stride(clover.Stride()),
          twisted(clover.Twisted()),
          mu2(clover.Mu2()),
          rho(clover.Rho()),
          bytes(clover.Bytes()),
          norm_bytes(clover.NormBytes()),
          backup_h(nullptr),
          backup_norm_h(nullptr)
        {
          if (clover.Order() != N) {
            errorQuda("Invalid clover order %d for FloatN (N=%d) accessor", clover.Order(), N);
          }
          this->clover = clover_ ? clover_ : (Float *)(clover.V(is_inverse));
          this->norm = norm_ ? norm_ : (norm_type *)(clover.Norm(is_inverse));
        }
 
        bool Twisted() const { return twisted; }
        real Mu2() const { return mu2; }
 
        __device__ __host__ inline clover_wrapper<real, Accessor> operator()(int x_cb, int parity, int chirality)
        {
          return clover_wrapper<real, Accessor>(*this, x_cb, parity, chirality);
        }
 
        __device__ __host__ inline const clover_wrapper<real, Accessor> operator()(
            int x_cb, int parity, int chirality) const
        {
          return clover_wrapper<real, Accessor>(const_cast<Accessor &>(*this), x_cb, parity, chirality);
        }
 
        __device__ __host__ inline void load(real v[block], int x, int parity, int chirality) const
        {
          norm_type nrm;
          if (isFixed<Float>::value) { nrm = vector_load<float>(norm, parity * norm_offset + chirality * stride + x); }
 
#pragma unroll
          for (int i=0; i<M; i++) {
            // first load from memory
            Vector vecTmp = vector_load<Vector>(clover, parity * offset + x + stride * (chirality * M + i));
            // second do scalar copy converting into register type
#pragma unroll
            for (int j = 0; j < N; j++) { copy_and_scale(v[i * N + j], reinterpret_cast<Float *>(&vecTmp)[j], nrm); }
          }
 
          if (add_rho) for (int i=0; i<6; i++) v[i] += rho;
        }
  
        __device__ __host__ inline void save(const real v[block], int x, int parity, int chirality)
        {
          real tmp[block];
 
          // find the norm of each chiral block
          if (isFixed<Float>::value) {
            norm_type scale = 0.0;
#pragma unroll
            for (int i = 0; i < block; i++) scale = fabsf((norm_type)v[i]) > scale ? fabsf((norm_type)v[i]) : scale;
            norm[parity*norm_offset + chirality*stride + x] = scale;
 
#ifdef __CUDA_ARCH__
            real scale_inv = __fdividef(fixedMaxValue<Float>::value, scale);
#else
            real scale_inv = fixedMaxValue<Float>::value / scale;
#endif
#pragma unroll
            for (int i = 0; i < block; i++) tmp[i] = v[i] * scale_inv;
          } else {
#pragma unroll
            for (int i = 0; i < block; i++) tmp[i] = v[i];
          }
 
#pragma unroll
          for (int i = 0; i < M; i++) {
            Vector vecTmp;
            // first do scalar copy converting into storage type
            for (int j = 0; j < N; j++) copy_scaled(reinterpret_cast<Float *>(&vecTmp)[j], tmp[i * N + j]);
            // second do vectorized copy into memory
            vector_store(clover, parity * offset + x + stride * (chirality * M + i), vecTmp);
          }
        }
 
        __device__ __host__ inline void load(real v[length], int x, int parity) const {
#pragma unroll
          for (int chirality = 0; chirality < 2; chirality++) load(&v[chirality * block], x, parity, chirality);
        }
 
        __device__ __host__ inline void save(const real v[length], int x, int parity) {
#pragma unroll
          for (int chirality = 0; chirality < 2; chirality++) save(&v[chirality * block], x, parity, chirality);
        }
 
        void save() {
          if (backup_h) errorQuda("Already allocated host backup");
          backup_h = safe_malloc(bytes);
          qudaMemcpy(backup_h, clover, bytes, cudaMemcpyDeviceToHost);
          if (norm_bytes) {
            backup_norm_h = safe_malloc(norm_bytes);
            qudaMemcpy(backup_norm_h, norm, norm_bytes, cudaMemcpyDeviceToHost);
          }
        }
 
        void load()
        {
          qudaMemcpy(clover, backup_h, bytes, cudaMemcpyHostToDevice);
          host_free(backup_h);
          backup_h = nullptr;
          if (norm_bytes) {
            qudaMemcpy(norm, backup_norm_h, norm_bytes, cudaMemcpyHostToDevice);
            host_free(backup_norm_h);
            backup_norm_h = nullptr;
          }
        }
 
        size_t Bytes() const
        {
          size_t bytes = length * sizeof(Float);
          if (isFixed<Float>::value) bytes += 2 * sizeof(norm_type);
          return bytes;
        }
      };
 
    template <typename real, int length> struct S { real v[length]; };
 
    template <typename Float, int length>
      struct QDPOrder {
        typedef typename mapper<Float>::type RegType;
        Float *clover;
        const int volumeCB;
        const int stride;
        const int offset;
 
        const bool twisted;
        const Float mu2;
 
      QDPOrder(const CloverField &clover, bool inverse, Float *clover_=0) 
      : volumeCB(clover.VolumeCB()), stride(volumeCB), offset(clover.Bytes()/(2*sizeof(Float))),
        twisted(clover.Twisted()), mu2(clover.Mu2()) {
        if (clover.Order() != QUDA_PACKED_CLOVER_ORDER) {
          errorQuda("Invalid clover order %d for this accessor", clover.Order());
        }
        this->clover = clover_ ? clover_ : (Float *)(clover.V(inverse));
      }
 
        bool  Twisted() const   {return twisted;}
        Float Mu2()     const   {return mu2;}
 
        __device__ __host__ inline void load(RegType v[length], int x, int parity) const {
          // factor of 0.5 comes from basis change
#if defined( __CUDA_ARCH__) && !defined(DISABLE_TROVE)
          typedef S<Float,length> structure;
          trove::coalesced_ptr<structure> clover_((structure*)clover);
          structure v_ = clover_[parity*volumeCB + x];
          for (int i=0; i<length; i++) v[i] = 0.5*(RegType)v_.v[i];
#else
          for (int i=0; i<length; i++) v[i] = 0.5*clover[parity*offset + x*length+i];
#endif
        }
  
        __device__ __host__ inline void save(const RegType v[length], int x, int parity) {
#if defined( __CUDA_ARCH__) && !defined(DISABLE_TROVE)
          typedef S<Float,length> structure;
          trove::coalesced_ptr<structure> clover_((structure*)clover);
          structure v_;
          for (int i=0; i<length; i++) v_.v[i] = 2.0*(Float)v[i];
          clover_[parity*volumeCB + x] = v_;
#else
          for (int i=0; i<length; i++) clover[parity*offset + x*length+i] = 2.0*v[i];
#endif
        }
 
        size_t Bytes() const { return length*sizeof(Float); }
      };
 
    template <typename Float, int length>
      struct QDPJITOrder {
        typedef typename mapper<Float>::type RegType;
        Float *diag;       
        Float *offdiag;   
        const int volumeCB;
        const int stride;
 
        const bool twisted;
        const Float mu2;
 
      QDPJITOrder(const CloverField &clover, bool inverse, Float *clover_=0) 
      : volumeCB(clover.VolumeCB()), stride(volumeCB), twisted(clover.Twisted()), mu2(clover.Mu2()) {
        if (clover.Order() != QUDA_QDPJIT_CLOVER_ORDER) {
          errorQuda("Invalid clover order %d for this accessor", clover.Order());
        }
        offdiag = clover_ ? ((Float **)clover_)[0] : ((Float **)clover.V(inverse))[0];
        diag = clover_ ? ((Float **)clover_)[1] : ((Float **)clover.V(inverse))[1];
      }
        
      bool  Twisted()   const   {return twisted;}
      Float Mu2()       const   {return mu2;}
 
        __device__ __host__ inline void load(RegType v[length], int x, int parity) const {
          // the factor of 0.5 comes from a basis change
          for (int chirality=0; chirality<2; chirality++) {
            // set diagonal elements
            for (int i=0; i<6; i++) {
              v[chirality*36 + i] = 0.5*diag[((i*2 + chirality)*2 + parity)*volumeCB + x];
            }
 
            // the off diagonal elements
            for (int i=0; i<30; i++) {
              int z = i%2;
              int off = i/2;
              const int idtab[15]={0,1,3,6,10,2,4,7,11,5,8,12,9,13,14};
              v[chirality*36 + 6 + i] = 0.5*offdiag[(((z*15 + idtab[off])*2 + chirality)*2 + parity)*volumeCB + x];
            }
 
          }
        }
  
        __device__ __host__ inline void save(const RegType v[length], int x, int parity) {
          // the factor of 2.0 comes from undoing the basis change
          for (int chirality=0; chirality<2; chirality++) {
            // set diagonal elements
            for (int i=0; i<6; i++) {
              diag[((i*2 + chirality)*2 + parity)*volumeCB + x] = 2.0*v[chirality*36 + i];
            }
 
            // the off diagonal elements
            for (int i=0; i<30; i++) {
              int z = i%2;
              int off = i/2;
              const int idtab[15]={0,1,3,6,10,2,4,7,11,5,8,12,9,13,14};
              offdiag[(((z*15 + idtab[off])*2 + chirality)*2 + parity)*volumeCB + x] = 2.0*v[chirality*36 + 6 + i];
            }
          }
        }
        
        size_t Bytes() const { return length*sizeof(Float); }
      };
      
 
    template <typename Float, int length>
      struct BQCDOrder {
        typedef typename mapper<Float>::type RegType;
        Float *clover[2];
        const int volumeCB;
        const int stride;
 
        const bool twisted;
        const Float mu2;
 
      BQCDOrder(const CloverField &clover, bool inverse, Float *clover_=0) 
      : volumeCB(clover.Stride()), stride(volumeCB), twisted(clover.Twisted()), mu2(clover.Mu2()) {
        if (clover.Order() != QUDA_BQCD_CLOVER_ORDER) {
          errorQuda("Invalid clover order %d for this accessor", clover.Order());
        }
        this->clover[0] = clover_ ? clover_ : (Float *)(clover.V(inverse));
        this->clover[1] = (Float *)((char *)this->clover[0] + clover.Bytes() / 2);
      }
 
 
        bool  Twisted() const   {return twisted;}
        Float Mu2()     const   {return mu2;}
 
        __device__ __host__ inline void load(RegType v[length], int x, int parity) const {
          int bq[36] = { 21, 32, 33, 0,  1, 20,                   // diagonal
                         28, 29, 30, 31, 6, 7,  14, 15, 22, 23,   // column 1  6
                         34, 35, 8, 9, 16, 17, 24, 25,            // column 2  16
                         10, 11, 18, 19, 26, 27,                  // column 3  24
                         2,  3,  4,  5,                           // column 4  30
                         12, 13};
          
          // flip the sign of the imaginary components
          int sign[36];
          for (int i=0; i<6; i++) sign[i] = 1;
          for (int i=6; i<36; i+=2) {
            if ( (i >= 10 && i<= 15) || (i >= 18 && i <= 29) )  { sign[i] = -1; sign[i+1] = -1; }
            else { sign[i] = 1; sign[i+1] = -1; }
          }
        
          const int M=length/2;
          for (int chirality=0; chirality<2; chirality++) 
            for (int i=0; i<M; i++) 
              v[chirality*M+i] = sign[i] * clover[parity][x*length+chirality*M+bq[i]];
        
        }
  
        // FIXME implement the save routine for BQCD ordered fields
        __device__ __host__ inline void save(RegType v[length], int x, int parity) {
 
        };
 
        size_t Bytes() const { return length*sizeof(Float); }
      };
 
  } // namespace clover
 
  // Use traits to reduce the template explosion
  template<typename Float,int N=72, bool add_rho=false> struct clover_mapper { };
 
  // double precision uses Float2
  template<int N, bool add_rho> struct clover_mapper<double,N,add_rho> { typedef clover::FloatNOrder<double, N, 2, add_rho> type; };
 
  // single precision uses Float4
  template<int N, bool add_rho> struct clover_mapper<float,N,add_rho> { typedef clover::FloatNOrder<float, N, 4, add_rho> type; };
 
  // half precision uses Float4
  template<int N, bool add_rho> struct clover_mapper<short,N,add_rho> { typedef clover::FloatNOrder<short, N, 4, add_rho> type; };
 
  // quarter precision uses Float4
  template <int N, bool add_rho> struct clover_mapper<int8_t, N, add_rho> {
    typedef clover::FloatNOrder<int8_t, N, 4, add_rho> type;
  };
 
} // namespace quda
 
#endif //_CLOVER_ORDER_H