Make PPCG usable from PW solver

Author: cheerly-pku

docs/advanced/input_files/input-main.md

@@ -977,7 +977,7 @@
 ### pw_diag_thr
 
 - **Type**: Real
-- **Description**: Only used when you use ks_solver = cg/dav/dav_subspace/bpcg. It indicates the threshold for the first electronic iteration, from the second iteration the pw_diag_thr will be updated automatically. For nscf calculations with planewave basis set, pw_diag_thr should be <= 1e-3.
+- **Description**: Only used when you use ks_solver = cg/dav/dav_subspace/bpcg/ppcg. It indicates the threshold for the first electronic iteration, from the second iteration the pw_diag_thr will be updated automatically. For nscf calculations with planewave basis set, pw_diag_thr should be <= 1e-3.
 - **Default**: 0.01
 
 ### diago_smooth_ethr
@@ -996,8 +996,8 @@
 ### pw_diag_nmax
 
 - **Type**: Integer
-- **Availability**: *basis_type==pw, ks_solver==cg/dav/dav_subspace/bpcg*
-- **Description**: Only useful when you use ks_solver = cg/dav/dav_subspace/bpcg. It indicates the maximal iteration number for cg/david/dav_subspace/bpcg method.
+- **Availability**: *basis_type==pw, ks_solver==cg/dav/dav_subspace/bpcg/ppcg*
+- **Description**: Only useful when you use ks_solver = cg/dav/dav_subspace/bpcg/ppcg. It indicates the maximal iteration number for cg/david/dav_subspace/bpcg/ppcg method.
 - **Default**: 50
 
 ### pw_diag_ndim
@@ -1112,6 +1112,7 @@
   - bpcg: The BPCG method, which is a block-parallel Conjugate Gradient (CG) method, typically exhibits higher acceleration in a GPU environment.
   - dav: The Davidson algorithm.
   - dav_subspace: The Davidson algorithm without orthogonalization operation, this method is the most recommended for efficiency. pw_diag_ndim can be set to 2 for this method.
+  - ppcg: The projection preconditioned conjugate-gradient method, currently available for CPU plane-wave calculations.
 
   For numerical atomic orbitals basis,
 

docs/parameters.yaml

@@ -521,6 +521,7 @@ parameters:
       * bpcg: The BPCG method, which is a block-parallel Conjugate Gradient (CG) method, typically exhibits higher acceleration in a GPU environment.
       * dav: The Davidson algorithm.
       * dav_subspace: The Davidson algorithm without orthogonalization operation, this method is the most recommended for efficiency. `pw_diag_ndim` can be set to 2 for this method.
+      * ppcg: The projection preconditioned conjugate-gradient method, currently available for CPU plane-wave calculations.
 
       For numerical atomic orbitals basis,
 
@@ -942,7 +943,7 @@ parameters:
     category: Plane wave related variables
     type: Real
     description: |
-      Only used when you use ks_solver = cg/dav/dav_subspace/bpcg. It indicates the threshold for the first electronic iteration, from the second iteration the pw_diag_thr will be updated automatically. For nscf calculations with planewave basis set, pw_diag_thr should be <= 1e-3.
+      Only used when you use ks_solver = cg/dav/dav_subspace/bpcg/ppcg. It indicates the threshold for the first electronic iteration, from the second iteration the pw_diag_thr will be updated automatically. For nscf calculations with planewave basis set, pw_diag_thr should be <= 1e-3.
     default_value: "0.01"
     unit: ""
     availability: ""
@@ -966,10 +967,10 @@ parameters:
     category: Plane wave related variables
     type: Integer
     description: |
-      Only useful when you use ks_solver = cg/dav/dav_subspace/bpcg. It indicates the maximal iteration number for cg/david/dav_subspace/bpcg method.
+      Only useful when you use ks_solver = cg/dav/dav_subspace/bpcg/ppcg. It indicates the maximal iteration number for cg/david/dav_subspace/bpcg/ppcg method.
     default_value: "50"
     unit: ""
-    availability: "basis_type==pw, ks_solver==cg/dav/dav_subspace/bpcg"
+    availability: "basis_type==pw, ks_solver==cg/dav/dav_subspace/bpcg/ppcg"
   - name: pw_diag_ndim
     category: Plane wave related variables
     type: Integer

source/source_hsolver/diago_ppcg.h

@@ -1,6 +1,8 @@
 #ifndef DIAGO_PPCG_H
 #define DIAGO_PPCG_H
 
+#include "source_base/module_device/types.h"
+
 #include <vector>
 #include <functional>
 #include <cmath>
@@ -17,21 +19,17 @@ namespace hsolver {
 // -----------------------------------------------------------------------------
 //
 // Supports two algorithmic strategies:
-//   BLOCK_SUBSPACE — block subspace diagonalization (File 1 approach).
 //   CONJUGATE_GRADIENT — band-by-band Polak-Ribiere CG with line minimization
 //     (File 2 approach).
+//   BLOCK_SUBSPACE — block subspace diagonalization (File 1 approach).
 //
-// The block-subspace strategy tends to be more robust near convergence;
-// conjugate-gradient is more memory efficient for large systems.
+// CONJUGATE_GRADIENT is the default because it is the tested production path.
+// BLOCK_SUBSPACE is kept as an explicit experimental strategy.
 // -----------------------------------------------------------------------------
 
 enum class PpcgStrategy { BLOCK_SUBSPACE, CONJUGATE_GRADIENT };
 
-// Device tags (extensible for GPU backends).
-namespace base_device {
-    struct DEVICE_CPU {};
-    struct DEVICE_GPU {};
-}
+namespace base_device = ::base_device;
 
 template <typename T, typename Device>
 class DiagoPPCG
@@ -54,7 +52,7 @@ class DiagoPPCG
               const int& sbsize,
               const int& rr_step,
               const bool gamma_g0_real,
-              const PpcgStrategy strategy = PpcgStrategy::BLOCK_SUBSPACE);
+              const PpcgStrategy strategy = PpcgStrategy::CONJUGATE_GRADIENT);
 
     // -------------------------------------------------------------------------
     // Main entry point

source/source_hsolver/hsolver_pw.cpp

@@ -11,18 +11,81 @@
 #include "source_hsolver/diago_cg.h"
 #include "source_hsolver/diago_dav_subspace.h"
 #include "source_hsolver/diago_david.h"
+#include "source_hsolver/diago_ppcg.h"
 #include "source_hsolver/diago_iter_assist.h"
 #include "source_io/module_parameter/parameter.h"
 #include "source_psi/psi.h"
 #include "source_estate/elecstate_tools.h"
 
 
 #include <algorithm>
+#include <type_traits>
 #include <vector>
 
 namespace hsolver
 {
 
+namespace
+{
+template <typename T, typename Device, typename Real, typename HPsiFunc, typename SPsiFunc>
+double run_ppcg_pw(const HPsiFunc& hpsi_func,
+                   const SPsiFunc& spsi_func,
+                   const int ld_psi,
+                   const int nband,
+                   const int dim,
+                   T* psi,
+                   Real* eigenvalue,
+                   const std::vector<double>& ethr_band,
+                   const Real* pre_condition,
+                   const double diag_thr,
+                   const int diag_iter_max,
+                   const int pw_diag_ndim,
+                   const bool gamma_only,
+                   std::true_type)
+{
+    const int sbsize = std::max(1, std::min(nband, pw_diag_ndim));
+    const int rr_step = std::max(1, pw_diag_ndim);
+
+    DiagoPPCG<T, Device> ppcg(static_cast<Real>(diag_thr),
+                              diag_iter_max,
+                              sbsize,
+                              rr_step,
+                              gamma_only,
+                              PpcgStrategy::CONJUGATE_GRADIENT);
+
+    return ppcg.diag(hpsi_func,
+                     spsi_func,
+                     ld_psi,
+                     nband,
+                     dim,
+                     psi,
+                     eigenvalue,
+                     ethr_band,
+                     pre_condition);
+}
+
+template <typename T, typename Device, typename Real, typename HPsiFunc, typename SPsiFunc>
+double run_ppcg_pw(const HPsiFunc&,
+                   const SPsiFunc&,
+                   const int,
+                   const int,
+                   const int,
+                   T*,
+                   Real*,
+                   const std::vector<double>&,
+                   const Real*,
+                   const double,
+                   const int,
+                   const int,
+                   const bool,
+                   std::false_type)
+{
+    ModuleBase::WARNING_QUIT("HSolverPW::hamiltSolvePsiK",
+                             "PPCG is currently implemented for CPU PW calculations only.");
+    return 0.0;
+}
+} // namespace
+
 template <typename T, typename Device>
 void HSolverPW<T, Device>::cal_smooth_ethr(const double& wk,
                                            const double* wg,
@@ -83,7 +146,7 @@ void HSolverPW<T, Device>::solve(hamilt::Hamilt<T, Device>* pHamilt,
     this->nproc_in_pool = nproc_in_pool_in;
 
     // report if the specified diagonalization method is not supported
-    const std::initializer_list<std::string> _methods = {"cg", "dav", "dav_subspace", "bpcg"};
+    const std::initializer_list<std::string> _methods = {"cg", "dav", "dav_subspace", "bpcg", "ppcg"};
     if (std::find(std::begin(_methods), std::end(_methods), this->method) == std::end(_methods))
     {
         ModuleBase::WARNING_QUIT("HSolverPW::solve", "This type of eigensolver is not supported!");
@@ -379,6 +442,24 @@ void HSolverPW<T, Device>::hamiltSolvePsiK(hamilt::Hamilt<T, Device>* hm,
                         ntry_max,
                         notconv_max));
     }
+    else if (this->method == "ppcg")
+    {
+        DiagoIterAssist<T, Device>::avg_iter += run_ppcg_pw<T, Device, Real>(
+            hpsi_func,
+            spsi_func,
+            psi.get_nbasis(),
+            psi.get_nbands(),
+            psi.get_current_ngk(),
+            psi.get_pointer(),
+            eigenvalue,
+            this->ethr_band,
+            pre_condition.data(),
+            this->diag_thr,
+            this->diag_iter_max,
+            PARAM.inp.pw_diag_ndim,
+            PARAM.globalv.gamma_only_pw,
+            std::is_same<Device, base_device::DEVICE_CPU>());
+    }
     ModuleBase::timer::end("HSolverPW", "solve_psik");
     return;
 }

source/source_hsolver/test/diago_ppcg_test.cpp

@@ -8,13 +8,9 @@
  *      Exact eigenvalues are the diagonal entries.  Simplest possible
  *      smoke test — should converge in very few iterations.
  *
- * Tests use the CONJUGATE_GRADIENT strategy which has a try/catch fallback
- * for LAPACK sygvd failures and is therefore more portable across different
- * LAPACK implementations.
- *
- * BLOCK_SUBSPACE strategy tests exist in git history but are disabled here
- * because they require a LAPACK with reliable dsygvd for small ill-conditioned
- * generalized eigenvalue problems.
+ * Tests primarily exercise the default CONJUGATE_GRADIENT strategy, with a
+ * BLOCK_SUBSPACE smoke test to keep the explicit experimental path finite on a
+ * small Hermitian problem.
  */
 
 #include "../diago_ppcg.h"
@@ -335,7 +331,7 @@ TEST_F(DiagoPPCG2x2Test, ConjugateGradient)
         << "2x2 CG: too many iterations";
 }
 
-TEST(DiagoPPCGComplexHermitianTest, ConjugateGradientKeepsImaginaryProjection)
+TEST(DiagoPPCGComplexHermitianTest, DefaultKeepsImaginaryProjection)
 {
     const int n_dim = 2;
     const int nband = 2;
@@ -362,8 +358,7 @@ TEST(DiagoPPCGComplexHermitianTest, ConjugateGradientKeepsImaginaryProjection)
         /* max_iter = */ 10,
         /* sbsize   = */ 2,
         /* rr_step  = */ 1,
-        /* gamma_g0 = */ false,
-        hsolver::PpcgStrategy::CONJUGATE_GRADIENT
+        /* gamma_g0 = */ false
     );
 
     auto h_op = [&H_mat, n_dim](T* in, T* out, int ld_in, int ncol) {
@@ -378,6 +373,49 @@ TEST(DiagoPPCGComplexHermitianTest, ConjugateGradientKeepsImaginaryProjection)
     EXPECT_NEAR(eval[1], 2.5 + delta, 1e-10);
 }
 
+TEST(DiagoPPCGComplexHermitianTest, BlockSubspaceSmokeNoNaN)
+{
+    const int n_dim = 2;
+    const int nband = 2;
+    const int ld = n_dim;
+
+    std::vector<T> H_mat(n_dim * n_dim, T(0));
+    H_mat[0 + 0 * n_dim] = T(2.0, 0.0);
+    H_mat[1 + 1 * n_dim] = T(3.0, 0.0);
+    H_mat[0 + 1 * n_dim] = T(0.0, 1.0);
+    H_mat[1 + 0 * n_dim] = T(0.0, -1.0);
+
+    std::vector<T> psi(ld * nband, T(0));
+    psi[0 + 0 * ld] = T(1.0, 0.0);
+    psi[1 + 1 * ld] = T(1.0, 0.0);
+
+    std::vector<Real> prec(n_dim, 2.0);
+    std::vector<double> ethr(nband, 1e-10);
+    std::vector<Real> eval(nband, 0.0);
+
+    hsolver::DiagoPPCG<T, hsolver::base_device::DEVICE_CPU> solver(
+        /* diag_thr = */ 1e-10,
+        /* max_iter = */ 8,
+        /* sbsize   = */ 2,
+        /* rr_step  = */ 1,
+        /* gamma_g0 = */ false,
+        hsolver::PpcgStrategy::BLOCK_SUBSPACE
+    );
+
+    auto h_op = [&H_mat, n_dim](T* in, T* out, int ld_in, int ncol) {
+        dense_h_multiply(H_mat.data(), n_dim, in, out, ld_in, ncol);
+    };
+
+    solver.diag(h_op, nullptr, ld, nband, n_dim,
+                psi.data(), eval.data(), ethr, prec.data());
+
+    const Real delta = std::sqrt(1.25);
+    for (int i = 0; i < nband; ++i)
+        EXPECT_TRUE(std::isfinite(eval[i])) << "BLOCK_SUBSPACE produced NaN/Inf";
+    EXPECT_NEAR(eval[0], 2.5 - delta, 1e-8);
+    EXPECT_NEAR(eval[1], 2.5 + delta, 1e-8);
+}
+
 // =============================================================================
 // Test fixture: degenerate eigenvalues
 // H = I + J  (identity plus all-ones), 4×4.

source/source_io/module_parameter/read_input_item_elec_stru.cpp

@@ -56,6 +56,7 @@ For plane-wave basis,
 * bpcg: The BPCG method, which is a block-parallel Conjugate Gradient (CG) method, typically exhibits higher acceleration in a GPU environment.
 * dav: The Davidson algorithm.
 * dav_subspace: The Davidson algorithm without orthogonalization operation, this method is the most recommended for efficiency. `pw_diag_ndim` can be set to 2 for this method.
+* ppcg: The projection preconditioned conjugate-gradient method, currently available for CPU plane-wave calculations.
 
 For numerical atomic orbitals basis,
 
@@ -131,7 +132,7 @@ Then the user has to correct the input file and restart the calculation.)";
         };
         item.check_value = [](const Input_Item& item, const Parameter& para) {
             const std::string& ks_solver = para.input.ks_solver;
-            const std::vector<std::string> pw_solvers = {"cg", "dav", "bpcg", "dav_subspace"};
+            const std::vector<std::string> pw_solvers = {"cg", "dav", "bpcg", "dav_subspace", "ppcg"};
             const std::vector<std::string> lcao_solvers = {
                 "genelpa",
                 "elpa",
@@ -1040,7 +1041,7 @@ Use case: When experimental or high-level theoretical results suggest that the S
         item.annotation = "threshold for eigenvalues is cg electron iterations";
         item.category = "Plane wave related variables";
         item.type = "Real";
-        item.description = "Only used when you use ks_solver = cg/dav/dav_subspace/bpcg. It indicates the threshold for the first electronic iteration, from the second iteration the pw_diag_thr will be updated automatically. For nscf calculations with planewave basis set, pw_diag_thr should be <= 1e-3.";
+        item.description = "Only used when you use ks_solver = cg/dav/dav_subspace/bpcg/ppcg. It indicates the threshold for the first electronic iteration, from the second iteration the pw_diag_thr will be updated automatically. For nscf calculations with planewave basis set, pw_diag_thr should be <= 1e-3.";
         item.default_value = "0.01";
         item.unit = "";
         item.availability = "";
@@ -1102,10 +1103,10 @@ Use case: When experimental or high-level theoretical results suggest that the S
         item.annotation = "max iteration number for cg";
         item.category = "Plane wave related variables";
         item.type = "Integer";
-        item.description = "Only useful when you use ks_solver = cg/dav/dav_subspace/bpcg. It indicates the maximal iteration number for cg/david/dav_subspace/bpcg method.";
+        item.description = "Only useful when you use ks_solver = cg/dav/dav_subspace/bpcg/ppcg. It indicates the maximal iteration number for cg/david/dav_subspace/bpcg/ppcg method.";
         item.default_value = "50";
         item.unit = "";
-        item.availability = "basis_type==pw, ks_solver==cg/dav/dav_subspace/bpcg";
+        item.availability = "basis_type==pw, ks_solver==cg/dav/dav_subspace/bpcg/ppcg";
         read_sync_int(input.pw_diag_nmax);
         this->add_item(item);
     }