invert_quda.h

Go to the documentation of this file.

#ifndef _INVERT_QUDA_H
#define _INVERT_QUDA_H

#include <quda.h>
#include <quda_internal.h>
#include <dirac_quda.h>
#include <color_spinor_field.h>
#include <vector>

namespace quda {

  struct SolverParam {
    QudaInverterType inv_type;

    QudaInverterType inv_type_precondition;


    void *preconditioner;

    void *deflation_op;

    QudaResidualType residual_type;

    QudaUseInitGuess use_init_guess;

    QudaComputeNullVector compute_null_vector;

    double delta;

    bool use_sloppy_partial_accumulator;

    int solution_accumulator_pipeline;

    int max_res_increase;

    int max_res_increase_total;

    int heavy_quark_check;

    int pipeline;

    double tol;

    double tol_restart;

    double tol_hq;

    bool compute_true_res;

    double true_res;

    double true_res_hq;

    int maxiter;

    int iter;

    QudaPrecision precision;

    QudaPrecision precision_sloppy;

    QudaPrecision precision_precondition;

    QudaPreserveSource preserve_source;

    int overlap_precondition;


    int num_src;

    // Multi-shift solver parameters

    int num_offset;

    double offset[QUDA_MAX_MULTI_SHIFT];

    double tol_offset[QUDA_MAX_MULTI_SHIFT];

    double tol_hq_offset[QUDA_MAX_MULTI_SHIFT];

    double true_res_offset[QUDA_MAX_MULTI_SHIFT];

    double iter_res_offset[QUDA_MAX_MULTI_SHIFT];

    double true_res_hq_offset[QUDA_MAX_MULTI_SHIFT];


    int Nsteps;

    int Nkrylov;

    int precondition_cycle;

    double tol_precondition;

    int maxiter_precondition;

    double omega;



    QudaSchwarzType schwarz_type;

    double secs;

    double gflops;

    // Incremental EigCG solver parameters
    QudaPrecision precision_ritz;//also search space precision

    int nev;//number of eigenvectors produced by EigCG
    int m;//Dimension of the search space
    int deflation_grid;
    int rhs_idx;

    int     eigcg_max_restarts;
    int     max_restart_num;
    double  inc_tol;
    double  eigenval_tol;

    QudaVerbosity verbosity_precondition; 

    bool is_preconditioner; 

    bool global_reduction; 

    QudaExtLibType extlib_type;

    SolverParam() : compute_null_vector(QUDA_COMPUTE_NULL_VECTOR_NO),
      compute_true_res(true), verbosity_precondition(QUDA_SILENT) { ; }

    SolverParam(const QudaInvertParam &param) : inv_type(param.inv_type),
      inv_type_precondition(param.inv_type_precondition), preconditioner(param.preconditioner), deflation_op(param.deflation_op),
      residual_type(param.residual_type), use_init_guess(param.use_init_guess),
      compute_null_vector(QUDA_COMPUTE_NULL_VECTOR_NO), delta(param.reliable_delta),
      use_sloppy_partial_accumulator(param.use_sloppy_partial_accumulator),
      solution_accumulator_pipeline(param.solution_accumulator_pipeline),
      max_res_increase(param.max_res_increase), max_res_increase_total(param.max_res_increase_total),
      heavy_quark_check(param.heavy_quark_check), pipeline(param.pipeline),
      tol(param.tol), tol_restart(param.tol_restart), tol_hq(param.tol_hq),
      compute_true_res(param.compute_true_res), true_res(param.true_res),
      true_res_hq(param.true_res_hq), maxiter(param.maxiter), iter(param.iter),
      precision(param.cuda_prec), precision_sloppy(param.cuda_prec_sloppy),
      precision_precondition(param.cuda_prec_precondition),
      preserve_source(param.preserve_source), num_src(param.num_src), num_offset(param.num_offset),
      Nsteps(param.Nsteps), Nkrylov(param.gcrNkrylov), precondition_cycle(param.precondition_cycle),
      tol_precondition(param.tol_precondition), maxiter_precondition(param.maxiter_precondition),
      omega(param.omega), schwarz_type(param.schwarz_type), secs(param.secs), gflops(param.gflops),
      precision_ritz(param.cuda_prec_ritz), nev(param.nev), m(param.max_search_dim),
      deflation_grid(param.deflation_grid), rhs_idx(0),
      eigcg_max_restarts(param.eigcg_max_restarts), max_restart_num(param.max_restart_num),
      inc_tol(param.inc_tol), eigenval_tol(param.eigenval_tol),
      verbosity_precondition(param.verbosity_precondition),
      is_preconditioner(false), global_reduction(true), extlib_type(param.extlib_type)
    {
      for (int i=0; i<num_offset; i++) {
  offset[i] = param.offset[i];
  tol_offset[i] = param.tol_offset[i];
  tol_hq_offset[i] = param.tol_hq_offset[i];
      }

      if(param.rhs_idx != 0 && (param.inv_type==QUDA_INC_EIGCG_INVERTER || param.inv_type==QUDA_GMRESDR_PROJ_INVERTER)){
        rhs_idx = param.rhs_idx;
      }
    }

    SolverParam(const SolverParam &param) : inv_type(param.inv_type),
      inv_type_precondition(param.inv_type_precondition), preconditioner(param.preconditioner), deflation_op(param.deflation_op),
      residual_type(param.residual_type), use_init_guess(param.use_init_guess),
      delta(param.delta), use_sloppy_partial_accumulator(param.use_sloppy_partial_accumulator),
      solution_accumulator_pipeline(param.solution_accumulator_pipeline),
      max_res_increase(param.max_res_increase), max_res_increase_total(param.max_res_increase_total),
      heavy_quark_check(param.heavy_quark_check), pipeline(param.pipeline),
      tol(param.tol), tol_restart(param.tol_restart), tol_hq(param.tol_hq),
      compute_true_res(param.compute_true_res), true_res(param.true_res),
      true_res_hq(param.true_res_hq), maxiter(param.maxiter), iter(param.iter),
      precision(param.precision), precision_sloppy(param.precision_sloppy),
      precision_precondition(param.precision_precondition),
      preserve_source(param.preserve_source), num_offset(param.num_offset),
      Nsteps(param.Nsteps), Nkrylov(param.Nkrylov), precondition_cycle(param.precondition_cycle),
      tol_precondition(param.tol_precondition), maxiter_precondition(param.maxiter_precondition),
      omega(param.omega), schwarz_type(param.schwarz_type), secs(param.secs), gflops(param.gflops),
      precision_ritz(param.precision_ritz), nev(param.nev), m(param.m),
      deflation_grid(param.deflation_grid), rhs_idx(0),
      eigcg_max_restarts(param.eigcg_max_restarts), max_restart_num(param.max_restart_num),
      inc_tol(param.inc_tol), eigenval_tol(param.eigenval_tol),
      verbosity_precondition(param.verbosity_precondition),
      is_preconditioner(param.is_preconditioner), global_reduction(param.global_reduction), extlib_type(param.extlib_type)
    {
      for (int i=0; i<num_offset; i++) {
  offset[i] = param.offset[i];
  tol_offset[i] = param.tol_offset[i];
  tol_hq_offset[i] = param.tol_hq_offset[i];
      }

      if((param.inv_type == QUDA_INC_EIGCG_INVERTER || param.inv_type == QUDA_EIGCG_INVERTER) && m % 16){//current hack for the magma library
        m = (m / 16) * 16 + 16;
        warningQuda("\nSwitched eigenvector search dimension to %d\n", m);
      }
      if(param.rhs_idx != 0 && (param.inv_type==QUDA_INC_EIGCG_INVERTER || param.inv_type==QUDA_GMRESDR_PROJ_INVERTER)){
        rhs_idx = param.rhs_idx;
      }
    }

    ~SolverParam() { }

    void updateInvertParam(QudaInvertParam &param, int offset=-1) {
      param.true_res = true_res;
      param.true_res_hq = true_res_hq;
      param.iter += iter;
      reduceDouble(gflops);
      param.gflops += gflops;
      param.secs += secs;
      if (offset >= 0) {
  param.true_res_offset[offset] = true_res_offset[offset];
        param.iter_res_offset[offset] = iter_res_offset[offset];
  param.true_res_hq_offset[offset] = true_res_hq_offset[offset];
      } else {
  for (int i=0; i<num_offset; i++) {
    param.true_res_offset[i] = true_res_offset[i];
          param.iter_res_offset[i] = iter_res_offset[i];
    param.true_res_hq_offset[i] = true_res_hq_offset[i];
  }
      }
      //for incremental eigCG:
      param.rhs_idx = rhs_idx;
    }

    void updateRhsIndex(QudaInvertParam &param) {
      //for incremental eigCG:
      rhs_idx = param.rhs_idx;
    }

  };

  class Solver {

  protected:
    SolverParam &param;
    TimeProfile &profile;

  public:
    Solver(SolverParam &param, TimeProfile &profile) : param(param), profile(profile) { ; }
    virtual ~Solver() { ; }

    virtual void operator()(ColorSpinorField &out, ColorSpinorField &in) = 0;

    virtual void solve(ColorSpinorField &out, ColorSpinorField &in);


    static Solver* create(SolverParam &param, DiracMatrix &mat, DiracMatrix &matSloppy,
        DiracMatrix &matPrecon, TimeProfile &profile);

    static double stopping(const double &tol, const double &b2, QudaResidualType residual_type);

    bool convergence(const double &r2, const double &hq2, const double &r2_tol,
         const double &hq_tol);

    bool convergenceHQ(const double &r2, const double &hq2, const double &r2_tol,
         const double &hq_tol);

    bool convergenceL2(const double &r2, const double &hq2, const double &r2_tol,
         const double &hq_tol);

    void PrintStats(const char*, int k, const double &r2, const double &b2, const double &hq2);

    void PrintSummary(const char *name, int k, const double &r2, const double &b2);

    virtual double flops() const { return 0; }
  };

  class CG : public Solver {

  private:
    const DiracMatrix &mat;
    const DiracMatrix &matSloppy;
    // pointers to fields to avoid multiple creation overhead
    ColorSpinorField *yp, *rp, *App, *tmpp;
    std::vector<ColorSpinorField*> p;
    bool init;

  public:
    CG(DiracMatrix &mat, DiracMatrix &matSloppy, SolverParam &param, TimeProfile &profile);
    virtual ~CG();

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

  class CGNE : public CG {

  private:
    DiracMMdag mmdag;
    DiracMMdag mmdagSloppy;
    ColorSpinorField *xp;
    bool init;

  public:
    CGNE(DiracMatrix &mat, DiracMatrix &matSloppy, SolverParam &param, TimeProfile &profile);
    virtual ~CGNE();

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

  class CGNR : public CG {

  private:
    DiracMdagM mdagm;
    DiracMdagM mdagmSloppy;
    ColorSpinorField *bp;
    bool init;

  public:
    CGNR(DiracMatrix &mat, DiracMatrix &matSloppy, SolverParam &param, TimeProfile &profile);
    virtual ~CGNR();

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



  class MPCG : public Solver {
    private:
      const DiracMatrix &mat;
      void computeMatrixPowers(cudaColorSpinorField out[], cudaColorSpinorField &in, int nvec);
      void computeMatrixPowers(std::vector<cudaColorSpinorField>& out, std::vector<cudaColorSpinorField>& in, int nsteps);


    public:
      MPCG(DiracMatrix &mat, SolverParam &param, TimeProfile &profile);
      virtual ~MPCG();

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



  class PreconCG : public Solver {
    private:
      const DiracMatrix &mat;
      const DiracMatrix &matSloppy;
      const DiracMatrix &matPrecon;

      Solver *K;
      SolverParam Kparam; // parameters for preconditioner solve

    public:
      PreconCG(DiracMatrix &mat, DiracMatrix &matSloppy, DiracMatrix &matPrecon, SolverParam &param, TimeProfile &profile);

      virtual ~PreconCG();

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


  class BiCGstab : public Solver {

  private:
    DiracMatrix &mat;
    const DiracMatrix &matSloppy;
    const DiracMatrix &matPrecon;

    // pointers to fields to avoid multiple creation overhead
    ColorSpinorField *yp, *rp, *pp, *vp, *tmpp, *tp;
    bool init;

  public:
    BiCGstab(DiracMatrix &mat, DiracMatrix &matSloppy, DiracMatrix &matPrecon,
       SolverParam &param, TimeProfile &profile);
    virtual ~BiCGstab();

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

  class SimpleBiCGstab : public Solver {

  private:
    DiracMatrix &mat;

    // pointers to fields to avoid multiple creation overhead
    cudaColorSpinorField *yp, *rp, *pp, *vp, *tmpp, *tp;
    bool init;

  public:
    SimpleBiCGstab(DiracMatrix &mat, SolverParam &param, TimeProfile &profile);
    virtual ~SimpleBiCGstab();

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

  class MPBiCGstab : public Solver {

  private:
    DiracMatrix &mat;
    void computeMatrixPowers(std::vector<cudaColorSpinorField>& pr, cudaColorSpinorField& p, cudaColorSpinorField& r, int nsteps);

  public:
    MPBiCGstab(DiracMatrix &mat, SolverParam &param, TimeProfile &profile);
    virtual ~MPBiCGstab();

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

  class BiCGstabL : public Solver {

  private:
    DiracMatrix &mat;
    const DiracMatrix &matSloppy;

    int nKrylov; // in the language of BiCGstabL, this is L.
    
    // Various coefficients and params needed on each iteration.
    Complex rho0, rho1, alpha, omega, beta; // Various coefficients for the BiCG part of BiCGstab-L. 
    Complex *gamma, *gamma_prime, *gamma_prime_prime; // Parameters for MR part of BiCGstab-L. (L+1) length.
    Complex **tau; // Parameters for MR part of BiCGstab-L. Tech. modified Gram-Schmidt coeffs. (L+1)x(L+1) length.
    double *sigma; // Parameters for MR part of BiCGstab-L. Tech. the normalization part of Gram-Scmidt. (L+1) length.
    
    // pointers to fields to avoid multiple creation overhead
    // full precision fields
    ColorSpinorField *r_fullp;   
    ColorSpinorField *yp;        
    // sloppy precision fields
    ColorSpinorField *tempp;     
    std::vector<ColorSpinorField*> r; // Current residual + intermediate residual values, along the MR.
    std::vector<ColorSpinorField*> u; // Search directions.
    
    // Saved, preallocated vectors. (may or may not get used depending on precision.)
    ColorSpinorField *x_sloppy_saved_p; 
    ColorSpinorField *r0_saved_p;       
    ColorSpinorField *r_sloppy_saved_p; 
    
    int reliable(double &rNorm, double &maxrx, double &maxrr, const double &r2, const double &delta);
    
    void computeTau(Complex **tau, double *sigma, std::vector<ColorSpinorField*> r, int begin, int size, int j);
    void updateR(Complex **tau, std::vector<ColorSpinorField*> r, int begin, int size, int j);
    void orthoDir(Complex **tau, double* sigma, std::vector<ColorSpinorField*> r, int j, int pipeline);
    
    void updateUend(Complex* gamma, std::vector<ColorSpinorField*> u, int nKrylov);
    void updateXRend(Complex* gamma, Complex* gamma_prime, Complex* gamma_prime_prime,
                                std::vector<ColorSpinorField*> r, ColorSpinorField& x, int nKrylov);
    
    bool init; 
    
    std::string solver_name; // holds BiCGstab-l, where 'l' literally equals nKrylov.

  public:
    BiCGstabL(DiracMatrix &mat, DiracMatrix &matSloppy, SolverParam &param, TimeProfile &profile);
    virtual ~BiCGstabL();

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

  class GCR : public Solver {

  private:
    const DiracMatrix &mat;
    const DiracMatrix &matSloppy;
    const DiracMatrix &matPrecon;

    Solver *K;
    SolverParam Kparam; // parameters for preconditioner solve

    int nKrylov;

    Complex *alpha;
    Complex **beta;
    double *gamma;

    bool init;

    ColorSpinorField *rp;       
    ColorSpinorField *yp;       
    ColorSpinorField *tmpp;     
    ColorSpinorField *x_sloppy; 
    ColorSpinorField *r_sloppy; 
    ColorSpinorField *r_pre;    
    ColorSpinorField *p_pre;    
    ColorSpinorField *rM;       

    std::vector<ColorSpinorField*> p;  // GCR direction vectors
    std::vector<ColorSpinorField*> Ap; // mat * direction vectors

  public:
    GCR(DiracMatrix &mat, DiracMatrix &matSloppy, DiracMatrix &matPrecon,
  SolverParam &param, TimeProfile &profile);

    GCR(DiracMatrix &mat, Solver &K, DiracMatrix &matSloppy, DiracMatrix &matPrecon,
  SolverParam &param, TimeProfile &profile);
    virtual ~GCR();

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

  class MR : public Solver {

  private:
    const DiracMatrix &mat;
    const DiracMatrix &matSloppy;
    ColorSpinorField *rp;
    ColorSpinorField *Arp;
    ColorSpinorField *tmpp;
    ColorSpinorField *yp;  //Holds initial guess if applicable
    bool init;
    bool allocate_r;
    bool allocate_y;

  public:
    MR(DiracMatrix &mat, DiracMatrix &matSloppy, SolverParam &param, TimeProfile &profile);
    virtual ~MR();

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

  // Steepest descent solver used as a preconditioner
  class SD : public Solver {
    private:
      const DiracMatrix &mat;
      cudaColorSpinorField *Ar;
      cudaColorSpinorField *r;
      cudaColorSpinorField *y;
      bool init;

    public:
      SD(DiracMatrix &mat, SolverParam &param, TimeProfile &profile);
      virtual ~SD();


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

  // Extended Steepest Descent solver used for overlapping DD preconditioning
  class XSD : public Solver {
    private:
      const DiracMatrix &mat;
      cudaColorSpinorField *xx;
      cudaColorSpinorField *bx;
      SD *sd; // extended sd is implemented using standard sd
      bool init;
      int R[4];

    public:
      XSD(DiracMatrix &mat, SolverParam &param, TimeProfile &profile);
      virtual ~XSD();

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


  class PreconditionedSolver : public Solver {

  private:
    Solver *solver;
    const Dirac &dirac;
    const char *prefix;

  public:
  PreconditionedSolver(Solver &solver, const Dirac &dirac, SolverParam &param, TimeProfile &profile, const char *prefix)
    : Solver(param, profile), solver(&solver), dirac(dirac), prefix(prefix) { }
    virtual ~PreconditionedSolver() { delete solver; }

    void operator()(ColorSpinorField &x, ColorSpinorField &b) {
      setOutputPrefix(prefix);

      QudaSolutionType solution_type = b.SiteSubset() == QUDA_FULL_SITE_SUBSET ? QUDA_MAT_SOLUTION : QUDA_MATPC_SOLUTION;

      ColorSpinorField *out=nullptr;
      ColorSpinorField *in=nullptr;

      dirac.prepare(in, out, x, b, solution_type);
      (*solver)(*out, *in);
      dirac.reconstruct(x, b, solution_type);

      setOutputPrefix("");
    }
  };


  class MultiShiftSolver {

  protected:
    SolverParam &param;
    TimeProfile &profile;

  public:
    MultiShiftSolver(SolverParam &param, TimeProfile &profile) :
    param(param), profile(profile) { ; }
    virtual ~MultiShiftSolver() { ; }

    virtual void operator()(std::vector<ColorSpinorField*> out, ColorSpinorField &in) = 0;
    bool convergence(const double *r2, const double *r2_tol, int n) const;
  };

  class MultiShiftCG : public MultiShiftSolver {

  protected:
    const DiracMatrix &mat;
    const DiracMatrix &matSloppy;

  public:
    MultiShiftCG(DiracMatrix &mat, DiracMatrix &matSloppy, SolverParam &param, TimeProfile &profile);
    virtual ~MultiShiftCG();

    void operator()(std::vector<ColorSpinorField*> out, ColorSpinorField &in);
  };



  class MinResExt {

  protected:
    const DiracMatrix &mat;
    bool orthogonal; 
    bool apply_mat; 
    TimeProfile &profile;

  public:
    MinResExt(DiracMatrix &mat, bool orthogonal, bool apply_mat, TimeProfile &profile);
    virtual ~MinResExt();

    void operator()(ColorSpinorField &x, ColorSpinorField &b,
        std::vector<std::pair<ColorSpinorField*,ColorSpinorField*> > basis);

    void operator()(ColorSpinorField &x, ColorSpinorField &b,
        std::vector<ColorSpinorField*> p,
        std::vector<ColorSpinorField*> q);
  };

  using ColorSpinorFieldSet = ColorSpinorField;

  //forward declaration
  class EigCGArgs;

  class IncEigCG : public Solver {

  private:
    DiracMatrix &mat;
    DiracMatrix &matSloppy;
    DiracMatrix &matPrecon;

    Solver *K;
    SolverParam Kparam; // parameters for preconditioner solve

    ColorSpinorFieldSet *Vm;  //eigCG search vectors  (spinor matrix of size eigen_vector_length x m)

    ColorSpinorField *rp;       
    ColorSpinorField *yp;       
    ColorSpinorField* p;  // conjugate vector
    ColorSpinorField* Ap; // mat * conjugate vector
    ColorSpinorField *tmpp;     
    ColorSpinorField* Az; // mat * conjugate vector from the previous iteration 
    ColorSpinorField *r_pre;    
    ColorSpinorField *p_pre;    

    EigCGArgs *eigcg_args;

    TimeProfile &profile;    //time profile for initCG solver

    bool init;

  public:

    IncEigCG(DiracMatrix &mat, DiracMatrix &matSloppy, DiracMatrix &matPrecon, SolverParam &param, TimeProfile &profile);

    virtual ~IncEigCG();

    void RestartVT(const double beta, const double rho);
    void UpdateVm(ColorSpinorField &res, double beta, double sqrtr2); 
    //EigCG solver:
    int eigCGsolve(ColorSpinorField &out, ColorSpinorField &in);
    //InitCG solver:
    int initCGsolve(ColorSpinorField &out, ColorSpinorField &in);
    //Incremental eigCG solver (for eigcg and initcg calls)
    void operator()(ColorSpinorField &out, ColorSpinorField &in);
  };

//forward declaration
 class GMResDRArgs;

 class GMResDR : public Solver {

  private:

    DiracMatrix &mat;
    DiracMatrix &matSloppy;
    DiracMatrix &matPrecon;

    Solver *K;
    SolverParam Kparam; // parameters for preconditioner solve

    ColorSpinorFieldSet *Vm;//arnoldi basis vectors, size (m+1)
    ColorSpinorFieldSet *Zm;//arnoldi basis vectors, size (m+1)

    ColorSpinorField *rp;       
    ColorSpinorField *yp;       
    ColorSpinorField *tmpp;     
    ColorSpinorField *r_sloppy; 
    ColorSpinorField *r_pre;    
    ColorSpinorField *p_pre;    

    TimeProfile &profile;    //time profile for initCG solver

    GMResDRArgs *gmresdr_args;

    bool init;

  public:

    GMResDR(DiracMatrix &mat, DiracMatrix &matSloppy, DiracMatrix &matPrecon, SolverParam &param, TimeProfile &profile);
    GMResDR(DiracMatrix &mat, Solver &K, DiracMatrix &matSloppy, DiracMatrix &matPrecon, SolverParam &param, TimeProfile &profile);

    virtual ~GMResDR();

    //GMRES-DR solver
    void operator()(ColorSpinorField &out, ColorSpinorField &in);
    //
    //GMRESDR method
    void RunDeflatedCycles (ColorSpinorField *out, ColorSpinorField *in, const double tol_threshold);
    //
    int FlexArnoldiProcedure (const int start_idx, const bool do_givens);

    void RestartVZH();

    void UpdateSolution(ColorSpinorField *x, ColorSpinorField *r, bool do_gels);

  };

} // namespace quda

#endif // _INVERT_QUDA_H