# Threads, BLAS & LAPACK

The first and easiest speedup uses **multiple cores of one machine** — shared-memory **threads** (the
[primer](parallel-computing-primer.md)'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()`:

```python
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):

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

Crucially — as promised in the [index](index.md) — **more threads change only the speed, not the answer**:

```python
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](parallel-computing-primer.md)). 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](../00-intro/installation-and-make-setup.md) (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

```python
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}
:label: ex-threads

1. 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?
2. On a 64-core node you run 8 MPI ranks and want to use the whole node without oversubscribing. What
   `nthread` do you pass?
3. Your `import qc` fails with `undefined symbol: omp_get_num_procs` on an MKL build. What is the fix?
:::

:::{solution} ex-threads
:class: dropdown

1. 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.
2. `nthread=8` — so `8 ranks × 8 threads = 64 = cores_per_node`, fully using the node with no oversubscription.
3. 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](mpi-and-interconnects.md).