Threads, BLAS & LAPACK

Threads, BLAS & LAPACK#

The first and easiest speedup uses multiple cores of one machine — shared-memory threads (the primer’s first model). This chapter shows the one knob you need, run(nthread=N), explains what qc-rs actually threads, and demystifies the BLAS library that does the heavy linear algebra underneath.

The one knob: `run(nthread=N)`#

Set the number of CPU cores a calculation may use with nthread= on .run():

import qc
water = "O 0 0 0.117; H 0 0.757 -0.469; H 0 -0.757 -0.469"

m = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom").scf(ref="r").run(nthread=8)

nthread is cores per (MPI) rank: it sets both the qc-rs kernel worker threads and the BLAS/LAPACK thread count together, with no environment-variable juggling. The count resolves as nthread argument > QC_THREADS environment variable > all available cores. QC_THREADS is the external fallback for when you cannot pass nthread (e.g. a script you do not control):

export QC_THREADS=8      # used when run(nthread=...) is not given

Crucially — as promised in the index — more threads change only the speed, not the answer:

e1 = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom").scf(ref="r").run(nthread=1).scf.energy
e4 = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom").scf(ref="r").run(nthread=4).scf.energy
# e1 == e4 == -76.026794   (bit-identical)

What qc-rs actually threads#

Not every step benefits equally (recall Amdahl). In qc-rs today:

Threaded: the KS-DFT exchange–correlation grid — the per-grid-block loop (LDA, GGA, and meta-GGA, for both RKS and UKS) runs on std::thread workers, each with its own scratch, reducing into the shared V_xc. This is the main thread-parallel region.
Threaded via BLAS: the dense linear algebra — the Fock diagonalization, the large matrix multiplies — runs in the multi-threaded BLAS/LAPACK library (below).
Still serial: the one- and two-electron integral assembly. On a small molecule this can dominate, so do not expect a small water job to speed up much with more cores — its bottleneck is not (yet) threaded.

A single big scratch buffer (the Workspace) is carved into per-thread lanes before the parallel loop, so the threads share one physical allocation with no per-iteration memory churn — the reason qc-rs threading is cheap.

The layer underneath: BLAS & LAPACK#

Almost every quantum-chemistry program leans on BLAS (Basic Linear Algebra Subprograms) and LAPACK for matrix multiplication and eigensolves — decades-tuned libraries that are the real engine of the dense algebra. qc-rs links one of two:

library	when
Intel MKL	Intel/AMD CPUs, especially HPC clusters (usually fastest there)
OpenBLAS	open-source default; portable, no Intel toolchain needed

You chose one at build time (the intel-mkl-system vs openblas-* feature). Both are multi-threaded, and run(nthread=) sets their thread count for you.

Important

MKL threading layer on Linux/HPC An MKL build needs the threading layer pinned, or you get an import error (undefined symbol: omp_get_num_procs) or silently-wrong threading. Set MKL_THREADING_LAYER=GNU (and the LP64 interface MKL_INTERFACE_LAYER=LP64) in your environment — make setup writes these into .vscode/my.env for you, and the shell env.sh exports them. This is the single most common HPC setup snag.

Oversubscription: the one trap#

Oversubscription happens when you ask for more threads than you have cores — for example running \(W\) qc-rs worker threads that each call a BLAS routine using \(M\) BLAS threads, giving \(W \times M\) threads fighting over the cores. That thrashes the CPU and runs slower, not faster. qc-rs avoids this by design (its threaded XC region uses a single-threaded GEMM, so there is no nested multi-threaded BLAS to oversubscribe), and because run(nthread=) sets one coherent count. The rule for you: when combining with MPI (next chapter), keep ranks_per_node × nthread ≤ cores_per_node.

Worked example: cores don’t change the answer#

import qc
water = "O 0 0 0.117; H 0 0.757 -0.469; H 0 -0.757 -0.469"

for n in (1, 2, 4):
    e = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom").scf(ref="r").run(nthread=n).scf.energy
    print(f"nthread={n}: E = {e:.6f}")
# nthread=1: E = -76.026794
# nthread=2: E = -76.026794
# nthread=4: E = -76.026794

The energy is identical; on a large enough DFT job the wall-clock would fall with nthread (measure it — report the absolute seconds, not just the ratio).

Exercise 17

You set run(nthread=16) on a water/cc-pVDZ RHF and it is no faster than nthread=2. Is qc-rs broken? What is the likely reason?
On a 64-core node you run 8 MPI ranks and want to use the whole node without oversubscribing. What nthread do you pass?
Your import qc fails with undefined symbol: omp_get_num_procs on an MKL build. What is the fix?

Solution to Exercise 17

Not broken. Water/cc-pVDZ is tiny and its bottleneck is the serial integral assembly, not the threaded XC/BLAS regions (Amdahl) — there is too little parallel work to benefit from 16 cores. Bigger molecules / DFT grids scale.
nthread=8 — so 8 ranks × 8 threads = 64 = cores_per_node, fully using the node with no oversubscription.
Set MKL_THREADING_LAYER=GNU (and MKL_INTERFACE_LAYER=LP64) in the environment; make setup puts these in .vscode/my.env / env.sh.

Threads scale you within one machine. To go beyond it — across many machines — you need MPI and the interconnect.