CULA (R14, CUDA 4.4) - 超音波流体屋のプログラム備忘録

CUDA 4.1, CULA Dense R14 版の(古い)情報です。
→ CULA Dense (R17, CUDA 5.5)
→ CULA Sparse (S5, CUDA 5.5)

CULA free editionならfloat, float complex限定、連立方程式・最小二乗法・LU分解あたりが使える。
ダウンロードは以下からID登録してから。
http://www.culatools.com/downloads/dense/
linuxは落としてきた.runを管理者権限でbashで実行。
今回はCULA Dense R14 free editoinを使用。
要CUDA 4.1だったので、CUDAも4.1にバージョンアップ。

ソース

CULA使い方の一例として、連立方程式 Ax = bを解く話。
んで今回はMatlab記法で A=ones(10 000, 10 000)-eye(10 000); b=ones(10 000, 1);
対角成分だけ0、後は1の行列。解は確か全部 1/9999。

#include <stdio.h>
#include <stdlib.h>
 
//cuda, cula関係のヘッダ
#include <cutil.h>
#include <cula.h>
#include <cula_lapack.h>
#include <cula_lapack_device.h>
 
int main(int argc, char** argv){
 
	CUT_DEVICE_INIT(argc, argv);
	culaStatus s; //失敗するとculaNoError以外が入るらしい
	s = culaInitialize();
 
	//ホスト領域の確保・値の代入
	float *a, *b;		
	int n=10000; 
	int m=10000;
	a = (float *)calloc(sizeof(float), m*n);
	b = (float *)calloc(sizeof(float), m);
	//やってみたらいわゆるColumn-major, CLapackと同じらしい
	int i, j;
	for(i=0; i<n; i++){
	for(j=0; j<m; j++){
		a[i+j*n]= (float) (i != j);	
	}
		b[i] = 1;	
	}
	//その他、行数・列数情報。CLapackと同じ。
	int ipiv[n];
	int nrhs = 1;
	int lda = n; 
	int ldb = n; 
 
	//時間計測用タイマー
	unsigned int timer = 0;
	CUT_SAFE_CALL( cutCreateTimer( &timer));
	CUT_SAFE_CALL( cutStartTimer( timer) );
 
	//連立方程式を解く。これ１行。
	s = culaSgesv(n, nrhs,  a, lda, ipiv, b, ldb);
 
	//タイマーストップ
	CUT_SAFE_CALL( cutStopTimer( timer));
	printf("Elapsed time : %f (ms)\n", cutGetTimerValue(timer));	 
 
	//エラー処理等したければ。 正常終了なら0が返る。 
	printf("culastate %d noerr:%d \n", s, culaNoError); 
 
	//結果表示
	for(j=0; j<m; j++) printf("%f ", b[j]);
 
	//終了処理
	free(a);
	free(b);
	culaShutdown();
	CUT_SAFE_CALL( cutDeleteTimer( timer));
	CUT_EXIT(argc, argv);
	return 0;
}
 

コンパイル

Makefileはこんな感じ。環境変数を工夫すればもっと短くてすむはず。

CC=nvcc
LIB=-lcutil_x86_64 -lcula_core -lcula_lapack
SRC=sample.cu
OUT=sample.out
LIB1=-L~/NVIDIA_GPU_Computing_SDK/C/lib
LIB2=-L/usr/local/cula/lib64
INC1=-I~/NVIDIA_GPU_Computing_SDK/CUDALibraries/common/inc
INC2=-I/usr/local/cula/include/
all:
	${CC} ${SRC} -o ${OUT} ${LIB} ${LIB1} ${LIB2} ${INC1} ${INC2}
 

計算時間

NVS 4200M のthinkpad T420sで13.4 sec。
CPU: Core i7-2640M 2.8GHz (smpとか無しで)
octaveで 109sec, clapackのdgesv_で496sec。
だいたい8倍くらいは速いかも。

あらかじめGPUに格納しときたいときは

culaDeviceSgesvほかdevice系の関数を使う。
行列Aとベクトルbに対応する領域、加えて関数に必要な配列ipivなども忘れずにデバイス上に確保してから、
同様の関数を採用すればよい。上記42行付近は以下の通りに書き換える。

{
 	double *da, *db;
 	int *dipiv;
 	CUDA_SAFE_CALL(cudaMalloc((void**) &da, sizeof(double)*n*n));
	CUDA_SAFE_CALL(cudaMalloc((void**) &db, sizeof(double)*n));
	CUDA_SAFE_CALL(cudaMalloc((void**) &dipiv, sizeof(int)*n))
 
	CUDA_SAFE_CALL(cudaMemcpy(da, a, sizeof(double)*n*n, cudaMemcpyHostToDevice));
	CUDA_SAFE_CALL(cudaMemcpy(db, b, sizeof(double)*n, cudaMemcpyHostToDevice));
 
	s = culaDeviceSgesv(n, nrhs,  da, lda, dipiv, db, ldb);
}
 

精度

単精度のせいか、大行列になると荒れる。
上記の例ですでにばらつきが大きい。
n=5000ぐらいが限界だろうか。
それ以上は倍精度対応版を買うしかないと思われる。

まあDense/Sparseともアカデミックなら$95、商用でも$395/年。
(買って落とした分はずっと使っていいけど、期限後のバージョンアップには更新がいるみたい)
年数十〜数百万ふんだくるどこぞの解析ソフトよかはお得じゃないかしら。

参考
(pdf)CULA Reference Guide
(pdf)CULA API Guide