Skip to content

fix: pin BLAS to one thread in PWM kernels to stop thread explosion - #37

Merged
aviezerl merged 2 commits into
masterfrom
fix/nested-blas-thread-explosion
Jul 13, 2026
Merged

fix: pin BLAS to one thread in PWM kernels to stop thread explosion#37
aviezerl merged 2 commits into
masterfrom
fix/nested-blas-thread-explosion

Conversation

@aviezerl

Copy link
Copy Markdown
Contributor

Problem

extract_pwm() (and compute_pwm() / compute_local_pwm() / calc_seq_pwm()) open thousands of OS threads and can fail outright on core-limited machines. Reported via iceqream::infer_trajectory_motifs(), whose default path runs extract_pwm().

Root cause is nested parallelism: these kernels run a small per-sequence BLAS dgemm (PWMWorker::compute_matrix_scores) inside an RcppParallel/TBB parallelFor over sequences. With a multi-threaded BLAS (MKL or OpenBLAS), each TBB worker spawns its own BLAS thread team, so one call opens n_threads × n_blas_threads threads:

setThreadOptions OS threads (128-core node, MKL)
16 1025 (64×16+1)
128 8193

At the set_parallel() default (~102 threads) that's ~6500 threads for a single extract_pwm. It oversubscribes badly and, on a cluster job with a thread/process (cgroup pids) limit, fails with a thread-creation error.

(A bare RcppParallel parallelFor opens exactly n_threads; the multiplier is entirely the inner threaded dgemm. OMP_NUM_THREADS=1 at launch also collapses it - confirming the diagnosis.)

Fix

New helper local_serial_blas() pins the BLAS to one thread for the duration of these kernels and restores the previous count on exit. The per-sequence dgemm is tiny and the parallelFor over sequences is the right parallelism granularity, so the inner BLAS should never be threaded here.

compute_pwm / compute_local_pwm / calc_seq_pwm:
    local_serial_blas()   # BLAS -> 1 thread, restored on return
    ...parallelFor...

Verification (MKL)

  • 16-thread extract_pwm: 1025 → 17 OS threads; compute_pwm / compute_local_pwm likewise.
  • Results identical (all.equal TRUE); throughput same or slightly better (the threaded dgemm was pure overhead at this size).
  • BLAS thread count restored after the call, so the rest of the session keeps threaded BLAS.
  • Full testthat suite passes.

Adds RhpcBLASctl to Imports (the portable cross-BLAS thread control - RhpcBLASctl::blas_set_num_threads). Bumps to 0.0.10.

https://claude.ai/code/session_01KfEW7AGsSzjhMqiY8htdnq

aviezerl added 2 commits June 22, 2026 16:23
compute_pwm(), compute_local_pwm() and calc_seq_pwm() (and thus extract_pwm())
run a per-sequence BLAS dgemm inside an RcppParallel/TBB parallelFor over
sequences. With a multi-threaded BLAS (MKL, OpenBLAS) every TBB worker spawns
its own BLAS thread team, so a single call opens n_threads * n_blas_threads OS
threads - e.g. 64x16 = 1025 at 16 threads on a 128-core node, and thousands at
the set_parallel() default. That oversubscribes the machine and, on a cluster
job with a thread/process (cgroup pids) limit, makes the call FAIL with a
thread-creation error (reported via iceqream::infer_trajectory_motifs ->
extract_pwm).

The per-sequence dgemm is tiny and the parallelFor over sequences is the right
granularity, so the inner BLAS must be serial. New helper local_serial_blas()
pins the BLAS to one thread for the duration of these kernels and restores the
previous count on exit (the rest of the session keeps its threaded BLAS).

Verified on MKL: a 16-thread extract_pwm drops from 1025 to 17 OS threads,
results identical, throughput same or slightly better. Adds RhpcBLASctl to
Imports (the only portable cross-BLAS thread control). Bumps to 0.0.10.

Claude-Session: https://claude.ai/code/session_01KfEW7AGsSzjhMqiY8htdnq
@aviezerl
aviezerl merged commit ca7d77d into master Jul 13, 2026
6 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant