color_spinor.h

Go to the documentation of this file.

#pragma once

#include <complex_quda.h>
#include <quda_matrix.h>

namespace quda {

  template<typename Float, typename T> struct colorspinor_wrapper;
  template<typename Float, typename T> struct colorspinor_ghost_wrapper;

  template <typename Float, int Nc, int Ns>
    struct ColorSpinor {

      complex<Float> data[Nc*Ns];

      __device__ __host__ inline ColorSpinor<Float, Nc, Ns>() {
#pragma unroll
  for (int i=0; i<Nc*Ns; i++) {
    data[i] = 0;
  }
      }

      __device__ __host__ inline ColorSpinor<Float, Nc, Ns>(const ColorSpinor<Float, Nc, Ns> &a) {
#pragma unroll
  for (int i=0; i<Nc*Ns; i++) {
    data[i] = a.data[i];
  }
      }

      __device__ __host__ inline ColorSpinor<Float, Nc, Ns>& operator=(const ColorSpinor<Float, Nc, Ns> &a) {
  if (this != &a) {
#pragma unroll
    for (int i=0; i<Nc*Ns; i++) {
      data[i] = a.data[i];
    }
  }
  return *this;
      }

      __device__ __host__ inline ColorSpinor<Float, Nc, Ns>& operator+=(const ColorSpinor<Float, Nc, Ns> &a) {
  if (this != &a) {
#pragma unroll
    for (int i=0; i<Nc*Ns; i++) {
      data[i] += a.data[i];
    }
  }
  return *this;
      }

      template<typename S>
      __device__ __host__ inline ColorSpinor<Float, Nc, Ns>(const colorspinor_wrapper<Float, S> &s);

      template<typename S>
      __device__ __host__ inline void operator=(const colorspinor_wrapper<Float, S> &s);

      template<typename S>
      __device__ __host__ inline ColorSpinor<Float, Nc, Ns>(const colorspinor_ghost_wrapper<Float, S> &s);

      template<typename S>
      __device__ __host__ inline void operator=(const colorspinor_ghost_wrapper<Float, S> &s);

      __device__ __host__ inline complex<Float>& operator()(int s, int c) { return data[s*Nc + c]; }

      __device__ __host__ inline const complex<Float>& operator()(int s, int c) const { return data[s*Nc + c]; }

      __device__ __host__ inline complex<Float>& operator()(int idx) { return data[idx]; }

      __device__ __host__ inline const complex<Float>& operator()(int idx) const { return data[idx]; }

    };

  template <typename Float, int Nc>
    struct ColorSpinor<Float, Nc, 4> {
    static const int Ns = 4;
    complex<Float> data[Nc*4];

    __device__ __host__ inline ColorSpinor<Float, Nc, 4>() {
#pragma unroll
      for (int i=0; i<Nc*Ns; i++) {
  data[i] = 0;
      }      
    }

    __device__ __host__ inline ColorSpinor<Float, Nc, 4>(const ColorSpinor<Float, Nc, 4> &a) {
#pragma unroll
      for (int i=0; i<Nc*Ns; i++) {
  data[i] = a.data[i];
      }      
    }

    __device__ __host__ inline ColorSpinor<Float, Nc, 4>& operator=(const ColorSpinor<Float, Nc, 4> &a) {
      if (this != &a) {
#pragma unroll
  for (int i=0; i<Nc*Ns; i++) {
    data[i] = a.data[i];
  }
      }
      return *this;
    }

    __device__ __host__ inline ColorSpinor<Float, Nc, 4>& operator+=(const ColorSpinor<Float, Nc, 4> &a) {
      if (this != &a) {
#pragma unroll
  for (int i=0; i<Nc*Ns; i++) {
    data[i] += a.data[i];
  }
      }
      return *this;
    }

    __device__ __host__ inline ColorSpinor<Float,Nc,4> gamma(int dim) {
      ColorSpinor<Float,Nc,4> a;
      complex<Float> j(0.0,1.0);

      switch (dim) {
      case 0: // x dimension
#pragma unroll
  for (int i=0; i<Nc; i++) {
    a(0,i) =  j*(*this)(3,i);
    a(1,i) =  j*(*this)(2,i);
    a(2,i) = -j*(*this)(1,i);
    a(3,i) = -j*(*this)(0,i);
  }
  break;
      case 1: // y dimension
#pragma unroll
  for (int i=0; i<Nc; i++) {
    a(0,i) =  (*this)(3,i);
    a(1,i) = -(*this)(2,i);
    a(2,i) = -(*this)(1,i);
    a(3,i) =  (*this)(0,i);
  }
  break;
      case 2: // z dimension
#pragma unroll
  for (int i=0; i<Nc; i++) {
    a(0,i) =  j*(*this)(2,i);
    a(1,i) = -j*(*this)(3,i);
    a(2,i) = -j*(*this)(0,i);
    a(3,i) =  j*(*this)(1,i);
  }
  break;
      case 3: // t dimension
#pragma unroll
  for (int i=0; i<Nc; i++) {
    a(0,i) =  (*this)(0,i);
    a(1,i) =  (*this)(1,i);
    a(2,i) = -(*this)(2,i);
    a(3,i) = -(*this)(3,i);
  }
  break;
      case 4: // gamma_5
#pragma unroll
  for (int i=0; i<Nc; i++) {
    a(0,i) = (*this)(2,i);
    a(1,i) = (*this)(3,i);
    a(2,i) = (*this)(0,i);
    a(3,i) = (*this)(1,i);
  }
  break;
      }

      return a;
    }

    __device__ __host__ inline ColorSpinor<Float,Nc,4> igamma(int dim) {
      ColorSpinor<Float,Nc,4> a;
      complex<Float> j(0.0,1.0);

      switch (dim) {
      case 0: // x dimension
#pragma unroll
  for (int i=0; i<Nc; i++) {
    a(0,i) = -(*this)(3,i);
    a(1,i) = -(*this)(2,i);
    a(2,i) =  (*this)(1,i);
    a(3,i) =  (*this)(0,i);
  }
  break;
      case 1: // y dimension
#pragma unroll
  for (int i=0; i<Nc; i++) {
    a(0,i) =  j*(*this)(3,i);
    a(1,i) = -j*(*this)(2,i);
    a(2,i) = -j*(*this)(1,i);
    a(3,i) =  j*(*this)(0,i);
  }
  break;
      case 2: // z dimension
#pragma unroll
  for (int i=0; i<Nc; i++) {
    a(0,i) = -(*this)(2,i);
    a(1,i) =  (*this)(3,i);
    a(2,i) =  (*this)(0,i);
    a(3,i) = -(*this)(1,i);
  }
  break;
      case 3: // t dimension
#pragma unroll
  for (int i=0; i<Nc; i++) {
    a(0,i) =  j*(*this)(0,i);
    a(1,i) =  j*(*this)(1,i);
    a(2,i) = -j*(*this)(2,i);
    a(3,i) = -j*(*this)(3,i);
  }
  break;
      case 4: // gamma_5
#pragma unroll
  for (int i=0; i<Nc; i++) {
    a(0,i) = complex<Float>(-(*this)(2,i).imag(), (*this)(2,i).real());
    a(1,i) = complex<Float>(-(*this)(3,i).imag(), (*this)(3,i).real());
    a(2,i) = complex<Float>(-(*this)(0,i).imag(), (*this)(0,i).real());
    a(3,i) = complex<Float>(-(*this)(1,i).imag(), (*this)(1,i).real());
  }
  break;
      }

      return a;
    }

    __device__ __host__ inline ColorSpinor<Float,Nc,2> chiral_project(int chirality) const {
      ColorSpinor<Float,Nc,2> proj;
#pragma unroll
      for (int s=0; s<Ns/2; s++) {
#pragma unroll
  for (int c=0; c<Nc; c++) {
    proj(s,c) = (*this)(chirality*Ns/2+s,c);
  }
      }
      return proj;
    }

    __device__ __host__ inline ColorSpinor<Float,Nc,2> project(int dim, int sign) {
      ColorSpinor<Float,Nc,2> proj;
      complex<Float> j(0.0,1.0);
      
      switch (dim) {
      case 0: // x dimension
  switch (sign) {
  case 1: // positive projector
#pragma unroll
    for (int i=0; i<Nc; i++) {
      proj(0,i) = (*this)(0,i) + j * (*this)(3,i);
      proj(1,i) = (*this)(1,i) + j * (*this)(2,i);
    }
    break;
  case -1: // negative projector
#pragma unroll
    for (int i=0; i<Nc; i++) {
      proj(0,i) = (*this)(0,i) - j * (*this)(3,i);
      proj(1,i) = (*this)(1,i) - j * (*this)(2,i);
    }
    break;
  }
  break;
      case 1: // y dimension
  switch (sign) {
  case 1: // positive projector
#pragma unroll
    for (int i=0; i<Nc; i++) {
      proj(0,i) = (*this)(0,i) + (*this)(3,i);
      proj(1,i) = (*this)(1,i) - (*this)(2,i);
    }
    break;
  case -1: // negative projector
#pragma unroll
    for (int i=0; i<Nc; i++) {
      proj(0,i) = (*this)(0,i) - (*this)(3,i);
      proj(1,i) = (*this)(1,i) + (*this)(2,i);
    }
    break;
  }
        break;
      case 2: // z dimension
  switch (sign) {
  case 1: // positive projector
#pragma unroll
    for (int i=0; i<Nc; i++) {
      proj(0,i) = (*this)(0,i) + j * (*this)(2,i);
      proj(1,i) = (*this)(1,i) - j * (*this)(3,i);
    }
    break;
  case -1: // negative projector
#pragma unroll
    for (int i=0; i<Nc; i++) {
      proj(0,i) = (*this)(0,i) - j * (*this)(2,i);
      proj(1,i) = (*this)(1,i) + j * (*this)(3,i);
    }
    break;
  }
  break;
      case 3: // t dimension
  switch (sign) {
  case 1: // positive projector
#pragma unroll
    for (int i=0; i<Nc; i++) {
      proj(0,i) = 2*(*this)(0,i);
      proj(1,i) = 2*(*this)(1,i);
    }
    break;
  case -1: // negative projector
#pragma unroll
    for (int i=0; i<Nc; i++) {
      proj(0,i) = 2*(*this)(2,i);
      proj(1,i) = 2*(*this)(3,i);
    }
    break;
  }
  break;
      }

      return proj;
    }

    __device__ __host__ inline ColorSpinor<Float,Nc,4> sigma(int mu, int nu) {
      ColorSpinor<Float,Nc,4> a;
      ColorSpinor<Float,Nc,4> &b = *this;
      complex<Float> j(0.0,1.0);

      switch(mu) {
      case 0:
  switch(nu) {
  case 1:
#pragma unroll
    for (int i=0; i<Nc; i++) {
      a(0,i) =  j*b(0,i);
      a(1,i) = -j*b(1,i);
      a(2,i) =  j*b(2,i);
      a(3,i) = -j*b(3,i);
    }
    break;
  case 2:
#pragma unroll
    for (int i=0; i<Nc; i++) {
      a(0,i) = -b(1,i);
      a(1,i) =  b(0,i);
      a(2,i) = -b(3,i);
      a(3,i) =  b(2,i);
    }
    break;
  case 3:
#pragma unroll
    for (int i=0; i<Nc; i++) {
      a(0,i) = -j*b(3,i);
      a(1,i) = -j*b(2,i);
      a(2,i) = -j*b(1,i);
      a(3,i) = -j*b(0,i);
    }
    break;
  }
  break;
      case 1:
  switch(nu) {
  case 0:
#pragma unroll
    for (int i=0; i<Nc; i++) {
      a(0,i) = -j*b(0,i);
      a(1,i) =  j*b(1,i);
      a(2,i) = -j*b(2,i);
      a(3,i) =  j*b(3,i);
    }
    break;
  case 2:
#pragma unroll
    for (int i=0; i<Nc; i++) {
      a(0,i) = j*b(1,i);
      a(1,i) = j*b(0,i);
      a(2,i) = j*b(3,i);
      a(3,i) = j*b(2,i);
    }
    break;
  case 3:
#pragma unroll
    for (int i=0; i<Nc; i++) {
      a(0,i) = -b(3,i);
      a(1,i) =  b(2,i);
      a(2,i) = -b(1,i);
      a(3,i) =  b(0,i);
    }
    break;
  }
  break;
      case 2:
  switch(nu) {
  case 0:
#pragma unroll
    for (int i=0; i<Nc; i++) {
      a(0,i) =  b(1,i);
      a(1,i) = -b(0,i);
      a(2,i) =  b(3,i);
      a(3,i) = -b(2,i);
    }
    break;
  case 1:
#pragma unroll
    for (int i=0; i<Nc; i++) {
      a(0,i) = -j*b(1,i);
      a(1,i) = -j*b(0,i);
      a(2,i) = -j*b(3,i);
      a(3,i) = -j*b(2,i);
    }
    break;
  case 3:
#pragma unroll
    for (int i=0; i<Nc; i++) {
      a(0,i) = -j*b(2,i);
      a(1,i) =  j*b(3,i);
      a(2,i) = -j*b(0,i);
      a(3,i) =  j*b(1,i);
    }
    break;
  }
  break;
      case 3:
  switch(nu) {
  case 0:
#pragma unroll
    for (int i=0; i<Nc; i++) {
      a(0,i) = j*b(3,i);
      a(1,i) = j*b(2,i);
      a(2,i) = j*b(1,i);
      a(3,i) = j*b(0,i);
    }
    break;
  case 1:
#pragma unroll
    for (int i=0; i<Nc; i++) {
      a(0,i) =  b(3,i);
      a(1,i) = -b(2,i);
      a(2,i) =  b(1,i);
      a(3,i) = -b(0,i);
    }
    break;
  case 2:
#pragma unroll
    for (int i=0; i<Nc; i++) {
      a(0,i) =  j*b(2,i);
      a(1,i) = -j*b(3,i);
      a(2,i) =  j*b(0,i);
      a(3,i) = -j*b(1,i);
    }
    break;
  }
  break;
      }

      return a;
    }


    __device__ __host__ inline complex<Float>& operator()(int s, int c) { return data[s*Nc + c]; }

    __device__ __host__ inline const complex<Float>& operator()(int s, int c) const { return data[s*Nc + c]; }

    __device__ __host__ inline complex<Float>& operator()(int idx) { return data[idx]; }

    __device__ __host__ inline const complex<Float>& operator()(int idx) const { return data[idx]; }

    template<typename S>
      __device__ __host__ inline ColorSpinor<Float, Nc, Ns>(const colorspinor_wrapper<Float, S> &s);

    template<typename S>
      __device__ __host__ inline void operator=(const colorspinor_wrapper<Float, S> &s);

    template<typename S>
      __device__ __host__ inline ColorSpinor<Float, Nc, Ns>(const colorspinor_ghost_wrapper<Float, S> &s);

    template<typename S>
      __device__ __host__ inline void operator=(const colorspinor_ghost_wrapper<Float, S> &s);

    __device__ __host__ inline void toNonRel() {
      ColorSpinor<Float,Nc,Ns> a;
#pragma unroll
      for (int c=0; c<Nc; c++) {
  a(0,c) =  (*this)(1,c)+(*this)(3,c);
  a(1,c) = -(*this)(2,c)-(*this)(0,c);
  a(2,c) = -(*this)(3,c)+(*this)(1,c);
  a(3,c) = -(*this)(0,c)+(*this)(2,c);
      }
      *this = a;
    }

    __device__ __host__ inline void toRel() {
      ColorSpinor<Float,Nc,Ns> a;
#pragma unroll
      for (int c=0; c<Nc; c++) {
  a(0,c) = -(*this)(1,c)-(*this)(3,c);
  a(1,c) =  (*this)(2,c)+(*this)(0,c);
  a(2,c) =  (*this)(3,c)-(*this)(1,c);
  a(3,c) =  (*this)(0,c)-(*this)(2,c);
      }
      *this = a;
    }

  __device__ __host__ void print() {
      for (int s=0; s<Ns; s++) {
  for (int c=0; c<Nc; c++) {
    printf("s=%d c=%d %e %e\n", s, c, data[s*Nc+c].real(), data[s*Nc+c].imag());
  }
      }
    }
  };

  template <typename Float, int Nc>
    struct ColorSpinor<Float, Nc, 2> {
    static const int Ns = 2;
    complex<Float> data[Nc*2];
    
    __device__ __host__ inline ColorSpinor<Float, Nc, 2>() {
#pragma unroll
      for (int i=0; i<Nc*Ns; i++) {
  data[i] = 0;
      }      
    }

    __device__ __host__ inline ColorSpinor<Float, Nc, 2>(const ColorSpinor<Float, Nc, 2> &a) {
#pragma unroll
      for (int i=0; i<Nc*Ns; i++) {
  data[i] = a.data[i];
      }      
    }


    __device__ __host__ inline ColorSpinor<Float, Nc, 2>& operator=(const ColorSpinor<Float, Nc, 2> &a) {
      if (this != &a) {
#pragma unroll
  for (int i=0; i<Nc*Ns; i++) {
    data[i] = a.data[i];
  }
      }
      return *this;
    }

    __device__ __host__ inline ColorSpinor<Float, Nc, 2>& operator+=(const ColorSpinor<Float, Nc, 2> &a) {
      if (this != &a) {
#pragma unroll
  for (int i=0; i<Nc*Ns; i++) {
    data[i] += a.data[i];
  }
      }
      return *this;
    }

    __device__ __host__ inline ColorSpinor<Float,Nc,4> chiral_reconstruct(int chirality) const {
      ColorSpinor<Float,Nc,4> recon;
#pragma unroll
      for (int s=0; s<Ns; s++) {
#pragma unroll
  for (int c=0; c<Nc; c++) {
    recon(chirality*Ns+s,c) = (*this)(s,c);
  }
      }
      return recon;
    }

      __device__ __host__ inline ColorSpinor<Float,Nc,4> reconstruct(int dim, int sign) {
      ColorSpinor<Float,Nc,4> recon;
      complex<Float> j(0.0,1.0);
      
      switch (dim) {
      case 0: // x dimension
  switch (sign) {
  case 1: // positive projector
#pragma unroll
    for (int i=0; i<Nc; i++) {
      recon(0,i) = (*this)(0,i);
      recon(1,i) = (*this)(1,i);
      recon(2,i) = -j*(*this)(1,i);
      recon(3,i) = -j*(*this)(0,i);
    }
    break;
  case -1: // negative projector
#pragma unroll
    for (int i=0; i<Nc; i++) {
      recon(0,i) = (*this)(0,i);
      recon(1,i) = (*this)(1,i);
      recon(2,i) = j*(*this)(1,i);
      recon(3,i) = j*(*this)(0,i);
    }
    break;
  }
  break;
      case 1: // y dimension
  switch (sign) {
  case 1: // positive projector
#pragma unroll
    for (int i=0; i<Nc; i++) {
      recon(0,i) = (*this)(0,i);
      recon(1,i) = (*this)(1,i);
      recon(2,i) = -(*this)(1,i);
      recon(3,i) = (*this)(0,i);
    }
    break;
  case -1: // negative projector
#pragma unroll
    for (int i=0; i<Nc; i++) {
      recon(0,i) = (*this)(0,i);
      recon(1,i) = (*this)(1,i);
      recon(2,i) = (*this)(1,i);
      recon(3,i) = -(*this)(0,i);
    }
    break;
  }
  break;      
      case 2: // z dimension
  switch (sign) {
  case 1: // positive projector
#pragma unroll
    for (int i=0; i<Nc; i++) {
      recon(0,i) = (*this)(0,i);
      recon(1,i) = (*this)(1,i);
      recon(2,i) = -j*(*this)(0,i);
      recon(3,i) = j*(*this)(1,i);
    }
    break;
  case -1: // negative projector
#pragma unroll
    for (int i=0; i<Nc; i++) {
      recon(0,i) = (*this)(0,i);
      recon(1,i) = (*this)(1,i);
      recon(2,i) = j*(*this)(0,i);
      recon(3,i) = -j*(*this)(1,i);
    }
    break;
  }
  break;
      case 3: // t dimension
  switch (sign) {
  case 1: // positive projector
#pragma unroll
    for (int i=0; i<Nc; i++) {
      recon(0,i) = (*this)(0,i);
      recon(1,i) = (*this)(1,i);
      recon(2,i) = 0;
      recon(3,i) = 0;
    }
    break;
  case -1: // negative projector
#pragma unroll
    for (int i=0; i<Nc; i++) {
      recon(0,i) = 0;
      recon(1,i) = 0;
      recon(2,i) = (*this)(0,i);
      recon(3,i) = (*this)(1,i);
    }
    break;
  }
  break;
      }
      return recon;
    }

    __device__ __host__ inline complex<Float>& operator()(int s, int c) { return data[s*Nc + c]; }

    __device__ __host__ inline const complex<Float>& operator()(int s, int c) const { return data[s*Nc + c]; }

    __device__ __host__ inline complex<Float>& operator()(int idx) { return data[idx]; }

    __device__ __host__ inline const complex<Float>& operator()(int idx) const { return data[idx]; }

    template<typename S>
      __device__ __host__ inline ColorSpinor<Float, Nc, Ns>(const colorspinor_wrapper<Float, S> &s);

    template<typename S>
      __device__ __host__ inline void operator=(const colorspinor_wrapper<Float, S> &s);

    template<typename S>
      __device__ __host__ inline ColorSpinor<Float, Nc, Ns>(const colorspinor_ghost_wrapper<Float, S> &s);

    template<typename S>
      __device__ __host__ inline void operator=(const colorspinor_ghost_wrapper<Float, S> &s);

    __device__ __host__ void print() {
      for (int s=0; s<Ns; s++) {
  for (int c=0; c<Nc; c++) {
    printf("s=%d c=%d %e %e\n", s, c, data[s*Nc+c].real(), data[s*Nc+c].imag());
  }
      }
    }
  };


  template<typename Float, int Nc, int Ns> __device__ __host__ inline
    Matrix<complex<Float>,Nc> outerProdSpinTrace(const ColorSpinor<Float,Nc,Ns> &a, const ColorSpinor<Float,Nc,Ns> &b) {

    Matrix<complex<Float>,Nc> out;

    // outer product over color
#pragma unroll
    for (int i=0; i<Nc; i++) {
#pragma unroll
      for(int j=0; j<Nc; j++) {
  // trace over spin (manual unroll for perf)
  out(j,i).real(                   a(0,j).real() * b(0,i).real() );
  out(j,i).real( out(j,i).real() + a(0,j).imag() * b(0,i).imag() );
  out(j,i).imag(                   a(0,j).imag() * b(0,i).real() );
  out(j,i).imag( out(j,i).imag() - a(0,j).real() * b(0,i).imag() );
  //out(j,i) = a(0,j) * conj(b(0,i));

#pragma unroll
  for (int s=1; s<Ns; s++) {
    out(j,i).real( out(j,i).real() + a(s,j).real() * b(s,i).real() );
    out(j,i).real( out(j,i).real() + a(s,j).imag() * b(s,i).imag() );
    out(j,i).imag( out(j,i).imag() + a(s,j).imag() * b(s,i).real() );
    out(j,i).imag( out(j,i).imag() - a(s,j).real() * b(s,i).imag() );
    // out(j,i) += a(s,j) * conj(b(s,i));
  }
      }
    }
    return out;
  }

  template<typename Float, int Nc, int Ns> __device__ __host__ inline
    ColorSpinor<Float,Nc,Ns> operator+(const ColorSpinor<Float,Nc,Ns> &x, const ColorSpinor<Float,Nc,Ns> &y) {

    ColorSpinor<Float,Nc,Ns> z;

#pragma unroll
    for (int i=0; i<Nc; i++) {
#pragma unroll
      for (int s=0; s<Ns; s++) {
  z.data[s*Nc + i] = x.data[s*Nc + i] + y.data[s*Nc + i];
      }
    }

    return z;
  }

  template<typename Float, int Nc, int Ns> __device__ __host__ inline
    ColorSpinor<Float,Nc,Ns> operator-(const ColorSpinor<Float,Nc,Ns> &x, const ColorSpinor<Float,Nc,Ns> &y) {

    ColorSpinor<Float,Nc,Ns> z;

#pragma unroll
    for (int i=0; i<Nc; i++) {
#pragma unroll
      for (int s=0; s<Ns; s++) {
  z.data[s*Nc + i] = x.data[s*Nc + i] - y.data[s*Nc + i];
      }
    }

    return z;
  }

  template<typename Float, int Nc, int Ns, typename S> __device__ __host__ inline
    ColorSpinor<Float,Nc,Ns> operator*(const S &a, const ColorSpinor<Float,Nc,Ns> &x) {

    ColorSpinor<Float,Nc,Ns> y;

#pragma unroll
    for (int i=0; i<Nc; i++) {
#pragma unroll
      for (int s=0; s<Ns; s++) {
  y.data[s*Nc + i] = a * x.data[s*Nc + i];
      }
    }

    return y;
  }

  template<typename Float, int Nc, int Ns> __device__ __host__ inline
    ColorSpinor<Float,Nc,Ns> operator*(const Matrix<complex<Float>,Nc> &A, const ColorSpinor<Float,Nc,Ns> &x) {

    ColorSpinor<Float,Nc,Ns> y;

#pragma unroll
    for (int i=0; i<Nc; i++) {
#pragma unroll
      for (int s=0; s<Ns; s++) {
  y.data[s*Nc + i].x  = A(i,0).real() * x.data[s*Nc + 0].real();
  y.data[s*Nc + i].x -= A(i,0).imag() * x.data[s*Nc + 0].imag();
  y.data[s*Nc + i].y  = A(i,0).real() * x.data[s*Nc + 0].imag();
  y.data[s*Nc + i].y += A(i,0).imag() * x.data[s*Nc + 0].real();
      }
#pragma unroll
      for (int j=1; j<Nc; j++) {
#pragma unroll
  for (int s=0; s<Ns; s++) {
    y.data[s*Nc + i].x += A(i,j).real() * x.data[s*Nc + j].real();
    y.data[s*Nc + i].x -= A(i,j).imag() * x.data[s*Nc + j].imag();
    y.data[s*Nc + i].y += A(i,j).real() * x.data[s*Nc + j].imag();
    y.data[s*Nc + i].y += A(i,j).imag() * x.data[s*Nc + j].real();
  }
      }
    }

    return y;
  }

  template<typename Float, int Nc, int Ns> __device__ __host__ inline
    ColorSpinor<Float,Nc,Ns> operator*(const HMatrix<Float,Nc*Ns> &A, const ColorSpinor<Float,Nc,Ns> &x) {

    ColorSpinor<Float,Nc,Ns> y;
    constexpr int N = Nc*Ns;

#pragma unroll
    for (int i=0; i<N; i++) {
      if (i==0) {
  y.data[i].x  = A(i,0).real() * x.data[0].real();
  y.data[i].y  = A(i,0).real() * x.data[0].imag();
      } else {
  y.data[i].x  = A(i,0).real() * x.data[0].real();
  y.data[i].x -= A(i,0).imag() * x.data[0].imag();
  y.data[i].y  = A(i,0).real() * x.data[0].imag();
  y.data[i].y += A(i,0).imag() * x.data[0].real();
      }
#pragma unroll
      for (int j=1; j<N; j++) {
  if (i==j) {
    y.data[i].x += A(i,j).real() * x.data[j].real();
    y.data[i].y += A(i,j).real() * x.data[j].imag();
  } else {
    y.data[i].x += A(i,j).real() * x.data[j].real();
    y.data[i].x -= A(i,j).imag() * x.data[j].imag();
    y.data[i].y += A(i,j).real() * x.data[j].imag();
    y.data[i].y += A(i,j).imag() * x.data[j].real();
  }
      }
    }

    return y;
  }

} // namespace quda