clover_deriv_quda.cu

Go to the documentation of this file.

#include <cstdio>
#include <cstdlib>
#include <cuda.h>
#include <tune_quda.h>
#include <gauge_field.h>
#include <cassert>

#include <jitify_helper.cuh>
#include <kernels/clover_deriv.cuh>

namespace quda {

#ifdef GPU_CLOVER_DIRAC

  template<typename Float, typename Arg>
  class CloverDerivative : public TunableVectorY {

  private:
    Arg arg;
    const GaugeField &meta;

#if defined(SHARED_ACCUMULATOR) && defined(SHARED_ARRAY)
    unsigned int sharedBytesPerThread() const { return 18*sizeof(Float) + 8; }
#elif defined(SHARED_ACCUMULATOR)
    unsigned int sharedBytesPerThread() const { return 18*sizeof(Float); }
#else
    unsigned int sharedBytesPerThread() const { return 0; }
#endif
    unsigned int sharedBytesPerBlock(const TuneParam &) const { return 0; }

    unsigned int minThreads() const { return arg.volumeCB; }
    bool tuneGridDim() const { return false; } // don't tune the grid dimension

  public:
    CloverDerivative(const Arg &arg, const GaugeField &meta) : TunableVectorY(2), arg(arg), meta(meta) {
      writeAuxString("threads=%d,prec=%lu,fstride=%d,gstride=%d,ostride=%d",
         arg.volumeCB,sizeof(Float),arg.force.stride,
         arg.gauge.stride,arg.oprod.stride);
#ifdef JITIFY
      create_jitify_program("kernels/clover_deriv.cuh");
#endif
    }
    virtual ~CloverDerivative() {}

    void apply(const cudaStream_t &stream){
      TuneParam tp = tuneLaunch(*this, getTuning(), getVerbosity());
#ifdef JITIFY
      using namespace jitify::reflection;
      jitify_error = program->kernel("quda::cloverDerivativeKernel")
                         .instantiate(Type<Float>(), Type<Arg>())
                         .configure(tp.grid, tp.block, tp.shared_bytes, stream)
                         .launch(arg);
#else
      cloverDerivativeKernel<Float><<<tp.grid,tp.block,tp.shared_bytes>>>(arg);
#endif
    } // apply

    bool advanceBlockDim(TuneParam &param) const {
      dim3 block = param.block;
      dim3 grid = param.grid;
      bool rtn = TunableVectorY::advanceBlockDim(param);
      param.block.z = block.z;
      param.grid.z = grid.z;

      if (!rtn) {
  if (param.block.z < 4) {
    param.block.z++;
    param.grid.z = (4 + param.block.z - 1) / param.block.z;
    rtn = true;
  } else {
    param.block.z = 1;
    param.grid.z = 4;
    rtn = false;
  }
      }
      return rtn;
    }

    void initTuneParam(TuneParam &param) const {
      TunableVectorY::initTuneParam(param);
      param.block.y = 1;
      param.block.z = 1;
      param.grid.y = 2;
      param.grid.z = 4;
    }

    void defaultTuneParam(TuneParam &param) const { initTuneParam(param); }

    // The force field is updated so we must preserve its initial state
    void preTune() { arg.force.save(); }
    void postTune() { arg.force.load(); }

    long long flops() const { return 16 * 198 * 3 * 4 * 2 * (long long)arg.volumeCB; }
    long long bytes() const { return ((8*arg.gauge.Bytes() + 4*arg.oprod.Bytes())*3 + 2*arg.force.Bytes()) * 4 * 2 * arg.volumeCB; }

    TuneKey tuneKey() const { return TuneKey(meta.VolString(), typeid(*this).name(), aux); }
  };

  template<typename Float>
  void cloverDerivative(cudaGaugeField &force,
      cudaGaugeField &gauge,
      cudaGaugeField &oprod,
      double coeff, int parity) {

    if (oprod.Reconstruct() != QUDA_RECONSTRUCT_NO) errorQuda("Force field does not support reconstruction");

    if (force.Order() != oprod.Order()) errorQuda("Force and Oprod orders must match");

    if (force.Reconstruct() != QUDA_RECONSTRUCT_NO) errorQuda("Force field does not support reconstruction");

    if (force.Order() == QUDA_FLOAT2_GAUGE_ORDER) {
      typedef gauge::FloatNOrder<Float, 18, 2, 18> F;
      typedef gauge::FloatNOrder<Float, 18, 2, 18> O;

      if (gauge.isNative()) {
  if (gauge.Reconstruct() == QUDA_RECONSTRUCT_NO) {
    typedef typename gauge_mapper<Float,QUDA_RECONSTRUCT_NO>::type G;
    typedef CloverDerivArg<Float,F,G,O> Arg;
    Arg arg(F(force), G(gauge), O(oprod), force.X(), oprod.X(), coeff, parity);
    CloverDerivative<Float, Arg> deriv(arg, gauge);
    deriv.apply(0);
#if 0
  } else if (gauge.Reconstruct() == QUDA_RECONSTRUCT_12) {
    typedef typename gauge_mapper<Float,QUDA_RECONSTRUCT_12>::type G;
    typedef CloverDerivArg<Float,F,G,O> Arg;
    Arg arg(F(force), G(gauge), O(oprod), force.X(), oprod.X(), coeff, parity);
    CloverDerivative<Float, Arg> deriv(arg, gauge);
    deriv.apply(0);
#endif
  } else {
    errorQuda("Reconstruction type %d not supported",gauge.Reconstruct());
  }
      } else {
  errorQuda("Gauge order %d not supported", gauge.Order());
      }
    } else {
      errorQuda("Force order %d not supported", force.Order());
    } // force / oprod order

    qudaDeviceSynchronize();
  }
#endif // GPU_CLOVER

  void cloverDerivative(
      cudaGaugeField &force, cudaGaugeField &gauge, cudaGaugeField &oprod, double coeff, QudaParity parity)
  {
#ifdef GPU_CLOVER_DIRAC
  assert(oprod.Geometry() == QUDA_TENSOR_GEOMETRY);
  assert(force.Geometry() == QUDA_VECTOR_GEOMETRY);

  for (int d=0; d<4; d++) {
    if (oprod.X()[d] != gauge.X()[d])
      errorQuda("Incompatible extended dimensions d=%d gauge=%d oprod=%d", d, gauge.X()[d], oprod.X()[d]);
  }

  int device_parity = (parity == QUDA_EVEN_PARITY) ? 0 : 1;

  if(force.Precision() == QUDA_DOUBLE_PRECISION){
    cloverDerivative<double>(force, gauge, oprod, coeff, device_parity);
#if 0
  } else if (force.Precision() == QUDA_SINGLE_PRECISION){
    cloverDerivative<float>(force, gauge, oprod, coeff, device_parity);
#endif
  } else {
    errorQuda("Precision %d not supported", force.Precision());
  }

  return;
#else
  errorQuda("Clover has not been built");
#endif
  }

} // namespace quda