BLAS とは
BLAS
は様々な人，企業が開発している:
・
Reference BLAS
基準として作られた
BLAS(
無料)
．速くない．・
Intel MKL
Intel
社が開発(
有料)
．CPU
毎に設計されて非常に速い．・
OpenBLAS
後藤和茂先生が作成した
GotoBLAS
が 引き継がれたBLAS
．速い(
無料)
．・
ATLAS (Automatically Tuned Linear Algebra Software)
自動でチューニングするBLAS
．CPU
に依存しない．46BLAS とは
BLAS :
関数名，サブルーチン名や引数，役割が同じで あるため，どの
BLAS
を使っても動作は同じ．しかし，
BLAS
によって速度が変わる!!
① 計算時間がかかるところを
BLAS
で作成② 使用者が
BLAS
を選択し，実行するコンピュータ毎に最適な選択し，高速に
!
47
BLAS とは
BLAS
は演算ごとにレベル分けされている: Level 1:
ベクトル
-
ベクトルの演算Level 2:
行列
-
ベクトルの演算Level 3:
行列
-
行列の演算48
Lapack とは
Lapack:
Linear Algebra PACKage
の略 線形代数ライブラリ．連立一次方程式や固有値問題などが解ける．
例えば
dgesv(An,Bn,An,ipiv,b,bn,info)
連立一次方程式のサブルーチンAx = b
の解
x
をb
に代入して出力される．49
Lapack とは
Lapack:
Lapack
はFortran
で記述されており，内部でBLAS
を用い ているためBLAS
を差し替えることでCPU
に依存した最 適化が可能!
世界中で利用されており，信頼性も高い
!!
現在はバージョン
3.5.0 Lapack
のリファレンス:
http://www.netlib.org/lapack/
50
Lapack とは
Lapack
は変数の型ごとに関数，サブルーチン名が変わる:
ge :
一般行列gb :
一般帯行列tr :
三角行列 など…
例
:
dgesv (
引き数) dgesv, dgbsv
51
連立一次方程式
☆文法
?gesv(n, nrhs, A, lda, ipiv, b, ldb, info) n : integer
型．A
の次元(n × n), b
の行数．nrhs : integer
型．b
の列数．A : ?
型のn × n
配列．lda : max(1,n)
b : ?
型のn × nrhs
配列ipiv : integer
型．n
次元の配列．info : integer
型．info = 0
なら正常info > 0
正則でない可能性info < 0 info
番目の値が不正 52Matlab
BLAS
やLapack
などのライブラリを使えば世界最高速のツールが使えるが
…
もっと手軽に世界最高のツールを利用したい
!!
Matlab を利用する !!
Matlab
はIntel MKL
をはじめとした世界の様々なライブラリを搭載したプログラミング言語．
そのため，世界最高峰の行列積や連立一次方程式 の近似解法が簡単に使用可能
!!
53
Matlab
端末上で
matlab
とコマンドをうつ．
と
Matlab
コマンドをうつ．Matlab が起動
54
>> n = 5
>> A = ones(n,n);
>> B = ones(n,n);
>> C = A*B;
Matlab
端末上で
matlab
とコマンドをうつ．
Matlab が起動
55
>> n = 5
>> A = ones(n,n);
>> B = ones(n,n);
>> C = A*B;
n
に5
を代入全要素
1
のn × n
行列A
とB
を作成 セミコロン(;)
があると非表示になる．
Matlab
端末上で
matlab
とコマンドをうつ．
Matlab が起動
56
>> n = 5
>> A = ones(n,n);
>> B = ones(n,n);
>> C = A*B;
A × B
の行列積を実行．Intel MKL
のBLAS
dgemm
が呼ばれる!!
Matlab
57
>> n = 1000, A = ones(n,n); B = ones(n,n);
>> tic; C = A*B; toc
Matlab
58
>> n = 1000, A = ones(n,n); B = ones(n,n);
>> tic; C = A*B; toc
一行で書くことも可 表示する場合はカンマで区切る
実行時間の計測
自作行列積と実行時間を比べよう !
Matlab
59
>> n = 10, A = rand(n,n); b = A*ones(n,1);
>> tic; x = A¥b, toc
Matlab
60
>> n = 10, A = rand(n,n); b = A*ones(n,1);
>> tic; x = A¥b, toc
n × n
の乱数行列A
を作成A
と全要素1
のベクトルの行列ベクトル積 連立一次方程式Ax = b
を満たすx
を求める．Lapack
の連立一次方程式を求める関数が呼び出される．
答えは
b = A*ones(n,1)
であるため，大体
x=ones(n,1)
になる．Matlab
61
>> n = 100, A = rand(n,n); A = A + A’;
>> tic; lambda = eig(A), toc
Matlab
62
>> n = 100, A = rand(n,n); A = A + A’;
>> tic; lambda = eig(A), toc
A’
はA
の転置行列．A+A’
は対称行列になる行列
A
の固有値問題A*x = lambda*x
を満た す固有値lambda
を求める関数．これも
Lapack
の固有値問題の関数が呼ばれ ている(
対称行列用)
．A
が対称行列のため，固有値はすべて実数 になる．おまけ : 精度保証付き数値計算
63
今まで利用してきた浮動小数点数とは簡単に言うと 超高速な近似計算である．そのため
を計算した場合，実に
10000000000
回以上の近似計算を 行う．このとき，正しい結果は得られているのであろうか?
とみると大体
1
に近いためあっていることが推測される．>> n = 3000, A = rand(n,n); b =A*ones(n,1);
>> tic; x = A¥b; toc
>> format long
>> x(1)