MILC batched deflation

include/quda_milc_interface.h

@@ -52,7 +52,52 @@ extern "C" {
   } QudaInvertArgs_t;
 
   /**
-   * Parameters related to deflated solvers.
+   * Parameters related to deflated linear solvers.
+   */
+  typedef struct {
+    size_t struct_size; /** Size of this struct in bytes. Used to check that host application and QUDA see the same struct size **/
+    double tol_restart;
+    QudaPrecision prec_eigensolver;
+    int poly_deg; /** Degree of the Chebyshev polynomial **/
+    double a_min; /** Range used in polynomial acceleration **/
+    double a_max;
+    QudaBoolean preserve_evals; /** Whether to preserve the evals or recompute them **/
+    int n_ev;                   /** Size of the eigenvector search space **/
+    int n_kr;                   /** Total size of Krylov space **/
+    int n_conv;                 /** Number of requested converged eigenvectors **/
+    int n_ev_deflate;           /** Number of requested converged eigenvectors to use in deflation **/
+    double tol;                 /** Tolerance on the least well known eigenvalue's residual **/
+    int max_restarts;           /** For IRLM/IRAM, quit after n restarts **/
+    int batched_rotate;       /** For the Ritz rotation, the maximal number of extra vectors the solver may allocate **/
+    int block_size;           /** For block method solvers, the block size **/
+    char vec_infile[256];     /** Filename prefix where to load the null-space vectors */
+    char vec_outfile[256];    /** Filename prefix for where to save the null-space vectors */
+    QudaParity vec_in_parity; /** Parity of the incoming eigenvectors **/
+    QudaPrecision save_prec;  /** The precision with which to save the vectors */
+    QudaBoolean partfile;     /** Whether to save eigenvectors in QIO singlefile or partfile format */
+    QudaBoolean io_parity_inflate; /** Whether to inflate single-parity eigen-vector I/O **/
+    QudaBoolean use_norm_op;
+    QudaBoolean use_pc;
+    QudaEigType eig_type;            /** Type of eigensolver algorithm to employ **/
+    QudaEigSpectrumType spectrum;    /** Which part of the spectrum to solve **/
+    double qr_tol;                   /** Tolerance on the QR iteration **/
+    QudaBoolean require_convergence; /** If true, the solver will error out if the convergence criteria are not met **/
+    int check_interval;              /** For IRLM/IRAM, check every nth restart **/
+    QudaBoolean use_dagger;          /** If use_dagger, use Mdag **/
+    QudaBoolean compute_gamma5;     /** Performs the \gamma_5 OP solve by post multiplying the eignvectors with \gamma_5
+                                       before computing the eigenvalues */
+    QudaBoolean compute_svd;        /** Performs an MdagM solve, then constructs the left and right SVD. **/
+    QudaBoolean use_eigen_qr;       /** Use Eigen routines to eigensolve the upper Hessenberg via QR **/
+    QudaBoolean use_poly_acc;       /** Use Polynomial Acceleration **/
+    QudaBoolean arpack_check;       /** In the test function, cross check the device result against ARPACK **/
+    char arpack_logfile[512];       /** For Arpack cross check, name of the Arpack logfile **/
+    int compute_evals_batch_size;   /** The batch size used when computing eigenvalues **/
+    QudaBoolean preserve_deflation; /** Whether to preserve the deflation space between solves **/
+
+  } QudaEigensolverArgs_t;
+
+  /**
+   * Parameters related to EigCG deflated solvers.
    */
 
   typedef struct {
@@ -163,6 +208,11 @@ extern "C" {
    */
   void qudaSetLayout(QudaLayout_t layout);
 
+  /**
+   * Clean up the QUDA deflation space.
+   */
+  void qudaCleanUpDeflationSpace();
+
   /**
    * Destroy the QUDA context.
    */
@@ -363,6 +413,32 @@ extern "C" {
 		  double* const final_rel_resid,
 		  int* num_iters);
 
+  /**
+   * Solve Ax=b with deflation for an improved staggered operator. All fields are fields
+   * passed and returned are host (CPU) field in MILC order.  This
+   * function requires that persistent gauge and clover fields have
+   * been created prior.  This interface is experimental.
+   *
+   * @param[in] external_precision Precision of host fields passed to QUDA (2 - double, 1 - single)
+   * @param[in] quda_precision Precision for QUDA to use (2 - double, 1 - single)
+   * @param[in] mass Fermion mass parameter
+   * @param[in] inv_args Struct setting some solver metadata
+   * @param[in] eig_args Struct setting some eigensolver metadata
+   * @param[in] target_residual Target residual
+   * @param[in] target_relative_residual Target Fermilab residual
+   * @param[in] milc_fatlink Fat-link field on the host
+   * @param[in] milc_longlink Long-link field on the host
+   * @param[in] source Right-hand side source field
+   * @param[out] solution Solution spinor field
+   * @param[in] final_residual True residual
+   * @param[in] final_relative_residual True Fermilab residual
+   * @param[in] num_iters Number of iterations taken
+   */
+  void qudaInvertDeflatable(int external_precision, int quda_precision, double mass, QudaInvertArgs_t inv_args,
+                            QudaEigensolverArgs_t eig_args, double target_residual, double target_fermilab_residual,
+                            const void *const milc_fatlink, const void *const milc_longlink, void *source,
+                            void *solution, double *const final_resid, double *const final_rel_resid, int *num_iters);
+
   /**
    * Prepare a staggered/HISQ multigrid solve with given fat and
    * long links. All fields passed are host (CPU) fields
@@ -455,6 +531,34 @@ extern "C" {
                       int* num_iters,
                       int num_src);
 
+  /**
+   * Solve Ax=b with deflation for an improved staggered operator with many right hand sides.
+   * All fields are fields passed and returned are host (CPU) field in MILC order.
+   * This function requires that persistent gauge and clover fields have
+   * been created prior.  This interface is experimental.
+   *
+   * @param[in] external_precision Precision of host fields passed to QUDA (2 - double, 1 - single)
+   * @param[in] quda_precision Precision for QUDA to use (2 - double, 1 - single)
+   * @param[in] mass Fermion mass parameter
+   * @param[in] inv_args Struct setting some solver metadata
+   * @param[in] eig_args Struct setting some eigensolver metadata
+   * @param[in] target_residual Target residual
+   * @param[in] target_relative_residual Target Fermilab residual
+   * @param[in] milc_fatlink Fat-link field on the host
+   * @param[in] milc_longlink Long-link field on the host
+   * @param[in] source array of right-hand side source fields
+   * @param[out] solution array of solution spinor fields
+   * @param[in] final_residual True residual
+   * @param[in] final_relative_residual True Fermilab residual
+   * @param[in] num_iters Number of iterations taken
+   * @param[in] num_src Number of source fields
+   */
+  void qudaInvertMsrcDeflatable(int external_precision, int quda_precision, double mass, QudaInvertArgs_t inv_args,
+                                QudaEigensolverArgs_t eig_args, double target_residual, double target_fermilab_residual,
+                                const void *const fatlink, const void *const longlink, void **sourceArray,
+                                void **solutionArray, double *const final_residual,
+                                double *const final_fermilab_residual, int *num_iters, int num_src);
+
   /**
    * Solve for multiple shifts (e.g., masses) using an improved
    * staggered operator.  All fields are fields passed and returned