multigrid.h

Go to the documentation of this file.

#pragma once
 
#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;
 
    int NblockOrtho;
 
    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;
 
    QudaFieldLocation setup_location;
 
    char filename[100];
 
    QudaTransferType transfer_type;
 
    bool use_mma;
 
    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]),
      transfer_type(param.transfer_type[level]),
      use_mma(param.use_mma == QUDA_BOOLEAN_TRUE)
    {
      // 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]),
      NblockOrtho(param.mg_global.n_block_ortho[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]),
      setup_location(param.mg_global.setup_location[level]),
      transfer_type(param.mg_global.transfer_type[level]),
      use_mma(param.use_mma)
    {
      // 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;
 
    cudaGaugeField *xInvKD;
 
    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();
 
    bool hermitian() { return false; };
 
    void reset(bool refresh=false);
 
    void dumpNullVectors() const;
 
    void createSmoother();
 
    void destroySmoother();
 
    void createCoarseDirac();
 
    void createCoarseSolver();
 
    void destroyCoarseSolver();
 
    void verify(bool recursively = false);
 
    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 mass, double mu, double mu_factor,
                QudaDiracType dirac, QudaMatPCType matpc);
 
  void StaggeredCoarseOp(GaugeField &Y, GaugeField &X, const Transfer &T, const cudaGaugeField &gauge,
                         const cudaGaugeField *XinvKD, double mass, 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 mass, double mu, double mu_factor,
                      QudaDiracType dirac, QudaMatPCType matpc, bool need_bidirectional, bool use_mma = false);
 
  void calculateYhat(GaugeField &Yhat, GaugeField &Xinv, const GaugeField &Y, const GaugeField &X, bool use_mma = false);
 
  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