# スレッド・BLAS・LAPACK

最初で最も易しい高速化は、**1 台のマシンの複数コア** — 共有メモリの**スレッド**（[入門](parallel-computing-primer.md)
の第 1 モデル）を使います。本章は必要な唯一のノブ `run(nthread=N)` を示し、qc-rs が実際に何をスレッド化するかを
説明し、その下で重い線形代数を担う **BLAS** ライブラリを解き明かします。

## 唯一のノブ:`run(nthread=N)`

計算が使ってよい CPU コア数を、`.run()` の **`nthread=`** で設定します:

```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` は**（MPI）ランクごとのコア数**です:qc-rs カーネルのワーカースレッドと BLAS/LAPACK のスレッド数を、
環境変数の調整なしに*両方*まとめて設定します。数は **`nthread` 引数 > `QC_THREADS` 環境変数 > 利用可能な全コア**
の順で解決します。`QC_THREADS` は `nthread` を渡せないとき（例:自分で制御できないスクリプト）の外部フォール
バックです:

```bash
export QC_THREADS=8      # run(nthread=...) が与えられないときに使われる
```

重要なのは — [index](index.md) で約束したとおり — **スレッドが増えても変わるのは速さだけで、答えではない**こと:

```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   （ビット単位で同一）
```

## qc-rs が実際にスレッド化するもの

すべてのステップが等しく恩恵を受けるわけではありません（[Amdahl](parallel-computing-primer.md) を思い出して）。
現在の qc-rs では:

- **スレッド化**:**KS-DFT 交換相関グリッド** — グリッドブロックごとのループ（LDA・GGA・meta-GGA、RKS・UKS
  とも）が `std::thread` ワーカーで走り、各々が自分のスクラッチを持ち、共有の `V_xc` に簡約します。主要な
  スレッド並列領域です。
- **BLAS 経由でスレッド化**:密な線形代数 — Fock の対角化、大きな行列積 — がマルチスレッド BLAS/LAPACK
  ライブラリ（下記）で走ります。
- **まだ逐次**:1 電子・2 電子の**積分組み立て**。小さな分子ではこれが支配しうるので、小さな水のジョブがコアを
  増やして大きく速くなるとは期待しないでください — そのボトルネックは（まだ）スレッド化されていません。

1 つの大きなスクラッチバッファ（`Workspace`）が並列ループの*前*にスレッドごとのレーンに切り分けられるので、
スレッドは 1 つの物理確保を反復ごとのメモリ churn なしに共有します — qc-rs のスレッド化が安価な理由です。

## 下層:BLAS と LAPACK

ほぼすべての量子化学プログラムは、行列積と固有値解法に **BLAS**（Basic Linear Algebra Subprograms）と **LAPACK**
に頼ります — 数十年チューニングされた、密な代数の真のエンジンです。qc-rs は次の 2 つのいずれかをリンクします:

| ライブラリ | いつ |
|---|---|
| **Intel MKL** | Intel/AMD CPU、特に HPC クラスタ（そこでは通常最速） |
| **OpenBLAS** | オープンソースの既定;移植性が高く、Intel ツールチェーン不要 |

[ビルド時](../00-intro/installation-and-make-setup.md)に 1 つを選びました（`intel-mkl-system` 対 `openblas-*`
機能）。どちらも**マルチスレッド**で、`run(nthread=)` がそのスレッド数を設定します。

:::{important} Linux/HPC での MKL スレッド層
MKL ビルドはスレッド層を固定する必要があり、さもないとインポートエラー（`undefined symbol:
omp_get_num_procs`）や、静かに誤ったスレッド化になります。環境で **`MKL_THREADING_LAYER=GNU`**（と LP64
インターフェース `MKL_INTERFACE_LAYER=LP64`）を設定してください — `make setup` がこれらを `.vscode/my.env` に
書き、シェルの `env.sh` がエクスポートします。HPC セットアップで最もよくある落とし穴です。
:::

## オーバーサブスクリプション:唯一の罠

**オーバーサブスクリプション**は、コア数より*多く*のスレッドを要求したときに起こります — 例えば $W$ 個の qc-rs
ワーカースレッドが各々 $M$ 個の BLAS スレッドを使うルーチンを呼ぶと、$W \times M$ 個のスレッドがコアを奪い合い
ます。CPU が thrash して、速くなるどころか*遅く*なります。qc-rs は設計上これを避けます（スレッド化された XC
領域は単一スレッド GEMM を使うので、オーバーサブスクライブする入れ子のマルチスレッド BLAS がない）し、
`run(nthread=)` が 1 つの整合した数を設定します。**あなた**への規則:MPI と組み合わせるとき（次章）、
`ノードあたりランク数 × nthread ≤ ノードあたりコア数` を保つこと。

## 総合例:コアは答えを変えない

```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
```

エネルギーは同一です;十分大きな DFT ジョブなら*実時間*は `nthread` とともに下がります（測って — 比だけでなく
絶対秒数を報告してください）。

:::{exercise}
:label: ex-threads

1. 水/cc-pVDZ の RHF に `run(nthread=16)` を設定したのに `nthread=2` と変わらない速さでした。qc-rs は壊れて
   いますか。考えられる理由は。
2. 64 コアのノードで 8 MPI ランクを走らせ、オーバーサブスクライブせずにノード全体を使いたいです。どの `nthread`
   を渡しますか。
3. MKL ビルドで `import qc` が `undefined symbol: omp_get_num_procs` で失敗します。対処は。
:::

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

1. 壊れていません。水/cc-pVDZ は小さく、そのボトルネックはスレッド化された XC/BLAS 領域ではなく**逐次の積分
   組み立て**（Amdahl）です — 16 コアで恩恵を受けるほどの並列仕事がありません。大きな分子 / DFT グリッドは
   スケールします。
2. `nthread=8` — `8 ランク × 8 スレッド = 64 = ノードあたりコア数`、オーバーサブスクライブなしでノードを完全に
   使います。
3. 環境で **`MKL_THREADING_LAYER=GNU`**（と `MKL_INTERFACE_LAYER=LP64`）を設定;`make setup` がこれらを
   `.vscode/my.env` / `env.sh` に入れます。
:::

スレッドは 1 台のマシン内でスケールさせます。それを超えて — 多数のマシンにわたって — 行くには [MPI と相互接続](mpi-and-interconnects.md)
が必要です。