multigrid.h

Go to the documentation of this file.

#ifndef _MG_QUDA_H
#define _MG_QUDA_H

#include <invert_quda.h>
#include <transfer.h>
#include <vector>
#include <complex_quda.h>

// at the moment double-precision multigrid is only enabled when debugging
#ifdef HOST_DEBUG
#define GPU_MULTIGRID_DOUBLE
#endif

namespace quda {

  // forward declarations
  class MG;
  class DiracCoarse;

  struct MGParam : SolverParam {

    QudaMultigridParam  &mg_global;

    int level;

    int Nlevel;

    int geoBlockSize[QUDA_MAX_DIM];

    int spinBlockSize;

    int Nvec;

    MG *coarse;

    MG *fine;

    std::vector<ColorSpinorField*> &B;

    int nu_pre;

    int nu_post;

    double smoother_tol;

    QudaMultigridCycleType cycle_type;

    QudaBoolean global_reduction;

    DiracMatrix *matResidual;

    DiracMatrix *matSmooth;

    DiracMatrix *matSmoothSloppy;

    QudaInverterType smoother;

    QudaSolutionType coarse_grid_solution_type;

    QudaSolveType smoother_solve_type;

    QudaFieldLocation location;

    char filename[100];

    MGParam(QudaMultigridParam &param,
      std::vector<ColorSpinorField*> &B,
      DiracMatrix *matResidual,
      DiracMatrix *matSmooth,
      DiracMatrix *matSmoothSloppy,
      int level=0) :
      SolverParam(*(param.invert_param)), 
      mg_global(param), 
      level(level),
      Nlevel(param.n_level),
      spinBlockSize(param.spin_block_size[level]),
      Nvec(param.n_vec[level]),
      B(B), 
      nu_pre(param.nu_pre[level]),
      nu_post(param.nu_post[level]),
      smoother_tol(param.smoother_tol[level]),
      cycle_type(param.cycle_type[level]),
      global_reduction(param.global_reduction[level]),
      matResidual(matResidual),
      matSmooth(matSmooth),
      matSmoothSloppy(matSmoothSloppy),
      smoother(param.smoother[level]),
      coarse_grid_solution_type(param.coarse_grid_solution_type[level]),
      smoother_solve_type(param.smoother_solve_type[level]),
      location(param.location[level])
      { 
  // set the block size
  for (int i=0; i<QUDA_MAX_DIM; i++) geoBlockSize[i] = param.geo_block_size[level][i];

  // set the smoother relaxation factor
  omega = param.omega[level];
      }

    MGParam(const MGParam &param, 
      std::vector<ColorSpinorField*> &B,
      DiracMatrix *matResidual,
      DiracMatrix *matSmooth,
      DiracMatrix *matSmoothSloppy,
      int level=0) :
      SolverParam(param),
      mg_global(param.mg_global),
      level(level),
      Nlevel(param.Nlevel),
      spinBlockSize(param.mg_global.spin_block_size[level]),
      Nvec(param.mg_global.n_vec[level]),
      coarse(param.coarse),
      fine(param.fine),
      B(B),
      nu_pre(param.mg_global.nu_pre[level]),
      nu_post(param.mg_global.nu_post[level]),
      smoother_tol(param.mg_global.smoother_tol[level]),
      cycle_type(param.mg_global.cycle_type[level]),
      global_reduction(param.mg_global.global_reduction[level]),
      matResidual(matResidual),
      matSmooth(matSmooth),
      matSmoothSloppy(matSmoothSloppy),
      smoother(param.mg_global.smoother[level]),
      coarse_grid_solution_type(param.mg_global.coarse_grid_solution_type[level]),
      smoother_solve_type(param.mg_global.smoother_solve_type[level]),
      location(param.mg_global.location[level])
      {
  // set the block size
  for (int i=0; i<QUDA_MAX_DIM; i++) geoBlockSize[i] = param.mg_global.geo_block_size[level][i];

  // set the smoother relaxation factor
  omega = param.mg_global.omega[level];
      }

  };

  class MG : public Solver {

  private:
    MGParam &param;

    Transfer *transfer;

    Solver *presmoother, *postsmoother;

    TimeProfile &profile_global;

    TimeProfile profile;

    char prefix[128];

    char coarse_prefix[128];

    MG *coarse;

    MG *fine;

    Solver *coarse_solver;

    MGParam *param_coarse;

    SolverParam *param_presmooth;

    SolverParam *param_postsmooth;

    SolverParam *param_coarse_solver;

    std::vector<ColorSpinorField*> *B;

    std::vector<ColorSpinorField*> *B_coarse;

    ColorSpinorField *r;

    ColorSpinorField *b_tilde;

    ColorSpinorField *r_coarse;

    ColorSpinorField *x_coarse;

    ColorSpinorField *tmp_coarse;

    Dirac *diracCoarseResidual;

    Dirac *diracCoarseSmoother;

    Dirac *diracCoarseSmootherSloppy;

    DiracMatrix *matCoarseResidual;

    DiracMatrix *matCoarseSmoother;

    DiracMatrix *matCoarseSmootherSloppy;

  public:
    MG(MGParam &param, TimeProfile &profile);

    virtual ~MG();

    void createSmoother();

    void destroySmoother();

    void reset();

    void verify();

    void operator()(ColorSpinorField &out, ColorSpinorField &in);

    void loadVectors(std::vector<ColorSpinorField*> &B);

    void saveVectors(std::vector<ColorSpinorField*> &B);

    void generateNullVectors(std::vector<ColorSpinorField*> B);

    double flops() const;

  };

  void ApplyCoarse(ColorSpinorField &out, const ColorSpinorField &inA, const ColorSpinorField &inB,
       const GaugeField &Y, const GaugeField &X, double kappa, int parity = QUDA_INVALID_PARITY,
       bool dslash=true, bool clover=true, bool dagger=false);

  void CoarseOp(GaugeField &Y, GaugeField &X, GaugeField &Xinv, GaugeField &Yhat, const Transfer &T,
    const cudaGaugeField &gauge, const cudaCloverField *clover,
    double kappa, double mu, double mu_factor, QudaDiracType dirac, QudaMatPCType matpc);

  void CoarseCoarseOp(GaugeField &Y, GaugeField &X, GaugeField &Xinv, GaugeField &Yhat, const Transfer &T,
          const GaugeField &gauge, const GaugeField &clover, const GaugeField &cloverInv,
          double kappa, double mu, double mu_factor, QudaDiracType dirac, QudaMatPCType matpc);

  struct multigrid_solver {
    Dirac *d;
    Dirac *dSmooth;
    Dirac *dSmoothSloppy;

    DiracM *m;
    DiracM *mSmooth;
    DiracM *mSmoothSloppy;

    std::vector<ColorSpinorField*> B;

    MGParam *mgParam;

    MG *mg;
    TimeProfile &profile;

    multigrid_solver(QudaMultigridParam &mg_param, TimeProfile &profile);

    virtual ~multigrid_solver()
    {
      profile.TPSTART(QUDA_PROFILE_FREE);
      if (mg) delete mg;

      if (mgParam) delete mgParam;

      for (unsigned int i=0; i<B.size(); i++) delete B[i];

      if (m) delete m;
      if (mSmooth) delete mSmooth;
      if (mSmoothSloppy) delete mSmoothSloppy;

      if (d) delete d;
      if (dSmooth) delete dSmooth;
      if (dSmoothSloppy && dSmoothSloppy != dSmooth) delete dSmoothSloppy;
      profile.TPSTOP(QUDA_PROFILE_FREE);
    }
  };

} // namespace quda

#endif // _MG_QUDA_H