ヤコビ法2011年1月8日 | ||||||||||||||||||||
ファイル内容ヤコビ法を GPU を使って実装してみましょう。 ベースの C プログラムjacobi.cu
void calcResidual(int n, const float *a, const float *b, const float *x, float *r)
{
int i, j;
for(i = 0; i < n; i++){
float v;
v = 0.;
for(j = 0; j < n; j++){
v += a[i*n + j]*x[j];
}
r[i] = v;
}
for(i = 0; i < n; i++) r[i] = b[i] - r[i];
}
int jacobi(int n, const float *a, const float *b, float *x, float tol,
int maxIter, float *resFinal, int *iterFinal)
{
...
bNorm = calcNorm(n, b);
if(bNorm == 0.) bNorm = 1.;
for(iter = 0; iter < maxIter; iter++){
for(i = 0; i < n; i++) x0[i] = x[i];
for(i = 0; i < n; i++){
int j;
float v;
v = 0.;
for(j = 0; j < n; j++){
if(i != j) v += a[i*n + j]*x0[j];
}
x[i] = ALPHA*(b[i] - v)/a[i*n + i] + (1. - ALPHA)*x0[i];
}
calcResidual(n, a, b, x, r);
rNorm = calcNorm(n, r);
res = rNorm/bNorm;
if(res < tol) break;
}
...
return 0;
}
上ではヤコビ法の部分だけ示しています。プログラムは以下のような書式の matrix.dat ファイルを読み込んで計算します。 matrix.dat 3 4. 1. 1. 2. 6. 1. 1. 2. 5. 4. -8. 7. はじめの数字は行列のサイズです。上の場合は 3×3 の行列を意味します。その下では 3×3 の行列の要素を記述しています。その下では右辺ベクトルの要素を記述しています。 CUDA 用のプログラムjacobix2.cu
float calcNorm(int n, const float *v)
{
float vNorm;
int i;
vNorm = 0.;
for(i = 0; i < n; i++) vNorm += v[i]*v[i];
vNorm = sqrt(vNorm);
return vNorm;
}
__global__ void subVector(int n, const float *a, const float*b, float *c)
{
int i = blockDim.x*blockIdx.x + threadIdx.x;
c[i] = a[i] - b[i];
}
__global__ void multMatrixVector(int n, const float *a, const float *b, float *c)
{
int i = blockDim.x*blockIdx.x + threadIdx.x;
int j;
float v;
v = 0.;
for(j = 0; j < n; j++){
v += a[i*n + j]*b[j];
}
c[i] = v;
}
void calcResidual(int n, const float *a, const float *b, const float *x, float *r)
{
int grid;
grid = (int)ceil((float)n/BLOCK_SIZE);
multMatrixVector<<<grid, BLOCK_SIZE>>>(n, a, x, r);
subVector<<<grid, BLOCK_SIZE>>>(n, b, r, r);
}
__global__ void jacobiInternal(int n, const float *a, const float *b, const float *x0, float *x)
{
int i = blockDim.x*blockIdx.x + threadIdx.x;
int j;
float v;
v = 0.;
for(j = 0; j < n; j++){
if(i != j) v += a[i*n + j]*x0[j];
}
x[i] = ALPHA*(b[i] - v)/a[i*n + i] + (1.f - ALPHA)*x0[i];
}
int jacobi(int n, const float *a, const float *b, float *x, float tol,
int maxIter, float *resFinal, int *iterFinal)
{
...
grid = (int)ceil((float)n/BLOCK_SIZE);
bNorm = calcNorm(n, b);
if(bNorm == 0.) bNorm = 1.;
for(iter = 0; iter < maxIter; iter++){
cudaMemcpy(x0_d, x_d, n*sizeof(float), cudaMemcpyDeviceToDevice);
jacobiInternal<<<grid, BLOCK_SIZE>>>(n, a_d, b_d, x0_d, x_d);
calcResidual(n, a_d, b_d, x_d, r_d);
cudaMemcpy(r, r_d, n*sizeof(float), cudaMemcpyDeviceToHost);
rNorm = calcNorm(n, r);
res = rNorm/bNorm;
if(res < tol) break;
}
...
return 0;
}
ヤコビ法の中心処理部分を jacobiInternal() としてカーネル関数にしています。スレッドは解ベクトル x の各要素に対応させています。また、残差計算 calcResidual() 用に行列とベクトルの乗算 multMatrixVector()、ベクトルの減算 subVector() をカーネル関数として定義しています。ここではグリッドサイズ (ブロック数) をつぎのように計算しています。 grid = (int)ceil((float)n/BLOCK_SIZE); これは行列サイズ n が BLOCK_SIZE で割り切れない場合にも対処するためのものです。 ガウス-ザイデル法ヤコビ法をガウス-ザイデル法に書き直すのは簡単なので、ついでにガウス-ザイデル法のプログラムも作りました。 gauss_seidel.cu
...
__global__ void gaussSeidelInternal(int n, const float *a, const float *b, float *x)
{
int i = blockDim.x*blockIdx.x + threadIdx.x;
int j;
float v;
v = 0.;
for(j = 0; j < n; j++){
if(i != j) v += a[i*n + j]*x[j];
}
x[i] = ALPHA*(b[i] - v)/a[i*n + i] + (1. - ALPHA)*x[i];
}
...
x の値を計算するときに前の値 x0 を使わないようにしただけです。 速度の比較上で見た 3 つのプログラムについて実行速度を比較します。入力データには OpenFOAM から取り出したサイズ 2500 の方程式を使います。収束のための許容値 (TOL) は 1e-3 とします。 マシンスペック
計算結果
注意: CUDA を使用するプログラムは、たぶん初期化のためだと思うが、たまに数秒余計に時間がかかることがある。 GPU 計算により計算速度が CPU 計算のおよそ 4 倍になりました。ガウス-ザイデル法については、ヤコビ法よりデータの独立性が弱まるために速度が遅くなるかもしれないと思いましたが、あまり変わりませんでした。 | ||||||||||||||||||||
| PENGUINITIS | ||||||||||||||||||||