メニーコアと格闘ブログ:CUDAを使ってプログラムとC言語プログラムの比較 その1

出典: トータル・ディスクロージャ・サイト(事実をありのままに)

CUDAプログラムは、C言語を拡張した形で記述します。

今回は簡単な数値計算を例に、C言語プログラムとCUDAプログラムの書式の違い、動作の違いを追いかけてみたいと思います。下で示す例題は、“GPGPU 勉強会”というホームページからの拝借です。

http://epa.scitec.kobe-u.ac.jp/~gpgpu/hiki/hiki.cgi?CUDA+Samples#1-BC-A1-B8-B5-C0-FE-B7-C1-B0-DC-CE-AE-B3-C8-BB-B6-CA-FD-C4-F8-BC-B0


[例題]

例題として用いる数値計算は、1次元線形移流拡散方程式と呼ばれるものです。

支配方程式と差分式は以下のように表わされます。

ファイル:CUDAを使ってプログラムとC言語プログラムの比較01.jpg

この差分式を解くプログラムについて、C言語CUDAを用いたプログラムを以下に示します。


[コンパイル方法]

  • C言語のコンパイル:

  $ gcc ad.c -o ad.exe -lm (-lmtで数学ライブラリをリンクします)

  • CUDAプログラムコンパイル:

  $ nvcc ad1.cu -o ad1.exe


[バグ?]

実はC言語のプログラム(ad.c)は、このままコンパイルすると正常動作しません。コンパイルは通りますが、実行結果に“nan”と出力されます。

46行目は

u[n] += du[n-1]*dt;

と、なるのではないかと思うのですが…(ad1.cuとの比較から類推しました)。



[C言語プログラム例: ad.c]

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#define NX 10000  //1次元要素数
#define DX 0.1f   //要素刻み幅

#define NT 10000  //計算時間
#define DT 1.0e-2 //時間刻み

#define NU 2.0e-1

void
set_initial(float *u, float *x, int nx, float dx)
{
  int n;

  for (n=0; n<nx; n++) {
    x[n] = dx * n;
  }
  for (n=1; n<=nx; n++) {
    u[n] = expf( -(x[n-1]-x[nx-1]/2)*(x[n-1]-x[nx/2])/(nx*DX/10) );
  }
  u[0] = u[nx];
  u[nx+1] = u[1];
}

void
output(FILE *file, float *u, float *x, int nx, float t)
{
  int n;
  fprintf(file, "# t = %f\n", t);
  for (n=0; n<nx; n+=(((nx-1)/100)+1))
    fprintf(file, "%f ", u[n]);
  fprintf(file, "\n");
}

void
integration(float *u, float *du, int nx, float nu, float dx, float dt)
{ // 移流拡散方程式を計算する関数
  int n;
  for(n=1; n<=nx; n++)
    du[n-1] = -(u[n+1]-u[n-1])/(2*dx) + nu * (u[n+1]-2*u[n]+u[n-1])/(dx*dx);
  for(n=1; n<=nx; n++)
    u[n] += du[n-1];
  u[0] = u[nx];
  u[nx+1] = u[1];
}

int
main()
{
  float x[NX];
  float u[NX+2];
  float du[NX];
  int n;
  float t;
  FILE *file;

  t = 0.0f;
  set_initial(u, x, NX, DX);

  file = fopen("result.dat", "w");

  fprintf(file, "# x: ");
  for(n=0; n<NX; n+=(((NX-1)/100)+1))
    fprintf(file, "%f ", x[n]);
  fprintf(file, "\n");

  output(file, &u[1], x, NX, t);

  for (n=1; n<=NT; n++) {
    t = DT * n;
    integration(u, du, NX, NU, DX, DT); //関数呼び出し
    if ((n%((NT-1)/100+1)) == 0)
      output(file, &u[1], x, NX, t);
  }

  fclose(file);

  return 0;
}


[CUDAプログラム例: ad1.cu]

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#define NX 10000
#define DX 0.1f

#define NT 100000
#define DT 1.0e-2

#define NU 2.0e-1

#define BLOCKDIM 192

void
set_initial(float *u, float *u_h, float *x, int nx, float dx)
{
  int n;

  for (n=0; n<nx; n++) {
    x[n] = dx * n;
  }

  for (n=1; n<=nx; n++) {
    u_h[n] = expf( -(x[n-1]-x[nx-1]/2)*(x[n-1]-x[nx/2])/(nx*dx/10) );
  }
  u_h[0] = u_h[nx];
  u_h[nx+1] = u_h[1];
  cudaMemcpy(u, u_h, sizeof(float)*(nx+2), cudaMemcpyHostToDevice); //③
}

void
output(FILE *file, float *u, float *u_h, float *x, int nx, float t)
{
  int n;

  cudaMemcpy(u_h, u, sizeof(float)*nx, cudaMemcpyDeviceToHost); //⑦
  fprintf(file, "# t = %f\n", t);
  for (n=0; n<nx; n+=(((nx-1)/100)+1))
    fprintf(file, "%f ", u_h[n]);
  fprintf(file, "\n");
}

__global__
void
integration_kernel0(float *u, float *du, int nx, float nu, float dx)
{ //デバイス側で処理される関数
  int n = BLOCKDIM * blockIdx.x + threadIdx.x + 1;
  float u0, u1, u2;
  if (n<=nx) {
    u0 = u[n-1];
    u1 = u[n];
    u2 = u[n+1];
    du[n-1] = -(u2-u0)/(2*dx) + nu * (u2-2*u1+u0)/(dx*dx);
  }
}  
__global__
void
integration_kernel1(float *u, float *du, int nx, float dt)
{ // デバイス側で処理される関数
  int n = blockDim.x * blockIdx.x + threadIdx.x;
  if (n<nx)
    u[n+1] += du[n] * dt;
}
void
integration(float *u, float *du, int nx, float nu, float dx, float dt)
{ //CPU側で処理される関数
  integration_kernel0<<<(nx-1)/BLOCKDIM+1, BLOCKDIM>>>(u, du, nx, nu, dx); //④ デバイス側 カーネル関数その1呼び出し
  integration_kernel1<<<(nx-1)/BLOCKDIM+1, BLOCKDIM>>>(u, du, nx, dt); //⑤ デバイス側 カーネル関数その2呼び出し
  cudaMemcpy(&u[0], &u[nx], sizeof(float), cudaMemcpyDeviceToDevice); //⑥ メモリ内容コピー
  cudaMemcpy(&u[nx+1], &u[1], sizeof(float), cudaMemcpyDeviceToDevice); //⑥ メモリ内容コピー
}

int
main()
{
  float x[NX];
  float *u_h;
  float *u, *du;
  int n;
  float t;
  FILE *file;

  cudaMalloc((void**)&u, sizeof(float)*(NX+2)); //① デバイス側 メモリ確保
  cudaMalloc((void**)&du, sizeof(float)*NX); //① デバイス側メモリ確保
  cudaMallocHost((void**)&u_h, sizeof(float)*(NX+2)); //② ホスト側メモリ確保?
  

  t = 0.0f;
  set_initial(u, u_h, x, NX, DX);

  file = fopen("result.dat", "w");

  fprintf(file, "# x: ");
  for (n=0; n<NX; n+=(((NX-1)/100)+1))
    fprintf(file, "%f ", x[n]);
  fprintf(file, "\n");

  output(file, &u[1], u_h, x, NX, t);

  for (n=1; n<=NT; n++) {
    t = DT * n;
    integration(u, du, NX, NU, DX, DT); // 関数呼び出し
    if ((n%((NT-1)/100+1)) == 0)
      output(file, &u[1], u_h, x, NX, t);
  }

  fclose(file);

  cudaFree(u);
  cudaFree(du);
  cudaFreeHost(u_h);

  return 0;
}



[実行結果]

$ time ./ad.exe 

real    0m5.720s
user    0m5.716s
sys     0m0.003s


$ time ./ad1.exe 

real    0m2.390s
user    0m0.085s
sys     0m2.234s


実行時間が半分になりました。


ナン 2009年8月20日 (木) 15:28 (JST)


目次

この記事へのコメントをお寄せください

  • サイトへの書き込みに差し支えございましたら トータルディスクロージャーサイトサポート係へメールをお送りください
  • トータル・ディスクロージャ・サイトに投稿された文章と画像は、すべてその著作権がHPCシステムズ株式会社に帰属し、HPCシステムズ株式会社が著作権を所有することに同意してください。
  • あなたの文章が他人によって自由に編集、配布されることを望まない場合は、投稿を控えてください。
  • コメントを書き込む場合は名前にひらがなを織り交ぜてください。
  • あなたの投稿する文章と画像はあなた自身によって書かれたものであるか、パブリック・ドメインかそれに類する自由なリソースからの複製であることを約束してください。あなたが著作権を保持していない作品を許諾なしに投稿してはいけません!

<comments hideform="false" />


Comments

oyakata said ...

CPU版に比べてGPU版の計算時間NTが10倍になってますけど、当然そろえた上での実行時間比較ですよね?

--oyakata 2010年2月2日 (火) 04:57 (UTC)

ヒロちゃん said ...

> CPU版に比べてGPU版の計算時間NTが10倍になってますけど、当然そろえた上での実行時間比較ですよね?

はい、NTをそろえた上で実行時間比較をしております。誤解をうむ表記をしてしまい申し訳ありませんでした。

念のため、別環境(CPU: Xeon E5540 @ 2.53GHz、GPU: Tesla C1060)ですが、NT を 100000 でそろえたところ、次となりました。 $ time ./ad.exe 16.077u 0.000s 0:16.08 99.9% 0+0k 0+88io 0pf+0w $ time ./ad1.exe 0.736u 4.573s 0:05.37 98.6% 0+0k 0+184io 0pf+0w


--ヒロちゃん 2010年2月5日 (金) 23:47 (UTC)

個人用ツール