Repository navigation
feat(hsolver): add focused PPCG solver support - #8264
cheerly-pku wants to merge 209 commits into
Conversation
Consider the previous contributions made by classmates, I'm only capable to make small difference without disrupting the entire program ---- like such a small "static".
…w_Small-Changes 2025PKUCourseHW5: Case: 1 - Change rank_seed_offset to static const
…ent) Add PPCG iterative diagonalization with two strategies: - CONJUGATE_GRADIENT: band-by-band Polak-Ribiere CG (verified working) - BLOCK_SUBSPACE: block subspace diagonalization Includes potrf retry fix: save/restore original matrix before applying diagonal shift, preventing accumulated shifts from corrupting the matrix. Test: 1D particle-in-a-box (n_dim=10), CG strategy matches exact eigenvalues with error 4.3e-12. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…ormalization Three fixes for numerical stability: 1. potrf: save/restore original matrix before diagonal shift retries, preventing accumulated shifts from corrupting the Cholesky factor. 2. sygvd/syevd: skip workspace query (lwork=-1) and allocate directly. The LAPACK replacement ignores workspace queries, causing the second call to operate on already-transformed data, corrupting eigenvalues. 3. Block subspace: add chol_qr + hpsi/spi recomputation after update_one_block and every rayleigh_ritz, keeping wavefunctions S-orthonormal and preventing numerical drift of H|psi> and S|psi>. Results (1D particle-in-a-box, S=I): - CG nband=1: error 4.3e-12 (unchanged, already working) - BLOCK_SUBSPACE nband=1: no longer NaN, converges (to wrong eigenvalue due to algorithmic limitation with S=I)
… RR steps 1. solve_small_generalized: save/restore M matrix before retry with shifts (prevents accumulation of shifts on sygvd-corrupted M) 2. BLOCK_SUBSPACE: add Krylov fallback for near-collinear p/w vectors When p is nearly parallel to w (cos^2 > 0.99), replace p with H·w to keep the 3-vector subspace [psi, w, p] full rank. This fixes NaN eigenvalues for nband>1 with S=I. 3. LAPACK: use standard workspace query (lwork=-1) pattern for syevd/sygvd More robust with real LAPACK implementations. 4. CG: add periodic Rayleigh-Ritz subspace rotation every rr_step iterations Corrects band ordering and eigenvalue estimates after band-by-band line minimization. Resets PR state after rotation.
The gamma_dot function returns only the real part of inner products, which is correct for Hermitian forms like <psi|H|psi> but wrong for projection coefficients where the imaginary part matters. In orth_gradient and project_against, the projection coefficient <psi_i | v> must use the full complex inner product to correctly remove the overlap. Using only Re(<psi_i | v>) leaves an imaginary component that corrupts the search direction, causing excited-state bands to converge to wrong eigenvalues. This fixes the CG strategy bands 1 and 2 converging to the highest eigenvalue (3.919) instead of the first excited states.
…PACE The chol_qr_active call after update_one_block re-orthonormalized psi but left the p vector in the old basis, creating an inconsistency. The p vector is constructed in update_one_block using the same subspace rotation as psi, so they start consistent. Adding chol_qr_active before the p vector is updated breaks this consistency.
This change was unrelated to the PPCG integration and should not have been included.
H and S are real symmetric operators whose eigenvectors are real. The previous complex random initialization produced complex off-diagonal elements in the H-gram matrix (max |Im| ~ 0.5 for nband=3), causing Re(<psi_i|H|psi_j>) != <psi_i|H|psi_j>. The gamma_dot function only returns the real part, so all subspace Gram matrices (built via gram()) computed wrong off-diagonals, leading to incorrect eigenvalues from sygvd. With real-only psi all inner products are real and gamma_dot is exact, so both BLOCK_SUBSPACE and CONJUGATE_GRADIENT strategies should now converge to the correct eigenvalues.
…tioning The 3-block subspace method builds a generalized eigenvalue problem with basis V = [psi, w, p] where p is constructed from the previous subspace eigenvectors (p_new += w_l * cw in update_one_block). This makes p a linear combination of the w vectors, causing the [w, p] block of the S-gram matrix M to become nearly rank-deficient. With nband=3 and sbsize=4 the 9x9 M matrix has condition number large enough that dsygvd produces negative eigenvalues for the positive-definite problem (observed: -0.26 at iter=2), and the eigenvalues diverge exponentially thereafter. Setting use_p=false reduces the subspace to [psi, w] (2-block), which is a preconditioned Davidson-like method. It converges robustly: the BLOCK_SUBSPACE test now passes in 57 ms with all 3 eigenvalues within 1e-8 of the exact values. The 3-block code path is preserved for future re-enablement once a more robust p-vector construction is implemented.
With rr_step=4, the non-RR iterations use Cholesky orthonormalization
which mixes bands through the upper-triangular U^{-1}, causing high-energy
bands to contaminate low-energy ones. This drives CG eigenvalues to the
spectrum maximum [3.31, 3.68, 3.92] instead of the correct lowest values
[0.081, 0.317, 0.690].
Using rr_step=1 forces Rayleigh-Ritz every iteration, which correctly
diagonalizes the subspace and preserves band ordering.
The orth_cholesky call before rayleigh_ritz mixes bands through the
upper-triangular U^{-1} factor, contaminating low-energy bands with
high-energy components. This drives CG eigenvalues to the spectrum
maximum instead of the minimum.
rayleigh_ritz solves the generalized eigenvalue problem K v = λ M v
via dsygvd, which correctly handles non-S-orthogonal bases. The
orth_cholesky is not needed and is actively harmful.
This makes the CG RR path consistent with BLOCK_SUBSPACE, which
calls rayleigh_ritz without prior orth_cholesky.
BLOCK_SUBSPACE starts with rayleigh_ritz (line 1085) which finds correct eigenvalues and rotates psi before the iteration loop. CG was using diagonal Rayleigh quotients instead — these are poor approximations for random initial guesses, producing wrong gradients that drive the band-by-band line_minimize toward high-energy eigenstates. With rr_step=1 (every-iteration RR), the CG loop itself is now correct, but without an initial RR the first line_minimize step already pushes psi in the wrong direction, and subsequent RR steps cannot fully recover.
Backup preserved at diago_ppcg_test.cpp.bak
The linear approximation α = -C/B drops the α² term from the Rayleigh quotient derivative dR/dα = 0. This picks one of the two stationary points (minimum or maximum) arbitrarily. For bands far from convergence it can select the MAXIMUM, driving ψ toward high-energy states instead of the desired lowest eigenvalues. Solve the full quadratic Aα² + Bα + C = 0, evaluate R(α) for both roots (and the linear guess), and pick the one with the lowest R. Also restore the CG unit test (rr_step=1, initial rayleigh_ritz).
…e use_p
Three changes to make both PPCG strategies correctly converge with rr_step=4:
1. CG non-RR path: After orth_cholesky, solve the nband x nband subspace
generalized eigenvalue problem instead of using diagonal Rayleigh quotients.
The upper-triangular U^{-1} from Cholesky mixes high-energy components into
low-energy bands, making diagonal RQs overestimate the eigenvalues. The
subspace solve gives correct Ritz values without rotating the states,
preserving Polak-Ribiere conjugate-direction accumulators.
2. BLOCK_SUBSPACE: Re-enable use_p=true (3-block [psi, w, p] subspace).
The Krylov fallback (replace p with H·w when p ~ w) was already in place
but dead because use_p was hardcoded to false. Now it activates on the
first iteration (p is zero-initialized) and whenever p becomes collinear
with w after update_one_block.
3. CG test: Change rr_step from 1 back to 4 so the non-RR Cholesky path
is exercised, validating the true Polak-Ribiere CG mechanism.
The 3-block [psi, w, p] subspace generalized eigenproblem becomes ill-conditioned when residuals are small (near convergence). The [w, p] Gram block shrinks, the M matrix approaches singularity, and dsygvd produces garbage eigenvectors that drive eigenvalues to catastrophic values (e.g., -137775 instead of 0.081). The p-bad H·w Krylov fallback fixes p~w collinearity but does not address the small-residual ill-conditioning, which is fundamental to the 3-block construction. Keep use_p=false for robust convergence.
The [w,p] block of the Gram matrix M shrinks as residuals converge, making M nearly singular and causing sygvd to produce garbage eigenvectors. Scaling w and p to unit S-norm keeps M well-conditioned (diagonal ~1) without changing the subspace — Ritz values are identical and Ritz vector coefficients cancel in update_one_block. This enables the full 3-block [psi,w,p] subspace (use_p=true) by addressing the fundamental ill-conditioning that the p-bad Krylov fallback alone could not handle.
The Krylov fallback (replace p with Hw when p~w) was flawed: when w is approximately an eigenvector (Hw ≈ λw), the replacement does not fix collinearity. After S-norm scaling, p ≈ w still, M_wp ≈ [1,1;1,1] is rank-1, and dsygvd fails. Instead, simply skip p for this iteration (use_p_now=false). update_one_block still produces a valid p for the next iteration from the w Ritz-vector contribution.
Add tests for: - 2x2 matrix (smallest non-trivial case) - Degenerate eigenvalues (H = I + J, multiplicity-3 degeneracy) - Larger 20x20 tridiagonal with 5 bands - Dense 8x8 matrix via Givens rotations (addresses full-matrix coverage) All use CONJUGATE_GRADIENT strategy which has sygvd fallback. BLOCK_SUBSPACE tests deferred due to dsygvd instability with some LAPACK builds. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Covers diagonal, tridiagonal, dense, pentadiagonal, degenerate, Neumann, S≠I, gamma_g0, single-band, all-band, many-band, bad preconditioner, tight threshold, scaled, gapped spectrum, rr_step=1, 1x1, and eigenvector quality checks. Adds QuickBenchmark (CI-friendly) and DISABLED_FullBenchmark.
Combine b5d5fdb and ad39c40 without rewriting either history. Retain upstream deletion of generated INPUT docs; updated C++ Input_Item metadata remains the documentation source. Verified build: abacus_pw_para plus focused hsolver targets. Focused CTest: 6/6, including 50 double, 4 float, 9 parameter and 6 PW tests and MPI np2/3/6 and two pools. Native and expanded integration issues remain under investigation; this local merge is not a release or full-validation claim.
Exercise kept counts 2, 3 and 7 in a rank-7 complex projected problem. Verify generalized residuals and overlap orthogonality for every valid output, including discarded directions consumed by low-root recovery. No production behavior or INPUT changes; no documentation regeneration is required. Verified freshly built MODULE_HSOLVER_ppcg: 51/51 pass, one benchmark remains disabled. No runtime solver change is included.
Resolve the stochastic HSolver test CMake conflict by keeping develop's formatting and the PPCG source registration. Preserve the original PR history and numerical implementation.
Retain PPCG tests in the configure-time fixture setup, register ppcg in the unified solver domain, and pass the MPI test executable explicitly. Add INPUT round-trip, metadata and basis-compatibility checks. Seven focused CTests pass; full build and native-case validation are recorded separately.
Reject the unsupported sdft/ppcg combination during unified INPUT validation, including zero deterministic bands, and cover both rejection and supported CG controls. Document this restriction and correct the guide to distinguish full Ritz cadence from independent H/S image refresh.
Guard the optional mallinfo2 API for supported glibc versions and report an unavailable metric elsewhere. Label post-diag retained heap separately from peak memory, and describe each solver's stopping and restart behavior accurately.
Share checked int-index and LAPACK workspace limits across CPU and GPU entry paths, reject rank-local invalid shapes collectively, and avoid overflow in block-count arithmetic. Add bounded serial and asymmetric MPI regressions. Strengthen the generalized fixture to check lowest analytic roots, fresh residuals and S-orthogonality, and split grouped test declarations.
Replace the all-one smoke thresholds with measured per-property limits: 1e-6 eV energy, 0.02 eV/Angstrom force aggregate and 0.5 kbar stress aggregate. Preserve existing inputs and references. Record RTX3090 CPU/GPU paired agreement and observed reference offsets from repeated two-rank and four-rank runs.
…curacy Preserve accurate iterative eigensolver updates and basis recovery, add the optional pw_diag_thr_max bound across SCF and fixed-cell relaxation cycles, and cover numerical and INPUT boundaries with focused tests. Retain the validated noncollinear density/Libxc fixes and full-precision PW output. Compare complete wavefunction and Wannier data with phase/gauge-aware checks. Record independent reference calibration and explicit accuracy settings in the affected case READMEs without relaxing acceptance thresholds. Local validation of this snapshot (GNU 12.3, OpenMPI 5.0.10, OMP=1): - Autotest.sh -n4 -o1 -j4 over 130 baseline and 130 PPCG-side calculations: all runtime exits zero; 121 substitutions: 64 pass, 57 warning, 0 fatal. - Nine unchanged/native controls pass. Two existing incomplete baseline reference records and remaining SCF/inner warnings are preserved. - Final relax_dftu_spin4 cross-checks: CG, Davidson-subspace and PPCG; MPI1/MPI4 reference checks pass at original tolerances. - HSolver/INPUT cap tests: 33/33. Final binary CLI controls: 7/7. - C++11 checks for cap changes and staged governance pass. INPUT help is generated from updated Input_Item registrations; generated documentation is not committed. Case READMEs describe reference provenance. The local report, thesis, presentation, and experiment archives are excluded.
Integrate origin/ppcg-v2-reopen without rewriting either history. Preserve its workspace guards, stochastic INPUT rejection, GPU calibration, benchmark portability, and empty-rank MPI support alongside the locally validated solver and reference changes. Retain both Davidson stability and ownership test harnesses. Apply empty-rank BLAS dimensions and zero projection contributions to the new orthogonalization paths. Keep complete shared-gauge Wannier comparisons and both INPUT descriptions. The merged tree is being rebuilt and retested before publication.
Set nscf_hse_same pw_diag_thr to 1e-12 with the default 50-iteration budget. The former threshold leaves a 1.05e-8 eV CG band error after the numerical update and crosses the stored 1e-8 eV comparison boundary. CG, Davidson-subspace and PPCG reproduce the existing band reference with the tighter input and pass the original Autotest energy/band checks at four MPI ranks. Four tighter controls independently agree within 2.4e-13 eV. Preserve references, comparison tolerances and the PPCG residual warning at the default budget. No production source changes.
The restart_spd fixture stops on eigenvalue changes but independently checks a residual. Its 1e-14 request stopped one expansion too early with OpenBLAS 0.3.20, yielding 8.14e-7 against the unchanged 2e-7 residual limit. Request 1e-15 for this fixture without changing production stopping rules, the 200-iteration budget, or any acceptance threshold. The exact test passes all 18 stability scenarios on OpenBLAS 0.3.20 and 0.3.33 (36 runs), plus all 47 focused CTests. The original CI failure is reproduced with 0.3.20; the corrected residual is 2.03e-11.
The shared Libxc uniform-gradient correction also changes this LCAO case. Its old charge cube failed CPU and GPU checks at the Fermi-energy header. Tighten SCF to 1e-10 and print 12-digit cube fields; generate references independently at 1e-12 while retaining all acceptance limits and the original 100-step budget. Validate all 3375 grid points and scalar checks with Libxc 5.1.7 and 7.0.0, explicit CPU/GPU runs and one/two MPI ranks. Eight original Autotest runs pass, plus a Libxc 5.1.7 CPU four-rank control. The 65-case GPU LCAO multik suite passes with the CI Libxc version.
Set pw_diag_thr=1e-12 to keep CG bands away from the original 1e-8 eV printed comparison boundary. Set pw_diag_nmax=100 so PPCG substitutions converge at the tighter residual request without the warnings seen at 50. Retain every reference and comparison tolerance; the 200-step control adds no improvement. The exact fixture passes all eight CPU/GPU, CG/PPCG, MPI 1/4 runs and 32 original property checks with zero inner residual warnings. Maximum band differences against the tighter Davidson-subspace control are 1.954e-11 eV for CG and 1.901e-13 eV for PPCG.
Tighten SCF and diagonalization settings for eight PW fixtures. Keep the original references for the three lattice cases; calibrate the five remaining scalar references with independent high-accuracy Davidson-subspace results cross-checked against CG and PPCG. Keep all acceptance tolerances, physical inputs, and the comparison pipeline unchanged. Fresh paired runs of all 135 active cases reduce PPCG reference warnings from 60 to 52, with 83 passes and zero fatal failures. All 270 calculations exit normally and no prior reference pass regresses. Forty focused CPU/GPU runs pass 220 property checks with no inner residual warning. Production source is unchanged. Full-sweep residual warnings and the two intentionally short SCF cases remain separately recorded.
The three LR force fixtures requested lr_thr=1e-2 while checking signed forces and their absolute-value sum to 1e-6. The shared Davidson changes alter this loose stopping path. An independent build of upstream c82d87d reproduces the old references with loose Davidson, while upstream and current LAPACK give identical printed signed force components. Request lr_thr=1e-16 and calibrate the scalar/signed references from upstream dense LAPACK. Retain Davidson, the existing retry policy, Z-CG coverage, ground-state inputs and all acceptance tolerances. Twelve final controls across OpenBLAS 0.3.20/0.3.33 and MPI1/4 pass; signed-force and sum differences stay below 2e-7. The original 27-case 08_RI suite passes 143 property checks, with the unrelated existing HF-ULR SCF nonconvergence recorded separately. Subsequent 09_DeePKS and 10_others suites also pass locally.
Merge 8a6df3b while retaining the PPCG, Davidson stability and CG accuracy checks. Follow the upstream test-to-unittests migration, keep the upstream CG tests and the existing accuracy harness in separate targets, and use the relocated linear-algebra tests. Preserve both sides of the PW case list. Validation: all selected targets rebuilt; 60 enabled CTests pass with each of OpenBLAS 0.3.20 and 0.3.33. The comparison benchmark remains disabled. Full CPU integration verification follows before publication.
The Read the Docs build for f151acf hit its 900-second limit while compiling ABACUS, before generating parameters or running Sphinx. The default untyped CMake build adds -O3 -g even though this executable is used only to export INPUT metadata. Set Debug with -O0 -g0 and disable native optimization for this documentation-only build. Keep assertions enabled and preserve the regular build defaults. Validation: the complete serial/no-LCAO install finished locally in 129.315 seconds with four build jobs. All 605 documented parameters, generated Markdown, and 73 HTML files match the merged Release build byte for byte. Strict HTML and LaTeX generation pass with warnings treated as errors. Remote RTD timing requires the next build.
Keep the PPCG implementation, explicit Rayleigh-Ritz cadence, focused unit tests, and dedicated integration fixtures on the current develop base. Remove unrelated solver, XC, charge-mixing, output, comparator, and reference calibrations from the old PR scope.
# Conflicts: # source/source_hsolver/CMakeLists.txt # source/source_hsolver/unittests/CMakeLists.txt # source/source_pw/module_stodft/test/CMakeLists.txt
Use the existing parallel reduction wrapper for trace writer selection and keep the production PPCG path on the conservative full Ritz cadence. Remove the unvalidated user-facing cadence parameter from the focused PR.
ABACUS-PR-BOT
left a comment
There was a problem hiding this comment.
PPCG's normal-return failure status needs to reach the PW caller, and concurrent test invocations need separate trace files. Static review at b37124c; builds, runtime behavior and reference values were not independently verified.
| this->pw_diag_ndim, | ||
| 1, | ||
| cpu_path); | ||
| DiagoIterAssist<T, Device>::avg_iter += ppcg_result.average_iterations; |
There was a problem hiding this comment.
When PPCG exhausts its iteration budget, diag() returns normally with converged=false, but this branch keeps only average_iterations. The returned eigenpairs then enter occupations/output without a convergence diagnostic, including NSCF. Please handle the returned status here and add a PW-dispatch regression with an exhausted budget.
| // process environments. With multiple pools, every rank requests the same | ||
| // path to verify that only the MPI world root writes it. | ||
| const char* trace_environment = "ABACUS_PPCG_RESIDUAL_TRACE"; | ||
| const std::string trace_path = testing::TempDir() + "ppcg_residual_trace_parallel_test.csv"; |
There was a problem hiding this comment.
Both the addtest entry and the shell-wrapper entry launch this binary, so parallel CTest can run two world roots against this same temporary filename. Each run removes, truncates and later deletes the other run's trace. Please use a unique per-invocation path or a shared CTest resource lock.
| * cg: The conjugate-gradient (CG) method. | ||
| * 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. | ||
| * dav_subspace: A block Davidson method. The CPU implementation orthogonalizes trial vectors in the overlap metric and filters linearly dependent expansion directions. `pw_diag_ndim` can be set to 2 for this method. |
There was a problem hiding this comment.
Non-blocking: The unchanged dav_subspace path only preconditions and normalizes each trial vector before appending the block; it does not perform the metric orthogonalization or dependent-direction filtering described here. Please keep this help text aligned with the implementation.
| } | ||
| } | ||
|
|
||
| // Bad preconditioner: use 1.0 instead of 2.0 |
There was a problem hiding this comment.
Non-blocking: Replacing 2I with I only rescales the residual, and scale_to_unit_snorm removes that factor before the subspace solve. Use a fixed nonuniform positive diagonal if this fixture is meant to exercise a poor preconditioner, and check convergence plus a fresh residual.
| @@ -0,0 +1,6 @@ | |||
| GPU smoke case for the PPCG host/device bridge. | |||
|
|
|||
| It uses the same GaAs input as the CPU PPCG case and checks the standard energy, | |||
There was a problem hiding this comment.
Non-blocking: This shares the GaAs structure, but its 2x2x2 k mesh differs from the CPU PPCG case's 1x1x2 mesh. Please say 'same structure' and identify the mesh so the different energy/force references are not mistaken for a CPU/GPU discrepancy.
| copy_block(subspace.sw_l, l, subspace.sbasis); | ||
| if (nblk >= 3) | ||
| { | ||
| copy_block(subspace.p_l, 2 * l, subspace.basis); |
There was a problem hiding this comment.
Non-blocking: Please name the repeated 2*l block offset before these calls, following AGENTS.md's rule for computed arguments. Apply the same cleanup to the new reduction counts and endpoint ternaries in the tests.
| { | ||
| // A single writer avoids file corruption when several MPI pools solve | ||
| // different k-points concurrently. | ||
| residual_trace.open(residual_trace_path); |
There was a problem hiding this comment.
Non-blocking: Please record the resolved full CSV path in the running log when this trace file is created or overwritten, through an explicit reporting interface. This keeps the new output consistent with the project's file-I/O convention.
Scope
This is a fresh PPCG-only successor to the closed #8042, started from
036b6414dand synchronized with currentdevelopat597deb7bf. It contains:The old PR's unrelated Libxc and charge-mixing changes, CG/Davidson changes, WFC/Wannier comparators, generic output changes, LR fixtures, and ordinary 01_PW reference calibrations are intentionally absent. The SCAN cutoff fixture 820 is also excluded because its reference was calibrated against the withdrawn Libxc behavior and needs an independent current-develop recalibration.
Verification
Fresh local build on
b37124cb7(started from036b6414d, synchronized with currentdevelop):Result: 5/5 tests passed. The MPI executable also passes with
--gtest_color=no, which previously exposed a positional-argument division-by-zero path. The input-domain and solver-label tests pass (2/2).The current
abacus_pw_parabinary ran 817, 818, and 819 with two MPI processes, one OpenMP/BLAS thread, and the existing acceptance limits; all three passed all 7 collected reference checks. CUDA execution was not available on this host, so the GPU fixture remains a CI/GPU-host check rather than a local numerical claim.The branch also fixes the previously identified PPCG workspace bound for root recovery (
[X,Q]can reach3*nband), removes test-directory trace files, replaces the newM_PImacros, centralizes MPI world-root detection in the parallel wrapper, and makes the parallel test ignore GoogleTest flags when parsing its optional pool count. All required GitHub checks now pass; the PR is mergeable and is waiting for maintainer review.Review notes
The implementation and tests are intentionally kept in one focused PR because the production dispatch, test registration, and correctness fixtures must be reviewed together. The excluded validation/calibration work can be proposed separately with its own current-develop evidence.