caop thrust#68
Open
jcphill wants to merge 13 commits into
Open
Conversation
…vector<size_t> Replace std::vector<std::vector<size_t>> with sbd::det_vector<size_t> for all CAOP basis/carryover parameters and locals across davidson.h, mult.h, mkham.h, carryover.h, sbdiag.h, lanczos.h, generalop.h (friend declarations), and main.cc. Comparator lambdas updated to const auto& to avoid vector temporaries on lower_bound. Smoke-validated on DeltaAI: Hubbard 8-site E=-1.645751311064591 matches reference exactly; TPB and GDB rebuild and pass unchanged.
Adds three optimizations to the CAOP method 0 mult kernel, collectively delivering ~22% wall-time reduction on a 1.6M-basis 30-orbital benchmark at 4 ranks (68.7s → 53.7s davidson time). Bitmask fast-rejection (generalop.h, mult.h): - GeneralOp::PrecomputeMasks() builds per-term m1_ (creation targets, must be occupied in bra) and m2_ (annihilation targets, must be unoccupied), stored as mutable det_vector<size_t>. - m2_ is masked against m1_ at construction to handle orbitals that appear in both creation and annihilation sets of the same term. - mult() calls PrecomputeMasks lazily on first invocation, then checks (~vb & m1_[n]) || (vb & m2_[n]) before entering the k-loop. Rejects ~94% of 2-body terms before any bit manipulation. - GeneralOp friend declaration updated to match new int bit_length. Popcnt sign factor (bit_manipulation.h): - bit_string_sign_factor replaces nested O(bit_length * det_words) loops with __builtin_popcountll: parity over full words before r, then partial word masked to bits below x. - Arguments x and r changed from size_t to int; loop variable k int. int bit_length throughout (mult.h, davidson.h, lanczos.h, sbdiag.h): - bit_length parameter changed from size_t to int in mult(), Davidson(), Lanczos(), and the local variable in caop::diag(). - Eliminates size_t/int mixed arithmetic in the inner loop (q/r/x now plain int division/modulo); static_cast<size_t> on fops_[k].q_ removed.
Adds include/sbd/caop/basic/mult_thrust.h with: - CaopBufferSlider<ElemT>: CPU-staged MPI slider for twk double-buffering - CaopMultKernel<ElemT>: BLOCK=128/SUBWARP=32/GROUPS=4 GPU functor with ballot-dispatch over operator terms, per-lane ket scratch in dynamic shared memory, fused word-descending binary search, and inline __popcll sign computation - CaopMultThrust<ElemT>: Init() precomputes tbs_seq for all tasks at Init time (bs is fixed across Davidson iterations), builds word-descending sorted operator tables with inversion-count reordering sign absorbed into d_coeff; run() does double-buffer twk pipeline with compute/copy overlap - Free mult() forwarding to CaopMultThrust (same signature as CPU) Edits include/sbd/caop/basic/mult.h to include mult_thrust.h and guard the CPU GeneralOp mult() body with #ifndef SBD_THRUST. Design documented in findings/caop-gpu-mult-design-2026-06-29.md. Co-Authored-By: Claude <noreply@anthropic.com>
CaopMultThrust::Init() accessed GeneralOp private members directly,
which compiles only if CaopMultThrust is a declared friend. Instead,
move all GeneralOp/ProductOp/CAOp member accesses into the free mult()
function (already declared friend in generalop.h), extract into
CaopRawTerm/CaopRawOp POD structs, and pass those to Init().
Also use a static CaopMultThrust<ElemT>* in mult() so Init() (which
does one blocking MpiSlide per task to precompute tbs_seq) runs only
on the first Davidson/Lanczos iteration, not on every H*v call.
Builds and links cleanly:
cmake -DSBD_GPU_BACKEND=thrust -DSBD_COMPLEX=ON
-DCMAKE_CXX_FLAGS_RELEASE="-O3" -> [100%] Built target caop_diag
…reuse after Init() no longer runs blocking MpiSlide loops for tbs precomputation. Instead, precompute_tbs_seq() runs on the first run() call (Davidson iteration 1 only), caching results in d_tbs_seq_ / n_kets_per_task_ on GPU. All subsequent run() calls use the cached GPU data directly with zero MPI or copy overhead for tbs. Changes: - Init(): stores bs_ (local bra det_vector) + Allreduce global_max_n_; no blocking MPI for tbs. - precompute_tbs_seq(): extracted from old Init; called once on first run(). - run(): guard on tbs_initialized_; no per-run Allreduce (uses global_max_n_); precomputed d_tbs_seq_[task] path unchanged after first call.
…e per word
compute_contribution() previously allocated G*SW*(ES+1) size_t words in shared
memory (vk_scr) as scratch for the ket being constructed. This was unnecessary:
the outer loop processes words descending (w = ES-1..0), so at word w all lower
words w2<w are still == my_bra[w2].
Replacement:
- cur_w: one size_t register per lane; tracks evolving state of word w
- n_occ_base: int computed once per word as sum of popcll(my_bra[w2]) for w2<w;
reused by all ops within the word (was recomputed per op against vk_scr)
- n_occ, n_occ_base: int (max value <= 64*elem_size, int is sufficient)
Shared memory savings: G*SW*(ES+1)*8 bytes = 4*32*2*8 = 2048 B/block for ES=1.
For ES=1 (30-orbital case) n_occ_base is always 0 and the loop is dead code.
…slide[0]==0 Move d_twk_[2], h_twk_[2], streams, and events from per-call stack to class members. Allocate on the first run() call alongside tbs_seq precompute; free in destructor. Eliminates one cudaMallocHost, two device_vector resizes, two stream creates, and four event creates per Davidson mult() call (~60 per solve). When slide_[0]==0 (the common case: task 0 is always the local rank), the previous code did twk_init = wk (full copy) then memcpy'd twk_init to h_twk_[0]. Now use a const pointer to wk.data() directly and memcpy once.
…ze>1 When ops span multiple words they are sorted word-DESC before building the GPU operator tables. The resulting permutation sign was baked into d_coeff unconditionally. For bosonic operators (sign=false) this sign is spurious — swapping bosonic operators leaves the coefficient unchanged — so the GPU gave the wrong energy for any bit_length that produces elem_size>1 (e.g. 30-orbital with bit_length=12 gives elem_size=3). The elem_size=1 case was unaffected because all ops land in a single word, so the sorted order is always identical to the original (inversions=0). Fix: gate the reordering sign on the sign flag. Also remove stray "TEST MpiAllreduce" debug print from mpi_utility_thrust.h. Validated: 1-rank smoke test with --bit_length 12 (elem_size=3) now gives 1.244264613763439 on GPU, matching the CPU reference (job 2584467).
Replace getline-per-line loop with a single ifs.read() into a flat buffer, then stride through at (total_bit_length+1) bytes per record and parse in-place with no intermediate string or vector allocations. Speedup on 48-site benchmark (4 ranks x 5M-line / 234 MB shards, NVMe Lustre, different node so cold): cold read: 44.9s -> 1.92s (23x) warm read: 4.8s -> 1.90s (2.5x) Cold improvement: one syscall vs ~30K 8KB reads on Lustre. Warm improvement: eliminates 5M std::string + 5M vector<size_t> allocations per rank per load. Energy unchanged: 13.00000000000032 across all 5 reps.
Throw std::runtime_error on any unexpected byte in the text reader.
Pre-compute n_occ_base as the full-bra popcount across all words before the w loop. At the top of each iteration subtract __popcll(cur_w) — which equals __popcll(my_bra[w]) before any ops touch cur_w — leaving the sum for words w2 < w needed by sign logic. Reduces popcount work from O(elem_size^2) to O(elem_size) per term. For elem_size=1 (30-orbital) no change. For larger elem_size the inner loop over w2 is eliminated entirely. Correctness verified: elem_size=2 (48-site, 4r, energy 13.0) and elem_size=3 (--bit_length 12, 1r, energy 1.244264613763439).
compute_contribution(): lo, hi, lo0, L, R, mid converted from size_t to int — GPU 32-bit arithmetic, avoids 64-bit index overhead. n_kets is already int so hi = n_kets needs no cast. Also: upper_bound now starts from lo0 (the lower_bound start) rather than the lower_bound result, giving more uniform memory access across threads. 48-site benchmark (4 ranks x 5M kets, 72 OMP, elem_size=2): baseline (5f00e1c): davidson ~16.1s n_occ_base (9750921) + lo0 + int: davidson ~14.2s (~12% improvement) int alone: neutral vs lo0-only (within noise)
Default SBD_CAOP_MULT_BLOCK_SIZE: 128 → 64 Default SBD_CAOP_MULT_MIN_BLOCKS_PER_SM: 4 → 32 bs64-mbpsm32: 32 regs/thread, 32 blocks/SM, 64 warps/SM (100% occupancy). Benchmark (48-site, 4 ranks x 5M kets, gh152): davidson 13.19s vs 14.22s baseline (mbpsm4/bs128) — 7.2% improvement. Tighter rep-to-rep variance than bs128 at the same occupancy level. Subwarp scan (sw8/sw16/sw32, all bs64-mbpsm32, gh152): sw8: 13.84s (+4.9% vs sw32) — 18 ballot iters/bra for 141-term ham sw16: 13.30s (+0.9% vs sw32) sw32: 13.19s (winner) — 5 ballot iters/bra; default unchanged Static_assert changes: SBD_CAOP_MULT_BLOCK_SIZE % 32 == 0 (must be multiple of full warp) 32 % SBD_CAOP_MULT_SUBWARP_SIZE == 0 (SW must divide 32) (removed: BLOCK_SIZE % SUBWARP_SIZE == 0 — implied by the two above)
Contributor
Author
|
This change includes #66. |
Contributor
|
I really like the direction of this PR. In particular, precomputing per-term masks (and, more generally, preparing execution-friendly data from the normal-ordered GeneralOp) makes a lot of sense to me. GeneralOp is mainly designed for operator algebra, so having a representation optimized for repeated application to a basis seems like a very natural optimization. One question about the lifetime of the cached driver: In our SQD workflow, sbd::diag may be called multiple times within the same process, and the Hamiltonian may change (for example during orbital optimization), or the basis may be expanded between calls. Since CaopMultThrust is stored in a function-local static, has this use case been considered/tested? I just wanted to make sure the cached data are correctly reinitialized when the Hamiltonian or basis changes. Also, I wonder if it would be worth adding an assertion to verify that the precomputed mask size always matches the basis elem_size, just to catch any future inconsistencies. |
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.
CPU binary:
diag-caop3-cpu-2026-07-07(42b2bf4), 72 OMP threads/rank.GPU binary:
diag-caop3-thrust-2026-07-07(42b2bf4), 1 GPU/rank.Workload: 48-site 20M-det,
--do_redist_basis 1 --iteration 6 --block 10