CLAPACKをVisual Studio 2010上で使う(32/64bit) (高速にLU分解がしたいので)　2012/05/18--

  • サーベイ
+ ATLASは，64bitでコンパイルするのが難しいらしい
+ GOTOBLAS + LAPACKは，Core i7(Sandy Bridge CPU)ではうまくいかない(?)。解決できるんだろうけど時間がかかりそう
+ OpenBLAS + LAPACKは，MinGW+Cygwin環境でコンパイルできた．が、headerでVisual Sutio 2010にはもう入っていないcomplex.hを要求されたりしてちょっと使い難い．

+よってここのCLAPACKを利用する．

  • *.h と *.libをダウンロード
+ここからf2c.hとclapack.hをもらってくる
+ここからコンパイル済みの*.libをもらってくる(自分でコンパイルした方が速そうだけど・・・上記のサーベイに丸1日かかってもう疲れた)
+Visual Studioのプロジェクトから見えるところに置いておく

  • ヘッダをインクルードする 
 extern "C"
{
#include "f2c.h"
#include "clapack.h"
}

#pragma comment(lib, "BLAS_nowrap.lib")
#pragma comment(lib, "libf2c.lib")
#pragma comment(lib, "clapack_nowrap.lib")

  • コードを書く (LU分解) 
void test()
{
　　　　integer N = 4;
　　　　integer INFO ; 
　　　　long iPiv[4 ];
　　　　double a[16]={
　　　　　　　1.0,　　5.0,　　9.0,　　4.0, 
　　　　　　　2.0,　11.0,　　7.0,　14.0, 
　　　　　　　3.0,　10.0,　　6.0,　15.0, 
　　　　　　13.0,　　2.0,　12.0,　　1.0　}; //行列はColumn majourなのに注意

　　　　dgetrf_( &N, &N, a, &N, iPiv, &INFO );

　　　　for (int y= 0; y< 4; ++y ){ 
　　　　　　　　for (int x= 0; x< 4; ++x ) printf( "a[%d][%d]=%f ",y,x, a[ 4 * x + y] );
　　　　　　　　printf( "\n ");
　　　　}

　　　　for (int y= 0; y< 4; ++y )　printf( "iPiv[%d]=%d \n", y, iPiv[y]);
}

  • 結果
\\aにLU分解結果が入る
a[0][0]=9.000000 a[0][1]=7.000000 　a[0][2]=　6.000000 a[0][3]=12.000000
a[1][0]=0.444444 a[1][1]=10.888889 a[1][2]=12.333333 a[1][3]=-4.333333
a[2][0]=0.555556 a[2][1]=0.653061 　a[2][2]=-1.387755 a[2][3]=-1.836735
a[3][0]=0.111111 a[3][1]=0.112245 　a[3][2]=-0.683824 a[3][3]=10.897059

//iPiv[i]には、i列目を分解する時に選択したピボット(行番号)が入る．
iPiv[0]=3
iPiv[1]=4
iPiv[2]=4
iPiv[3]=4

//これはつまり
1 2 3 4
↓(iPiv[0]=3)
3 2 1 4
↓(iPiv[1]=4)
3 4 1 2
↓(iPiv[2]=4)
3 4 2 1
という順序に列が並んでいるという事

  • 結局

このサイトのコンパイル済みの.libで，3000×3000程度の密行列のLU分解に30秒強の時間がかかった．
ニューメリカルレシピに乗ってたコードをOpenMPを使ってちょっとチューニングした自分のコードの方が速くて残念．
OpenMPの製作者ががんばって開発を進めているようなので、そちらに期待しつつ、自分のLU分解を使うことにする．

- (追記)BundlerのソースコードにCLapack + BLASをVisual Studioのプロジェクト化したものがあった。これを利用したら速いバージョンのCLAPACKが得られるかも。