C
CUDA
GPU

CUDAによる素数リスト生成プログラム

More than 1 year has passed since last update.

環境

CUDA 9.0
NVIDIA GeForce GTX960
core i7 5820K
DDR4 32Gbyte

ソース

エラトステネスの篩は後述する。
まずは1から引数指定の数まで、順に除算して素数判定する方法で実装する。

prime_cuda.cu
#include <stdio.h>
#include <stdlib.h>
#include <cuda.h>
#include <helper_cuda.h>
#include <helper_functions.h>
#include <time.h>

#define Grid_x 1024
#define Grid_y 1024
#define Grid_z 1
#define Block_x 16
#define Block_y 8
#define Block_z 1


__global__ void thread_num(unsigned long long int  *device_result , unsigned long long int cycle);

/* timer */
int timer(void){
  time_t now = time(NULL);
  struct tm *pnow = localtime(&now);
  char buff[128]="";
  sprintf(buff,"%d:%d:%d",pnow->tm_hour,pnow->tm_min,pnow->tm_sec);
  printf("%s\n",buff);
  return 0;
}

int main(int argc, char** argv){

  FILE *outputfile;
  unsigned long long int num=0;
  unsigned long long int cycle = 0; 

  if (argv[1] == NULL) {
    printf("err\n");
    exit(1);
  }
  num=atoll(argv[1]);

  outputfile = fopen("./prime_num.txt", "w");
  if (outputfile == NULL) {
    printf("cannot open\n");
    exit(1);
  }
  timer();

  /*ブロックサイズとグリッドサイズの設定*/
  dim3 grid(Grid_x,Grid_y,Grid_z);
  dim3 block(Block_x,Block_y,Block_z);

  /*ホスト側の変数設定*/
  unsigned long long int  thread_size = (Grid_x * Grid_y * Grid_z) * (Block_x * Block_y * Block_z);
  unsigned long long int  *host_result;
  host_result = (unsigned long long int *)malloc(thread_size * sizeof(unsigned long long int));

  cycle = num / thread_size;

  /*デバイス側の変数設定*/
  unsigned long long int  *device_result;

  for (unsigned long long int i = 0 ; i < cycle + 1 ; i++){

  /*デバイスメモリ領域の確保*/
  checkCudaErrors(cudaMalloc((void**)&device_result, thread_size * sizeof(unsigned long long int)));

  /*タイマーを作成して計測開始*/
  cudaEvent_t start;
  cudaEvent_t stop;
  checkCudaErrors(cudaEventCreate(&start));
  checkCudaErrors(cudaEventCreate(&stop));
  checkCudaErrors(cudaEventRecord(start, NULL));
  printf("Range : %llu - %llu\n",thread_size*i , thread_size*(i+1)-1);

  /*カーネルの起動*/
  thread_num<<<grid , block>>>(device_result, i);
  cudaThreadSynchronize();

  /*タイマーを停止し実行時間を表示*/
  checkCudaErrors(cudaEventRecord(stop, NULL));
  checkCudaErrors(cudaEventSynchronize(stop));
  float msecTotal = 0.0f;
  checkCudaErrors(cudaEventElapsedTime(&msecTotal, start, stop));
  printf("Processing time: %f (msec)\n", msecTotal);

  /*再度タイマー開始*/
  checkCudaErrors(cudaEventRecord(start, NULL));

  /*結果の領域確保とデバイス側からのメモリ転送*/
  host_result = (unsigned long long int *)malloc(thread_size * sizeof(unsigned long long int));
  checkCudaErrors(cudaMemcpy(host_result, device_result, thread_size * sizeof(unsigned long long int ) , cudaMemcpyDeviceToHost));

  /*タイマーを停止し実行時間を表示*/
  checkCudaErrors(cudaEventRecord(stop, NULL));
  checkCudaErrors(cudaEventSynchronize(stop));
  checkCudaErrors(cudaEventElapsedTime(&msecTotal, start, stop));
  printf("Memory copy time: %f (msec)\n", msecTotal);
  printf("Now Writing...\n");

  for( unsigned long long int  l = 0; l < thread_size; l++ ){ 
    if (host_result[l] != 0){
      fprintf(outputfile,"%llu\n", host_result[l]);
    }
  }
  /*ホスト・デバイスメモリの開放*/
  free(host_result);
  checkCudaErrors(cudaFree(device_result));
  }

  fclose(outputfile);
  timer();

  /*終了処理*/
  cudaThreadExit();
  exit(0);
}

__global__ void thread_num(unsigned long long int  *device_result , unsigned long long int cycle){
  /*スレッドIDの割り当て*/
  /* メモ : 
     blockDim = block size , 
     threadIdx = 0~blockDim-1 , 
     blockIdx = 0~grid size-1 , 
     max thread = blockDim * max blockIdx + max threadIdx 
  */
  unsigned long long int  thread_idx = threadIdx.x+blockDim.x*blockIdx.x;
  unsigned long long int  thread_idy = threadIdx.y+blockDim.y*blockIdx.y;
  unsigned long long int  thread_idz = threadIdx.z+blockDim.z*blockIdx.z;
  unsigned long long int thread_id = ( blockDim.x * (Grid_x - 1) + blockDim.x ) * ( blockDim.y * (Grid_y - 1) + blockDim.y ) * thread_idz + ( blockDim.x * (Grid_x - 1) + blockDim.x ) * thread_idy + thread_idx ;
  unsigned long long int dev = 0 ;
  int flag = 0;

  /*素数判定*/
  if ( thread_id == 1 ){
    device_result[thread_id] = 0;
  }else if ( thread_id == 2 ){
    device_result[thread_id] = 2;
  }else if ( thread_id % 2 == 0 ){ 
    device_result[thread_id] = 0;
  }else{
    dev = 3;
    while ( ( dev * dev ) <= thread_id + cycle * (Grid_x * Grid_y * Grid_z) * (Block_x * Block_y * Block_z)){
      if ( (thread_id + cycle * (Grid_x * Grid_y * Grid_z) * (Block_x * Block_y * Block_z) ) % dev == 0 ){ 
    flag=1;
    break;
      }
      dev += 2;
    }
    if (flag == 0){
      device_result[thread_id] = thread_id + cycle * (Grid_x * Grid_y * Grid_z) * (Block_x * Block_y * Block_z);
    }else if (flag == 1){
      device_result[thread_id] = 0;
    }
  }
}

説明

一度の計算で、unsigned long long int型で1024x1024x16x8個分の領域をグラフィックボード側のメモリ空間に用意。
スレッドごとに自スレッドの番号が素数かどうかを判定する。
1サイクルで各スレッドは1つの数の素数判定をすることになる。

引数で指定した数に達するまで1サイクル1024x1024x16x8個の計算を繰り返す。

※テクスチャメモリなどの大容量のメモリを使いこなせると1サイクルで全計算をできそうだが、まだ初心者のためメモリの使い方は今後の課題。

ビルドは下記のとおり。

nvcc -m64 prime_cuda.cu -o prime_cuda.exe -O3

実行例

500000000までの範囲で素数判定を行う。
1024x1024x16x8の倍数で500000000を超える最小数まで計算される。
この場合は0から536870911までとなる。
素数リストはprime_num.txtに書き込まれる。

./prime_cuda.exe 500000000

Range : 0 - 134217727
Processing time: 21444.285156 (msec)
Memory copy time: 175.316101 (msec)
Now Writing...
Range : 134217728 - 268435455
Processing time: 38330.570312 (msec)
Memory copy time: 174.686310 (msec)
Now Writing...
Range : 268435456 - 402653183
Processing time: 48694.289062 (msec)
Memory copy time: 174.285660 (msec)
Now Writing...
Range : 402653184 - 536870911
Processing time: 58622.562500 (msec)
Memory copy time: 174.144485 (msec)
Now Writing...

素数リスト

2
3
5
7
...
(略)
...
536870849
536870869
536870879
536870909

実行時間

下図のように4200000000を超えたあたりで計算時間が増加している。
これはint(32bit)とlong long int(64bit)の境界と考えられる。
※64bitでビルドしているはずなのだが。

prime.png

Cとの比較

prime_cuda.cuと同じ方法で素数リストを作成する。
openMPを用いて、マルチスレッド処理を行う。

prime_openmp.c
#include <stdio.h>
#include <stdlib.h>
#include <omp.h>
#include <math.h>

int main (int argc , char* argv[])
{
  long long num = 0 ;
  int *str;
  int *strcp;

  num = atoll(argv[1]);
  str = (int *)malloc(sizeof(int) * num);
  strcp = (int *)malloc(sizeof(int) * num);
  for (long long i = 0 ; i < num ; i++)
    {
      str[i]=i;
      strcp[i]=0;
    }
#pragma omp parallel for schedule(dynamic)
  for (long long i = 0 ; i < num ; i++)
    {
      if (i < 2)
        {
      str[i]=0;
        }
      else if (i > 2)
    {
      for (long long j = 2 ; j < (long long)pow((double)i,0.5) + 1; j++)
        {
          strcp[i] = str[j];
          if (strcp[i] != 0)
        {
          if((str[i] % strcp[i]) == 0)
            {
              str[i] = 0 ;
              break;
                    }
                }
            }
        }
    }
  for (long long i = 0 ; i < num ; i++)
    {
      if(str[i] != 0)
    {
      printf("%lld\n",str[i]);
        }
    }
  free(str);
  free(strcp);
  return 0;
}
gcc -fopenmp prime_openmp.c -o prime_openmp.out -std=c99 -O3 -lm

比較結果

134217728(CUDAでの計算の1サイクル)までの範囲で比較した結果、
CUDAでは21秒、
Cでは10分程度かかった。
openMPとCUDAとの並列度の差が表れている。

エラトステネスの篩との比較

Cでエラトステネスの篩を実装し、CUDAと速度比較をする。
ここでもopenMPを使う。

prime_eratosthenes.c
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <omp.h>

int main(int argc , char *argv[])
{
  long long i, j;
  long long num = 0;
  long long sq = 0;
  long long *prime;

  num = atoll(argv[1]);
  sq=(long long)pow((double)num , 0.5);
  prime = (long long *)malloc(sizeof(long long) * num);


  for (i = 0; i < num; i++)
    {
      prime[i] = 1;  
    }
  prime[0] = 0;  

  #pragma omp parallel for private(j) schedule(dynamic)
  for (i = 1; i < sq; i++) {
    if (prime[i] == 1) 
      {
    for (j = (i+1); (i+1) * j <= num; j++)
      {
        prime[(i+1) * j - 1] = 0;
      }
      }
  }

  for (i = 0; i < num; i++)
    if (prime[i] == 1)  
      printf("%d\n", i + 1);

  free(prime);

  return (0);
}

同じく134217728までの範囲で実行時間は1秒程度。
CUDAよりも速い。

CUDA版エラトステネスの篩

最後にCUDAでエラトステネスの篩を実装し、
Cで実装したエラトステネスの篩の実行時間(1秒)
と比較する。

prime_cuda_eratosthenes.cu
#include <stdio.h>
#include <stdlib.h>
#include <cuda.h>
#include <helper_cuda.h>
#include <helper_functions.h>
#include <time.h>
#include <math.h>
#include <omp.h>

#define Grid_x 256
#define Grid_y 512
#define Grid_z 1
#define Block_x 32
#define Block_y 32
#define Block_z 1


__global__ void thread_num(unsigned int  *device_result);

/* timer */
int timer(void){
  time_t now = time(NULL);
  struct tm *pnow = localtime(&now);
  char buff[128]="";
  sprintf(buff,"%d:%d:%d",pnow->tm_hour,pnow->tm_min,pnow->tm_sec);
  printf("%s\n",buff);
  return 0;
}

int main(int argc, char** argv){

  FILE *outputfile;

  outputfile = fopen("./prime_num.txt", "w");
  timer();

  /*ブロックサイズとグリッドサイズの設定*/
  dim3 grid(Grid_x,Grid_y,Grid_z);
  dim3 block(Block_x,Block_y,Block_z);

  /*ホスト側の変数設定*/
  unsigned int  mem_size = (Grid_x * Grid_y * Grid_z) * (Block_x * Block_y * Block_z);
  unsigned int  *host_result;

  host_result = (unsigned int *)malloc(mem_size * sizeof(unsigned int));
#pragma omp parallel for 
  for (unsigned int i = 0 ; i < mem_size ; i++)
    {
      host_result[i] = i ;
    }

  /*デバイス側の変数設定*/
  unsigned int  *device_result;

  /*デバイスメモリ領域の確保*/
  checkCudaErrors(cudaMalloc((void**)&device_result, mem_size * sizeof(unsigned int)));

  /*ホスト側からのメモリ転送*/
  checkCudaErrors(cudaMemcpy(device_result, host_result, mem_size * sizeof(unsigned int) , cudaMemcpyHostToDevice));

  /*タイマーを作成して計測開始*/
  cudaEvent_t start;
  cudaEvent_t stop;
  checkCudaErrors(cudaEventCreate(&start));
  checkCudaErrors(cudaEventCreate(&stop));
  checkCudaErrors(cudaEventRecord(start, NULL));
  printf("Range : 0 - %u\n",mem_size);

  /*カーネルの起動*/
  thread_num<<<grid , block>>>(device_result);
  cudaThreadSynchronize();

  /*タイマーを停止し実行時間を表示*/
  checkCudaErrors(cudaEventRecord(stop, NULL));
  checkCudaErrors(cudaEventSynchronize(stop));
  float msecTotal = 0.0f;
  checkCudaErrors(cudaEventElapsedTime(&msecTotal, start, stop));
  printf("Processing time: %f (msec)\n", msecTotal);

  /*再度タイマー開始*/
  checkCudaErrors(cudaEventRecord(start, NULL));

  /*結果の領域確保とデバイス側からのメモリ転送*/
  checkCudaErrors(cudaMemcpy(host_result, device_result, mem_size * sizeof(unsigned int) , cudaMemcpyDeviceToHost));

  /*タイマーを停止し実行時間を表示*/
  checkCudaErrors(cudaEventRecord(stop, NULL));
  checkCudaErrors(cudaEventSynchronize(stop));
  checkCudaErrors(cudaEventElapsedTime(&msecTotal, start, stop));
  printf("Memory copy time: %f (msec)\n", msecTotal);

  for( unsigned int  l = 0; l < mem_size; l++ )
  { 
    if (host_result[l] != 0)
    {
      fprintf(outputfile, "%u\n", host_result[l]);
    }
  }
  fclose(outputfile);
  /*ホスト・デバイスメモリの開放*/
  free(host_result);
  checkCudaErrors(cudaFree(device_result));
  timer();

  /*終了処理*/
  cudaThreadExit();
  exit(0);
}

__global__ void thread_num(unsigned int  *device_result)
{
  /*スレッドIDの割り当て*/
  /* メモ : 
     blockDim = block size , 
     threadIdx = 0~blockDim-1 , 
     blockIdx = 0~grid size-1 , 
     max thread = blockDim * max blockIdx + max threadIdx 
  */
  unsigned int  thread_idx = threadIdx.x+blockDim.x*blockIdx.x;
  unsigned int  thread_idy = threadIdx.y+blockDim.y*blockIdx.y;
  unsigned int  thread_idz = threadIdx.z+blockDim.z*blockIdx.z;
  unsigned int thread_id = ( blockDim.x * (Grid_x - 1) + blockDim.x ) * ( blockDim.y * (Grid_y - 1) + blockDim.y ) * thread_idz + ( blockDim.x * (Grid_x - 1) + blockDim.x ) * thread_idy + thread_idx ;

  unsigned int mem_size = (Grid_x * Grid_y * Grid_z) * (Block_x * Block_y * Block_z);

  /*素数判定*/

  if (thread_id < 2)
    {
      device_result[thread_id] = 0;
    }
  else
    {
      for (unsigned int i = 2 ; i < (mem_size/thread_id) ; i++)
        {
           device_result[thread_id * i] = 0;
        }
    }
}

Range : 0 - 134217728
Processing time: 2741.332764 (msec)
Memory copy time: 41.863617 (msec)

計算時間は2.7秒程度。
CPUでの演算の方が速い。
エラトステネスの篩は並列度の効果が薄くなる。
そのためコア単体の性能が大きく影響していると考えられる。