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;

    std::vector<ColorSpinorField *> evecs;

    std::vector<Complex> evals;

    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;

    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 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; 

    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(false),
      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),
      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),
      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),
      mg_instance(false),
      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),
      deflate(false),
      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),
      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),
      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),
      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;
    }

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

  };

  class Solver {

  protected:
    SolverParam &param;
    TimeProfile &profile;
    int node_parity;

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

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

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

    static Solver* create(SolverParam &param, DiracMatrix &mat, DiracMatrix &matSloppy,
        DiracMatrix &matPrecon, 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);

    EigenSolver *eig_solve;
    bool deflate_init = false;
    std::vector<ColorSpinorField *> defl_tmp1;
    std::vector<ColorSpinorField *> defl_tmp2;

    void constructDeflationSpace(const ColorSpinorField &meta, const DiracMatrix &mat, bool svd);

    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, *rnewp, *pp, *App, *tmpp, *tmp2p, *tmp3p, *rSloppyp, *xSloppyp;
    std::vector<ColorSpinorField*> p;
    bool init;

  public:
    CG(DiracMatrix &mat, DiracMatrix &matSloppy, 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);
  };



  class CG3 : public Solver {

  private:
    const DiracMatrix &mat;
    const DiracMatrix &matSloppy;
    // 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(DiracMatrix &mat, DiracMatrix &matSloppy, SolverParam &param, TimeProfile &profile);
    virtual ~CG3();

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



  class CG3NE : public Solver {

  private:
    const DiracMatrix &mat;
    const DiracMatrix &matSloppy;
    DiracDagger matDagSloppy;
    // pointers to fields to avoid multiple creation overhead
    ColorSpinorField *yp, *rp, *AdagrSp, *AAdagrSp, *rSp, *xSp, *xS_oldp, *tmpSp, *rS_oldp;
    bool init;

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

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

  class CGNE : public CG {

  private:
    DiracMMdag mmdag;
    DiracMMdag mmdagSloppy;
    ColorSpinorField *xp;
    ColorSpinorField *yp;
    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 *y_sloppy; 
    ColorSpinorField *r_sloppy; 

    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 *r_sloppy;
    ColorSpinorField *Arp;
    ColorSpinorField *tmpp;
    ColorSpinorField *tmp_sloppy;
    ColorSpinorField *x_sloppy;
    bool init;

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

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

  class CACG : public Solver {

  private:
    const DiracMatrix &mat;
    const DiracMatrix &matSloppy;
    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(DiracMatrix &mat, DiracMatrix &matSloppy, SolverParam &param, TimeProfile &profile);
    virtual ~CACG();

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

  class CACGNE : public CACG {

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

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

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

  class CACGNR : public CACG {

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

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

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

  class CAGCR : public Solver {

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

    // 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(DiracMatrix &mat, DiracMatrix &matSloppy, SolverParam &param, TimeProfile &profile);
    virtual ~CAGCR();

    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) {
      pushOutputPrefix(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);

      popOutputPrefix();
    }
  };

  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*>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(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:
    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 increment(ColorSpinorField &V, int nev);

    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