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を使っているとき、
gfortran -O3 main.f90 -lblas
これで大方通るはず。
通らない時は、BLASが置いてあるディレクトリにpathが通っていない。$ ldconfig -pとかで、一応候補にあるかどうか確認してみる。なぜだか、ライブラリのpathを知っていて、そこにはLD_LIBRARY_PATHを通していない場合は以下のようにして、リンクできる:
gfortran -O3 main.f90 -L/path/to/be/specified/ -lblas
自分でビルドしたBLASをリンク
BLAS (Basic Linear Algebra Subprograms)からコードを取ってきて自分でビルドして、リンクできる。ライブラリというよりも、分割コンパイルしたオブジェクトを自分のコードにくっつけると思った方が、感覚としては近いかも。自分が管理者ではない計算機の、どこにBLASがあるのかわからないとき、ひとまず自分でビルドしたものをリンクすれば、動作確認はできる。
例として、blas-3.8.0.tgzをダウンロードしてきたとする。
展開
ダウンロードしたコードを展開して、適当なディレクトリ(例えば'/home/user/sources/`とか)にうつす。
tar -zxvf blas-3.8.0.tgz
mv ./BLAS-3.8.0 /home/user/sources/
ビルド
展開したディレクトリに移って、
make
でおしまい。Makefileはmake.incをインクルードしていて、コンパイルのオプションとか出来たオブジェクトの名前とかを変えられる。デフォルトのままmakeすると静的ライブラリblas_LINUX.aが出来上がる。
所望の場所に配置
自分で使う事を念頭に、管理者権限の必要のない場所に置く。例えば/home/user/opt/BLAS-3.8.0/libを作って、その中に置く。
マイBLASをリンク
gfortran -O3 main.f90 /home/user/opt/BLAS-3.8.0/lib/blas_LINUX.a
でコンパイルすると、晴れてマイBLASがリンクできる。
OpenBLAS
OpenBLAS : An optimized BLAS library
大変、速いらしい。
展開＆ビルド
tar -zxvfで展開。普通にmakeでビルドが通る。数分ビルドにはかかるかも。
標準出力の内容をを見る限り、BLAS, CBLAS, LAPACK, LAPACKEがビルドされているみたい。
ビルド後にmake PREFIX=/path/to/be/installed installで、ライブラリ群が所望の場所にインストールされる。
リンクは
gfortran -O3 main.f90 -L/home/user/opt/OpenBLAS-0.3.7/lib -lopenblas
か
gfortran -O3 main.f90 /home/user/opt/OpenBLAS-0.3.7/lib/libopenblas.a -lpthread
で通る。
## トラブルシューティング
Pathの通しそびれ
./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をみると
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を追加すると、動く。
静的ライブラリをリンクするとライブラリ不足で怒られる
誤って
gfortran -O3 main.f90 /home/user/opt/OpenBLAS-0.3.7/lib/libopenblas.a
でコンパイルすると
/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のオプションをつけて
gfortran -fopenmp -O3 main.f90 /home/user/opt/OpenBLAS-0.3.7/lib/libopenblas.a
とすると一応コンパイルは通る。~~ただ、a.outを実行するとSegmentation fault (core dump)で落ちる。~~と思っていたら、知らないうちに動いていた。(なんでだ)
足りないライブラリをリンクすればよいので、
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)
アーキテクチャごとに最適化したBLASをインストールできる。
展開＆ビルド
Basic Steps of an ATLAS installを参考に。
tar xf atlas3.10.3.tar.bz2で展開。
出来たディレクトリATLASをATLAS3.10.3に変更
このディレクトリの中にディレクトリBLDdirを作成。(ソースのあるディレクトリでconfigureしようとすると、怒られる。)
CPUのthrottlingを止めるのはだるかったので、-cripple-atlas-performanceをつけてそのまま./configure.
../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)の文字列が見えたが、どうも共有ライブラリ関係がうまくいっていない様子(というか、そもそもビルドのところ共有ライブラリが出来ていないようにみえる)。
ひとまず、静的ライブラリのみ使うことにする。
リンクは
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からダウンロード。ダウンロードにあたって、名前やらメールアドレス、所属の登録が必要。登録したメールアドレスに*.licが届く。(が、このライセンスファイルを使わずインストールできてしまった)
展開＆インストール
tar -zxvf l_mkl_2020.0.166.tgzで展開。
出来たl_mkl_2020.0.166に入って、
sh ./install.sh
で、指示に従って入力していく。すると、/home/user/intelが作られて、この下にライブラリ群がインストールされている。
リンク
自分の環境に合わせて、Intel® Math Kernel Library Link Line Advisor | Intel® Softwareでリンク。
インテルコア、静的リンク、32-bit 整数、スレッド並列なしだと
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の速度比較を行ってみる。
他の環境で、調べていただいたものがコメント欄にあるので、そちらも参照していただくと良いと思います。
使うコード
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回、倍精度実数の行列積を計算して、平均した計算速度が出力される
出力例
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まで出るかもしれないので、一応この値は超えていない。
あと、演算自体は同じなので、ライブラリ間の速度の比は間違えてない。