CUDAで行列演算：行列式

CUDAで行列演算 その4 行列式†

行列の加減算や乗算などは要素ごとに計算を分解することが用意で，GPUでの並列計算に適しているといえます．しかし，行列式の計算は処理が直列につながるようなものが多く，あまり計算の分解が容易でないため，GPUでの並列化には向いていないと思いますが，一応実装してみます．

ガウスの消去法†

一般的な行列式の計算方法は，第j列についての余因子展開を用いて，

となります．Mijは余因子行列です． これをプログラムで実装する場合，再帰的な処理を用いて行うと簡単に実装できます． しかし，CUDAカーネルでは再帰的な関数呼び出しは許可されていません． そのため，ガウスの消去法を用いてみます．

ガウスの消去法は逆行列を求めるための手法で，前進消去によって上三角行列の算出し，後退代入で逆行列を求める直接解法です．この前進消去によって求められた上三角行列の対角成分の積が行列式の値になります．行列A(n×n)の要素をaijとすると，

をk=1〜n-1, i=k+1〜n, j=k〜n について繰り返すことで上三角行列が得られます．

CUDAでの実装†

メモリ確保とデータ転送†

行列の定義はCUDAで行列演算：乗算と同じものを用います．今回は入力行列は1つです．

Matrix hA, hU;
	hA.width = hA.height = 16;
	hU.width = hU.height = hA.width;

	hA.elements = new float[hA.width*hA.height];
	hU.elements = new float[hU.width*hU.height];

	RandomInit(hA.elements, hA.width*hA.height);

hUは結果の上三角行列を格納するためのものです． RandomInitはCUDAで行列演算：加減算のものと同じで，テスト用に乱数を行列の要素に代入する関数です．

デバイス側のメモリ確保とデータ転送も同様に

Matrix dA;
	dA.width = hA.width;	dA.height = hA.height;

	int size;
	// デバイスメモリの確保とホストからの転送
	size = dA.width*dA.height*sizeof(float);
	cutilSafeCall(cudaMalloc((void**)&dA.elements, size));
	cutilSafeCall(cudaMemcpy(dA.elements, hA.elements, size, cudaMemcpyHostToDevice));

カーネル†

カーネルの実行時にシェアードメモリの大きさを指定します． これは，アクセスの多い k行目のn要素をシェアードメモリに格納するためです． 今回は行列の大きさをスレッド数(THREAD_NUM=512)までとしているため， 1ブロックのシェアードメモリにn要素全てを格納していますが， これを超える場合はもう少し工夫が必要であると思います．

dim3 block(THREAD_NUM, 1, 1);
	dim3 grid((dA.height+block.x-1)/block.x, 1, 1);

	matrixGE<<< grid, block, sizeof(float)*dA.width >>>(dA);

	// カーネル実行エラーのチェック
	cutilCheckMsg("Kernel execution failed");

	// GPUスレッドが終わるのを待つ
	cutilSafeCall(cudaThreadSynchronize());

カーネル名はmatrixGEとしています． 並列に計算できるi=k+1〜nで分解するため，カーネルは最大n個必要です．

カーネルの実装は，

__global__ 
void matrixGE(Matrix M)
{
	extern __shared__ float Mk[];

	int i = blockIdx.x*blockDim.x+threadIdx.x;

	if(i < M.height){
		int n = M.width;

		for(int k = 0; k < n-1; ++k){
			Mk[i] = M.elements[k*n+i];
			__syncthreads();

			if(i >= k+1){
				float tmp = M.elements[i*n+k]/Mk[k];

				for(int j = k; j < n; ++j){
					M.elements[i*n+j] -= tmp*Mk[j];
				}
			}

			__syncthreads();
		}

	}
}

シェアードメモリの動的な確保については，シェアードメモリを参照． スレッドごとにループ内をほとんどスルーする場合が多く，あまり効率がよくない気がしますがとりあえずの実装ということで．

ホストへの結果の転送と後片付け†

結果をホストに転送し，確保したメモリを解放します．

// デバイスからホストへ結果を転送
	size = dA.width*dA.height*sizeof(float);
	cutilSafeCall(cudaMemcpy(hU.elements, dA.elements, size, cudaMemcpyDeviceToHost));

	// デバイスメモリ解放
	cutilSafeCall(cudaFree(dA.elements));

そして得られた上三角行列から

float detA = 1.0;
	for(int i = 0; i < hU.width; ++i){
		detA *= hU.elements[i*hU.width+i];
	}

行列式の値 detA が得られます．

計算時間†

CPU用にもガウスの消去法を用いた行列式の計算プログラムを実装し， 計算時間と結果を比較しました．カーネル実装上の制限から計測に使った行列の大きさはn=256としました．

時間の単位は ms(ミリ秒)です．一回計測しただけです． やはりあまりよくないです．ただもう少し最適化できそうですが．