pythonでlinpackしてみる

# linpackをpythonで書いてみる

## 概要

過去何度かスーパーコンピュータの性能を示すbenchmarkの一つlinpackを手元で試してみようという記事を書いた。

* [AWS-GPUとスパコンを比較する方法(P3インスタンス編)-スパコン用ベンチマークソフトを動かす](https://qiita.com/shibacow/items/459949bb3bc81a58526a)
* [AWS-GPUとスパコンを比較する方法(P2インスタンス編)-スパコン用ベンチマークソフトを動かしてみる](https://qiita.com/shibacow/items/ab6ed97dd94578cd1d0e)
* [AWS-GPUとスパコンを比較する方法-スパコン用ベンチマークソフトを動かしてみる](https://qiita.com/shibacow/items/aa83cfc81ed493a4c25d)
* [AWSとスパコンを比較する方法-スパコン用ベンチマークソフトを動かしてみる](https://qiita.com/shibacow/items/08490ac51aeb0d131d2e)
* [お手元のマシンとスパコンを比較する方法-スパコン用ベンチマークソフトを動かしてみる](https://qiita.com/shibacow/items/753e1d46882152e46757)

なんどもlinpackと書いているが、実際のところ、linpackが何をしているかよく理解していない。
そのため、pythonでlinpackを書いてみて、何をしているか把握しようと考えた。

## linpackとは

[linpack](https://ja.wikipedia.org/wiki/LINPACK)は要は、n×n の線型方程式系 Ax = b を解く（xを求める）プログラムらしい。
アルゴリズム的な話を[こちら](https://www.slideshare.net/kawa0810/linpack)の資料を参考にした。
pythonで実装がないかと調べた見たところ、
[資料](http://www.hpcchallenge.org/presentations/sc2007/hpcc-awards2-2007-star-p.pdf)が見つかった。

## pythonで解いてみる。

[numpy](http://www.numpy.org/)で解く場合要は、

```py
import numpy as np
n = 100
A = np.random.rand(n, n)
b = np.random.rand(n, 1)
x = np.linalg.solve(A, b)
```

と解いているということらしい。

上で紹介した[この資料](http://www.hpcchallenge.org/presentations/sc2007/hpcc-awards2-2007-star-p.pdf)を参考にプログラムを書いてみた。

プログラムは、こちらの [github](https://github.com/shibacow/hpl_test/blob/master/hpl.py)にあげている。

ちゃんと意味はわかっておらず、ベンチマーク本体は、`np.linalg.solve(A, b)` で下の、`r0 = cn.linalg.norm(r, cn.inf)` は検算なのだろう。

手元のマシンで、nの数値を、10000にして、大体20Gflops/Secくらい出たので、それっぽいかなと感じた。

## cuda(cupy) を使う

GPUを使うと早いらしいので、NVIDAのGPUを使って演算出来ないか試してみた。pythonには、CUDAをnumpyライクに扱う [CuPy](https://cupy.chainer.org/) がありそちらを使って演算してみる。手元にcudaの動くマシンは無かったので、GCPのGCEインスタンスのGPU(P100)をアタッチして試した。

`cupy` は `numpy` と互換性があるようで、以下のように書き換えれば動いた


```python
import cupy as cu
n = 100
A = cu.random.rand(n, n)
b = cu.random.rand(n, 1)
x = cu.linalg.solve(A, b)

```

`nvida-smi` を見るときちんと動いているようだった。

15000次元で、テストしたところ。

`Linpack benchmark 532.0809  GFlops/Sec` とあり、 **532Gflops** 出たようだ（数値は怪しいのであまり鵜呑みにしないでください）。数値は怪しいが、少なくとも、numpyで動かしたときよりは早かった。

また、numpyでも試した。わかったことは、numpyもcpuがあれば並列で動くようだ。

## cuda(cupy) を単精度で使う

上の、cupyを使った際に、`a = cn.random.rand(n, n).astype(cn.float)` を指定した。numpyは、floatを指定すると、倍精度(float64)が指定されるそうだ。cudaなどのgpuはfloat32の方が性能が出るので、`a = cn.random.rand(n, n).astype(cn.float32);` と`float32`を明示して単精度で演算してみた。
そうすると、25000次元で解くことが出来て、計算速度も、 **1110.7489  GFlops/Sec** と1Tflopsを超えた(数値的に怪しいので鵜呑みにしないでください。ただ、上記プログラムより早かったです)。

ちなみに、numpy(CPU)()でやる場合、8CPUで `Linpack benchmark 81.1896  GFlops/Sec` で **81.1896GFlops** であった。大体、13倍くらい早い。

## まとめ

 * スパコンで使われるlinpackが一体何をしているか分からなかったので、pythonで書いてどういうものか見てみた。
 * numpyだと Ax=b を解くのは、 `x = np.linalg.solve(A, b)` とすれば良いようだ。これで、スパコンがやっている演算と **基本的** には同じ事ができる。
 * ついでにcupyも試してみた。GPUを使うと早くなるようだ。


## 追記

### NVida A100で試す。

GCPでは、A100もあるので、Nvidia A100で試す。SPOT VMで立ち上がるかも調べた。シンガポールでは、 `a2-highgpu-1g` が立ち上がるので、`--provisioning-model=SPOT` をつけて立ち上げる。1時間3ドルくらいかかるが、SPOTVMのため、40分0.58ドルくらいだった。

上記の hpl.pyを試す。単精度で試した。

```
INFO:root:performance=15341.996324664775 verified=True r0=0.48521143198013306 r1=0.003971508238464594 r2=0.0008978209225460887
INFO:root:Linpack benchmark 15341.9963  GFlops/Sec
```
なので、**15.3** TFlops出ているようだ。次元 45000 で試した。

https://qiita.com/yukoba/items/10d0ba3fb1d19a6ab6a5

を見ると、A100のFP32は19.5TFlopsのようなので、大体75%ぐらいの性能は出ているようだ。 FP64の最大値も、19.5TFlops だ。マックス性能が出る前に、GPUのメモリ容量を超えたと思われる。

ちなみに倍精度で試すと、 **10.6** TFlops程度であった。 次元 30000 で試した。

nvidia-smiを見ると、GPUを使っているようだった。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/26986/715f2d3d-9990-affd-f3c0-2d6316216f48.png)


### NVida T4 で試す。

同様に、NVida T4 でも試した。こちらもSPOT VMで起動。1時間半くらい使い、0.09ドル程度であった。FP32 で、 **2.681** TFlos,FP64で、 **0.238** TFlopsであった。

```
INFO:root:Linpack benchmark dimensions=30000  gpu memory=13734 loop=4 type=fp32 2681.5822  GFlops/Sec
INFO:root:Linpack benchmark dimensions=20000  gpu memory=12208 loop=2 type=fp64 238.7384  GFlops/Sec
```