N体問題
今回の内容
天体の運動方程式
天体運動のGPU実装
最適化による性能変化
計算の種類
画像処理，差分法
空間に固定された観測点を配置
観測点上で物理量がどのように
変化するかを追跡
Euler型
多粒子の運動
観測点を配置せず，観測点が粒
子と共に移動
Lagrange型
観測点
（固定）
観測点
（粒子と共に移動）
N体問題
万有引力によって運動する天体群のシミュレーション
天体を質点と考え，衝突を無視
天体は自身以外の天体から受ける万有引力のみで運動
天体
iの加速度
bN
i
j
j
ji
ji
j
i
m
,
1
3
r
r
g
a
r
ji
: 2天体
j, iの距離
g : 重力加速度
m
j
: 天体の質量
N
b
: 天体の総数
時間微分の離散化
天体
iの位置rと速度v
時間微分の離散化
テイラー展開を利用
i
i
dt
d
v
r
i
i
dt
d
a
v
2 22 3 33!
3
!
2
)
(
)
(
dt
d
Δt
dt
d
Δt
dt
d
Δt
t
Δt
t
r
r
r
r
r
v
r
r
r
r
r
2 22 3 33
!
3
!
2
1
)
(
)
(
dt
d
Δt
dt
d
Δt
Δt
Δt
t
Δt
t
dt
d
時間微分の離散化
離散時刻を上付き添字
nで表現
位置
r，速度vの時間発展
時刻
nでの天体位置r，速度vが既知
加速度
a
n
を求めることで時刻
n+1の天体位置を予測
n
t
r
r
(
)
1
)
(
t
Δt
r
n
r
Δt
n
i
n
i
n
i
r
v
r
1
Δt
n
i
n
i
n
i
v
a
v
1
プログラムの流れ
プログラムは2重ループ
全ペア相互作用(all pair interaction)
全ての天体が全ての天体のデータを利用
一度読み込んだデータの再利用が重要
計算量は天体個数の2乗に比例
for 全ての天体について{
for 全ての天体について{
天体間の距離を計算
if(同じ天体でなければ){
万有引力を計算
加速度aの値を更新
}
}
}
プログラムの流れ
プログラムは2重ループ
全ペア相互作用(all pair interaction)
全ての天体が全ての天体のデータを利用
一度読み込んだデータの再利用が重要
計算量は天体個数の2乗に比例
for 全ての天体について{
for 全ての天体について{
天体間の距離を計算
if(同じ天体でなければ){
万有引力を計算
加速度aの値を更新
}
}
}
プログラムの流れ
プログラムは2重ループ
全ペア相互作用(all pair interaction)
全ての天体が全ての天体のデータを利用
一度読み込んだデータの再利用が重要
計算量は天体個数の2乗に比例
for 全ての天体について{
for 全ての天体について{
天体間の距離を計算
if(同じ天体でなければ){
万有引力を計算
加速度aの値を更新
}
}
}
プログラムの流れ
プログラムは2重ループ
全ペア相互作用(all pair interaction)
全ての天体が全ての天体のデータを利用
一度読み込んだデータの再利用が重要
計算量は天体個数の2乗に比例
for 全ての天体について{
for 全ての天体について{
天体間の距離を計算
if(同じ天体でなければ){
万有引力を計算
加速度aの値を更新
}
}
}
プログラムの流れ
プログラムは2重ループ
全ペア相互作用(all pair interaction)
全ての天体が全ての天体のデータを利用
一度読み込んだデータの再利用が重要
計算量は天体個数の2乗に比例
for 全ての天体について{
for 全ての天体について{
天体間の距離を計算
if(同じ天体でなければ){
万有引力を計算
加速度aの値を更新
}
}
}
プログラムの流れ
万有引力は同一天体には影響しない
同一天体間の万有引力の計算は不可能（ゼロ除算）
ifで区別
if文を排除し，ゼロ除算を回避するため，微小値を付加
: 軟化因子
同一天体なら分子が0となるので計算結果に影響しない
bN
i
j
j
ji
ji
j
n
i
m
,
1
2
/
3
2
2
r
r
g
a
//重力加速度が1になるように規格化 #include<stdio.h> #include<stdlib.h> #include<math.h> #define N (4096) #define dt 0.001 #define Soften (1e‐6) //プロトタイプ宣言 void initial(…); void kernel(…); void integrate(…); int main(void){ float *x,*y,*z,*m; float *vx,*vy,*vz; float *ax,*ay,*az; x = (float *)malloc(N*sizeof(float)); y = (float *)malloc(N*sizeof(float)); z = (float *)malloc(N*sizeof(float)); m = (float *)malloc(N*sizeof(float)); vx = (float *)malloc(N*sizeof(float)); vy = (float *)malloc(N*sizeof(float)); vz = (float *)malloc(N*sizeof(float)); ax = (float *)malloc(N*sizeof(float)); ay = (float *)malloc(N*sizeof(float)); az = (float *)malloc(N*sizeof(float)); //初期値設定 initial(x,y,z,vx,vy,vz,m); //for(){ //本来なら必要な回数だけ繰り返す //加速度の計算 kernel(x,y,z,vx,vy,vz,ax,ay,az,m); //時間積分 integrate(x,y,z,vx,vy,vz,ax,ay,az); //} return 0; }
CPUプログラム
nbody.c
//初期値の設定
void initial(float *x,float *y,float *z,float *vx,float *vy,float *vz,float *m){ int i; //乱数で配置を決定 //x,y座標，x,y方向速度が‐1~1の範囲に収まるように決定 srand(N); for(i=0;i<N;i++){ x[i] = (float)rand()/RAND_MAX*2.0 ‐ 1.0; y[i] = (float)rand()/RAND_MAX*2.0 ‐ 1.0; z[i] = 0.0f; m[i] = 1.0f; vx[i] = (float)rand()/RAND_MAX*2.0 ‐ 1.0; vy[i] = (float)rand()/RAND_MAX*2.0 ‐ 1.0; vz[i] = 0.0f; } }
CPUプログラム
nbody.c
//時間積分 void integrate(float *x, float *y, float *z, float *vx, float *vy, float *vz, float *ax, float *ay, float *az){ int i; //Euler法で位置と速度を積分 //必ず位置の積分を先に実行 for(i=0;i<N;i++){ x[i] = x[i] + dt*vx[i]; y[i] = y[i] + dt*vy[i]; z[i] = z[i] + dt*vz[i]; vx[i] = vx[i] + dt*ax[i]; vy[i] = vy[i] + dt*ay[i]; vz[i] = vz[i] + dt*az[i]; } }
CPUプログラム
nbody.c
//加速度の計算 void kernel(float *x, float *y, float *z, float *vx, float *vy, float *vz, float *ax, float *ay, float *az, float *m){ int i,j; float rx,ry,rz; float dist2, dist6, invDist3,s; for(i=0;i<N;i++){ ax[i] = 0.0f; ay[i] = 0.0f; az[i] = 0.0f; } for(i=0;i<N;i++){ for(j=0;j<N;j++){ //if(i==j)continue; rx=x[j]‐x[i]; ry=y[j]‐y[i]; rz=z[j]‐z[i]; //2天体間の距離を計算 dist2 = rx*rx + ry*ry + rz*rz + Soften;//軟化パラメータ // m/r^3 の計算 dist6 = dist2*dist2*dist2; invDist3 = 1.0/sqrt(dist6); s = m[j]*invDist3; //天体jによる加速度を加算 ax[i] = ax[i] + rx*s; ay[i] = ay[i] + ry*s; az[i] = az[i] + rz*s; } } }
CPUプログラム
nbody.c
#include<stdio.h> #include<stdlib.h> #include<math.h> #define N (64*64) #define dt 0.001 #define Soften (1e‐6) #include "kernel0.cu"//カーネルを切替 int main(void){ //GPUのメモリ上に確保 float *x,*y,*z,*m; float *vx,*vy,*vz; float *ax,*ay,*az; //CPUのメモリ上に確保 初期設定用 float *host_x,*host_y,*host_z,*host_m; float *host_vx,*host_vy,*host_vz; cudaMalloc((void **)&x, (N*sizeof(float))); cudaMalloc((void **)&y, (N*sizeof(float))); cudaMalloc((void **)&z, (N*sizeof(float))); cudaMalloc((void **)&m, (N*sizeof(float))); cudaMalloc((void **)&vx, (N*sizeof(float))); cudaMalloc((void **)&vy, (N*sizeof(float))); cudaMalloc((void **)&vz, (N*sizeof(float))); cudaMalloc((void **)&ax, (N*sizeof(float))); cudaMalloc((void **)&ay, (N*sizeof(float))); cudaMalloc((void **)&az, (N*sizeof(float))); host_x = (float *)malloc(N*sizeof(float)); //以下，host_y等も確保 //初期設定はCPUと同じ initial(host_x,host_y,host_z, host_vx,host_vy,host_vz,host_m); cudaMemcpy(x, host_x, N*sizeof(float), cudaMemcpyHostToDevice); //以下，host_y等もGPUへコピー
GPUプログラム（1スレッド版）
nbody.cu
//for(){ //本来は必要な回数だけ繰り返す kernel<<<NB,NT>>> (x,y,z,vx,vy,vz,ax,ay,az,m); integrate<<<NB,NT>>> (x,y,z,vx,vy,vz,ax,ay,az); //必要な結果をCPUへコピー //} free(host_x); free(host_y); free(host_z); free(host_m); free(host_vx); free(host_vy); free(host_vz); cudaFree(x); cudaFree(y); cudaFree(z); cudaFree(m); cudaFree(vx); cudaFree(vy); cudaFree(vz); cudaFree(ax); cudaFree(ay); cudaFree(az); return 0; }
GPUプログラム（1スレッド版）
nbody.cu
#define NT 1 #define NB 1 //加速度の計算 __global__ void kernel(float *x,float *y, float *z,float *vx, float *vy, float *vz,float *ax, float *ay, float *az,float *m){ int i,j; float rx,ry,rz; float dist2, dist6, invDist3,s; for(i=0;i<N;i++){ ax[i] = 0.0f; ay[i] = 0.0f; az[i] = 0.0f; } for(i=0;i<N;i++){ for(j=0;j<N;j++){ if(i==j)continue; rx=x[j]‐x[i]; ry=y[j]‐y[i]; rz=z[j]‐z[i]; //2天体間の距離を計算 dist2 = rx*rx + ry*ry + rz*rz; // m/r^3 の計算 dist6 = dist2*dist2*dist2; invDist3 = 1.0/sqrt(dist6); s = m[j]*invDist3; //天体jによる加速度を加算 ax[i] = ax[i] + rx*s; ay[i] = ay[i] + ry*s; az[i] = az[i] + rz*s; } } }
GPUプログラム（1スレッド版）
kernel0.cu
//時間積分 __global__ void integrate(float *x, float *y, float *z, float *vx, float *vy, float *vz, float *ax, float *ay, float *az){ int i; for(i=0;i<N;i++){ x[i] = x[i] + dt*vx[i]; y[i] = y[i] + dt*vy[i]; z[i] = z[i] + dt*vz[i]; vx[i] = vx[i] + dt*ax[i]; vy[i] = vy[i] + dt*ay[i]; vz[i] = vz[i] + dt*az[i]; } }
GPUプログラム（1スレッド版）
kernel0.cu
実行時間
天体の個数 N = 4096
カーネル
実行時間 [ms]
#define NT 1 #define NB 1 //加速度の計算 __global__ void kernel(float *x,float *y, float *z,float *vx, float *vy, float *vz,float *ax, float *ay, float *az,float *m){ int i,j; float rx,ry,rz; float dist2, dist6, invDist3,s; for(i=0;i<N;i++){ ax[i] = 0.0f; ay[i] = 0.0f; az[i] = 0.0f; } for(i=0;i<N;i++){ for(j=0;j<N;j++){ rx=x[j]‐x[i]; ry=y[j]‐y[i]; rz=z[j]‐z[i]; //2天体間の距離を計算 dist2 = rx*rx + ry*ry + rz*rz +Soften; // m/r^3 の計算 dist6 = dist2*dist2*dist2; invDist3 = 1.0/sqrt(dist6); s = m[j]*invDist3; //天体jによる加速度を加算 ax[i] = ax[i] + rx*s; ay[i] = ay[i] + ry*s; az[i] = az[i] + rz*s; } } }
GPUプログラム（1スレッド版，if消去）
kernel1.cu
実行時間
天体の個数 N = 4096
カーネル
実行時間 [ms]
1スレッド実行
33×10
31スレッドが天体一つの加速度を計算
天体のデータを1次元配列で保持
加速度を求めたい天体iのforループ
相互作用の計算のために参照される天体jのループ
x[i]
x[j]
加速度を求めたい
天体
相互作用の計算の
ために参照される
天体
・・・
GPUによる並列処理の方針
ベクトル和と本質的には同じ
・・・
・・・
・・・
c[i]
a[i]
b[i]
＋
for(i=0;i<N;i++){
c[i]=a[i]+b[i];
}
＋
＋
＋
＋
＋
GPUによる並列処理の方針
1スレッドが一つの配列添字を計算
・・・
・・・
・・・
c[i]
a[i]
b[i]
i=blockIdx.x*blockDim.x+threadIdx.x;
c[i]=a[i]+b[i];
＋
＋
＋
＋
＋
＋
GPUによる並列処理の方針
加速度の計算
・・・
x[j]
for(i=0;i<N;i++){
for(j=0;j<N;j++){
加速度を積算
}
}
・・・
x[i]
加速度を求めたい
天体
相互作用の計算の
ために参照される
jのループ
GPUによる並列処理の方針
加速度の計算
・・・
x[j]
for(i=0;i<N;i++){
for(j=0;j<N;j++){
加速度を積算
}
}
・・・
x[i]
加速度を求めたい
天体
相互作用の計算の
ために参照される
jのループ
GPUによる並列処理の方針
加速度の計算
・・・
x[j]
for(i=0;i<N;i++){
for(j=0;j<N;j++){
加速度を積算
}
}
・・・
x[i]
加速度を求めたい
天体
相互作用の計算の
ために参照される
jのループ
GPUによる並列処理の方針
加速度の計算
・・・
x[j]
for(i=0;i<N;i++){
for(j=0;j<N;j++){
加速度を積算
}
}
・・・
x[i]
加速度を求めたい
天体
相互作用の計算の
ために参照される
jのループ
GPUによる並列処理の方針
加速度の計算
・・・
x[j]
for(i=0;i<N;i++){
for(j=0;j<N;j++){
加速度を積算
}
}
・・・
x[i]
加速度を求めたい
天体
相互作用の計算の
ために参照される
jのループ
GPUによる並列処理の方針
加速度の計算
・・・
x[j]
for(i=0;i<N;i++){
for(j=0;j<N;j++){
加速度を積算
}
}
・・・
x[i]
加速度を求めたい
天体
相互作用の計算の
ために参照される
jのループ
GPUによる並列処理の方針
ベクトル和と同様に1スレッドが一つの天体を計算
i=blockIdx.x*blockDim.x+threadIdx.x;
for(j=0;j<N;j++){
加速度を積算
}
・・・
x[j]
・・・
x[i]
加速度を求めたい
天体
相互作用の計算の
ために参照される
GPUによる並列処理の方針
ベクトル和と同様に1スレッドが一つの天体を計算
i=blockIdx.x*blockDim.x+threadIdx.x;
for(j=0;j<N;j++){
加速度を積算
}
・・・
x[j]
・・・
x[i]
加速度を求めたい
天体
相互作用の計算の
ために参照される
GPUによる並列処理の方針
ベクトル和と同様に1スレッドが一つの天体を計算
i=blockIdx.x*blockDim.x+threadIdx.x;
for(j=0;j<N;j++){
加速度を積算
}
・・・
x[j]
・・・
x[i]
加速度を求めたい
天体
相互作用の計算の
ために参照される
GPUによる並列処理の方針
ベクトル和と同様に1スレッドが一つの天体を計算
i=blockIdx.x*blockDim.x+threadIdx.x;
for(j=0;j<N;j++){
加速度を積算
}
・・・
x[j]
・・・
x[i]
加速度を求めたい
天体
相互作用の計算の
ために参照される
GPUによる並列処理の方針
ベクトル和と同様に1スレッドが一つの天体を計算
i=blockIdx.x*blockDim.x+threadIdx.x;
for(j=0;j<N;j++){
加速度を積算
}
・・・
x[j]
・・・
x[i]
加速度を求めたい
天体
相互作用の計算の
ために参照される
GPUによる並列処理の方針
ベクトル和と同様に1スレッドが一つの天体を計算
i=blockIdx.x*blockDim.x+threadIdx.x;
for(j=0;j<N;j++){
加速度を積算
}
・・・
x[j]
・・・
x[i]
加速度を求めたい
天体
相互作用の計算の
ために参照される
#define NT 256 #define NB (N/NT) //加速度の計算 __global__ void kernel(float *x,float *y, float *z,float *vx, float *vy, float *vz,float *ax, float *ay, float *az,float *m){ int i,j; float rx,ry,rz; float dist2, dist6, invDist3,s; i = blockDim.x*blockIdx.x+threadIdx.x; for(j=0;j<N;j++){ rx=x[j]‐x[i]; ry=y[j]‐y[i]; rz=z[j]‐z[i]; //2天体間の距離を計算 dist2 = rx*rx + ry*ry + rz*rz + Soften;//軟化パラメータ // m/r^3 の計算 dist6 = dist2*dist2*dist2; invDist3 = 1.0/sqrt(dist6); s = m[j]*invDist3; //天体jによる加速度を加算 ax[i] = ax[i] + rx*s; ay[i] = ay[i] + ry*s; az[i] = az[i] + rz*s; } }
GPUプログラム（1スレッドが天体一つを計算）
kernel2.cu
//時間積分
__global__ void integrate(float *x, float *y, float *z,
float *vx, float *vy, float *vz, float *ax, float *ay, float *az){
int i=blockIdx.x*blockDim.x + threadIdx.x;
x[i] = x[i] + dt*vx[i]; y[i] = y[i] + dt*vy[i]; z[i] = z[i] + dt*vz[i]; vx[i] = vx[i] + dt*ax[i]; vy[i] = vy[i] + dt*ay[i]; vz[i] = vz[i] + dt*az[i]; }
GPUプログラム（1スレッドが天体一つを計算）
kernel2.cu
実行時間
天体の個数 N = 4096
スレッド数
NT = 256
カーネル
実行時間 [ms]
1スレッド実行
33×10
31スレッド実行（ifの消去）
29×10
31スレッド1天体
7.63
#define NT 256 #define NB (N/NT) //加速度の計算 __global__ void kernel(float *x,float *y, float *z,float *vx, float *vy, float *vz,float *ax, float *ay, float *az,float *m){ int i,j; float rx,ry,rz; float dist2, dist6, invDist3,s; //加速度のデータをレジスタに置く float r_ax, r_ay, r_az; i = blockDim.x*blockIdx.x+threadIdx.x;
r_ax = r_ay = r_az = 0.0f; for(j=0;j<N;j++){ rx=x[j]‐x[i]; ry=y[j]‐y[i]; rz=z[j]‐z[i]; //2天体間の距離を計算 dist2 = rx*rx + ry*ry + rz*rz + Soften;//軟化パラメータ // m/r^3 の計算 dist6 = dist2*dist2*dist2; invDist3 = 1.0/sqrt(dist6); s = m[j]*invDist3; //天体jによる加速度を加算 r_ax = r_ax + rx*s; r_ay = r_ay + ry*s; r_az = r_az + rz*s; } ax[i] = r_ax; ay[i] = r_ay; az[i] = r_az; }
GPUプログラム（レジスタ利用）
kernel3.cu
実行時間
天体の個数 N = 4096
スレッド数
NT = 256
カーネル
実行時間 [ms]
1スレッド実行
33×10
31スレッド実行（ifの消去）
29×10
31スレッド1天体
7.63
レジスタ利用
4.49
GPUプログラム（共有メモリによる再利用）
x[i]
x[j]
・・・
共有メモリ
blockIdx.x=0
blockIdx.x=1
i= 0 1 2 3 4 5 6 7
threadIdx.x= 0 1 2 3 0 1 2 3
x[0] x[1] x[2] x[3] x[0] x[1] x[2] x[3]
j= 0 1 2 3 4 5 6 7
0+threadIdx.x
0+threadIdx.x
GPUプログラム（共有メモリによる再利用）
x[i]
x[j]
・・・
共有メモリ
blockIdx.x=0
blockIdx.x=1
i=0~3の天体と
j=0~3の天体で
加速度を計算
i=4~7の天体と
j=0~3の天体で
加速度を計算
threadIdx.x= 0 1 2 3 0 1 2 3
x[0] x[1] x[2] x[3] x[0] x[1] x[2] x[3]
j= 0 1 2 3 4 5 6 7
i= 0 1 2 3 4 5 6 7
GPUプログラム（共有メモリによる再利用）
x[i]
x[j]
・・・
共有メモリ
blockIdx.x=0
blockIdx.x=1
i=0~3の天体と
j=0~3の天体で
加速度を計算
i=4~7の天体と
j=0~3の天体で
加速度を計算
threadIdx.x= 0 1 2 3 0 1 2 3
x[0] x[1] x[2] x[3] x[0] x[1] x[2] x[3]
j= 0 1 2 3 4 5 6 7
i= 0 1 2 3 4 5 6 7
GPUプログラム（共有メモリによる再利用）
x[i]
x[j]
・・・
共有メモリ
blockIdx.x=0
blockIdx.x=1
i=0~3の天体と
j=0~3の天体で
加速度を計算
i=4~7の天体と
j=0~3の天体で
加速度を計算
threadIdx.x= 0 1 2 3 0 1 2 3
x[0] x[1] x[2] x[3] x[0] x[1] x[2] x[3]
j= 0 1 2 3 4 5 6 7
i= 0 1 2 3 4 5 6 7
GPUプログラム（共有メモリによる再利用）
x[i]
x[j]
・・・
共有メモリ
blockIdx.x=0
blockIdx.x=1
i=0~3の天体と
j=0~3の天体で
加速度を計算
i=4~7の天体と
j=0~3の天体で
加速度を計算
threadIdx.x= 0 1 2 3 0 1 2 3
x[0] x[1] x[2] x[3] x[0] x[1] x[2] x[3]
j= 0 1 2 3 4 5 6 7
i= 0 1 2 3 4 5 6 7
GPUプログラム（共有メモリによる再利用）
x[i]
x[j]
・・・
共有メモリ
blockIdx.x=0
blockIdx.x=1
i= 0 1 2 3 4 5 6 7
threadIdx.x= 0 1 2 3 0 1 2 3
x[4] x[5] x[6] x[7] x[4] x[5] x[6] x[7]
j= 0 1 2 3 4 5 6 7
4(=NT)+threadIdx.x
4(=NT)+threadIdx.x
GPUプログラム（共有メモリによる再利用）
x[i]
x[j]
・・・
共有メモリ
blockIdx.x=0
blockIdx.x=1
i=0~3の天体と
j=4~7の天体で
加速度を計算
i=4~7の天体と
j=4~7の天体で
加速度を計算
threadIdx.x= 0 1 2 3 0 1 2 3
j= 0 1 2 3 4 5 6 7
i= 0 1 2 3 4 5 6 7
x[4] x[5] x[6] x[7] x[4] x[5] x[6] x[7]
GPUプログラム（共有メモリによる再利用）
x[i]
x[j]
・・・
共有メモリ
blockIdx.x=0
blockIdx.x=1
i=0~3の天体と
j=4~7の天体で
加速度を計算
i=4~7の天体と
j=4~7の天体で
加速度を計算
threadIdx.x= 0 1 2 3 0 1 2 3
j= 0 1 2 3 4 5 6 7
i= 0 1 2 3 4 5 6 7
x[4] x[5] x[6] x[7] x[4] x[5] x[6] x[7]
GPUプログラム（共有メモリによる再利用）
x[i]
x[j]
・・・
共有メモリ
blockIdx.x=0
blockIdx.x=1
i=0~3の天体と
j=4~7の天体で
加速度を計算
i=4~7の天体と
j=4~7の天体で
加速度を計算
threadIdx.x= 0 1 2 3 0 1 2 3
j= 0 1 2 3 4 5 6 7
i= 0 1 2 3 4 5 6 7
x[4] x[5] x[6] x[7] x[4] x[5] x[6] x[7]
GPUプログラム（共有メモリによる再利用）
x[i]
x[j]
・・・
共有メモリ
blockIdx.x=0
blockIdx.x=1
i=0~3の天体と
j=4~7の天体で
加速度を計算
i=4~7の天体と
j=4~7の天体で
加速度を計算
threadIdx.x= 0 1 2 3 0 1 2 3
j= 0 1 2 3 4 5 6 7
i= 0 1 2 3 4 5 6 7
x[4] x[5] x[6] x[7] x[4] x[5] x[6] x[7]
#define NT 256 #define NB (N/NT) //加速度の計算 __global__ void kernel(…){ int i,j; float rx,ry,rz; float dist2, dist6, invDist3,s; //加速度と座標のデータをレジスタに置く float r_ax,r_ay,r_az,r_x,r_y,r_z; __shared__ float s_x[NT],s_y[NT],s_z[NT],s_m[NT]; i = blockDim.x*blockIdx.x+threadIdx.x; r_ax = r_ay = r_az = 0.0f;
r_x = x[i]; r_y = y[i]; r_z = z[i]; for(j=0;j<N;j+=NT){ s_x[threadIdx.x] = x[j+threadIdx.x]; s_y[threadIdx.x] = y[j+threadIdx.x]; s_z[threadIdx.x] = z[j+threadIdx.x]; s_m[threadIdx.x] = m[j+threadIdx.x]; __syncthreads(); for(js = 0;js<NT;js++){ rx=s_x[js]‐r_x; ry=s_y[js]‐r_y; rz=s_z[js]‐r_z; //2天体間の距離を計算 dist2 = rx*rx + ry*ry + rz*rz + Soften;//軟化パラメータ // m/r^3 の計算 dist6 = dist2*dist2*dist2; invDist3 = 1.0/sqrt(dist6); s = s_m[js]*invDist3; //天体jによる加速度を加算 r_ax = r_ax + rx*s; r_ay = r_ay + ry*s; r_az = r_az + rz*s; } } ax[i] = r_ax; ay[i] = r_ay; az[i] = r_az; }
GPUプログラム（共有メモリによる再利用）
kernel4.cu
実行時間
天体の個数 N = 4096
スレッド数
NT = 256
カーネル
実行時間 [ms]
1スレッド実行
33×10
31スレッド実行（ifの消去）
29×10
31スレッド1天体
7.63
レジスタ利用
4.49
共有メモリ利用
4.27
#define NT 256 #define NB (N/NT) //加速度の計算 __global__ void kernel(…){ int i,j; float rx,ry,rz; float dist2, dist6, invDist3,s; //加速度と座標のデータをレジスタに置く float r_ax,r_ay,r_az,r_x,r_y,r_z; __shared__ float s_x[NT],s_y[NT],s_z[NT],s_m[NT]; i = blockDim.x*blockIdx.x+threadIdx.x; r_ax = r_ay = r_az = 0.0f;
r_x = x[i]; r_y = y[i]; r_z = z[i]; for(j=0;j<N;j+=NT){ s_x[threadIdx.x] = x[j+threadIdx.x]; s_y[threadIdx.x] = y[j+threadIdx.x]; s_z[threadIdx.x] = z[j+threadIdx.x]; s_m[threadIdx.x] = m[j+threadIdx.x]; #pragma unroll for(js = 0;js<NT;js++){ rx=s_x[js]‐r_x; ry=s_y[js]‐r_y; rz=s_z[js]‐r_z; //2天体間の距離を計算 dist2 = rx*rx + ry*ry + rz*rz + Soften;//軟化パラメータ // m/r^3 の計算 dist6 = dist2*dist2*dist2; invDist3 = 1.0/sqrt(dist6); s = s_m[js]*invDist3; //天体jによる加速度を加算 r_ax = r_ax + rx*s; r_ay = r_ay + ry*s; r_az = r_az + rz*s; } } ax[i] = r_ax; ay[i] = r_ay; az[i] = r_az; }
