Repository navigation
feat(hsolver): add focused PPCG solver support - #8264
Open
cheerly-pku wants to merge 215 commits into
Open
cheerly-pku wants to merge 215 commits into
cheerly-pku wants to merge 215 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.
…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
ABACUS-PR-BOT
left a comment
Collaborator
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.
Propagate PPCG convergence diagnostics to the PW caller, stop uncertified NSCF results, and test both SCF and NSCF dispatch. Isolate concurrent MPI trace files and report resolved paths through an explicit running-log stream. Correct the poor-preconditioner regression and scoped help/fixture descriptions.
Author
|
Fresh verification report for the current head Key results:
All remote checks for this head now pass, including CUDA build/test, GNU/Intel builds, Test, governance, pre-commit, PyTest, abacuslite, and documentation. The seven earlier inline bot findings each have a separate reply and are resolved; no new unanswered inline finding was present during this update. |
ABACUS-PR-BOT
left a comment
ABACUS-PR-BOT
left a comment
Collaborator
There was a problem hiding this comment.
Two issues need attention: PPCG executes diagnostic work with tracing disabled, and fresh PW NSCF runs ignore the configured iteration budget. One non-blocking benchmark wording correction is also included. Static source review only; no builds, numerical runs or reported CI results were independently verified.
mohanchen
reviewed
Nov 5, 2026
Collaborator
There was a problem hiding this comment.
it's not necessary to have such a large test file
This branch has not been deployed
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Adds PPCG to PW Hamiltonian diagonalization, with residual-based convergence checks, CPU/MPI tests, and native CPU/GPU fixtures. This is the focused successor to closed #8042; the unrelated Libxc, charge-mixing, CG/Davidson, output-comparator, LR, and ordinary PW fixture changes remain outside this PR. SCAN fixture 820 remains excluded pending independent calibration.
Validated head:
9de3fb7072a224c5938a982c6f97eeb597703178, based on develop at597deb7bf. All numbers below were measured again on 2026-11-05 using this head. The executable's startup banner was checked against the commit after refreshing CMake's configure-time version information.Review changes
is_open() && good()before enabling trace-only work. A regression compares diagnostic-helper call counts with tracing unset, file opening unsuccessful, and tracing enabled. It reproduced the defect before the fix.pw_diag_nmaxthrough the real PPCG NSCF parameter setup. Before the fix, the handoff test requested 77 but received 30. Controls preserve the other solvers' NSCF setup and all solvers' SCF setup.avg_itermeans average subspace iterations per band.Fresh correctness checks
GNU 12.3, Release, Intel Xeon Gold 6242, OpenMPI 4.1.2, OpenMP enabled, LCAO/CUDA disabled, ELPA and OpenBLAS. The existing
build-reviewcache was reused. Runtime variables wereOMP_NUM_THREADS=1,OPENBLAS_NUM_THREADS=1,MKL_NUM_THREADS=1; residual tracing was unset except in trace-specific checks.tests/01_PW, wasbash ../integrate/Autotest.sh -a <absolute-path-to-abacus_pw_para> -n 2 -o 1 -j 1 -f CASES_CPU.txt; that selection contained only these three fixtures.pw_diag_nmax=1exited as expected at iteration 1;pw_diag_nmax=77completed successfully, with the last k-point trace reaching iteration 36. This exercises a budget beyond the previous static value of 30.abacus_pw_para --versionreportedv3.11.0-beta12;-h pw_diag_nmaxand--check-inputin a valid case directory succeeded.git diff --checkand the diff-oriented governance check passed. Existing OpenMPI/hwloc startup warnings remain in the logs.The native final energies (817/818/819) were
-4862.3309704994053391,-227.6796325580887697,-1882.1779472344287569eV.Fresh efficiency measurements
Each case/solver used a clean copy of the same committed fixture; only
ks_solverdiffered within each group. Runs were sequential, trace disabled, usingmpirun --bind-to core --map-by core -np 2 abacus_pw_para. There was one warm-up and five measured runs per solver, with rotated solver order. Wall time surrounds the MPI subprocess;solve_psikis the ABACUS solver timer.solve_psikmedian [min, max] sOf the 60 measured runs, 55 converged. The five excluded samples are all 819/CG, which reached the fixture's 50-step SCF limit. PPCG is slower than the Davidson variants on these small regression cases. Large-system performance remains to be measured.
A separate 12-run correctness comparison used the same overrides for every solver:
scf_thr_type=2,scf_thr=1e-10,scf_nmax=150,pw_diag_thr=1e-10,pw_diag_nmax=200. Davidson, Davidson-subspace and PPCG converged in all three cases. CG reached the 150-step SCF limit in all three and is excluded from the numerical comparison. Across the six converged Davidson/PPCG comparisons, the largest differences were 2.728E-12 eV in total energy, 2.060E-8 eV/Å in a force component, and 5.424E-7 kbar in a stress component. These runs are used only for numerical agreement.The local artifact directory is
/mnt/data/Cheerly/ppcg-bot-followup-20261105/:independent-validation.md/.jsonaudit the raw logs, 84 INPUTs, 336 other tracked files and pseudopotentials; protocols record the command, threads, commit and executable hash. The fresh native and NSCF records are innative/andnscf-cli/. Executable SHA-256:27666bef7165a9086d6f13d063a39d25ab80d40959bfaa731574f8430e6fc8ac. The earlier batch with a stale configure-time version banner is archived underbefore-version-refresh/and supplies none of the table's numbers.Remote checks and documentation
The new push starts a separate CI run for
9de3fb707; its status is shown in this PR's checks. The local checks above cover the CPU build and runtime paths. CUDA and the other toolchains rely on their corresponding CI jobs.No docs needed for these follow-up fixes: they restore the documented iteration-budget and opt-in trace behavior without changing the INPUT schema, defaults or acceptance tolerances. The PR retains its PPCG help metadata and solver documentation. No new global dependency or default interface argument is introduced by these follow-ups.