invert_quda.h

Go to the documentation of this file.

#pragma once
 
#include <quda.h>
#include <quda_internal.h>
#include <dirac_quda.h>
#include <color_spinor_field.h>
#include <qio_field.h>
#include <eigensolve_quda.h>
#include <vector>
#include <memory>
 
namespace quda {
 
  struct SolverParam {
 
    QudaInverterType inv_type;
 
    QudaInverterType inv_type_precondition;
 
    void *preconditioner;
 
    void *deflation_op;
 
    QudaResidualType residual_type;
 
    bool deflate;
 
    QudaEigParam eig_param;
 
    QudaUseInitGuess use_init_guess;
 
    QudaComputeNullVector compute_null_vector;
 
    double delta;
 
    bool use_alternative_reliable;
 
    bool use_sloppy_partial_accumulator;
 
    int solution_accumulator_pipeline;
 
    int max_res_increase;
 
    int max_res_increase_total;
 
    int max_hq_res_increase;
 
    int max_hq_res_restart_total;
 
    int heavy_quark_check;
 
    int pipeline;
 
    double tol;
 
    double tol_restart;
 
    double tol_hq;
 
    bool compute_true_res;
 
    bool sloppy_converge;
 
    double true_res;
 
    double true_res_hq;
 
    int maxiter;
 
    int iter;
 
    QudaPrecision precision;
 
    QudaPrecision precision_sloppy;
 
    QudaPrecision precision_refinement_sloppy;
 
    QudaPrecision precision_precondition;
 
    QudaPrecision precision_eigensolver;
 
    QudaPreserveSource preserve_source;
 
    bool return_residual;
 
    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;
 
    QudaCABasis ca_basis;
 
    double ca_lambda_min;
 
    double ca_lambda_max; // -1 -> power iter generate
 
    QudaSchwarzType schwarz_type;
 
    double secs;
 
    double gflops;
 
    // Incremental EigCG solver parameters
    QudaPrecision precision_ritz;//also search space precision
 
    int n_ev; // 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; 
 
    bool mg_instance;
 
    QudaExtLibType extlib_type;
 
    SolverParam() :
      compute_null_vector(QUDA_COMPUTE_NULL_VECTOR_NO),
      compute_true_res(true),
      sloppy_converge(false),
      verbosity_precondition(QUDA_SILENT),
      mg_instance(false)
    {
      ;
    }
 
    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),
      deflate(param.eig_param != 0),
      use_init_guess(param.use_init_guess),
      compute_null_vector(QUDA_COMPUTE_NULL_VECTOR_NO),
      delta(param.reliable_delta),
      use_alternative_reliable(param.use_alternative_reliable),
      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),
      max_hq_res_increase(param.max_hq_res_increase),
      max_hq_res_restart_total(param.max_hq_res_restart_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),
      sloppy_converge(false),
      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_refinement_sloppy(param.cuda_prec_refinement_sloppy),
      precision_precondition(param.cuda_prec_precondition),
      precision_eigensolver(param.cuda_prec_eigensolver),
      preserve_source(param.preserve_source),
      return_residual(preserve_source == QUDA_PRESERVE_SOURCE_NO ? true : false),
      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),
      ca_basis(param.ca_basis),
      ca_lambda_min(param.ca_lambda_min),
      ca_lambda_max(param.ca_lambda_max),
      schwarz_type(param.schwarz_type),
      secs(param.secs),
      gflops(param.gflops),
      precision_ritz(param.cuda_prec_ritz),
      n_ev(param.n_ev),
      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),
      mg_instance(false),
      extlib_type(param.extlib_type)
    {
      if (deflate) { eig_param = *(static_cast<QudaEigParam *>(param.eig_param)); }
      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),
      deflate(param.deflate),
      eig_param(param.eig_param),
      use_init_guess(param.use_init_guess),
      compute_null_vector(param.compute_null_vector),
      delta(param.delta),
      use_alternative_reliable(param.use_alternative_reliable),
      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),
      sloppy_converge(param.sloppy_converge),
      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_refinement_sloppy(param.precision_refinement_sloppy),
      precision_precondition(param.precision_precondition),
      precision_eigensolver(param.precision_eigensolver),
      preserve_source(param.preserve_source),
      return_residual(param.return_residual),
      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),
      ca_basis(param.ca_basis),
      ca_lambda_min(param.ca_lambda_min),
      ca_lambda_max(param.ca_lambda_max),
      schwarz_type(param.schwarz_type),
      secs(param.secs),
      gflops(param.gflops),
      precision_ritz(param.precision_ritz),
      n_ev(param.n_ev),
      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),
      mg_instance(param.mg_instance),
      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;
 
      param.ca_lambda_min = ca_lambda_min;
      param.ca_lambda_max = ca_lambda_max;
 
      if (deflate) *static_cast<QudaEigParam *>(param.eig_param) = eig_param;
    }
 
    void updateRhsIndex(QudaInvertParam &param) {
      //for incremental eigCG:
      rhs_idx = param.rhs_idx;
    }
 
  };
 
  class Solver {
 
  protected:
    const DiracMatrix &mat;
    const DiracMatrix &matSloppy;
    const DiracMatrix &matPrecon;
    const DiracMatrix &matEig;
 
    SolverParam &param;
    TimeProfile &profile;
    int node_parity;
    EigenSolver *eig_solve; 
    bool deflate_init;      
    bool deflate_compute;   
    bool recompute_evals;   
    std::vector<ColorSpinorField *> evecs; 
    std::vector<Complex> evals;            
  public:
    Solver(const DiracMatrix &mat, const DiracMatrix &matSloppy, const DiracMatrix &matPrecon,
           const DiracMatrix &matEig, SolverParam &param, TimeProfile &profile);
    virtual ~Solver();
 
    virtual void operator()(ColorSpinorField &out, ColorSpinorField &in) = 0;
 
    virtual void blocksolve(ColorSpinorField &out, ColorSpinorField &in);
 
    const DiracMatrix &M() { return mat; }
    const DiracMatrix &Msloppy() { return matSloppy; }
    const DiracMatrix &Mprecon() { return matPrecon; }
    const DiracMatrix &Meig() { return matEig; }
 
    virtual bool hermitian() = 0;
 
    static Solver *create(SolverParam &param, const DiracMatrix &mat, const DiracMatrix &matSloppy,
                          const DiracMatrix &matPrecon, const DiracMatrix &matEig, TimeProfile &profile);
 
    static double stopping(double tol, double b2, QudaResidualType residual_type);
 
    bool convergence(double r2, double hq2, double r2_tol, double hq_tol);
 
    bool convergenceHQ(double r2, double hq2, double r2_tol, double hq_tol);
 
    bool convergenceL2(double r2, double hq2, double r2_tol, double hq_tol);
 
    void PrintStats(const char *name, int k, double r2, double b2, double hq2);
 
    void PrintSummary(const char *name, int k, double r2, double b2, double r2_tol, double hq_tol);
 
    double precisionEpsilon(QudaPrecision prec = QUDA_INVALID_PRECISION) const;
 
    void constructDeflationSpace(const ColorSpinorField &meta, const DiracMatrix &mat);
 
    void destroyDeflationSpace();
 
    void extendSVDDeflationSpace();
 
    void injectDeflationSpace(std::vector<ColorSpinorField *> &defl_space);
 
    void extractDeflationSpace(std::vector<ColorSpinorField *> &defl_space);
 
    int deflationSpaceSize() const { return (int)evecs.size(); };
 
    void setDeflateCompute(bool flag) { deflate_compute = flag; };
 
    void setRecomputeEvals(bool flag) { recompute_evals = flag; };
 
    virtual double flops() const { return 0; }
  };
 
  class CG : public Solver {
 
  private:
    // pointers to fields to avoid multiple creation overhead
    ColorSpinorField *yp, *rp, *rnewp, *pp, *App, *tmpp, *tmp2p, *tmp3p, *rSloppyp, *xSloppyp;
    std::vector<ColorSpinorField*> p;
    bool init;
 
  public:
    CG(const DiracMatrix &mat, const DiracMatrix &matSloppy, const DiracMatrix &matPrecon, const DiracMatrix &matEig,
       SolverParam &param, TimeProfile &profile);
    virtual ~CG();
    void operator()(ColorSpinorField &out, ColorSpinorField &in){
      (*this)(out, in, nullptr, 0.0);
    };
 
    void operator()(ColorSpinorField &out, ColorSpinorField &in, ColorSpinorField *p_init, double r2_old_init);
 
    void blocksolve(ColorSpinorField& out, ColorSpinorField& in);
 
    virtual bool hermitian() { return true; } 
  };
 
  class CGNE : public CG
  {
 
  private:
    DiracMMdag mmdag;
    DiracMMdag mmdagSloppy;
    DiracMMdag mmdagPrecon;
    DiracMMdag mmdagEig;
    ColorSpinorField *xp;
    ColorSpinorField *yp;
    bool init;
 
  public:
    CGNE(const DiracMatrix &mat, const DiracMatrix &matSloppy, const DiracMatrix &matPrecon, const DiracMatrix &matEig,
         SolverParam &param, TimeProfile &profile);
    virtual ~CGNE();
 
    void operator()(ColorSpinorField &out, ColorSpinorField &in);
 
    virtual bool hermitian() { return false; } 
  };
 
  class CGNR : public CG
  {
 
  private:
    DiracMdagM mdagm;
    DiracMdagM mdagmSloppy;
    DiracMdagM mdagmPrecon;
    DiracMdagM mdagmEig;
    ColorSpinorField *bp;
    bool init;
 
  public:
    CGNR(const DiracMatrix &mat, const DiracMatrix &matSloppy, const DiracMatrix &matPrecon, const DiracMatrix &matEig,
         SolverParam &param, TimeProfile &profile);
    virtual ~CGNR();
 
    void operator()(ColorSpinorField &out, ColorSpinorField &in);
 
    virtual bool hermitian() { return false; } 
  };
 
  class CG3 : public Solver
  {
 
  private:
    // pointers to fields to avoid multiple creation overhead
    ColorSpinorField *yp, *rp, *tmpp, *ArSp, *rSp, *xSp, *xS_oldp, *tmpSp, *rS_oldp, *tmp2Sp;
    bool init;
 
  public:
    CG3(const DiracMatrix &mat, const DiracMatrix &matSloppy, const DiracMatrix &matPrecon, SolverParam &param,
        TimeProfile &profile);
    virtual ~CG3();
 
    void operator()(ColorSpinorField &out, ColorSpinorField &in);
 
    virtual bool hermitian() { return true; } 
  };
 
  class CG3NE : public CG3
  {
 
  private:
    DiracMMdag mmdag;
    DiracMMdag mmdagSloppy;
    DiracMMdag mmdagPrecon;
    ColorSpinorField *xp;
    ColorSpinorField *yp;
    bool init;
 
  public:
    CG3NE(const DiracMatrix &mat, const DiracMatrix &matSloppy, const DiracMatrix &matPrecon, SolverParam &param,
          TimeProfile &profile);
    virtual ~CG3NE();
 
    void operator()(ColorSpinorField &out, ColorSpinorField &in);
 
    virtual bool hermitian() { return false; } 
  };
 
  class CG3NR : public CG3
  {
 
  private:
    DiracMdagM mdagm;
    DiracMdagM mdagmSloppy;
    DiracMdagM mdagmPrecon;
    ColorSpinorField *bp;
    bool init;
 
  public:
    CG3NR(const DiracMatrix &mat, const DiracMatrix &matSloppy, const DiracMatrix &matPrecon, SolverParam &param,
          TimeProfile &profile);
    virtual ~CG3NR();
 
    void operator()(ColorSpinorField &out, ColorSpinorField &in);
 
    virtual bool hermitian() { return false; } 
  };
 
  class MPCG : public Solver {
    private:
      void computeMatrixPowers(cudaColorSpinorField out[], cudaColorSpinorField &in, int nvec);
      void computeMatrixPowers(std::vector<cudaColorSpinorField>& out, std::vector<cudaColorSpinorField>& in, int nsteps);
 
    public:
      MPCG(const DiracMatrix &mat, SolverParam &param, TimeProfile &profile);
      virtual ~MPCG();
 
      void operator()(ColorSpinorField &out, ColorSpinorField &in);
      virtual bool hermitian() { return true; } 
  };
 
 
  class PreconCG : public Solver {
    private:
      Solver *K;
      SolverParam Kparam; // parameters for preconditioner solve
 
    public:
      PreconCG(const DiracMatrix &mat, const DiracMatrix &matSloppy, const DiracMatrix &matPrecon,
               const DiracMatrix &matEig, SolverParam &param, TimeProfile &profile);
 
      virtual ~PreconCG();
 
      void operator()(ColorSpinorField &out, ColorSpinorField &in);
      virtual bool hermitian() { return true; } 
  };
 
 
  class BiCGstab : public Solver {
 
  private:
    // pointers to fields to avoid multiple creation overhead
    ColorSpinorField *yp, *rp, *pp, *vp, *tmpp, *tp;
    bool init;
 
  public:
    BiCGstab(const DiracMatrix &mat, const DiracMatrix &matSloppy, const DiracMatrix &matPrecon,
             const DiracMatrix &matEig, SolverParam &param, TimeProfile &profile);
    virtual ~BiCGstab();
 
    void operator()(ColorSpinorField &out, ColorSpinorField &in);
 
    virtual bool hermitian() { return false; } 
  };
 
  class SimpleBiCGstab : public Solver {
 
  private:
 
    // pointers to fields to avoid multiple creation overhead
    cudaColorSpinorField *yp, *rp, *pp, *vp, *tmpp, *tp;
    bool init;
 
  public:
    SimpleBiCGstab(const DiracMatrix &mat, SolverParam &param, TimeProfile &profile);
    virtual ~SimpleBiCGstab();
 
    void operator()(ColorSpinorField &out, ColorSpinorField &in);
 
    virtual bool hermitian() { return false; } 
  };
 
  class MPBiCGstab : public Solver {
 
  private:
    void computeMatrixPowers(std::vector<cudaColorSpinorField>& pr, cudaColorSpinorField& p, cudaColorSpinorField& r, int nsteps);
 
  public:
    MPBiCGstab(const DiracMatrix &mat, SolverParam &param, TimeProfile &profile);
    virtual ~MPBiCGstab();
 
    void operator()(ColorSpinorField &out, ColorSpinorField &in);
 
    virtual bool hermitian() { return false; } 
  };
 
  class BiCGstabL : public Solver {
 
  private:
    int n_krylov; // 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 n_krylov);
    void updateXRend(Complex *gamma, Complex *gamma_prime, Complex *gamma_prime_prime,
                     std::vector<ColorSpinorField *> r, ColorSpinorField &x, int n_krylov);
 
    bool init;
 
    std::string solver_name; // holds BiCGstab-l, where 'l' literally equals n_krylov.
 
  public:
    BiCGstabL(const DiracMatrix &mat, const DiracMatrix &matSloppy, SolverParam &param, TimeProfile &profile);
    virtual ~BiCGstabL();
 
    void operator()(ColorSpinorField &out, ColorSpinorField &in);
 
    virtual bool hermitian() { return false; } 
  };
 
  class GCR : public Solver {
 
  private:
    const DiracMdagM matMdagM; // used by the eigensolver
 
    Solver *K;
    SolverParam Kparam; // parameters for preconditioner solve
 
    int n_krylov;
 
    Complex *alpha;
    Complex **beta;
    double *gamma;
 
    bool init;
 
    ColorSpinorField *rp;       
    ColorSpinorField *tmpp;     
    ColorSpinorField *tmp_sloppy; 
    ColorSpinorField *r_sloppy; 
 
    std::vector<ColorSpinorField*> p;  // GCR direction vectors
    std::vector<ColorSpinorField*> Ap; // mat * direction vectors
 
  public:
    GCR(const DiracMatrix &mat, const DiracMatrix &matSloppy, const DiracMatrix &matPrecon, const DiracMatrix &matEig,
        SolverParam &param, TimeProfile &profile);
 
    GCR(const DiracMatrix &mat, Solver &K, const DiracMatrix &matSloppy, const DiracMatrix &matPrecon,
        const DiracMatrix &matEig, SolverParam &param, TimeProfile &profile);
    virtual ~GCR();
 
    void operator()(ColorSpinorField &out, ColorSpinorField &in);
 
    virtual bool hermitian() { return false; } 
  };
 
  class MR : public Solver {
 
  private:
    ColorSpinorField *rp;
    ColorSpinorField *r_sloppy;
    ColorSpinorField *Arp;
    ColorSpinorField *tmpp;
    ColorSpinorField *tmp_sloppy;
    ColorSpinorField *x_sloppy;
    bool init;
 
  public:
    MR(const DiracMatrix &mat, const DiracMatrix &matSloppy, SolverParam &param, TimeProfile &profile);
    virtual ~MR();
 
    void operator()(ColorSpinorField &out, ColorSpinorField &in);
 
    virtual bool hermitian() { return false; } 
  };
 
  class CACG : public Solver {
 
  private:
    bool init;
 
    bool lambda_init;
    QudaCABasis basis;
 
    Complex *Q_AQandg; // Fused inner product matrix
    Complex *Q_AS; // inner product matrix
    Complex *alpha; // QAQ^{-1} g
    Complex *beta; // QAQ^{-1} QpolyS
 
    ColorSpinorField *rp;
    ColorSpinorField *tmpp;
    ColorSpinorField *tmpp2;
    ColorSpinorField *tmp_sloppy;
    ColorSpinorField *tmp_sloppy2;
 
    std::vector<ColorSpinorField*> S;  // residual vectors
    std::vector<ColorSpinorField*> AS; // mat * residual vectors. Can be replaced by a single temporary.
    std::vector<ColorSpinorField*> Q;  // CG direction vectors
    std::vector<ColorSpinorField*> Qtmp; // CG direction vectors for pointer swap
    std::vector<ColorSpinorField*> AQ; // mat * CG direction vectors.
                                       // it's possible to avoid carrying these
                                       // around, but there's a stability penalty,
                                       // and computing QAQ becomes a pain (though
                                       // it does let you fuse the reductions...)
 
    void create(ColorSpinorField &b);
 
    void compute_alpha();
 
    void compute_beta();
 
    int reliable(double &rNorm,  double &maxrr, int &rUpdate, const double &r2, const double &delta);
 
  public:
    CACG(const DiracMatrix &mat, const DiracMatrix &matSloppy, const DiracMatrix &matPrecon, const DiracMatrix &matEig,
         SolverParam &param, TimeProfile &profile);
    virtual ~CACG();
 
    void operator()(ColorSpinorField &out, ColorSpinorField &in);
 
    virtual bool hermitian() { return true; } 
  };
 
  class CACGNE : public CACG {
 
  private:
    DiracMMdag mmdag;
    DiracMMdag mmdagSloppy;
    DiracMMdag mmdagPrecon;
    DiracMMdag mmdagEig;
    ColorSpinorField *xp;
    ColorSpinorField *yp;
    bool init;
 
  public:
    CACGNE(const DiracMatrix &mat, const DiracMatrix &matSloppy, const DiracMatrix &matPrecon,
           const DiracMatrix &matEig, SolverParam &param, TimeProfile &profile);
    virtual ~CACGNE();
 
    void operator()(ColorSpinorField &out, ColorSpinorField &in);
 
    virtual bool hermitian() { return false; } 
  };
 
  class CACGNR : public CACG {
 
  private:
    DiracMdagM mdagm;
    DiracMdagM mdagmSloppy;
    DiracMdagM mdagmPrecon;
    DiracMdagM mdagmEig;
    ColorSpinorField *bp;
    bool init;
 
  public:
    CACGNR(const DiracMatrix &mat, const DiracMatrix &matSloppy, const DiracMatrix &matPrecon,
           const DiracMatrix &matEig, SolverParam &param, TimeProfile &profile);
    virtual ~CACGNR();
 
    void operator()(ColorSpinorField &out, ColorSpinorField &in);
 
    virtual bool hermitian() { return false; } 
  };
 
  class CAGCR : public Solver {
 
  private:
    const DiracMdagM matMdagM; // used by the eigensolver
    bool init;
    const bool use_source; // whether we can reuse the source vector
 
    // Basis. Currently anything except POWER_BASIS causes a warning
    // then swap to POWER_BASIS.
    QudaCABasis basis;
 
    Complex *alpha; // Solution coefficient vectors
 
    ColorSpinorField *rp;
    ColorSpinorField *tmpp;
    ColorSpinorField *tmp_sloppy;
 
    std::vector<ColorSpinorField*> p;  // GCR direction vectors
    std::vector<ColorSpinorField*> q;  // mat * direction vectors
 
    void create(ColorSpinorField &b);
 
    void solve(Complex *psi_, std::vector<ColorSpinorField*> &q, ColorSpinorField &b);
 
public:
  CAGCR(const DiracMatrix &mat, const DiracMatrix &matSloppy, const DiracMatrix &matPrecon, const DiracMatrix &matEig,
        SolverParam &param, TimeProfile &profile);
  virtual ~CAGCR();
 
  void operator()(ColorSpinorField &out, ColorSpinorField &in);
 
  virtual bool hermitian() { return false; } 
  };
 
  // Steepest descent solver used as a preconditioner
  class SD : public Solver {
    private:
      cudaColorSpinorField *Ar;
      cudaColorSpinorField *r;
      cudaColorSpinorField *y;
      bool init;
 
    public:
      SD(const DiracMatrix &mat, SolverParam &param, TimeProfile &profile);
      virtual ~SD();
 
      void operator()(ColorSpinorField &out, ColorSpinorField &in);
 
      virtual bool hermitian() { return false; } 
  };
 
  // Extended Steepest Descent solver used for overlapping DD preconditioning
  class XSD : public Solver
  {
  private:
    cudaColorSpinorField *xx;
    cudaColorSpinorField *bx;
    SD *sd; // extended sd is implemented using standard sd
    bool init;
    int R[4];
 
  public:
    XSD(const DiracMatrix &mat, SolverParam &param, TimeProfile &profile);
    virtual ~XSD();
 
    void operator()(ColorSpinorField &out, ColorSpinorField &in);
 
    virtual bool hermitian() { return false; } 
  };
 
  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(solver.M(), solver.Msloppy(), solver.Mprecon(), solver.Meig(), param, profile),
    solver(&solver),
    dirac(dirac),
    prefix(prefix)
  {
  }
 
    virtual ~PreconditionedSolver() { delete solver; }
 
    void operator()(ColorSpinorField &x, ColorSpinorField &b) {
      pushOutputPrefix(prefix);
 
      QudaSolutionType solution_type = b.SiteSubset() == QUDA_FULL_SITE_SUBSET ? QUDA_MAT_SOLUTION : QUDA_MATPC_SOLUTION;
 
      ColorSpinorField *out=nullptr;
      ColorSpinorField *in=nullptr;
 
      if (dirac.hasSpecialMG()) {
        dirac.prepareSpecialMG(in, out, x, b, solution_type);
      } else {
        dirac.prepare(in, out, x, b, solution_type);
      }
      (*solver)(*out, *in);
      if (dirac.hasSpecialMG()) {
        dirac.reconstructSpecialMG(x, b, solution_type);
      } else {
        dirac.reconstruct(x, b, solution_type);
      }
 
      popOutputPrefix();
    }
 
    Solver &ExposeSolver() const { return *solver; }
 
    virtual bool hermitian() { return solver->hermitian(); } 
  };
 
  class MultiShiftSolver {
 
  protected:
    const DiracMatrix &mat;
    const DiracMatrix &matSloppy;
    SolverParam &param;
    TimeProfile &profile;
 
  public:
    MultiShiftSolver(const DiracMatrix &mat, const DiracMatrix &matSloppy, SolverParam &param, TimeProfile &profile) :
      mat(mat),
      matSloppy(matSloppy),
      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 {
 
  public:
    MultiShiftCG(const DiracMatrix &mat, const DiracMatrix &matSloppy, SolverParam &param, TimeProfile &profile);
    virtual ~MultiShiftCG();
    void operator()(std::vector<ColorSpinorField*>x, ColorSpinorField &b, std::vector<ColorSpinorField*> &p, double* r2_old_array );
 
    void operator()(std::vector<ColorSpinorField*> out, ColorSpinorField &in){
      std::unique_ptr<double[]> r2_old(new double[QUDA_MAX_MULTI_SHIFT]);
      std::vector<ColorSpinorField*> p;
 
      (*this)(out, in, p, r2_old.get());
 
      for (auto &pp : p) delete pp;
    }
 
  };
 
 
  class MinResExt {
 
  protected:
    const DiracMatrix &mat;
    bool orthogonal; 
    bool apply_mat; 
    bool hermitian; 
    TimeProfile &profile;
 
    void solve(Complex *psi_, std::vector<ColorSpinorField*> &p,
               std::vector<ColorSpinorField*> &q, ColorSpinorField &b, bool hermitian);
 
  public:
    MinResExt(const DiracMatrix &mat, bool orthogonal, bool apply_mat, bool hermitian, 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:
    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(const DiracMatrix &mat, const DiracMatrix &matSloppy, const DiracMatrix &matPrecon, SolverParam &param,
           TimeProfile &profile);
 
  virtual ~IncEigCG();
 
  void increment(ColorSpinorField &V, int n_ev);
 
  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);
 
  bool hermitian() { return true; } // EigCG is only for Hermitian systems
  };
 
//forward declaration
 class GMResDRArgs;
 
 class GMResDR : public Solver {
 
  private:
    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(const DiracMatrix &mat, const DiracMatrix &matSloppy, const DiracMatrix &matPrecon, SolverParam &param,
            TimeProfile &profile);
    GMResDR(const DiracMatrix &mat, Solver &K, const DiracMatrix &matSloppy, const 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);
 
    bool hermitian() { return false; } // GMRESDR for any linear system
 };
 
 struct deflation_space : public Object {
   bool svd;                              
   std::vector<ColorSpinorField *> evecs; 
   std::vector<Complex> evals;            
 };
 
} // namespace quda