GPUプログラミング・基礎編
東京工業大学学術国際情報センター
この資料は，本選参加者に対する事前資料として配布するものです．
できれば目を通しておいて下さい．ただし，わからなくても心配無用です．本選初日に，この資料を使って説明会をしますので．もちろん，CUDA プログラミングを事前に練習しておく必要もありません．

ＧＰＵプログラミング入門：行列積計算を例題に行列積とは：行列×行列の計算
AB

ここでは n x n の正方行列を対象とする

C j行

×nCij

i行

n

計算の大まかな流れ：
for (i = 0; i < n; i++) { // Cの第i行に注目
for (j = 0; j < n; j++) { // Cの第j列に注目
for (k = 0; k < n; k++) { // 総和計算のループ
c[i][j] += a[i][k]*b[k][j];} } }GPUの計算速度の威力
n=2048の場合の計算時間をTSUBAME2.5上で測ってみた• CPU版 (C言語で書かれたmm.c)
Æ 10.1秒かかった• GPU版 (「CUDA」で書かれたmm_gpu.cu)→ 0.32秒。CPU版より約31倍も速い！！

なぜ？ GPU は１つの中に 2688 個もの演算装置（計算コア）がいる・・・・・・

ただ，2688 倍
はちょっと無理1 つ 1 つはCPU よりかなり劣る

これらが上手く分業して計算できれば，すごく速くなる上手な分業化が必要！
各スレッド（≒計算コア）が Cij を受け持つと．．．
A

B

C j行

×n

Ciji行

n

計算の大まかな流れ：
for (i = 0; i < n; i++) { // Cの第i行に注目
for (j = 0; j < n; j++) { // Cの第j列に注目
for (k = 0; k < n; k++) { // 総和計算のループ
c[i][j] += a[i][k]*b[k][j];各スレッドでの計算は} } }これだけになる・・・・・・

これが実際のプログラム
正確には mm_gpu.cu というプログラム内の GPU 用関数__global__ void mm_gpu(double *A, double *B, double *C, int n){
int i, j, k;
i = blockIdx.y * blockDim.y + threadIdx.y;j = blockIdx.x * blockDim.x + threadIdx.x;// 自分の背番号から担当する (i,j) を決める
if (i >= n || j >= n) return; // 行列からはみ出す部分は計算しない
for (k = 0; k < n; k++) { // 総和を計算する部分
C[i*n+j] += A[i*n+k] * B[k*n+j]; // C[i][k] += A[i][k] * B[k][j]に相当}}

まずは，これを理解しよう！

1. CUDA プログラムの概要
CUDA プログラムとは，CPU+GPU を動かすためのプログラム．C 言語の拡張版と思えばよい．
以下が，行列積のための CUDA プログラムmm_gpu.cu である．
// A, Bの内容を設定し、Cをゼロクリア(略)#include <stdio.h>#include <stdlib.h>// A, B, Cのためにデバイスメモリを確保
cudaMalloc((void**)&DA, sizeof(double)*n*n);cudaMalloc((void**)&DB, sizeof(double)*n*n);cudaMalloc((void**)&DC, sizeof(double)*n*n);__global__ void mm_gpu(double *A, double *B, double *C, int n){
int i, j, k;
i = blockIdx.y * blockDim.y + threadIdx.y;j = blockIdx.x * blockDim.x + threadIdx.x;// 自分の背番号から担当する (i,j) を決める
if (i >= n || j >= n) return; // 行列からはみ出す部分は計算しない// A, B, Cの内容を、ホストメモリからデバイスメモリへコピー
cudaMemcpy(DA, A, sizeof(double)*n*n, cudaMemcpyHostToDevice);
cudaMemcpy(DB, B, sizeof(double)*n*n, cudaMemcpyHostToDevice);
cudaMemcpy(DC, C, sizeof(double)*n*n, cudaMemcpyHostToDevice);for (k = 0; k < n; k++) { // 総和を計算する部分
C[i*n+j] += A[i*n+k] * B[k*n+j]; // C[i][k] += A[i][k] * B[k][j]に相当}// GPUカーネル関数を呼び出す!! 約n*n個のスレッドを使う
mm_gpu<<<dim3((n+BS‐1)/BS, ((n+BS‐1)/BS)), dim3(BS, BS)>>>(DA, DB, DC, n);}
int main(int argc, char *argv[]){
int i, j, n;double *A, *B, *C;double *DA, *DA, *DC;n = atoi(argv[1]); // 行列の大きさ// A, B, Cのためにホストメモリを確保A = (double *)malloc(sizeof(double)*n*n);
B = (double *)malloc(sizeof(double)*n*n);
C = (double *)malloc(sizeof(double)*n*n);// 結果のCを、デバイスメモリからホストメモリへコピー
cudaMemcpy(C, DC, sizeof(double)*n*n, cudaMemcpyDeviceToHost);// Cを出力などに利用(略)return 0;}※TSUBAMEにある本物mm_gpu.cuは時間測定なども含む

1. CUDA プログラムの概要
CUDA プログラムの構成
ホスト関数+ GPUカーネル関数

二種類の関数の定義が xxx.cu ファイル内に混ざっている• ホスト関数– CPU上で実行される関数– ほぼ通常のC言語。main関数から処理がはじまる– GPUに対してデータ転送やGPUカーネル関数呼び出しを実行• GPUカーネル関数– GPU上で実行される関数– ホストプログラムから呼び出されて実行– (単にカーネル関数と呼ぶ場合も)

xxx.cuファイルは、nvccコマンドでコンパイルする

1. CUDA プログラムの概要
CPUGPU•

main()から開始し、当初はCPUだけ動く– GPUは、CPUから処理を依頼された時だ
け動く– GPU上で動く関数 = GPUカーネル関数

call•

CPUとGPU間ではメモリは別々– CPU から GPU にデータを転送して実行させるのが基本

return•GPU上では多数のスレッド達が動く．
スレッド達の間では同じメモリが見える(GPU上の)スレッドたち8

1. CUDA プログラムの概要
CUDAプログラムの流れの例
CPU上
GPU上
int main() {GPUカーネル関数の印

GPU側メモリにデータ用領域を確保
入力データをGPUへ転送__global__ void kernel_func() {カーネル関数をreturn; 実行

GPUカーネル関数を呼び出し出力をCPU側メモリへ転送}

}

入力

入力

出力

出力

CPU側メモリ（ホストメモリ）

GPU側メモリ（デバイスメモリ）

この２種類のメモリの区別は忘れずに！

2. mm_gpu.cu の詳細 (CPU側)
int main(int argc, char *argv[]){
int i, j, n;double *A, *B, *C; // ホストメモリ用のポインタ
double *DA, *DA, *DC; // デバイスメモリ用のポインタ
n = atoi(argv[1]); // 行列の大きさ// A, B, Cのためにホストメモリを確保
A = (double *)malloc(sizeof(double)*n*n);
B = (double *)malloc(sizeof(double)*n*n);
C = (double *)malloc(sizeof(double)*n*n);// A, Bの内容を設定し、Cをゼロクリア
for (i = 0; i < n; i++) {for (j = 0; j < n; j++) {A[i*n+j] = 3.0; // ここでのA, Bの内容はいい加減
B[i*n+j] = 0.1; //に決めた．（あくまで例なので）
C[i*n+j] = 0.0;
} }実行はmain関数から、CPU上で始まる← これは後で使う
行列A, B, Cのためのメモリ領域をmallocで準備
これらはホストメモリ上に確保される
このプログラムでは、二次元配列でなく「サイズn*nの一次元配列」にした一次元になっているため、たとえばCijをアクセスした
いときには、C[i*n+j]とする

2. mm_gpu.cu の詳細(CPU側)// A, B, Cのためにデバイスメモリを確保
cudaMalloc((void**)&DA, sizeof(double)*n*n);cudaMalloc((void**)&DB, sizeof(double)*n*n);cudaMalloc((void**)&DC, sizeof(double)*n*n);
結果のポインタはここに入る
これを実行すると、DA, DB, DCはそれぞれ、「sizeof(double)*n*n」バイトの大きさの、デバイスメモリ上の領域を指す

ADA

BDBCDC
ホストメモリデバイスメモリ

注意! ホスト関数(CPU側)からは、DA, DB, DCの中身には触れない。
DA[i*n+j] = … などしてはダメ！

2. mm_gpu.cu の詳細(CPU側)// A, B, Cの内容を、ホストメモリからデバイスメモリへコピー
cudaMemcpy(DA, A, sizeof(double)*n*n, cudaMemcpyHostToDevice);
cudaMemcpy(DB, B, sizeof(double)*n*n, cudaMemcpyHostToDevice);
cudaMemcpy(DC, C, sizeof(double)*n*n, cudaMemcpyHostToDevice);ADA

B

DBCDC

デバイスメモリ
ホストメモリcudaMemcpyの引数は4つ
1. コピー先
2. コピー元
3. 何バイトコピーしたいか
ここまでで、
4. コピーの種類
やっとGPUが動ける
cudaMemcpyHostToDevice: ホスト→デバイスのとき準備ができた！
cudaMemcpyDeviceToHost: デバイスÆホストのとき

2. mm_gpu.cu の詳細(CPU側)
いよいよGPU上で動く関数を呼び出そう !!
// GPUカーネル関数を呼び出す!! 約n*n個のスレッドを使う// このプログラムではBSは16mm_gpu <<<dim3((n+BS‐1)/BS, ((n+BS‐1)/BS)), dim3(BS, BS)>>>

呼びたい関数名(DA, DB, DC, n);

渡したい引数いくつのスレッドで関数を実行するか

この例はちょっと難しいのですが・・・
結局は、約(n * n)個のスレッドを使います
それぞれのスレッドには、二次元の背番号をつけますnBSBSn

2. mm_gpu.cu の詳細 (GPU側)
GPU上で、超多数のスレッドがmm_gpu関数を実行開始する
GPUカーネル関数の印
自分のY方向の背番号を計算
自分のX方向の背番号を計算__global__ void mm_gpu(double *A, double *B, double *C, int n){
デバイスメモリ上のint i, j, k;
ポインタが来ているはず
i = blockIdx.y * blockDim.y + threadIdx.y;j = blockIdx.x * blockDim.x + threadIdx.x;// 自分の背番号から担当する (i,j) を決める
if (i >= n || j >= n) return; // 行列からはみ出す部分は計算しない

僕の背番号は
i=0, j=0
僕の背番号は
i=5, j=3

2. mm_gpu.cu の詳細(GPU側)
背番号にしたがって、Cijを一個だけ計算する
for (k = 0; k < n; k++) { // 総和を計算する部分
C[i*n+j] += A[i*n+k] * B[k*n+j];}}一次元配列なのでこの表記

AB

C j行×nCijn
参考：
CPUの場合は一人で全体を計算するので3重ループだった=∑

i行
×for (i = 0; i < n; i++) { // Cの第i行に注目
for (j = 0; j < n; j++) { // Cの第j列に注目
for (k = 0; k < n; k++) { // 総和計算のループ
c[i][j] += a[i][k]*b[k][j];} } }2. mm_gpu.cu の詳細(CPU側)
GPUカーネル関数からmain関数へ戻ってきた// 結果のCを、デバイスメモリからホストメモリへコピー
cudaMemcpy(C, DC, sizeof(double)*n*n, cudaMemcpyDeviceToHost);// Cを出力などに利用(略)}デバイスからホスト行きなのでここでは、Cだけをコピーすることにした。A, Bはいらないので

ADA

B

DBCDC
ホストメモリデバイスメモリ3. 次のステップへ「動く」プログラムを作った後は、「速い」プログラムにしたい
ポイントは色々ある<<<○, ○ >>>のスレッド数の調整によって速度が変わる
GPU上で、隣のスレッドどうしは、同時に配列の近い場所にアクセスすると速い (コアレスドアクセス)
GPU上でif文で分岐すると非効率になることがある
不要なcudaMemcpyは削ったほうがよい
詳しくは、次のステップの資料で！