Refactor/expansion of MG setup methods

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

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

@@ -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;
@@ -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),
@@ -1694,4 +1694,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

@@ -231,7 +231,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;
@@ -394,9 +394,8 @@ namespace quda {
 
     /**
        @brief Load the null space vectors in from file
-       @param B Loaded null-space vectors (pre-allocated)
     */
-    void loadVectors(std::vector<ColorSpinorField *> &B);
+    void loadVectors();
 
     /**
        @brief Save the null space vectors in from file
@@ -405,22 +404,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 orthonormalize_vectors(const std::vector<ColorSpinorField*>& vecs, int count = -1) const;
 
     /**
        @brief Return the total flops done on this and all coarser levels.
@@ -440,6 +476,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));
+    }
+
   };
 
   /**

include/quda.h

@@ -638,6 +638,9 @@ extern "C" {
     /** Maximum number of iterations for refreshing the null-space vectors */
     int setup_maxiter_refresh[QUDA_MAX_MG_LEVEL];
 
+    /** Maximum number of iterations for inverse iteration polishing of null-space vectors */
+    int setup_maxiter_inverse_iterations_polish[QUDA_MAX_MG_LEVEL];
+
     /** Basis to use for CA solver setup */
     QudaCABasis setup_ca_basis[QUDA_MAX_MG_LEVEL];
 
@@ -650,8 +653,29 @@ extern "C" {
     /** Maximum eigenvalue for Chebyshev CA basis */
     double setup_ca_lambda_max[QUDA_MAX_MG_LEVEL];
 
+    /** Number of random starting vectors for Chebyshev filter null space generation */
+    int filter_startup_vectors[QUDA_MAX_MG_LEVEL];
+
+    /** Number of iterations for initial Chebyshev filter */
+    int filter_startup_iterations[QUDA_MAX_MG_LEVEL];
+
+    /** Frequency of rescales during initial filtering which helps avoid overflow */
+    int filter_startup_rescale_frequency[QUDA_MAX_MG_LEVEL];
+
+    /** Number of iterations between null vectors generated from each starting vector */
+    int filter_iterations_between_vectors[QUDA_MAX_MG_LEVEL];
+
+    /** Conservative estimate of largest eigenvalue of operator used for Chebyshev filter setup */
+    double filter_lambda_min[QUDA_MAX_MG_LEVEL];
+
+    /** Lower bound of eigenvalues that are not enhanced by the initial Chebyshev filter */
+    double filter_lambda_max[QUDA_MAX_MG_LEVEL];
+
     /** Null-space type to use in the setup phase */
-    QudaSetupType setup_type;
+    QudaNullVectorSetupType setup_type[QUDA_MAX_MG_LEVEL];
+
+    /** Null-space type to use to "fill in" extra vectors after restricting some */
+    QudaNullVectorSetupType setup_restrict_remaining_type[QUDA_MAX_MG_LEVEL];
 
     /** Pre orthonormalize vectors in the setup phase */
     QudaBoolean pre_orthonormalize;
@@ -714,7 +738,7 @@ extern "C" {
     int smoother_schwarz_cycle[QUDA_MAX_MG_LEVEL];
 
     /** The type of residual to send to the next coarse grid, and thus the
-	type of solution to receive back from this coarse grid */
+        type of solution to receive back from this coarse grid */
     QudaSolutionType coarse_grid_solution_type[QUDA_MAX_MG_LEVEL];
 
     /** The type of smoother solve to do on each grid (e/o preconditioning or not)*/
@@ -740,12 +764,6 @@ extern "C" {
         memory */
     QudaBoolean setup_minimize_memory;
 
-    /** Whether to compute the null vectors or reload them */
-    QudaComputeNullVector compute_null_vector;
-
-    /** Whether to generate on all levels or just on level 0 */
-    QudaBoolean generate_all_levels;
-
     /** Whether to run the verification checks once set up is complete */
     QudaBoolean run_verify;
 

lib/check_params.h

@@ -760,12 +760,6 @@ void printQudaMultigridParam(QudaMultigridParam *param) {
   int n_level = param->n_level;
 #endif
 
-#ifdef INIT_PARAM
-  P(setup_type, QUDA_NULL_VECTOR_SETUP);
-#else
-  P(setup_type, QUDA_INVALID_SETUP_TYPE);
-#endif
-
 #ifdef INIT_PARAM
   P(pre_orthonormalize, QUDA_BOOLEAN_FALSE);
 #else
@@ -779,6 +773,14 @@ void printQudaMultigridParam(QudaMultigridParam *param) {
 #endif
 
   for (int i=0; i<n_level; i++) {
+#ifdef INIT_PARAM
+    P(setup_type[i], QUDA_SETUP_NULL_VECTOR_INVERSE_ITERATIONS);
+    P(setup_restrict_remaining_type[i], QUDA_SETUP_NULL_VECTOR_INVERSE_ITERATIONS);
+#else
+    P(setup_type[i], QUDA_SETUP_NULL_VECTOR_INVALID);
+    P(setup_restrict_remaining_type[i], QUDA_SETUP_NULL_VECTOR_INVALID);
+#endif
+
 #ifdef INIT_PARAM
     P(verbosity[i], QUDA_SILENT);
 #else
@@ -803,10 +805,12 @@ void printQudaMultigridParam(QudaMultigridParam *param) {
     P(setup_tol[i], 5e-6);
     P(setup_maxiter[i], 500);
     P(setup_maxiter_refresh[i], 0);
+    P(setup_maxiter_inverse_iterations_polish[i], 0);
 #else
     P(setup_tol[i], INVALID_DOUBLE);
     P(setup_maxiter[i], INVALID_INT);
     P(setup_maxiter_refresh[i], INVALID_INT);
+    P(setup_maxiter_inverse_iterations_polish[i], INVALID_INT);
 #endif
 
 #ifdef INIT_PARAM
@@ -821,6 +825,23 @@ void printQudaMultigridParam(QudaMultigridParam *param) {
     P(setup_ca_lambda_max[i], INVALID_DOUBLE);
 #endif
 
+#ifdef INIT_PARAM
+    P(filter_startup_vectors[i], 1);
+    P(filter_startup_iterations[i], 1000);
+    P(filter_startup_rescale_frequency[i], 50);
+    P(filter_iterations_between_vectors[i], 150);
+    P(filter_lambda_min[i], 1.0);
+    P(filter_lambda_max[i], -1.0);
+#else
+    P(filter_startup_vectors[i], INVALID_INT);
+    P(filter_startup_iterations[i], INVALID_INT);
+    P(filter_startup_rescale_frequency[i], INVALID_INT);
+    P(filter_iterations_between_vectors[i], INVALID_INT);
+    P(filter_lambda_min[i], INVALID_DOUBLE);
+    P(filter_lambda_max[i], INVALID_DOUBLE);
+#endif
+
+
 #ifdef INIT_PARAM
     P(n_block_ortho[i], 1);
     P(block_ortho_two_pass[i], QUDA_BOOLEAN_TRUE);
@@ -924,20 +945,6 @@ void printQudaMultigridParam(QudaMultigridParam *param) {
   P(setup_minimize_memory, QUDA_BOOLEAN_INVALID);
 #endif
 
-  P(compute_null_vector, QUDA_COMPUTE_NULL_VECTOR_INVALID);
-  P(generate_all_levels, QUDA_BOOLEAN_INVALID);
-
-#ifdef CHECK_PARAM
-  // if only doing top-level null-space generation, check that n_vec
-  // is equal on all levels
-  if (param->generate_all_levels == QUDA_BOOLEAN_FALSE && param->compute_null_vector == QUDA_COMPUTE_NULL_VECTOR_YES) {
-    for (int i=1; i<n_level-1; i++)
-      if (param->n_vec[0] != param->n_vec[i])
-	errorQuda("n_vec %d != %d must be equal on all levels if generate_all_levels == false",
-		  param->n_vec[0], param->n_vec[i]);
-  }
-#endif
-
   P(run_verify, QUDA_BOOLEAN_INVALID);
 
 #ifdef INIT_PARAM

lib/eigensolve_quda.cpp

@@ -337,8 +337,8 @@ namespace quda
     // C_1 is the current 'out' vector.
 
     // Clone 'in' to two temporary vectors.
-    auto tmp1 = std::make_unique<ColorSpinorField>(in);
-    auto tmp2 = std::make_unique<ColorSpinorField>(out);
+    auto tmp_1 = std::make_unique<ColorSpinorField>(in);
+    auto tmp_2 = std::make_unique<ColorSpinorField>(out);
 
     // Using Chebyshev polynomial recursion relation,
     // C_{m+1}(x) = 2*x*C_{m} - C_{m-1}
@@ -355,14 +355,14 @@ namespace quda
 
       // FIXME - we could introduce a fused matVec + blas kernel here, eliminating one temporary
       // mat*C_{m}(x)
-      matVec(mat, out, *tmp2);
+      matVec(mat, out, *tmp_2);
 
-      blas::axpbypczw(d3, *tmp1, d2, *tmp2, d1, out, *tmp1);
-      std::swap(tmp1, tmp2);
+      blas::axpbypczw(d3, *tmp_1, d2, *tmp_2, d1, out, *tmp_1);
+      std::swap(tmp_1, tmp_2);
 
       sigma_old = sigma;
     }
-    blas::copy(out, *tmp2);
+    blas::copy(out, *tmp_2);
 
     // Save Chebyshev tuning
     saveTuneCache();