FortranコードにBLASをリンクする

Fortranで書いたコードにBLASをリンクしたい。
前提として、fortranのコンパイラ(gfortranとか)があることを想定。念頭に置いているのは、Linux様環境。
実環境: Kubuntu 19.04, Intel(R) Core(TM) i5-8259U (147.2 GFLOPs)
コンパイラ: gfortran GCC 8.3.0 (Ubuntu 8.3.0-6ubuntu1), apt-get で入れたもの

# Linuxのパッケージ管理で入れたBLASをリンク

`apt-get`なり、`yum`なり、`zypper`で入れたBLASをリンクする。gfortranを使っているとき、

```bash
gfortran -O3 main.f90 -lblas
```
これで大方通るはず。

通らない時は、BLASが置いてあるディレクトリにpathが通っていない。`$ ldconfig -p`とかで、一応候補にあるかどうか確認してみる。なぜだか、ライブラリのpathを知っていて、そこには`LD_LIBRARY_PATH`を通していない場合は以下のようにして、リンクできる:

```bash
gfortran -O3 main.f90 -L/path/to/be/specified/ -lblas
```

# 自分でビルドしたBLASをリンク
[ BLAS (Basic Linear Algebra Subprograms)](http://www.netlib.org/blas/)からコードを取ってきて自分でビルドして、リンクできる。ライブラリというよりも、分割コンパイルしたオブジェクトを自分のコードにくっつけると思った方が、感覚としては近いかも。自分が管理者ではない計算機の、どこにBLASがあるのかわからないとき、ひとまず自分でビルドしたものをリンクすれば、動作確認はできる。
例として、blas-3.8.0.tgzをダウンロードしてきたとする。

## 展開
ダウンロードしたコードを展開して、適当なディレクトリ(例えば'/home/user/sources/`とか)にうつす。

```bash
tar -zxvf blas-3.8.0.tgz
mv ./BLAS-3.8.0 /home/user/sources/
```

## ビルド
展開したディレクトリに移って、

```bash
make 
```
でおしまい。`Makefile`は`make.inc`をインクルードしていて、コンパイルのオプションとか出来たオブジェクトの名前とかを変えられる。デフォルトのままmakeすると静的ライブラリ`blas_LINUX.a`が出来上がる。

## 所望の場所に配置
自分で使う事を念頭に、管理者権限の必要のない場所に置く。例えば`/home/user/opt/BLAS-3.8.0/lib`を作って、その中に置く。

## マイBLASをリンク

```bash
gfortran -O3 main.f90 /home/user/opt/BLAS-3.8.0/lib/blas_LINUX.a
```
でコンパイルすると、晴れてマイBLASがリンクできる。

# OpenBLAS
[OpenBLAS : An optimized BLAS library](https://www.openblas.net/)
大変、速いらしい。

## 展開＆ビルド
`tar -zxvf`で展開。普通に`make`でビルドが通る。数分ビルドにはかかるかも。
標準出力の内容をを見る限り、BLAS, CBLAS, LAPACK, LAPACKEがビルドされているみたい。
ビルド後に`make PREFIX=/path/to/be/installed install`で、ライブラリ群が所望の場所にインストールされる。

リンクは

```bash
gfortran -O3 main.f90 -L/home/user/opt/OpenBLAS-0.3.7/lib -lopenblas
```
か

```bash
gfortran -O3 main.f90 /home/user/opt/OpenBLAS-0.3.7/lib/libopenblas.a -lpthread
```
で通る。

##　トラブルシューティング
### Pathの通しそびれ
```bash
./a.out: error while loading shared libraries: libopenblas.so.0: cannot open shared object file: No such file or directory
``` 
共有ライブラリなので、pathが通っていなくても、コンパイルは通る。しかし、pathが通っていないと、実行時に`libopenblas.so.0`を探し出せなくてコケる。例えば、`ldd ./a.out`をみると

```bash
        linux-vdso.so.1 (0x00007ffd9e2d3000)
        libopenblas.so.0 => not found
        libgfortran.so.5 => /lib/x86_64-linux-gnu/libgfortran.so.5 (0x00007f75be493000)
        libc.so.6 => /lib/x86_64-linux-gnu/libc.so.6 (0x00007f75be2a8000)
        libquadmath.so.0 => /lib/x86_64-linux-gnu/libquadmath.so.0 (0x00007f75be25e000)
        libm.so.6 => /lib/x86_64-linux-gnu/libm.so.6 (0x00007f75be110000)
        libgcc_s.so.1 => /lib/x86_64-linux-gnu/libgcc_s.so.1 (0x00007f75be0f6000)
        /lib64/ld-linux-x86-64.so.2 (0x00007f75bed1d000)
```
となって、`libopenblas`が見つかっていないことがわかる。`export LD_LIBRARY_PATH=/home/user/opt/OpenBLAS-0.3.7/lib:$LD_LIBRARY_PATH`でpathを追加すると、動く。

### 静的ライブラリをリンクするとライブラリ不足で怒られる
誤って

```bash
gfortran -O3 main.f90　/home/user/opt/OpenBLAS-0.3.7/lib/libopenblas.a
```
でコンパイルすると

```bash
/usr/bin/ld: /home/shinohara/opt/OpenBLAS-0.3.7/lib/libopenblas.a(blas_server.o): in function `blas_thread_init':
blas_server.c:(.text+0x40a): undefined reference to `pthread_create'
/usr/bin/ld: /home/shinohara/opt/OpenBLAS-0.3.7/lib/libopenblas.a(blas_server.o): in function `goto_set_num_threads':
blas_server.c:(.text+0x982): undefined reference to `pthread_create'
/usr/bin/ld: /home/shinohara/opt/OpenBLAS-0.3.7/lib/libopenblas.a(blas_server.o): in function `blas_thread_shutdown_':
blas_server.c:(.text+0xbcb): undefined reference to `pthread_join'
collect2: error: ld returned 1 exit status
```
と出てくるので、その場しのぎでopenmpのオプションをつけて

```bash
gfortran -fopenmp -O3 main.f90　/home/user/opt/OpenBLAS-0.3.7/lib/libopenblas.a
```
とすると一応コンパイルは通る。~~ただ、a.outを実行するとSegmentation fault (core dump)で落ちる。~~と思っていたら、知らないうちに動いていた。(なんでだ)

足りないライブラリをリンクすればよいので、

```bash
gfortran -O3 main.f90　/home/user/opt/OpenBLAS-0.3.7/lib/libopenblas.a -lpthread
```
でも、可。

### 思いのほか遅い
OpenBLASはスレッド並列化されていて、何も指定しないとシステムの最大スレッド並列数で実行される。これで遅く見えている可能性があるので、`export OMP_NUM_THREADS=1`等で、適切な値にするとよい。
というか、スレッド並列だと、`cpu_time`は正しく時間を計測してくれていない気がする。

# ATLAS
[Automatically Tuned Linear Algebra Software (ATLAS)](http://math-atlas.sourceforge.net/)
アーキテクチャごとに最適化したBLASをインストールできる。

## 展開＆ビルド
[Basic Steps of an ATLAS install](http://math-atlas.sourceforge.net/atlas_install/node6.html)を参考に。
`tar xf atlas3.10.3.tar.bz2`で展開。
出来たディレクトリ`ATLAS`を`ATLAS3.10.3`に変更
このディレクトリの中にディレクトリ`BLDdir`を作成。(ソースのあるディレクトリでconfigureしようとすると、怒られる。)
CPUのthrottlingを止めるのはだるかったので、`-cripple-atlas-performance`をつけてそのまま./configure.

```bash
../configure --prefix=/home/user/opt/ATLAS3.10.3 --cripple-atlas-performance
make build #結構時間かかる 数時間とか(自分のケースでは約5時間)
make check
make ptchek
make time
make install
```
何か、make installの過程で`make[1]: [Make.top:673 install_lib] Error 1 (ignored)`の文字列が見えたが、どうも共有ライブラリ関係がうまくいっていない様子(というか、そもそもビルドのところ共有ライブラリが出来ていないようにみえる)。
ひとまず、静的ライブラリのみ使うことにする。

リンクは

```bash
gfortran -O3 main.f90 /home/user/opt/ATLAS3.10.3/lib/libf77blas.a /home/user/opt/ATLAS3.10.3/lib/libatlas.a 
```
で通る。この`libf77blas.a`と`libatlas.a`の両方を入れないと怒られる。

# MKL
言わずと知れた、Intel謹製のライブラリ。BLAS, LAPACK, FFTのみならず、他にもいろいろな数値計算ライブラリが同梱。
[Intel® Math Kernel Library (Intel® MKL) | Intel® Software](https://software.intel.com/en-us/mkl)からダウンロード。ダウンロードにあたって、名前やらメールアドレス、所属の登録が必要。登録したメールアドレスに`*.lic`が届く。(が、このライセンスファイルを使わずインストールできてしまった)

## 展開＆インストール
`tar -zxvf l_mkl_2020.0.166.tgz`で展開。
出来た`l_mkl_2020.0.166`に入って、

```bash
sh ./install.sh
```
で、指示に従って入力していく。すると、`/home/user/intel`が作られて、この下にライブラリ群がインストールされている。

## リンク
自分の環境に合わせて、[Intel® Math Kernel Library Link Line Advisor | Intel® Software](https://software.intel.com/en-us/articles/intel-mkl-link-line-advisor/)でリンク。

インテルコア、静的リンク、32-bit 整数、スレッド並列なしだと

```bash:compilse.sh
MKLROOT=/home/user/intel/compilers_and_libraries_2020.0.166/linux/mkl
OPT=" -m64 -I${MKLROOT}/include "
LN=" -Wl,--start-group ${MKLROOT}/lib/intel64/libmkl_gf_lp64.a ${MKLROOT}/lib/intel64/libmkl_sequential.a ${MKLROOT}/lib/intel64/libmkl_core.a -Wl,--end-group -lpthread -lm -ldl "
gfortran -O3 ${OPT} main.f90 ${LN}

```
でコンパイル。

# 速度を比較してみよう
色々なBLASを入れてみたので、1000次元の行列でのDGEMMの速度比較を行ってみる。
**他の環境で、調べていただいたものがコメント欄にあるので、そちらも参照していただくと良いと思います。**

## 使うコード

```fortran:main.f90
Program BLAS_check
  Implicit None
  Integer(4), Parameter :: N = 1000, Niter = 10
  Integer(4) :: i,j,k,l
  Integer(4) :: iter
  Real(8) :: A(N,N), B(N,N), C(N,N)
  Real(8) :: ts, te
  Real(8) :: GFLOPs
  Character(256) char

  GFLOPs = dble(2*N**3 + 2*N**2)/1.0e9
  Print *, 'Total operations: ', GFLOPs ,'GFLOPs'
!Initialization                                                                                                                                                                                                                                                                                                                                                                    
  Call Random_Number(A)
  Call Random_Number(B)
  A = A - 0.5d0
  B = B - 0.5d0

!Hand-writing MM product                                                                                                                                                                                                                                                                                                                                                           
  char = "Hand-writting MM product"
  C = 0.d0
  Call Cpu_Time(ts)
  Do iter = 1, Niter
    C = 2.0d0*Hand_matmul(N,A,B) + 0.0d0*C
  End Do
  Call Cpu_Time(te)
  Call Write_time(char,ts,te,GFLOPs,Niter)

!BLAS                                                                                                                                                                                                                                                                                                                                                                              
  char = "Linked BLAS"
  C = 0.d0
  Call Cpu_Time(ts)
  Do iter = 1, Niter
    Call DGEMM('N','N',N,N,N,2.0d0,A,N,B,N,0.d0,C,N)
  End Do
  Call Cpu_Time(te)
  Call Write_time(char,ts,te,GFLOPs,Niter)

!Matmul                                                                                                                                                                                                                                                                                                                                                                            
  char = "The matmul function"
  C = 0.d0
  Call Cpu_Time(ts)
  Do iter = 1, Niter
    C = 2.0d0*Matmul(A,B) + 0.0d0*C
  End Do
  Call Cpu_Time(te)
  Call Write_time(char,ts,te,GFLOPs,Niter)

Contains
!==                                                                                                                                                                                                                                                                                                                                                                                
  Function Hand_matmul(N,A,B) result(C)
    Implicit none
    Integer(4), intent(in) :: N
    Real(8), intent(in) :: A(1:N,1:N), B(1:N,1:N)
    Real(8) :: C(1:N,1:N)
    Integer(4) :: i,j,k

    C(1:N,1:N) = 0.d0
    Do i = 1, N
      Do j = 1,N
        Do k = 1,N
          C(i,j) = C(i,j) + A(i,k)*B(k,j)
        End Do
      End Do
    End Do
    Return
  End Function Hand_matmul
!==                                                                                                                                                                                                                                                                                                                                                                                
  Subroutine Write_time(char,ts,te,GFLOPs,Niter)
    Implicit none
    Character(256), intent(in) :: char
    Real(8), intent(in) :: ts, te
    Real(8), intent(in) :: GFLOPs
    Integer(4), intent(in) :: Niter

    Print *,trim(char)
    Print *,'Elapse time =', (te-ts)/dble(Niter), ' sec'
    Print *,GFLOPs/(te-ts)*dble(Niter), ' GFLOPs/s'

  Return
  End Subroutine Write_time
!==                                                                                                                                                                                                                                                                                                                                                                                
End Program BLAS_Check
```
10回、倍精度実数の行列積を計算して、平均した計算速度が出力される
### 出力例

```bash
 Total operations:    2.0019999999999998      GFLOPs
 Hand-writting MM product
 Elapse time =   1.9024758000000002       sec
   1.0523129913137395       GFLOPs/s
 Linked BLAS
 Elapse time =  0.72820680000000026       sec
   2.7492190405253005       GFLOPs/s
 The matmul function
 Elapse time =   9.3727200000000011E-002  sec
   21.359861384955483       GFLOPs/s
```

## 結果

| BLASの種類                 |               速度 |    速度向上(パッケージ管理アプリを基準)    |
|:---------------------------|------------------:|:------------------:|
| パッケージ管理アプリのBLAS  |      2.84 GFLOPs/s|        1           |
| 自分でビルドしたBLAS        |      6.52 GFLOPs/s|       x 2.30       |
| OpenBLAS                   |      43.1 GFLOPs/s|       x 15.2       |
| ATLAS                      |      32.2 GFLOPs/s|       x 11.3       |
| MKL                        |      53.2 GFLOPs/s|       x 18.7       |
|---------------------------|------------------|------------------|
| 手で書いた行列行列積        |      1.18 GFLOPs/s|       x 0.42       |
| matmul関数                 |      22.1 GFLOPs/s|       x 7.78       |

上の五つがライブラリをリンクした際のBLASの測定値。下の二つは、コードをべた書きして行列行列積を書いたときの計測時間と、Fortranの組み込み関数matmulを使ったときの計算時間。
計算は全て`export OMP_NUM_THREADS=1`として、評価した(特にOpenBLASがスレッド並列化されているので)。MKLはsequential版をリンクしている。
MKLをリンクした時が一番速い。OpenBLAS, ATLASも中々速い。自分でnetlibからコードを取ってきてリンクすると、一応パッケージ管理アプリでいれたものよりは倍くらい速くなった。
あと、意外とmatmulが優秀。

### 注記
Intel(R) Core(TM) i5-8259Uは4コアなので、単一コアの理論性能は36.8 GFLOPS/sのはずなのだけれど、OpenBLASとMKLはこれを超えている。何か計算量の評価(FLOPsの値)とか間違っていないか、ちょっと不安。
一応、最大クロック3.8 GHzまで上がったとすると、60.8 GFLOPs/sまで出るかもしれないので、一応この値は超えていない。
あと、演算自体は同じなので、ライブラリ間の速度の比は間違えてない。