hoge diary - November 5, 2007

まあ，わざわざ書くまでもないことかも知れないのですが．

3x3くらいの逆行列だったらライブラリ使うより自分で計算した方が速かった，という話です．

対戦相手は CLAPACK + CPPLapack の組み合わせ．

使った CLAPACK は Version 3.0 で，ATLAS は使っていません．CPPLapack は 2007/10/28 に svn checkout したものです (2005_03_25 版だと下のコンパイルエラーが出るので，当時の最新版をチェックアウトして使いました．そしてごめんなさい．チェックアウトしたリビジョンを忘れてしまいました)．

dgbmatrix-misc.hpp(176) : error C2061: 構文エラー : 識別子 'A'
dgbmatrix-misc.hpp(178) : error C2440: '=' : 'double **(__cdecl *)(void)' から 'double **' に変換できません。この変換が可能なコンテキストはありません。

話が脱線しますが，2005_03_25版で，上のエラーが出ている行を見ると，下のように書いてあります．

double** A_darray(A.Darray);

これは double** の値を初期化するつもりで書いていますが，コンパイラは関数宣言に見えるものは全て関数宣言として扱うので，これを関数宣言として扱おうとしてコンパイルエラー，というわけです．

最新版では次のように書いてありますので問題ありません．

double** A_darray = A.Darray;

上の CPPLapack のエラーに関する話は，Belution.com VC++ フォーラムのどこかに書いてあったものです．ですが，残念なことに今 Belution.com のサーバがダウンしているようで全くフォーラムを見られません．

閑話休題．まずは CLAPACK+CPPLapack を用いた場合の逆行列計算関数(Inverse1)のソースを載せます．

CPPL::dgematrix Inverse1( const CPPL::dgematrix& mat )
{
  return CPPL::i(mat);
}

まあ，これは CPPLapack の関数を呼んでいるだけです (しかし... 逆行列の関数名が i とは...)．

続いて，手計算バージョン(Inverse2)のソース．

CPPL::dgematrix Inverse2( const CPPL::dgematrix& mat )
{
  CPPL::dgematrix ret(3, 3);

  double a = mat(0, 0);
  double b = mat(0, 1);
  double c = mat(0, 2);
  double d = mat(1, 0);
  double e = mat(1, 1);
  double f = mat(1, 2);
  double g = mat(2, 0);
  double h = mat(2, 1);
  double i = mat(2, 2);

  double det = (a*e*i + b*f*g + c*d*h) - (a*f*h + b*d*i + c*e*g);

  ret(0, 0) = (e*i-f*h)/det;
  ret(0, 1) = (h*c-i*b)/det;
  ret(0, 2) = (b*f-c*e)/det;
  ret(1, 0) = (f*g-d*i)/det;
  ret(1, 1) = (i*a-g*c)/det;
  ret(1, 2) = (c*d-a*f)/det;
  ret(2, 0) = (d*h-e*g)/det;
  ret(2, 1) = (g*b-h*a)/det;
  ret(2, 2) = (a*e-b*d)/det;

  return ret;
}

それぞれを 1000000 回ずつ呼び出して，実行時間を計測しました．WinXP(SP2)+VC8(Express)です．また，時間計測には高分解能パフォーマンスカウンタ(QueryPerformanceCounter() API)を使いました．

その時間計測結果がこちら．

関数名	時間[秒]
Inverse1	1.403501
Inverse2	0.400147

ええ，3x3くらいなら自前で計算した方が速い，というわけです．

そして最後にオチ．「逆行列を100万回計算してもたった 1.4 秒しか掛かってないので，こんなところを速くしても全体は速くならないっ!」(アムダールの法則)

以上．おあとがよろしいようで．