CUDAで行列演算：行列式
をテンプレートにして作成
[
トップ
|
新規
|
一覧
|
検索
|
最終更新
|
ヘルプ
]
開始行:
*CUDAで行列演算 その4 行列式 [#he65117c]
行列の加減算や乗算などは要素ごとに計算を分解することが用...
**ガウスの消去法 [#j0864ff0]
一般的な行列式の計算方法は，第j列についての余因子展開を用...
#ref(eq_det1.jpg,nolink)
となります．M&size(10){ij};は余因子行列です．
これをプログラムで実装する場合，再帰的な処理を用いて行う...
しかし，CUDAカーネルでは再帰的な関数呼び出しは許可されて...
そのため，ガウスの消去法を用いてみます．
ガウスの消去法は逆行列を求めるための手法で，前進消去によ...
#ref(eq_gauss.jpg,nolink)
をk=1〜n-1, i=k+1〜n, j=k〜n について繰り返すことで上三角...
**CUDAでの実装 [#jf63c9fd]
***メモリ確保とデータ転送 [#kc6b68fc]
行列の定義は[[CUDAで行列演算：乗算]]と同じものを用います...
#code(C){{
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で行列演算：加減算]]のものと同じで，テ...
デバイス側のメモリ確保とデータ転送も同様に
#code(C){{
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,...
}}
***カーネル [#e0254807]
カーネルの実行時にシェアードメモリの大きさを指定します．
これは，アクセスの多い k行目のn要素をシェアードメモリに格...
今回は行列の大きさをスレッド数(THREAD_NUM=512)までとして...
1ブロックのシェアードメモリにn要素全てを格納していますが，
これを超える場合はもう少し工夫が必要であると思います．
#code(C){{
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個...
カーネルの実装は，
#code(C){{
__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();
}
}
}
}}
シェアードメモリの動的な確保については，[[シェアードメモ...
スレッドごとにループ内をほとんどスルーする場合が多く，あ...
**ホストへの結果の転送と後片付け [#idfeb3fb]
結果をホストに転送し，確保したメモリを解放します．
#code(C){{
// デバイスからホストへ結果を転送
size = dA.width*dA.height*sizeof(float);
cutilSafeCall(cudaMemcpy(hU.elements, dA.elements, size,...
// デバイスメモリ解放
cutilSafeCall(cudaFree(dA.elements));
}}
そして得られた上三角行列から
#code(C){{
float detA = 1.0;
for(int i = 0; i < hU.width; ++i){
detA *= hU.elements[i*hU.width+i];
}
}}
行列式の値 detA が得られます．
**計算時間 [#dba92d65]
CPU用にもガウスの消去法を用いた行列式の計算プログラムを実...
計算時間と結果を比較しました．カーネル実装上の制限から計...
|GPU|28.04|
|CPU| 4.60|
時間の単位は ms(ミリ秒)です．一回計測しただけです．
やはりあまりよくないです．ただもう少し最適化できそうです...
終了行:
*CUDAで行列演算 その4 行列式 [#he65117c]
行列の加減算や乗算などは要素ごとに計算を分解することが用...
**ガウスの消去法 [#j0864ff0]
一般的な行列式の計算方法は，第j列についての余因子展開を用...
#ref(eq_det1.jpg,nolink)
となります．M&size(10){ij};は余因子行列です．
これをプログラムで実装する場合，再帰的な処理を用いて行う...
しかし，CUDAカーネルでは再帰的な関数呼び出しは許可されて...
そのため，ガウスの消去法を用いてみます．
ガウスの消去法は逆行列を求めるための手法で，前進消去によ...
#ref(eq_gauss.jpg,nolink)
をk=1〜n-1, i=k+1〜n, j=k〜n について繰り返すことで上三角...
**CUDAでの実装 [#jf63c9fd]
***メモリ確保とデータ転送 [#kc6b68fc]
行列の定義は[[CUDAで行列演算：乗算]]と同じものを用います...
#code(C){{
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で行列演算：加減算]]のものと同じで，テ...
デバイス側のメモリ確保とデータ転送も同様に
#code(C){{
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,...
}}
***カーネル [#e0254807]
カーネルの実行時にシェアードメモリの大きさを指定します．
これは，アクセスの多い k行目のn要素をシェアードメモリに格...
今回は行列の大きさをスレッド数(THREAD_NUM=512)までとして...
1ブロックのシェアードメモリにn要素全てを格納していますが，
これを超える場合はもう少し工夫が必要であると思います．
#code(C){{
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個...
カーネルの実装は，
#code(C){{
__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();
}
}
}
}}
シェアードメモリの動的な確保については，[[シェアードメモ...
スレッドごとにループ内をほとんどスルーする場合が多く，あ...
**ホストへの結果の転送と後片付け [#idfeb3fb]
結果をホストに転送し，確保したメモリを解放します．
#code(C){{
// デバイスからホストへ結果を転送
size = dA.width*dA.height*sizeof(float);
cutilSafeCall(cudaMemcpy(hU.elements, dA.elements, size,...
// デバイスメモリ解放
cutilSafeCall(cudaFree(dA.elements));
}}
そして得られた上三角行列から
#code(C){{
float detA = 1.0;
for(int i = 0; i < hU.width; ++i){
detA *= hU.elements[i*hU.width+i];
}
}}
行列式の値 detA が得られます．
**計算時間 [#dba92d65]
CPU用にもガウスの消去法を用いた行列式の計算プログラムを実...
計算時間と結果を比較しました．カーネル実装上の制限から計...
|GPU|28.04|
|CPU| 4.60|
時間の単位は ms(ミリ秒)です．一回計測しただけです．
やはりあまりよくないです．ただもう少し最適化できそうです...
ページ名:
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
