Skip to content

Refactor hsolver: DiagoIterAssist takes H|psi>/S|psi> functors - #7951

Open
Critsium-xy wants to merge 3 commits into
deepmodeling:developfrom
Critsium-xy:refactor/hsolver-iter-assist-closures
Open

Refactor hsolver: DiagoIterAssist takes H|psi>/S|psi> functors#7951
Critsium-xy wants to merge 3 commits into
deepmodeling:developfrom
Critsium-xy:refactor/hsolver-iter-assist-closures

Conversation

@Critsium-xy

Copy link
Copy Markdown
Collaborator

Reminder

  • I have read AGENTS.md and docs/developers_guide/agent_governance.md.
  • I have linked an issue or explained why this PR does not need one.
  • I have added adequate unit tests and/or case tests, or explained why not.
  • I have listed the exact verification commands run and their results.
  • I have described user-visible behavior changes, including INPUT parameter changes.
  • I have explained core-module impact for ESolver, HSolver, ElecState, Hamilt, Operator, Psi, or other source/ changes.
  • I have requested any needed governance exception below.

Linked Issue

No issue. Third step of making source_hsolver a self-contained numerical
module, after #7920 and #7948. Behaviour-neutral.

There is a standing request for exactly this change in the tree:
source_lcao/module_lr/hsolver_lrtd.hpp:122 reads

diagH_subspace needs refactor: replace Hamilt* with hpsi_func — or I
cannot use is_subspace=true as my HamiltLR does not inherit Hamilt.

This PR does the Hamilt*hpsi_func part. Actually enabling the LR path is
left to whoever owns that module, since HamiltLR needs functors of its own.

Unit Tests and/or Case Tests for my changes

Commands run (Linux, gcc, cmake 3.31 + ninja, 15 cores). Base commit
9cbcc0b55 and this branch were built and tested in the same build tree.
ENABLE_ELPA=ON is needed for full coverage of the module.

cmake -B build -G Ninja -DBUILD_TESTING=ON -DENABLE_LCAO=ON -DENABLE_MPI=ON \
      -DENABLE_OPENMP=ON -DENABLE_ELPA=ON
cmake --build build -j15 -- -k 0
OMP_NUM_THREADS=1 ctest --test-dir build -E '^(0[1-9]|1[0-9])_'
python tools/03_code_analysis/agent_governance_check.py --staged

Result summary

check base 9cbcc0b55 this branch
full build, ELPA on 3328/3328, 0 errors 3328/3328, 0 errors
unit tests 91% passed, 29 failed of 339 91% passed, 29 failed of 339
governance check 0 blockers

Failing sets are identical test-for-test (diff of the sorted lists is
empty); they are this sandbox's pre-existing failures.

The unit suite does not actually exercise these three functions, so that
result on its own is weak evidence. The CG unit tests pass
DiagoIterAssist<T>::need_subspace (which is false) into DiagoCG, so their
subspace_func — the one that calls diag_subspace — is never invoked. I
therefore drove each of the three converted entry points through integration
cases instead, and checked the total energy against the committed result.ref:

case reaches nspin branch etot (eV) result.ref Δ
01_PW/022_PW_CG diag_subspace ×14 1 -198.2238296277164 -198.2238296207179 7e-9
01_PW/036_PW_AF diag_subspace ×60 2 -5866.197297502092 -5866.197297502891 8e-10
10_others/01_NP_KP_sp diag_subspace_init ×64 1 -723.7838346685444 -723.78383467 2e-9
01_PW/scf_deltaspin4 cal_hs_subspace 4 -6369.1982731684611 -6369.198273168575 1e-10
17_DS_DFTU/14_PW_DS_S4_XYZ cal_hs_subspace 4 -6369.198950976948 -6369.19895097706 1e-10

Every case reproduces its reference to the precision the reference is stored at.
The call counts are read out of each run's own timing table, so the first three
rows are direct evidence that the converted code ran. cal_hs_subspace has no
timer of its own, but it is what fills the h_k/s_k subspace matrices that
diag_responce (×148 and ×156 in those two runs) then diagonalizes — a broken
substitution there could not land within 1e-10 of the reference.

The last two rows matter most: they are nspin 4, i.e. npol == 2, which is
the case the CPU branch's standing comment warns about and the one the unit
tests never reach.

Checks not run, with reason

  • GPU (base_device::DEVICE_GPU) branches of diag_subspace_init and of
    cal_mw_from_lambda.cpp: no GPU runtime on the machine available to me. They
    are compiled by CI's CUDA job, and the substitution there is the same
    forwarding shim as the CPU branches, with the dimensions passed through
    unchanged (which is the reason the functors take them explicitly).
  • 01_PW/BUG_SCF_DSPIN would have been an nspin 4 + cg case, but it exits
    non-zero on the base commit as well — it is a known-broken case, as its name
    says.

What's changed?

DiagoIterAssist was the last algorithm class in source_hsolver still
holding a hamilt::Hamilt<T, Device>*. It used it for exactly two things,
pHamilt->ops->hPsi(...) and pHamilt->sPsi(...), so its three
Hamiltonian-taking entry points — diag_subspace, diag_subspace_init and
cal_hs_subspace — now take two functors instead:

using HPsiFunc = std::function<
    void(T* psi_in, T* hpsi_out, const int ld_psi, const int current_nbasis, const int nvec)>;
using SPsiFunc = std::function<
    void(const T* psi_in, T* spsi_out, const int nrow, const int npw, const int nbands)>;

DiagoCG, DiagoDavid, DiagoDavSubspace and DiagoBPCG have taken closures
for a while, and HSolverPW::hamiltSolvePsiK already builds them; this brings
the last one into line.

Why these functors carry their dimensions instead of capturing them. The
existing 4-argument HPsiFunc in DiagoCG and friends works because those
solvers drive one fixed wavefunction layout, so the caller can capture a single
cur_nbasis. DiagoIterAssist drives three:

psi handed to hPsi sPsi(nrow, npw, …)
diag_subspace, cal_hs_subspace the caller's own psi (dmax, dmin, nstart)
diag_subspace_init, GPU branch 1-band staging buffer (dmin, dmin, 1)
diag_subspace_init, CPU branch nstart-band staging buffer (psi_nc, psi_nc, nstart)

sPsi's npw genuinely differs between the two branches of one function, so a
captured constant cannot reproduce it, and the CPU branch carries a standing
comment that current_nbasis must be npw without npol or Nonlocal::act's
gemm K stops matching vkb's row count. Passing the dimensions through keeps
every caller's functor a plain forwarder and makes the substitution checkable by
reading it, which matters because the SOC and GPU paths are not covered by the
unit tests available to me.

Why the substitution is equivalent. Where the old code handed
Operator::hPsi a psi::Psi plus a Range, the functor rebuilds a one-k
wrapper over the same pointer. That is the same thing because:

  • psi.get_pointer() returns psi_current, which fix_k sets to
    psi + ik * nbands * nbasis — exactly the pointer
    to_range(Range(1, current_k, 0, nstart-1)) computes;
  • Psi::get_npol() reads PARAM.inp.nspin globally rather than a per-object
    member, so the wrapper reports the same npol (and hence the same
    nbands * npol inside Operator::hPsi) as the original object;
  • get_nbasis() and get_current_nbas() are passed through explicitly, per the
    table above.

The ops == nullptr early exit in diag_subspace_init becomes an
empty-functor check. The two callers that can reach it (psi_prepare.cpp,
hsolver_lcaopw.cpp) build the functor only when ops is allocated and
otherwise leave it default-constructed, so the warning and the
copy-psi-to-evc fallback fire under exactly the same condition as before.

One honest caveat about that: psi_prepare.cpp builds a single functor and
uses it for both diag_subspace_init and diag_subspace, and diag_subspace
has never had a null-ops guard. So in the (unsupported, already-broken) state
where ops is null and that branch is taken, the failure changes shape — it
used to dereference a null pointer, and now throws std::bad_function_call.
Neither is a working state; I did not want to invent a new guard here.

before #7920 #7948 this PR
source_hsolver files including source_hamilt 16 6 5
#include "source_hamilt/..." lines there 17 6 5

The five that remain are hsolver_lcao.h, hsolver_lcaopw.{h,cpp} and
hsolver_pw.{h,cpp} — the HSolver* façades, which legitimately own a
Hamilt* because they run the k-loop. Moving those out of source_hsolver
into the physics modules they belong to is the next and last step, not this one.

Not in this PR: DiagoIterAssist's five mutable statics (PW_DIAG_THR,
PW_DIAG_NMAX, avg_iter, need_subspace, SCF_ITER). They are worth
removing — they are the AGENTS.md rule 2 pattern — but they create no
source_hamilt dependency, and only one algorithm actually reads one of them
(diago_bpcg.cpp:302 reads SCF_ITER); the rest is parameter plumbing between
setup_diago_params_* and the HSolverPW constructor, which already has
matching instance members. That is a separate, self-contained change and mixing
it in here would double the diff and the risk.

Governance Notes

  • INPUT/docs changes: none. No Input_Item, no user-visible behavior, no
    output format is touched.
  • Core module impact: HSolver, and Hamilt as seen from HSolver. No
    hamilt:: class is modified — DiagoIterAssist simply stops asking for one.
    Psi and Operator are untouched. Callers in source_psi and
    source_lcao/module_deltaspin gain a forwarding functor each; the values that
    reach Operator::hPsi and Hamilt::sPsi are unchanged.
  • Exceptions requested: none.
  • No default arguments were added to keep old call sites compiling; all six
    call sites are updated in this PR (AGENTS.md rule 5). The two pre-existing
    defaults on diag_subspace (n_band = 0, is_S_orthogonal = false) are
    kept as they were.
  • The forwarding functors are written out at each of the six call sites rather
    than behind a shared helper. A helper would have to live in source_hamilt
    to avoid re-introducing the dependency this PR removes; that seemed like the
    wrong thing to add here, and the natural moment to reconsider is after the
    façades move.

🤖 Generated with Claude Code

DiagoIterAssist was the last algorithm class in source_hsolver holding a
hamilt::Hamilt<T, Device>*. It used it for exactly two things --
`pHamilt->ops->hPsi(...)` and `pHamilt->sPsi(...)` -- so diag_subspace,
diag_subspace_init and cal_hs_subspace now take two functors instead:

    HPsiFunc = void(T* psi_in, T* hpsi_out, int ld_psi, int current_nbasis, int nvec)
    SPsiFunc = void(const T* psi_in, T* spsi_out, int nrow, int npw, int nbands)

DiagoCG, DiagoDavid, DiagoDavSubspace and DiagoBPCG have taken closures for
a while; this makes the last one consistent with them.

Both functors carry their dimensions explicitly rather than letting the
caller capture one fixed set, because DiagoIterAssist drives three
different wavefunction layouts: the caller's own psi in diag_subspace and
cal_hs_subspace, a one-band staging buffer in the GPU branch of
diag_subspace_init, and an nstart-band one in its CPU branch. Those
branches pass different values for sPsi's nrow/npw, and the CPU branch
carries a standing note that current_nbasis must be npw *without* npol or
Nonlocal::act's gemm goes wrong. Passing the dimensions through keeps every
caller's functor a plain forwarder and makes the substitution checkable by
inspection, which matters because the SOC and GPU paths are not covered by
the unit tests available to me.

The `pHamilt->ops == nullptr` early exit in diag_subspace_init becomes an
empty-functor check; the two callers that can hit it build the functor only
when ops is allocated.

source_hsolver -> source_hamilt includes: 6 files / 6 lines -> 5 files /
5 lines. The five that remain are the HSolver* facades, which legitimately
own a Hamilt pointer.

No behaviour change in any supported state. One failure mode changes shape:
psi_prepare.cpp builds one functor for both diag_subspace_init and
diag_subspace, and diag_subspace has never had a null-ops guard, so with a
null ops on that branch the old code dereferenced a null pointer and the new
one throws std::bad_function_call. Neither is a working state.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
@mohanchen mohanchen added the Refactor Refactor ABACUS codes label Sep 12, 2026
Comment thread source/source_hsolver/test/diago_cg_float_test.cpp Outdated
@mohanchen mohanchen added the Diago Issues related to diagonalizaiton methods label Sep 12, 2026
Critsium-xy and others added 2 commits September 12, 2026 13:34
Address review feedback on deepmodeling#7951: the H|psi>/S|psi> lambdas are now
declared as DiagoIterAssist::HPsiFunc / SPsiFunc (and DiagoCG::SubspaceFunc
for the CG subspace wrapper), and the Psi wrappers are declared with their
concrete types, so the reader no longer has to infer them.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01ShMLqMWfjYUnad69m5Lcgq
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

Diago Issues related to diagonalizaiton methods Refactor Refactor ABACUS codes

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants