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

namespace quda
{

  class EigenSolver
  {

protected:
    QudaEigParam *eig_param;
    TimeProfile &profile;

    // Problem parameters
    //------------------
    int nEv;          
    int nKr;          
    int nConv;        
    double tol;       
    bool reverse;     
    char spectrum[3]; 
    // Algorithm variables
    //--------------------
    bool converged;
    int restart_iter;
    int max_restarts;
    int check_interval;
    int iter;
    int iter_converged;
    int iter_locked;
    int iter_keep;
    int num_converged;
    int num_locked;
    int num_keep;

    double *residua;

    // Device side vector workspace
    std::vector<ColorSpinorField *> r;
    std::vector<ColorSpinorField *> d_vecs_tmp;

    ColorSpinorField *tmp1;
    ColorSpinorField *tmp2;

    Complex *Qmat;

public:
    EigenSolver(QudaEigParam *eig_param, TimeProfile &profile);

    virtual ~EigenSolver();

    virtual void operator()(std::vector<ColorSpinorField *> &kSpace, std::vector<Complex> &evals) = 0;

    static EigenSolver *create(QudaEigParam *eig_param, const DiracMatrix &mat, TimeProfile &profile);

    void matVec(const DiracMatrix &mat, ColorSpinorField &out, const ColorSpinorField &in);

    void chebyOp(const DiracMatrix &mat, ColorSpinorField &out, const ColorSpinorField &in);

    Complex blockOrthogonalize(std::vector<ColorSpinorField *> v, std::vector<ColorSpinorField *> r, int j);

    void deflate(std::vector<ColorSpinorField *> vec_defl, std::vector<ColorSpinorField *> vec,
                 std::vector<ColorSpinorField *> evecs, std::vector<Complex> evals);

    void deflateSVD(std::vector<ColorSpinorField *> vec_defl, std::vector<ColorSpinorField *> vec,
                    std::vector<ColorSpinorField *> evecs, std::vector<Complex> evals);

    void computeSVD(const DiracMatrix &mat, std::vector<ColorSpinorField *> &evecs, std::vector<Complex> &evals);

    void computeEvals(const DiracMatrix &mat, std::vector<ColorSpinorField *> &evecs, std::vector<Complex> &evals, int k);

    static void loadVectors(std::vector<ColorSpinorField *> &eig_vecs, std::string file);

    static void saveVectors(const std::vector<ColorSpinorField *> &eig_vecs, std::string file);

    void loadFromFile(const DiracMatrix &mat, std::vector<ColorSpinorField *> &eig_vecs, std::vector<Complex> &evals);
  };

  class TRLM : public EigenSolver
  {

public:
    const DiracMatrix &mat;
    TRLM(QudaEigParam *eig_param, const DiracMatrix &mat, TimeProfile &profile);

    virtual ~TRLM();

    // Variable size matrix
    std::vector<double> ritz_mat;

    // Tridiagonal/Arrow matrix, fixed size.
    double *alpha;
    double *beta;

    // Used to clone vectors and resize arrays.
    ColorSpinorParam csParam;

    void operator()(std::vector<ColorSpinorField *> &kSpace, std::vector<Complex> &evals);

    void lanczosStep(std::vector<ColorSpinorField *> v, int j);

    void reorder(std::vector<ColorSpinorField *> &kSpace);

    void eigensolveFromArrowMat(int nLocked, int arror_pos);

    void computeKeptRitz(std::vector<ColorSpinorField *> &kSpace);

  };

  void arpack_solve(std::vector<ColorSpinorField *> &h_evecs, std::vector<Complex> &h_evals, const DiracMatrix &mat,
                    QudaEigParam *eig_param, TimeProfile &profile);

} // namespace quda