ＬＡＰＡＣＫとは？

LAPACK(Linear Algebra PACKage)は，netlib で公開されている線形計算ライブラリです．各ルーチンはFORTRAN 77で記述されています．

LAPACK には，約300種類，各精度を合計すると約1,100本のルーチンが用意されています．主な機能は以下の通りです:
LAPACK は，一般に問題を解くための機能を提供するドライバルーチン(driverroutines)， 個々の問題を解くための機能を提供する計算ルーチン(comutationalroutines)，および， 補助的な計算や共通に使用される手続きを提供する補助ルーチン(auxiliary routines)で構成されています．

LAPACK は，内部で BLAS (Basic Linear Algebra Subprograms) ライブラリを呼び出します． LAPACK での実際の計算は BLAS を用いて行われています．

BLAS には様々な種類がありますが，各CPUベンダーから最適化されたBLASが提供されています．たとえば，Intel Pentium の Linux 用 BLAS はここから手にはいります．最適化された BLAS を使うだけで，演算速度が飛躍的に効率化します．

上記 BLAS の他に，各システムに最適化した BLAS を生成する ATLAS というソフトもあります．

CLAPACK とは？

CLAPACK とは，FORTRAN77 で記述されている LAPACK を f2c（FORTRAN のソースを C のソースに変換するソフト） を用いることにより C のソースから利用できるようにしたライブラリです．

CLAPACK のインストール

ここからclapack.tgz をダウンロードする。
適当なディレクトリで

% tar zxvf clapack.tgz
% cd CLAPACK
% cp -f INSTALL/make.inc.LINUX ./make.inc
% make f2clib
% make blaslib
% cd SRC
% make
% cd ../TESTING/MATGEN/
% make
% cd ../../

ルートになって

cp lapack_LINUX.a /usr/local/lib/liblapack.a
cp blas_LINUX.a /usr/local/lib/libblas.a
cp F2CLIBS/libF77.a /usr/local/lib
/sbin/ldconfig

これでlapackのライブラリを使えるようになる。

これを sample.c として保存し、

gcc sample.c -llapack -lblas -lF77

としてコンパイル実行できれば成功！
答えは全て1.0になるはず。

性能評価

CLAPACK の性能評価結果を以下に示します．比較対象としては，Ax = b という連立一次方程式を A を 1000 x 1000 の full matrix として

1) ガウスの消去法を用いて計算する
2) lapack の dgesv 関数を用いて計算する

の２通りで計算しました．

使用したマシンは

CPU : PentiumIII 700MHz
MEM : 320 MB

です．

計算に使用したソースを以下からDLできます．

random.c - a.dat(行列A), b.dat(ベクトルb)を生成する
gaus.c - ガウスの消去法
lapack.c - CLAPACK の dgesv 関数使用

使い方としては，random.c をまず実行し，a.dat, b.dat を生成してから，gaus.c or lapack.c を実行してください．こちらも答えはすべて 1.0 になるはずです。

計算結果

ガウスの消去法 : １８秒
CLAPACK : １４秒

30 % の性能アップです．

上記性能比較は最適化されていない BLAS(CLAPACK をコンパイルした時生成される)を用いた計算なので，最適化された BLAS を用いればより効率的な計算が行えると思います．

注意点

CLAPCK は f2c を用いて FORTRAN のソースを（ある意味強引に）C ソースに変換したものであり，すべて FORTRAN の規則がベースとなります．その為，以下のような注意が必要です．

1) 関数名の後に「_」がつく．

これは f2c の慣習で，例えば FORTRAN で「DGESV」という関数は「dgesv_」という関数となります．

2) 変数をポインタで関数へ渡す必要がある．

FORTRAN の規則により，関数を呼び出すときすべての引数をポインタとして渡さなければなりません．このため，関数の引数に直接数値を指定することはできません．たとえば，

FORTRAN:
call dgetrf(5, 5, A, 5, ipiv, info)

C:
N = M = LDA = 5;
dgetrf_(&M, &N, A, &LDA, ipiv, &info);

3) 文字列引数はとれない．

関数の引数が文字列の場合でも，その文字列の最初の文字だけが意味を持ちます．たとえば，

FORTRAN:
call dpotrf('Upper', n, a, lda, info)

C:
char s='U';
dpotrf_(&s, &n, a, &lda, &info);

4) 二次元配列を用いるときには要注意！！

FORTRAN と C では二次元配列の定義に違いがあります．

DOUBLE PRECISION A(LDA, N)

と宣言された二次元 FORTRAN 配列は LDA x N の倍語のメモリを取り列方式に従って格納されます．つまり，列内の要素は連続し，行内の要素はLDA 倍語の幅だけ離れます．しかし C 言語では二次元配列は行方式に従います．さらに，二次元の C 配列の行は必ずしも連続である必要はありません．

double A[LDA][N];

は長さ N の行への LDA ポインタを持ちます．これらのポインタは原則的にはメモリのどこにあってもいいものです．このような二次元 C 配列を CLAPACK ルーチンに引き渡すと必ずといっていいほど間違った結果をもたらします．

対応策として，大きさ LDA x N 倍語の一次元 C 配列を用いなければなりません．（あるいは同じ量の空間を malloc します）
CLAPACK が期待する配列を得るためには，次の方法を使うのがよいでしょう．

double*A;
A = malloc(LDA*N*sizeof(double));

一例として Ax = b の行列計算を一次元配列を用いて計算することを考えます．
A が 3x3 の以下のような行列であるとします．

A =
1 2 3
4 5 6
7 8 9

この場合，CLAPACK へ渡す一次配列へ格納する順番としては 1 4 7 2 5 8 3 6 9 の順番に格納していく必要があります．

参考文献

「Linux 数値計算ツール」
大石 進一
コロナ社


早稲田大学大石先生のHP
LAPACKの公式ホームページ
LAPACK FAQ

戻る