Add PPCG iterative diagonalization solver to source/source_hsolver/.#7404
Open
Silver-Moon-Over-Snow wants to merge 14 commits into
Open
Add PPCG iterative diagonalization solver to source/source_hsolver/.#7404Silver-Moon-Over-Snow wants to merge 14 commits into
Silver-Moon-Over-Snow wants to merge 14 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>
Silver-Moon-Over-Snow
changed the title
…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.
mohanchen
reviewed
May 31, 2026
| void Stochastic_WF<T, Device>::init_sto_orbitals(const int seed_in) | ||
| { | ||
| const unsigned int rank_seed_offset = 10000; | ||
| static const unsigned int rank_seed_offset = 10000; |
Collaborator
There was a problem hiding this comment.
What is the rationale for modifying this variable? It does not appear related to PPCG.
There was a problem hiding this comment.
I'm sorry, this is a part of the previous assignment and I forgot to handle that. I will delete it later. thank you.
…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.
Silver-Moon-Over-Snow
force-pushed
the
ppcg
branch
from
a3d481b to
3713712
Compare
…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.
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
2 participants
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.
Summary
Add PPCG iterative diagonalization solver to
source/source_hsolver/.Two strategies
CONJUGATE_GRADIENTBLOCK_SUBSPACEChanges
diago_ppcg.h— Template class headerdiago_ppcg.cpp— Full implementation with LAPACK bindingstest/diago_ppcg_test.cpp— GTest unit test (1D particle-in-a-box)CMakeLists.txt/test/CMakeLists.txt— Build integrationBug fix
potrfretry: saves original matrix before applying diagonal shift, restores before each retry — prevents accumulated shifts from corrupting the matrix.Test result (CG)
1D particle-in-a-box (n=10): error = 4.3e-12 ✅
Known limitation
BLOCK_SUBSPACEstrategy has numerical instability with S=I, unrelated to this fix.