N 体問題 長岡技術科学大学電気電子情報工学専攻出川智啓

Transcription

1 N 体問題 長岡技術科学大学電気電子情報工学専攻出川智啓

2 今回の内容 天体の運動方程式 天体運動の GPU 実装 最適化による性能変化 #pragma unroll 855

3 計算の種類 画像処理, 差分法 空間に固定された観測点を配置 観測点 ( 固定 ) 観測点上で物理量がどのように変化するかを追跡 Euler 型 多粒子の運動 観測点を配置せず, 観測点が粒子と共に移動 Lagrange 型 観測点 ( 粒子と共に移動 ) 856

4 N 体問題 万有引力によって運動する天体群のシミュレーション 天体を質点と考え, 衝突を無視 天体は自身以外の天体から受ける万有引力のみで運動 天体 i の加速度 a N b i m j j 1, j i g r r ji 3 ji r ji : 2 天体 j, i の距離 g : 重力加速度 m j : 天体の質量 N b : 天体の総数 857

5 時間微分の離散化 天体 i の位置 r と速度 v dri dt v i dvi a dt i 時間微分の離散化 テイラー展開を利用 r( t dr Δt d r Δt d r Δt) r( t) Δt 2 3 dt 2! dt 3! dt dr dt r( t Δt) Δt r( t) 1 Δt Δt 2! 2 d 2 dt r 2 Δt 3! 3 d 3 dt r 3 v 858

6 時間微分の離散化 離散時刻を上付き添字 n で表現 r( t) r n r( t Δt) r n 1 位置 r, 速度 v の時間発展 時刻 n での天体位置 r, 速度 v が既知 加速度 a n を求めることで時刻 n+1 の天体位置を予測 r n 1 i r n i v n i Δt v n 1 i v n i a n i Δt 859

7 プログラムの流れ プログラムは 2 重ループ 全ペア相互作用 (all pair interaction) 全ての天体が全ての天体のデータを利用 一度読み込んだデータの再利用が重要 計算量は天体個数の 2 乗に比例 for 全ての天体について { for 全ての天体について { 天体間の距離を計算 if( 同じ天体でなければ ){ 万有引力を計算加速度 a の値を更新 860

8 プログラムの流れ プログラムは 2 重ループ 全ペア相互作用 (all pair interaction) 全ての天体が全ての天体のデータを利用 一度読み込んだデータの再利用が重要 計算量は天体個数の 2 乗に比例 for 全ての天体について { for 全ての天体について { 天体間の距離を計算 if( 同じ天体でなければ ){ 万有引力を計算加速度 a の値を更新 861

9 プログラムの流れ プログラムは 2 重ループ 全ペア相互作用 (all pair interaction) 全ての天体が全ての天体のデータを利用 一度読み込んだデータの再利用が重要 計算量は天体個数の 2 乗に比例 for 全ての天体について { for 全ての天体について { 天体間の距離を計算 if( 同じ天体でなければ ){ 万有引力を計算加速度 a の値を更新 862

10 プログラムの流れ プログラムは 2 重ループ 全ペア相互作用 (all pair interaction) 全ての天体が全ての天体のデータを利用 一度読み込んだデータの再利用が重要 計算量は天体個数の 2 乗に比例 for 全ての天体について { for 全ての天体について { 天体間の距離を計算 if( 同じ天体でなければ ){ 万有引力を計算加速度 a の値を更新 863

11 プログラムの流れ プログラムは 2 重ループ 全ペア相互作用 (all pair interaction) 全ての天体が全ての天体のデータを利用 一度読み込んだデータの再利用が重要 計算量は天体個数の 2 乗に比例 for 全ての天体について { for 全ての天体について { 天体間の距離を計算 if( 同じ天体でなければ ){ 万有引力を計算加速度 a の値を更新 864

12 プログラムの流れ 万有引力は同一天体には影響しない 同一天体間の万有引力の計算は不可能 ( ゼロ除算 ) ifで区別 if 文を排除し, ゼロ除算を回避するため, 微小値を付加 : 軟化因子 同一天体なら分子が 0 となるので計算結果に影響しない a n i N b j 1, j i m j g r ji 2 r ji 2 3/ 2 865

13 CPU プログラム // 重力加速度が 1 になるように規格化 #include<stdio.h> #include<stdlib.h> #include<math.h> #define N (4096) #define dt #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; nbody.c 866

14 CPU プログラム // 初期値の設定 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* ; y[i] = (float)rand()/rand_max* ; z[i] = 0.0f; m[i] = 1.0f; vx[i] = (float)rand()/rand_max* ; vy[i] = (float)rand()/rand_max* ; vz[i] = 0.0f; nbody.c 867

15 CPU プログラム // 時間積分 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]; nbody.c 868

16 CPU プログラム // 加速度の計算 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; 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; 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; // 天体 j による加速度を加算 ax[i] = ax[i] + rx*s; ay[i] = ay[i] + ry*s; az[i] = az[i] + rz*s; rx=x[j] x[i]; ry=y[j] y[i]; nbody.c 869

17 GPU プログラム (1 スレッド版 ) #include<stdio.h> #include<stdlib.h> #include<math.h> #define N (64*64) #define dt #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 へコピー nbody.cu 870

18 GPU プログラム (1 スレッド版 ) //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 へコピー // cudafree(vy); cudafree(vz); cudafree(ax); cudafree(ay); cudafree(az); return 0; 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); nbody.cu 871

19 GPU プログラム (1 スレッド版 ) #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; 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; for(i=0;i<n;i++){ for(j=0;j<n;j++){ if(i==j)continue; kernel0.cu 872

20 GPU プログラム (1 スレッド版 ) // 時間積分 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]; kernel0.cu 873

21 実行時間 天体の個数 N = 4096 カーネル 実行時間 [ms] 1 スレッド実行

22 GPU プログラム (1 スレッド版,if 消去 ) #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; 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; for(i=0;i<n;i++){ for(j=0;j<n;j++){ kernel1.cu 875

23 実行時間 天体の個数 N = 4096 カーネル 実行時間 [ms] 1 スレッド実行 スレッド実行 (if の消去 )

24 1 スレッドが天体一つの加速度を計算 天体のデータを1 次元配列で保持 加速度を求めたい天体 iのforループ 相互作用の計算のために参照される天体 j のループ 加速度を求めたい天体 x[i] 相互作用の計算のために参照される天体 x[j] 877

25 GPU による並列処理の方針 ベクトル和と本質的には同じ for(i=0;i<n;i++){ c[i]=a[i]+b[i]; a[i] b[i] c[i] 878

26 GPU による並列処理の方針 1 スレッドが一つの配列添字を計算 i=blockidx.x*blockdim.x+threadidx.x; c[i]=a[i]+b[i]; a[i] b[i] c[i] 879

27 GPU による並列処理の方針 加速度の計算 for(i=0;i<n;i++){ for(j=0;j<n;j++){ 加速度を積算 加速度を求めたい x[i] 天体 相互作用の計算のために参照される天体 x[j] j のループ 880

28 GPU による並列処理の方針 加速度の計算 for(i=0;i<n;i++){ for(j=0;j<n;j++){ 加速度を積算 加速度を求めたい x[i] 天体 相互作用の計算のために参照される天体 x[j] j のループ 881

29 GPU による並列処理の方針 加速度の計算 for(i=0;i<n;i++){ for(j=0;j<n;j++){ 加速度を積算 加速度を求めたい x[i] 天体 相互作用の計算のために参照される天体 x[j] j のループ 882

30 GPU による並列処理の方針 加速度の計算 for(i=0;i<n;i++){ for(j=0;j<n;j++){ 加速度を積算 加速度を求めたい x[i] 天体 相互作用の計算のために参照される天体 x[j] j のループ 883

31 GPU による並列処理の方針 加速度の計算 for(i=0;i<n;i++){ for(j=0;j<n;j++){ 加速度を積算 加速度を求めたい x[i] 天体 相互作用の計算のために参照される天体 x[j] j のループ 884

32 GPU による並列処理の方針 加速度の計算 for(i=0;i<n;i++){ for(j=0;j<n;j++){ 加速度を積算 加速度を求めたい x[i] 天体 相互作用の計算のために参照される天体 x[j] j のループ 885

33 GPU による並列処理の方針 ベクトル和と同様に 1 スレッドが一つの天体を計算 i=blockidx.x*blockdim.x+threadidx.x; for(j=0;j<n;j++){ 加速度を積算 加速度を求めたい x[i] 天体 相互作用の計算のために参照される天体 x[j] 886

34 GPU による並列処理の方針 ベクトル和と同様に 1 スレッドが一つの天体を計算 i=blockidx.x*blockdim.x+threadidx.x; for(j=0;j<n;j++){ 加速度を積算 加速度を求めたい x[i] 天体 相互作用の計算のために参照される天体 x[j] 887

35 GPU による並列処理の方針 ベクトル和と同様に 1 スレッドが一つの天体を計算 i=blockidx.x*blockdim.x+threadidx.x; for(j=0;j<n;j++){ 加速度を積算 加速度を求めたい x[i] 天体 相互作用の計算のために参照される天体 x[j] 888

36 GPU による並列処理の方針 ベクトル和と同様に 1 スレッドが一つの天体を計算 i=blockidx.x*blockdim.x+threadidx.x; for(j=0;j<n;j++){ 加速度を積算 加速度を求めたい x[i] 天体 相互作用の計算のために参照される天体 x[j] 889

37 GPU による並列処理の方針 ベクトル和と同様に 1 スレッドが一つの天体を計算 i=blockidx.x*blockdim.x+threadidx.x; for(j=0;j<n;j++){ 加速度を積算 加速度を求めたい x[i] 天体 相互作用の計算のために参照される天体 x[j] 890

38 GPU による並列処理の方針 ベクトル和と同様に 1 スレッドが一つの天体を計算 i=blockidx.x*blockdim.x+threadidx.x; for(j=0;j<n;j++){ 加速度を積算 加速度を求めたい x[i] 天体 相互作用の計算のために参照される天体 x[j] 891

39 GPU プログラム (1 スレッドが天体一つを計算 ) #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; kernel2.cu 892

40 GPU プログラム (1 スレッドが天体一つを計算 ) // 時間積分 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]; kernel2.cu 893

41 実行時間 天体の個数 N = 4096 スレッド数 NT = 256 カーネル 実行時間 [ms] 1 スレッド実行 スレッド実行 (if の消去 ) スレッド 1 天体

42 GPU プログラム ( レジスタ利用 ) #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; //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; 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]; ax[i] = r_ax; ay[i] = r_ay; az[i] = r_az; kernel3.cu 895

43 実行時間 天体の個数 N = 4096 スレッド数 NT = 256 カーネル 実行時間 [ms] 1スレッド実行 スレッド実行 (ifの消去) スレッド1 天体 7.63 レジスタ利用

44 GPU プログラム ( 共有メモリによる再利用 ) blockidx.x=0 blockidx.x=1 threadidx.x= i= x[i] x[0] x[1] x[2] x[3] x[0] x[1] x[2] x[3] 共有メモリ 0+threadIdx.x 0+threadIdx.x j= x[j] 897

45 GPU プログラム ( 共有メモリによる再利用 ) blockidx.x=0 blockidx.x=1 threadidx.x= i=0~3 の天体と j=0~3 の天体で加速度を計算 i= x[i] i=4~7 の天体と j=0~3 の天体で加速度を計算 x[0] x[1] x[2] x[3] x[0] x[1] x[2] x[3] 共有メモリ j= x[j] 898

46 GPU プログラム ( 共有メモリによる再利用 ) blockidx.x=0 blockidx.x=1 threadidx.x= i=0~3 の天体と j=0~3 の天体で加速度を計算 i= x[i] i=4~7 の天体と j=0~3 の天体で加速度を計算 x[0] x[1] x[2] x[3] x[0] x[1] x[2] x[3] 共有メモリ j= x[j] 899

47 GPU プログラム ( 共有メモリによる再利用 ) blockidx.x=0 blockidx.x=1 threadidx.x= i=0~3 の天体と j=0~3 の天体で加速度を計算 i= x[i] i=4~7 の天体と j=0~3 の天体で加速度を計算 x[0] x[1] x[2] x[3] x[0] x[1] x[2] x[3] 共有メモリ j= x[j] 900

48 GPU プログラム ( 共有メモリによる再利用 ) blockidx.x=0 blockidx.x=1 threadidx.x= i=0~3 の天体と j=0~3 の天体で加速度を計算 i= x[i] i=4~7 の天体と j=0~3 の天体で加速度を計算 x[0] x[1] x[2] x[3] x[0] x[1] x[2] x[3] 共有メモリ j= x[j] 901

49 GPU プログラム ( 共有メモリによる再利用 ) blockidx.x=0 blockidx.x=1 threadidx.x= i= x[i] x[4] x[5] x[6] x[7] x[4] x[5] x[6] x[7] 共有メモリ 4(=NT)+threadIdx.x 4(=NT)+threadIdx.x j= x[j] 902

50 GPU プログラム ( 共有メモリによる再利用 ) blockidx.x=0 blockidx.x=1 threadidx.x= i=0~3 の天体と j=4~7 の天体で加速度を計算 i= x[i] i=4~7 の天体と j=4~7 の天体で加速度を計算 x[4] x[5] x[6] x[7] x[4] x[5] x[6] x[7] 共有メモリ j= x[j] 903

51 GPU プログラム ( 共有メモリによる再利用 ) blockidx.x=0 blockidx.x=1 threadidx.x= i=0~3 の天体と j=4~7 の天体で加速度を計算 i= x[i] i=4~7 の天体と j=4~7 の天体で加速度を計算 x[4] x[5] x[6] x[7] x[4] x[5] x[6] x[7] 共有メモリ j= x[j] 904

52 GPU プログラム ( 共有メモリによる再利用 ) blockidx.x=0 blockidx.x=1 threadidx.x= i=0~3 の天体と j=4~7 の天体で加速度を計算 i= x[i] i=4~7 の天体と j=4~7 の天体で加速度を計算 x[4] x[5] x[6] x[7] x[4] x[5] x[6] x[7] 共有メモリ j= x[j] 905

53 GPU プログラム ( 共有メモリによる再利用 ) blockidx.x=0 blockidx.x=1 threadidx.x= i=0~3 の天体と j=4~7 の天体で加速度を計算 i= x[i] i=4~7 の天体と j=4~7 の天体で加速度を計算 x[4] x[5] x[6] x[7] x[4] x[5] x[6] x[7] 共有メモリ j= x[j] 906

54 GPU プログラム ( 共有メモリによる再利用 ) #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; kernel4.cu 907

55 実行時間 天体の個数 N = 4096 スレッド数 NT = 256 カーネル 実行時間 [ms] 1スレッド実行 スレッド実行 (ifの消去) スレッド1 天体 7.63 レジスタ利用 4.49 共有メモリ利用

56 GPU プログラム ( ループアンロール ) #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(); #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; kernel5.cu 909

57 実行時間 天体の個数 N = 4096 スレッド数 NT = 256 カーネル 実行時間 [ms] 1スレッド実行 スレッド実行 (ifの消去) スレッド1 天体 7.63 レジスタ利用 4.49 共有メモリ利用 4.27 共有メモリ利用 + ループアンロール

58 実行時間 天体の個数 N = スレッド数 NT = 256 カーネル 実行時間 [ms] 1スレッド1 天体 916 レジスタ利用 638 共有メモリ利用 646 共有メモリ利用 + ループアンロール 499 共有メモリを単純に利用しただけでは遅くなる 共有メモリを利用することは高速化に有効 共有メモリを利用するために追加される処理 が遅くなる要因 911