簡単なCUDAプログラミングを試してみます。

【入門】CUDAプログラミングを試してみる

CUDAプログラミングの入門として、簡単なプログラムの実装例と、その実行例をご紹介していきます。挙動を確認することで、CUDAプログラミングへの理解を深められたらと思います。

0.環境構築

今回のプログラムを動かすためには、以下をインストールする必要があります。
  • gcc
  • cmake
  • nvidia-cuda-toolkit
  • NVIDIA Nsight Systems) ←プロファイルするためのツール

(Dockerを使う場合)

########## Pull ##########
FROM nvidia/cuda:11.1.1-base-ubuntu20.04
########## Non-interactive ##########
ENV DEBIAN_FRONTEND=noninteractive
########## CUDA ##########
RUN apt-get update && \
	apt-get install -y \
		build-essential \
		cmake \
		nvidia-cuda-toolkit && \
	apt-get remove -y gcc && \
	ln -s /usr/bin/gcc-8 /usr/bin/gcc && \
	ln -s /usr/bin/g++-8 /usr/bin/g++ && \
	ln -s /usr/bin/gcc-8 /usr/bin/cc && \
	ln -s /usr/bin/g++-8 /usr/bin/c++ && \
	cd ~/ && \
	wget https://developer.download.nvidia.com/devtools/repos/ubuntu2004/amd64/nsight-systems-2022.2.1_2022.2.1.31-1_amd64.deb && \
	apt-get install -y ./nsight-systems-2022.2.1_2022.2.1.31-1_amd64.deb

1.Hello World!をprintf

まずは、GPUを使って、並列で「Hello World!」と表示させてみます。

実装

#include <stdio.h>

__global__ void helloWorld()
{
	printf("Hello World!\n");
}

int main()
{
	const size_t num_blocks = 1;
	const size_t num_threads_per_block = 3;

	helloWorld<<<num_blocks, num_threads_per_block>>>();
	cudaDeviceSynchronize();
}

ビルド

$ nvcc ../src/hello_world.cu -o hello_world
なお、他のプログラムも同様にビルドできるため、これ以降ではビルド方法の記載を省略します。

実行

$ ./hello_world
なお、他のプログラムも同様に実行できるため、これ以降では実行方法の記載を省略します。

出力

Hello World!
Hello World!
Hello World!

NG例

以下のようにcudaDeviceSynchronize()を書かないと、同期される前にプログラムが終了してしまいます。
#include <stdio.h>

__global__ void helloWorld()
{
	printf("Hello World!\n");
}

int main()
{
	const size_t num_blocks = 1;
	const size_t num_threads_per_block = 3;

	helloWorld<<<num_blocks, num_threads_per_block>>>();
	// cudaDeviceSynchronize();
}
実行結果は、以下のように何も出力されません。
 

2.インデックスをprintf

並列プロセスの各インデックスを表示させてみます。

実装

#include <stdio.h>

__global__ void printIndicies()
{
	int index = threadIdx.x + blockIdx.x * blockDim.x;
	printf("index = %d\n", index);
}

int main()
{
	const size_t num_blocks = 2;
	const size_t num_threads_per_block = 3;

	printIndicies<<<num_blocks, num_threads_per_block>>>();
	cudaDeviceSynchronize();
}

出力

index = 0
index = 1
index = 2
index = 3
index = 4
index = 5

3.配列の各要素をprintf

配列の各要素にアクセスして値を表示させてみます。

実装

#include <stdio.h>

void printArrayCPU(int* arr, size_t num_elements)
{
	for(size_t i = 0; i < num_elements; i++){
		printf("arr[%d] = %d\n", i, arr[i]);
	}
}

__global__ void printArray(int* arr, size_t num_elements)
{
	int index = threadIdx.x + blockIdx.x * blockDim.x;
	printf("arr[%d] = %d\n", index, arr[index]);
}

int main()
{
	const size_t num_blocks = 2;
	const size_t num_threads_per_block = 3;
	const size_t num_elements = num_blocks * num_threads_per_block;

	int* arr;
	size_t bytes = num_elements * sizeof(int);
	cudaMallocManaged(&arr, bytes);

	for(size_t i = 0; i < num_elements; i++)   arr[i] = i;

	printf("----- CPU -----\n");
	printArrayCPU(arr, num_elements);

	printf("----- GPU -----\n");
	printArray<<<num_blocks, num_threads_per_block>>>(arr, num_elements);
	cudaDeviceSynchronize();
}

出力

----- CPU -----
arr[0] = 0
arr[1] = 1
arr[2] = 2
arr[3] = 3
arr[4] = 4
arr[5] = 5
----- GPU -----
arr[0] = 0
arr[1] = 1
arr[2] = 2
arr[3] = 3
arr[4] = 4
arr[5] = 5

NG例

以下のように、cudaMallocManaged()でメモリを確保しておかないと、GPUでそのメモリにアクセスできません。
#include <stdio.h>

void printArrayCPU(int* arr, size_t num_elements)
{
	for(size_t i = 0; i < num_elements; i++){
		printf("arr[%d] = %d\n", i, arr[i]);
	}
}

__global__ void printArray(int* arr, size_t num_elements)
{
	int index = threadIdx.x + blockIdx.x * blockDim.x;
	printf("arr[%d] = %d\n", index, arr[index]);
}

int main()
{
	const size_t num_blocks = 2;
	const size_t num_threads_per_block = 3;
	const size_t num_elements = num_blocks * num_threads_per_block;

	int* arr;
	size_t bytes = num_elements * sizeof(int);
	arr = (int *)malloc(bytes);

	for(size_t i = 0; i < num_elements; i++)   arr[i] = i;

	printf("----- CPU -----\n");
	printArrayCPU(arr, num_elements);

	printf("----- GPU -----\n");
	printArray<<<num_blocks, num_threads_per_block>>>(arr, num_elements);
	cudaDeviceSynchronize();
}
実行結果は、以下のように、GPUプロセスは出力されません。
----- CPU -----
arr[0] = 0
arr[1] = 1
arr[2] = 2
arr[3] = 3
arr[4] = 4
arr[5] = 5
----- GPU -----

4.Grid-stride loopで配列の各要素をprintf

Grid-stride loopというテクニックで、配列の各要素にアクセスして値を表示させてみます。

実装

#include <stdio.h>

__global__ void printArray(int* arr, size_t num_elements)
{
	int index = threadIdx.x + blockIdx.x * blockDim.x;
	int grid_stride = gridDim.x * blockDim.x;

	for(int i = index; i < num_elements; i += grid_stride){
		printf("arr[%d] = %d\n", i, arr[i]);
	}
}

int main()
{
	const size_t num_blocks = 2;
	const size_t num_threads_per_block = 3;
	const size_t num_elements = 10;

	int* arr;
	size_t bytes = num_elements * sizeof(int);
	cudaMallocManaged(&arr, bytes);

	for(size_t i = 0; i < num_elements; i++)   arr[i] = i;

	printArray<<<num_blocks, num_threads_per_block>>>(arr, num_elements);
	cudaDeviceSynchronize();
}

出力

arr[0] = 0
arr[1] = 1
arr[2] = 2
arr[3] = 3
arr[4] = 4
arr[5] = 5
arr[6] = 6
arr[7] = 7
arr[8] = 8
arr[9] = 9

5.デバイス情報をprintf

デバイス情報を表示させてみます。

実装

#include <stdio.h>

__global__ void helloWorld()
{
	printf("Hello World!\n");
}

int main()
{
	int device_id;
	cudaGetDevice(&device_id);

	cudaDeviceProp props;
	cudaGetDeviceProperties(&props, device_id);

	printf("Device ID: %d\n", device_id);
	printf("Number of SMs: %d\n", props.multiProcessorCount);
	printf("Compute Capability Major: %d\n", props.major);
	printf("Compute Capability Minor: %d\n", props.minor);
	printf("Warp Size: %d\n", props.warpSize);
}

出力

Device ID: 0
Number of SMs: 28
Compute Capability Major: 8
Compute Capability Minor: 6
Warp Size: 32

6.グリッドサイズの決め方

グリッドサイズは、SMの数の倍の数が良いとされています。

実装

#include <stdio.h>

__global__ void hoge()
{
	int hoge = 1 + 1;
}

int main()
{
	/*device query*/
	int device_id;
	cudaGetDevice(&device_id);
	int num_sm;
	cudaDeviceGetAttribute(&num_sm, cudaDevAttrMultiProcessorCount, device_id);

	printf("num_sm = %d\n", num_sm);

	/*Grid sizes that are multiples of the number of available SMs can increase performance*/
	const size_t num_blocks = 32 * num_sm;
	const size_t num_threads_per_block = 256;

	hoge<<<num_blocks, num_threads_per_block>>>();
	cudaDeviceSynchronize();
}

出力

num_sm = 28

7.プリフェッチ

プリフェッチを記述しておくことで、CPU-GPU間のデータコピーを素早く実行できます。プリフェッチを記述しないと、データの格納先を探索する時間がかかってしまうのです。

実装

#include <stdio.h>

void printArrayCPU(int* arr, size_t num_elements)
{
	for(size_t i = 0; i < num_elements; i++){
		printf("arr[%d] = %d\n", i, arr[i]);
	}
}

__global__ void printArray(int* arr, size_t num_elements)
{
	int index = threadIdx.x + blockIdx.x * blockDim.x;
	printf("arr[%d] = %d\n", index, arr[index]);
}

int main()
{
	const size_t num_blocks = 2;
	const size_t num_threads_per_block = 3;
	const size_t num_elements = num_blocks * num_threads_per_block;

	int* arr;
	size_t bytes = num_elements * sizeof(int);
	cudaMallocManaged(&arr, bytes);

	for(size_t i = 0; i < num_elements; i++)   arr[i] = i;

	int device_id;
	cudaGetDevice(&device_id);
	cudaMemPrefetchAsync(arr, bytes, device_id);

	printf("----- GPU -----\n");
	printArray<<<num_blocks, num_threads_per_block>>>(arr, num_elements);
	cudaDeviceSynchronize();

	printf("----- CPU -----\n");
	cudaMemPrefetchAsync(arr, bytes, cudaCpuDeviceId);
	printArrayCPU(arr, num_elements);
}

出力

----- GPU -----
arr[0] = 0
arr[1] = 1
arr[2] = 2
arr[3] = 3
arr[4] = 4
arr[5] = 5
----- CPU -----
arr[0] = 0
arr[1] = 1
arr[2] = 2
arr[3] = 3
arr[4] = 4
arr[5] = 5

プロファイル

実行ファイルをプロファイルしてみると、各処理の時間を見ることができます。プリフェッチすることで、CPUとGPU間の転送を早くできます。
$ nsys profile --stats=true ~/cuda_cpp_tutorial/build/print_array_with_prefetch

----- GPU -----
arr[0] = 0
arr[1] = 1
arr[2] = 2
arr[3] = 3
arr[4] = 4
arr[5] = 5
----- CPU -----
arr[0] = 0
arr[1] = 1
arr[2] = 2
arr[3] = 3
arr[4] = 4
arr[5] = 5
Generating '/tmp/nsys-report-d472.qdstrm'
[1/8] [========================100%] report1.nsys-rep
[2/8] [========================100%] report1.sqlite
[3/8] Executing 'nvtxsum' stats report
SKIPPED: /root/cuda_cpp_tutorial/build/report1.sqlite does not contain NV Tools Extension (NVTX) data.
[4/8] Executing 'osrtsum' stats report
Operating System Runtime API Statistics:
 Time (%)  Total Time (ns)  Num Calls   Avg (ns)    Med (ns)   Min (ns)  Max (ns)  StdDev (ns)       Name     
 --------  ---------------  ---------  ----------  ----------  --------  --------  -----------  --------------
     70.2        231463209         17  13615482.9  10060535.0      3250  68769206   17963544.0  poll          
     14.2         46707197        454    102879.3     11345.0      1010  13064068     743131.1  ioctl         
      9.6         31819233         13   2447633.3     26640.0     10670  20429791    6081350.0  sem_timedwait 
      5.5         18114990         30    603833.0      3790.0      1210  17910430    3268710.9  fopen         
      0.2           752665         27     27876.5      3110.0      2150    466242      88468.7  mmap64        
      0.1           389781         44      8858.7      8125.5      3300     25340       4368.9  open64        
      0.1           176071          5     35214.2     42150.0     18800     45871      11728.4  pthread_create
      0.0           111910         18      6217.2      5510.0      1270     28480       6326.0  mmap          
      0.0            41450          1     41450.0     41450.0     41450     41450          0.0  fgets         
      0.0            32010          6      5335.0      5040.0      2450      7900       2152.6  open          
      0.0            24480          7      3497.1      3190.0      2390      5690       1155.1  munmap        
      0.0            23310         11      2119.1      2250.0      1090      3960        893.8  write         
      0.0            20590          9      2287.8      1250.0      1070      6800       2079.4  fcntl         
      0.0            19100          9      2122.2      1710.0      1080      3980       1060.0  fclose        
      0.0            17040          8      2130.0      1935.0      1350      3320        663.7  read          
      0.0            15720          2      7860.0      7860.0      3500     12220       6166.0  socket        
      0.0            13250          2      6625.0      6625.0      5820      7430       1138.4  fread         
      0.0            11230          2      5615.0      5615.0      1420      9810       5932.6  fwrite        
      0.0             9390          1      9390.0      9390.0      9390      9390          0.0  connect       
      0.0             8610          2      4305.0      4305.0      1040      7570       4617.4  fflush        
      0.0             7480          1      7480.0      7480.0      7480      7480          0.0  pipe2         
      0.0             2380          1      2380.0      2380.0      2380      2380          0.0  bind          
[5/8] Executing 'cudaapisum' stats report
CUDA API Statistics:
 Time (%)  Total Time (ns)  Num Calls   Avg (ns)     Med (ns)    Min (ns)   Max (ns)   StdDev (ns)          Name         
 --------  ---------------  ---------  -----------  -----------  ---------  ---------  -----------  ---------------------
     99.8        190180906          1  190180906.0  190180906.0  190180906  190180906          0.0  cudaMallocManaged    
      0.1           216510          2     108255.0     108255.0      21790     194720     122280.0  cudaMemPrefetchAsync 
      0.0            56570          1      56570.0      56570.0      56570      56570          0.0  cudaDeviceSynchronize
      0.0            46641          1      46641.0      46641.0      46641      46641          0.0  cudaLaunchKernel     
[6/8] Executing 'gpukernsum' stats report
CUDA Kernel Statistics:
 Time (%)  Total Time (ns)  Instances  Avg (ns)  Med (ns)  Min (ns)  Max (ns)  StdDev (ns)                Name              
 --------  ---------------  ---------  --------  --------  --------  --------  -----------  --------------------------------
    100.0            43168          1   43168.0   43168.0     43168     43168          0.0  printArray(int *, unsigned long)
[7/8] Executing 'gpumemtimesum' stats report
CUDA Memory Operation Statistics (by time):
 Time (%)  Total Time (ns)  Count  Avg (ns)  Med (ns)  Min (ns)  Max (ns)  StdDev (ns)              Operation            
 --------  ---------------  -----  --------  --------  --------  --------  -----------  ---------------------------------
     61.1             1856      1    1856.0    1856.0      1856      1856          0.0  [CUDA Unified Memory memcpy HtoD]
     38.9             1184      1    1184.0    1184.0      1184      1184          0.0  [CUDA Unified Memory memcpy DtoH]
[8/8] Executing 'gpumemsizesum' stats report
CUDA Memory Operation Statistics (by size):
 Total (MB)  Count  Avg (MB)  Med (MB)  Min (MB)  Max (MB)  StdDev (MB)              Operation            
 ----------  -----  --------  --------  --------  --------  -----------  ---------------------------------
      0.004      1     0.004     0.004     0.004     0.004        0.000  [CUDA Unified Memory memcpy DtoH]
      0.004      1     0.004     0.004     0.004     0.004        0.000  [CUDA Unified Memory memcpy HtoD]
Generated:
    /root/cuda_cpp_tutorial/build/report1.nsys-rep
    /root/cuda_cpp_tutorial/build/report1.sqlite

8.エラーチェック

CUDAプログラミングにおけるエラーチェックをご紹介します。

実装

#include <stdio.h>

__global__ void printIndicies()
{
	int index = threadIdx.x + blockIdx.x * blockDim.x;
	printf("index = %d\n", index);
}

int main()
{
	const size_t num_blocks = 1;
	const size_t num_threads_per_block = 1025; // > 1024

	printIndicies<<<num_blocks, num_threads_per_block>>>();
	cudaDeviceSynchronize();

	cudaError_t err = cudaGetLastError();
	if(err != cudaSuccess)  printf("Error: %s\n", cudaGetErrorString(err));
}

出力

Error: invalid configuration argument

NG例

以下のように、エラーチェックを実装しないと、実行エラーに気づくことができません。
#include <stdio.h>

__global__ void printIndicies()
{
	int index = threadIdx.x + blockIdx.x * blockDim.x;
	printf("index = %d\n", index);
}

int main()
{
	const size_t num_blocks = 1;
	const size_t num_threads_per_block = 1025; // > 1024

	printIndicies<<<num_blocks, num_threads_per_block>>>();
	cudaDeviceSynchronize();
}
実行結果は、以下のように何も出力されません。
 

さいごに

CUDAプログラミングの基礎的な挙動を確認しました。 参考になれば幸いです。
Ad.