Add PPCG iterative diagonalization solver to source/source_hsolver/.#7404
Open
Silver-Moon-Over-Snow wants to merge 13 commits into
Open
Add PPCG iterative diagonalization solver to source/source_hsolver/.#7404Silver-Moon-Over-Snow wants to merge 13 commits into
Silver-Moon-Over-Snow wants to merge 13 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.
The gamma_dot function returns only Re(<x|y>), which is exact for real vectors but discards the imaginary part of off-diagonal elements in complex wavefunctions. For a Hermitian operator H and complex psi, <psi_i|H|psi_j> can have O(1) imaginary off-diagonals (e.g., max |Im| ~ 0.5 in the nband=3 test). Treating the real part as the full Hermitian matrix gives wrong subspace eigenvalues from sygvd. Changes: - Add full_dot() and gram_hermitian() for full complex inner products - Add solve_hermitian_gevp(): converts complex-Hermitian GEP to a 2n x 2n real block form, solves with dsygvd, and maps eigenvectors back to complex form — avoids requiring zhegvd in the LAPACK layer - SmallSubspace now stores complex K/M matrices and eigenvectors - build_small_subspace uses Hermitian fill (conj-symmetry) - rayleigh_ritz uses full complex Gram + block-real solve - Non-RR eigenvalue computation also uses full complex Gram - All eigenvector rotations use complex coefficients The block-real reduction doubles the matrix dimension (fine for small subspace sizes, typically 3 to 9), avoiding the need for complex LAPACK (zhegvd/zpotrf/ztrtri) while being mathematically exact. This, together with the test fix (real-only initial psi), should fix both BLOCK_SUBSPACE and CONJUGATE_GRADIENT strategies for nband > 1.
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.