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>

namespace quda {

  struct SolverParam {
    QudaInverterType inv_type;

    QudaInverterType inv_type_precondition;
    
    QudaResidualType residual_type;
    
    QudaUseInitGuess use_init_guess;       

    double delta;           

    bool use_sloppy_partial_accumulator;

    int max_res_increase;

    int max_res_increase_total;

    int pipeline;

    double tol;

    double tol_restart;             

    double tol_hq;          

    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;

    // 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 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;
    
    SolverParam(QudaInvertParam &param) : inv_type(param.inv_type), 
      inv_type_precondition(param.inv_type_precondition), 
      residual_type(param.residual_type), use_init_guess(param.use_init_guess),
      delta(param.reliable_delta), use_sloppy_partial_accumulator(param.use_sloppy_partial_accumulator), 
      max_res_increase(param.max_res_increase), max_res_increase_total(param.max_res_increase_total), pipeline(param.pipeline), tol(param.tol), tol_restart(param.tol_restart), tol_hq(param.tol_hq), 
      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_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) 
    { 
      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_EIGCG_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;
      param.gflops = (param.gflops*param.secs + gflops*secs) / (param.secs + secs);
      param.secs += secs;
      if (offset >= 0) {
        param.true_res_offset[offset] = true_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.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()(cudaColorSpinorField &out, cudaColorSpinorField &in) = 0;

    // solver factory
    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);

  };

  class CG : public Solver {

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

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

    void operator()(cudaColorSpinorField &out, cudaColorSpinorField &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()(cudaColorSpinorField &out, cudaColorSpinorField &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()(cudaColorSpinorField &out, cudaColorSpinorField &in);
  };


  class BiCGstab : public Solver {

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

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

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

    void operator()(cudaColorSpinorField &out, cudaColorSpinorField &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()(cudaColorSpinorField &out, cudaColorSpinorField &in);
  };
  
  class MPBiCGstab : public Solver {

  private:
    DiracMatrix &mat;

    // pointers to fields to avoid multiple creation overhead
    cudaColorSpinorField *yp, *rp, *pp, *vp, *tmpp, *tp;
    bool init;
    void computeMatrixPowers(std::vector<cudaColorSpinorField>& pr, cudaColorSpinorField& p, cudaColorSpinorField& r, int nsteps);

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

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



  class GCR : public Solver {

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

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

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

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

  class MR : public Solver {

  private:
    const DiracMatrix &mat;
    cudaColorSpinorField *rp;
    cudaColorSpinorField *Arp;
    cudaColorSpinorField *tmpp;
    bool init;
    bool allocate_r;

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

    void operator()(cudaColorSpinorField &out, cudaColorSpinorField &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()(cudaColorSpinorField &out, cudaColorSpinorField &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()(cudaColorSpinorField &out, cudaColorSpinorField &in);
  };


  // multigrid solver
  class alphaSA : public Solver {

  protected:
    const DiracMatrix &mat;

  public:
    alphaSA(DiracMatrix &mat, SolverParam &param, TimeProfile &profile);
    virtual ~alphaSA() { ; }

    void operator()(cudaColorSpinorField **out, cudaColorSpinorField &in);
  };

  class MultiShiftSolver {

  protected:
    SolverParam &param;
    TimeProfile &profile;

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

    virtual void operator()(cudaColorSpinorField **out, cudaColorSpinorField &in) = 0;
  };

  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()(cudaColorSpinorField **out, cudaColorSpinorField &in);
  };

  class MinResExt {

  protected:
    const DiracMatrix &mat;
    TimeProfile &profile;

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

    void operator()(cudaColorSpinorField &x, cudaColorSpinorField &b, cudaColorSpinorField **p,
                    cudaColorSpinorField **q, int N);
  };

  class DeflatedSolver {

  protected:
    SolverParam &param;
    TimeProfile &profile;

    //WARNING: eigcg_precision may not coinside with param.precision and param.precision_sloppy (both used for the initCG).
    //
    QudaPrecision eigcg_precision;//may be double or single.

  public:
    DeflatedSolver(SolverParam &param, TimeProfile &profile) : param(param), profile(profile) 
    { 
       eigcg_precision = param.precision_sloppy;//for mixed presicion use param.precision_sloppy 
    }

    virtual ~DeflatedSolver() { ; }

    virtual void operator()(cudaColorSpinorField *out, cudaColorSpinorField *in) = 0;

//    virtual void Deflate(cudaColorSpinorField &out, cudaColorSpinorField &in) = 0;//extrenal method (not implemented yet)
    virtual void StoreRitzVecs(void *host_buffer, double *inv_eigenvals, const int *X, QudaInvertParam *inv_par, const int nev, bool cleanResources = false) = 0;//extrenal method

    virtual void CleanResources() = 0;

    // solver factory
    static DeflatedSolver* create(SolverParam &param, DiracMatrix &mat, DiracMatrix &matSloppy, DiracMatrix &matCGSloppy, DiracMatrix &matDeflate, TimeProfile &profile);

    bool convergence(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);

  };

  struct DeflationParam;//Forward declaration

  class IncEigCG : public DeflatedSolver {

  private:
    DiracMatrix &mat;
    DiracMatrix &matSloppy;

    DiracMatrix &matCGSloppy;

    const DiracMatrix &matDefl;

    QudaPrecision search_space_prec;
    cudaColorSpinorField *Vm;  //search vectors  (spinor matrix of size eigen_vector_length x m)

    SolverParam initCGparam; // parameters for initCG solver
    TimeProfile &profile;    //time profile for initCG solver

    bool eigcg_alloc;
    bool use_eigcg;

  public:

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

    virtual ~IncEigCG();

    //EigCG solver
    int EigCG(cudaColorSpinorField &out, cudaColorSpinorField &in);

    //Incremental eigCG solver (for eigcg and initcg calls)
    void operator()(cudaColorSpinorField *out, cudaColorSpinorField *in);

    //Compute  u dH^{-1} u^{dagger}b: 
    void DeflateSpinor(cudaColorSpinorField &out, cudaColorSpinorField &in, DeflationParam *param, bool set2zero = true);
    //
    void DeflateSpinorReduced(cudaColorSpinorField &out, cudaColorSpinorField &in, DeflationParam *param, bool set2zero = true);

    //Deflation space management
    void CreateDeflationSpace(cudaColorSpinorField &eigcgSpinor, DeflationParam *&param);

    //extend projection matrix:
    //compute Q' = DiracM Q, (here U = [V, Q] - total Ritz set)
    //construct H-matrix components with Q'^{dag} Q', V^{dag} Q' and Q'^{dag} V
    //extend H-matrix with the components
    void ExpandDeflationSpace(DeflationParam *param, const int new_nev);
    //
    void DeleteDeflationSpace(DeflationParam *&param);
    //
    void DeleteEigCGSearchSpace();
    //
    void SaveEigCGRitzVecs(DeflationParam *param, bool cleanResources = false);
    //
    void StoreRitzVecs(void *host_buf, double *inv_eigenvals, const int *X, QudaInvertParam *inv_par, const int nev, bool cleanResources = false);
    //
    void CleanResources(); 

    void LoadEigenvectors(DeflationParam *param, int max_nevs, double tol = 1e-3);

    void ReportEigenvalueAccuracy(DeflationParam *param, int nevs_to_print);

  };

} // namespace quda

#endif // _INVERT_QUDA_H