Refactor/expansion of MG setup methods

include/blas_quda.h

@@ -202,8 +202,7 @@ namespace quda {
     }
 
     /**
-       @brief Compute the block "caxpy" with over the s
-et of
+       @brief Compute the block "caxpy" with over the set of
        ColorSpinorFields.  E.g., it computes
 
        y = x * a + y

include/dirac_quda.h

@@ -1300,7 +1300,7 @@ namespace quda {
        @brief Return the one-hop field for staggered operators for MG setup
 
        @return Gauge field
-   */
+    */
     virtual cudaGaugeField *getStaggeredShortLinkField() const { return gauge; }
 
     /**

include/enum_quda.h

@@ -197,12 +197,12 @@ typedef enum QudaResidualType_s {
   QUDA_INVALID_RESIDUAL = QUDA_INVALID_ENUM
 } QudaResidualType;
 
-// Which basis to use for CA algorithms
-typedef enum QudaCABasis_s {
+// Which basis to use for polynomials, CA solver basis
+typedef enum QudaPolynomialBasis_s {
   QUDA_POWER_BASIS,
   QUDA_CHEBYSHEV_BASIS,
   QUDA_INVALID_BASIS = QUDA_INVALID_ENUM
-} QudaCABasis;
+} QudaPolynomialBasis;
 
 // Whether the preconditioned matrix is (1-k^2 Deo Doe) or (1-k^2 Doe Deo)
 //
@@ -446,17 +446,15 @@ typedef enum QudaDeflatedGuess_s {
   QUDA_DEFLATED_GUESS_INVALID = QUDA_INVALID_ENUM
 } QudaDeflatedGuess;
 
-typedef enum QudaComputeNullVector_s {
-  QUDA_COMPUTE_NULL_VECTOR_NO,
-  QUDA_COMPUTE_NULL_VECTOR_YES,
-  QUDA_COMPUTE_NULL_VECTOR_INVALID = QUDA_INVALID_ENUM
-} QudaComputeNullVector;
-
-typedef enum QudaSetupType_s {
-  QUDA_NULL_VECTOR_SETUP,
-  QUDA_TEST_VECTOR_SETUP,
-  QUDA_INVALID_SETUP_TYPE = QUDA_INVALID_ENUM
-} QudaSetupType;
+typedef enum QudaNullVectorSetupType_s {
+  QUDA_SETUP_NULL_VECTOR_INVERSE_ITERATIONS,
+  QUDA_SETUP_NULL_VECTOR_CHEBYSHEV_FILTER,
+  QUDA_SETUP_NULL_VECTOR_EIGENVECTORS,
+  QUDA_SETUP_NULL_VECTOR_TEST_VECTORS,
+  QUDA_SETUP_NULL_VECTOR_RESTRICT_FINE,
+  QUDA_SETUP_NULL_VECTOR_FREE_FIELD,
+  QUDA_SETUP_NULL_VECTOR_INVALID = QUDA_INVALID_ENUM
+} QudaNullVectorSetupType;
 
 typedef enum QudaTransferType_s {
   QUDA_TRANSFER_AGGREGATE,

include/enum_quda_fortran.h

@@ -181,7 +181,7 @@
 #define QUDA_HEAVY_QUARK_RESIDUAL 4
 #define QUDA_INVALID_RESIDUAL QUDA_INVALID_ENUM
 
-#define QudaCABasis integer(4)
+#define QudaPolynomialBasis integer(4)
 #define QUDA_POWER_BASIS 0
 #define QUDA_CHEBYSHEV_BASIS 1
 #define QUDA_INVALID_BASIS QUDA_INVALID_ENUM
@@ -406,15 +406,14 @@
 #define QUDA_USE_INIT_GUESS_YES 1
 #define QUDA_USE_INIT_GUESS_INVALID QUDA_INVALID_ENUM
 
-#define QudaComputeNullVector integer(4)
-#define QUDA_COMPUTE_NULL_VECTOR_NO  0 
-#define QUDA_COMPUTE_NULL_VECTOR_YES 1
-#define QUDA_COMPUTE_NULL_VECTOR_INVALID QUDA_INVALID_ENUM
-
-#define QudaSetupType integer(4)
-#define QUDA_NULL_VECTOR_SETUP 0
-#define QUDA_TEST_VECTOR_SETUP 1
-#define QUDA_INVALID_SETUP_TYPE QUDA_INVALID_ENUM
+#define QudaNullVectorSetupType integer(4)
+#define QUDA_SETUP_NULL_VECTOR_INVERSE_ITERATIONS  0 
+#define QUDA_SETUP_NULL_VECTOR_CHEBYSHEV_FILTER 1
+#define QUDA_SETUP_NULL_VECTOR_EIGENVECTOES 2
+#define QUDA_SETUP_NULL_VECTOR_TEST_VECTORS 3
+#define QUDA_SETUP_NULL_VECTOR_RESTRICT_FINE 4
+#define QUDA_SETUP_NULL_VECTOR_FREE_FIELD 5
+#define QUDA_SETUP_NULL_VECTOR_INVALID QUDA_INVALID_ENUM
 
 #define QudaTransferType integer(4)
 #define QUDA_TRANSFER_AGGREGATE 0

include/invert_quda.h

@@ -60,8 +60,8 @@ namespace quda {
     /**< Whether to use an initial guess in the solver or not */
     QudaUseInitGuess use_init_guess;
 
-    /**< Whether to solve linear system with zero RHS */
-    QudaComputeNullVector compute_null_vector;
+    /**< Whether or not to allow a zero RHS solve for near-null vector generation */
+    bool compute_null_vector;
 
     /**< Reliable update tolerance */
     double delta;
@@ -202,7 +202,7 @@ namespace quda {
     double omega;
 
     /** Basis for CA algorithms */
-    QudaCABasis ca_basis;
+    QudaPolynomialBasis ca_basis;
 
     /** Minimum eigenvalue for Chebyshev CA basis */
     double ca_lambda_min;
@@ -211,7 +211,7 @@ namespace quda {
     double ca_lambda_max; // -1 -> power iter generate
 
     /** Basis for CA algorithms in a preconditioner */
-    QudaCABasis ca_basis_precondition;
+    QudaPolynomialBasis ca_basis_precondition;
 
     /** Minimum eigenvalue for Chebyshev CA basis in a preconditioner */
     double ca_lambda_min_precondition;
@@ -269,7 +269,7 @@ namespace quda {
        Default constructor
      */
     SolverParam() :
-      compute_null_vector(QUDA_COMPUTE_NULL_VECTOR_NO),
+      compute_null_vector(false),
       compute_true_res(true),
       sloppy_converge(false),
       verbosity_precondition(QUDA_SILENT),
@@ -291,7 +291,7 @@ namespace quda {
       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),
+      compute_null_vector(false),
       delta(param.reliable_delta),
       use_alternative_reliable(param.use_alternative_reliable),
       use_sloppy_partial_accumulator(param.use_sloppy_partial_accumulator),
@@ -708,40 +708,6 @@ namespace quda {
     */
     void setRecomputeEvals(bool flag) { recompute_evals = flag; };
 
-    /**
-       @brief Compute power iterations on a Dirac matrix
-       @param[in] diracm Dirac matrix used for power iterations
-       @param[in] start Starting rhs for power iterations; value preserved unless it aliases tempvec1 or tempvec2
-       @param[in,out] tempvec1 Temporary vector used for power iterations (FIXME: can become a reference when std::swap
-       can be used on ColorSpinorField)
-       @param[in,out] tempvec2 Temporary vector used for power iterations (FIXME: can become a reference when std::swap
-       can be used on ColorSpinorField)
-       @param[in] niter Total number of power iteration iterations
-       @param[in] normalize_freq Frequency with which intermediate vector gets normalized
-       @param[in] args Parameter pack of ColorSpinorFields used as temporary passed to Dirac
-       @return Norm of final power iteration result
-    */
-    template <typename... Args>
-    static double performPowerIterations(const DiracMatrix &diracm, const ColorSpinorField &start,
-                                         ColorSpinorField &tempvec1, ColorSpinorField &tempvec2, int niter,
-                                         int normalize_freq, Args &&...args);
-
-    /**
-       @brief Generate a Krylov space in a given basis
-       @param[in] diracm Dirac matrix used to generate the Krylov space
-       @param[out] Ap dirac matrix times the Krylov basis vectors
-       @param[in,out] p Krylov basis vectors; assumes p[0] is in place
-       @param[in] n_krylov Size of krylov space
-       @param[in] basis Basis type
-       @param[in] m_map Slope mapping for Chebyshev basis; ignored for power basis
-       @param[in] b_map Intercept mapping for Chebyshev basis; ignored for power basis
-       @param[in] args Parameter pack of ColorSpinorFields used as temporary passed to Dirac
-    */
-    template <typename... Args>
-    static void computeCAKrylovSpace(const DiracMatrix &diracm, std::vector<ColorSpinorField> &Ap,
-                                     std::vector<ColorSpinorField> &p, int n_krylov, QudaCABasis basis, double m_map,
-                                     double b_map, Args &&...args);
-
     /**
      * @brief Return flops
      * @return flops expended by this operator
@@ -1182,7 +1148,7 @@ namespace quda {
     bool init;
 
     bool lambda_init;
-    QudaCABasis basis;
+    QudaPolynomialBasis basis;
 
     std::vector<double> Q_AQandg; // Fused inner product matrix
     std::vector<double> Q_AS;     // inner product matrix
@@ -1320,7 +1286,7 @@ namespace quda {
     bool init;
 
     bool lambda_init;  // whether or not lambda_max has been initialized
-    QudaCABasis basis; // CA basis
+    QudaPolynomialBasis basis; // CA basis
 
     std::vector<Complex> alpha; // Solution coefficient vectors
 
@@ -1694,4 +1660,10 @@ namespace quda {
  */
  bool is_ca_solver(QudaInverterType type);
 
+ /**
+   @brief Returns if a solver supports deflation or not
+   @return true if solver supports deflation, false otherwise
+ */
+ bool is_deflatable_solver(QudaInverterType type);
+
 } // namespace quda

include/multigrid.h

@@ -286,7 +286,7 @@ namespace quda {
     SolverParam *param_coarse_solver;
 
     /** The coarse-grid representation of the null space vectors */
-    std::vector<ColorSpinorField*> *B_coarse;
+    std::vector<ColorSpinorField*> B_coarse;
 
     /** Residual vector */
     ColorSpinorField *r;
@@ -449,7 +449,7 @@ namespace quda {
 
     /**
        @brief Load the null space vectors in from file
-       @param B Loaded null-space vectors (pre-allocated)
+       @param B Load null-space vectors to here
     */
     void loadVectors(std::vector<ColorSpinorField *> &B);
 
@@ -460,22 +460,59 @@ namespace quda {
     void saveVectors(const std::vector<ColorSpinorField *> &B) const;
 
     /**
-       @brief Generate the null-space vectors
+       @brief Initialize a set of vectors to random noise
+       @param B set of vectors
+    */
+    void initializeRandomVectors(const std::vector<ColorSpinorField*> &B) const;
+
+    /**
+       @brief Wrapping function that manages logic for constructing near-null vectors
+       @param refresh whether or not we're only refreshing near null vectors
+    */
+    void constructNearNulls(bool refresh = false);
+
+    /**
+       @brief Generate the null-space vectors via inverse iterations
        @param B Generated null-space vectors
-       @param refresh Whether we refreshing pre-exising vectors or starting afresh
+       @param iterations if non-zero, override the max number of iterations
     */
-    void generateNullVectors(std::vector<ColorSpinorField*> &B, bool refresh=false);
+    void generateInverseIterations(std::vector<ColorSpinorField*> &B, int iterations = 0);
+
+    /**
+       @brief Generate the null-space vectors via a Chebyshev filter; this approach is
+              based on arXiv:2103.05034, P. Boyle and A. Yamaguchi.
+       @param B Generated null-space vectors
+    */
+    void generateChebyshevFilter(std::vector<ColorSpinorField*> &B);
 
     /**
        @brief Generate lowest eigenvectors
+       @param B Generated null-space vectors
+       @param post_restrict whether or not we're generating extra eigenvectors after restricting
+              some fine eigenvectors; if true we override the requested n_conv in the multigrid
+              eigenvalue struct
     */
-    void generateEigenVectors();
+    void generateEigenvectors(std::vector<ColorSpinorField*> &B, bool post_restrict = false);
+
+    /**
+       @brief Generate near-null vectors via restricting finer near-nulls, generating extras if need be
+       @param refresh Whether or not we're only refreshing extra near nulls
+    */
+    void generateRestrictedVectors(bool refresh = false);
 
     /**
        @brief Build free-field null-space vectors
-       @param B Free-field null-space vectors
+       @param B Generated null-space vectors
+    */
+    void generateFreeVectors(std::vector<ColorSpinorField*> &B);
+
+    /**
+      @brief Orthonormalize a vector of ColorSpinorField, erroring out if
+             the fields are not sufficiently linearly independent
+      @param vecs vector of ColorSpinorFields to normalize
+      @param count number of near-null vectors to orthonormalize, default all
     */
-    void buildFreeVectors(std::vector<ColorSpinorField*> &B);
+    void orthonormalizeVectors(const std::vector<ColorSpinorField*>& vecs, int count = -1) const;
 
     /**
        @brief Return the total flops done on this and all coarser levels.
@@ -495,6 +532,14 @@ namespace quda {
 
       return (param.level == 0 || kd_nearnull_gen);
     }
+
+    /**
+      @brief Return if we're on the coarsest grid right now
+    */
+    bool is_coarsest_grid() const {
+      return (param.level == (param.Nlevel - 1));
+    }
+
   };
 
   /**