multigrid.h

Go to the documentation of this file.

#pragma once

#include <invert_quda.h>
//#include <eigensolve_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;

    int NblockOrtho;

    MG *coarse;

    MG *fine;

    std::vector<ColorSpinorField*> &B;

    std::vector<Complex> evals;

    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;

    QudaFieldLocation setup_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]),
      NblockOrtho(param.n_block_ortho[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]),
      setup_location(param.setup_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, std::vector<Complex> evals,
            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]),
      NblockOrtho(param.mg_global.n_block_ortho[level]),
      coarse(param.coarse),
      fine(param.fine),
      B(B),
      evals(evals),
      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]),
      setup_location(param.mg_global.setup_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;

    bool resetTransfer;

    Solver *presmoother, *postsmoother;

    TimeProfile &profile_global;

    TimeProfile profile;

    char prefix[128];

    char coarse_prefix[128];

    MG *coarse;

    Solver *coarse_solver;

    MGParam *param_coarse;

    SolverParam *param_presmooth;

    SolverParam *param_postsmooth;

    SolverParam *param_coarse_solver;

    std::vector<ColorSpinorField*> *B_coarse;

    ColorSpinorField *r;

    ColorSpinorField *b_tilde;

    ColorSpinorField *r_coarse;

    ColorSpinorField *x_coarse;

    ColorSpinorField *tmp_coarse;

    ColorSpinorField *tmp2_coarse;

    const Dirac *diracResidual;

    const Dirac *diracSmoother;

    const Dirac *diracSmootherSloppy;

    Dirac *diracCoarseResidual;

    Dirac *diracCoarseSmoother;

    Dirac *diracCoarseSmootherSloppy;

    DiracMatrix *matCoarseResidual;

    DiracMatrix *matCoarseSmoother;

    DiracMatrix *matCoarseSmootherSloppy;

    RNG *rng;

    void pushLevel(int level) const;

    void popLevel(int level) const;

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

    virtual ~MG();

    void reset(bool refresh=false);

    void dumpNullVectors() const;

    void createSmoother();

    void destroySmoother();

    void createCoarseDirac();

    void createCoarseSolver();

    void destroyCoarseSolver();

    void verify();

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

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

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

    void generateNullVectors(std::vector<ColorSpinorField*> &B, bool refresh=false);

    void generateEigenVectors();

    void buildFreeVectors(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, const int *commDim=0,
                   QudaPrecision halo_precision=QUDA_INVALID_PRECISION);

  void CoarseOp(GaugeField &Y, GaugeField &X, 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, const Transfer &T, const GaugeField &gauge,
                      const GaugeField &clover, const GaugeField &cloverInv, double kappa, double mu, double mu_factor,
                      QudaDiracType dirac, QudaMatPCType matpc, bool need_bidirectional);

  void calculateYhat(GaugeField &Yhat, GaugeField &Xinv, const GaugeField &Y, const GaugeField &X);

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

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

    std::vector<ColorSpinorField*> B;
    std::vector<Complex> evals;

    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