2011年3月3日

GPGPUハンズオンプログラミング演習

株式会社クロスアビリティ
rkoga@x‐ability.jp
3 Mar 2011

Copyright (C) 2011 X-Ability Co.,Ltd.
All rights reserved.

講師： 古賀良太／古川祐貴
取締役•

計算化学ソルバー XA‐CHEM‐SUITEの開発

コンサルティングパートナー•

並列ソフトウェアの開発・ビルド、サーバ販売•

3 Mar 2011

ソフトウェア代理店

Copyright (C) 2011 X-Ability Co.,Ltd.
All rights reserved.

会社紹介
• 社名
株式会社クロスアビリティ (X‐Ability Co.,Ltd) ※役員３名のみ• 業務内容計算科学関連ソフトウェアの開発、販売
フィールドルータの開発、販売• ビジネスモデル
産学連携によるプロダクト開発、ベンダとの連携による販売• 設立
2008年1月• 主な製品XA‐CHEM‐SUITE (XA‐CUDA‐QM etc.), Field Router
3 Mar 2011

Copyright (C) 2011 X-Ability Co.,Ltd.
All rights reserved.

前座
• 対象
C/C++は一通り理解しているがGPGPUは初めての方※対象でない方は“できる限り”個別に対応します• 実践前提初プログラミング言語の学習で座学は意味がない
少し書いてみてから各種ツールの有用性がわかる• 絵は基本的に使わない
ただしグラフィック処理の出力は別（今回はなし）• 参考書籍に書いてあることは出来る限り書かない
GPGPUの障壁をさげることが本講習会の目的
コーディング＆コンパイルの体験と導入Tipsが重要
3 Mar 2011

Copyright (C) 2011 X-Ability Co.,Ltd.
All rights reserved.

• 参考書籍• 参考URL– CUDA Programming Guide– CUDA Occupancy Calculator
http://developer.download.nvidia.com/compute/cuda/CUDA_Occupancy_calculator.xls

3 Mar 2011

Copyright (C) 2011 X-Ability Co.,Ltd.
All rights reserved.

GPGPUとは？
• GPUを用いて汎用計算（科学技術計算など）を行うことをさす。
• 基本PCI Express x16バスにGPUを挿すだけだが、、、– デバイスドライバ（無料）のインストール要– 容量の大きい電源に交換要な場合が多い（電気食い）– マルチGPUの場合、host‐device転送バスが同一バンド幅かチェック要• １つの命令で多数複数スレッドの同時演算が可能– 理論的には相当に高速だが（Fermi coreは価格性能比でCore i7の10倍）、
Host(CPU+mem+HDD)‐GPU転送コストがかかる + 特別な言語(CUDA)でコーディングしなおす必要がある。

• 高速だが容量が小さいon‐chipメモリと低速だが容量が大きい（といっても数GB）external memoryによる構成• 倍精度演算が単精度演算よりダイブ遅い– 倍精度ユニットが高速になったが、ソフトウェアが活用してない場合がある
3 Mar 2011

Copyright (C) 2011 X-Ability Co.,Ltd.
All rights reserved.

CUDAプログラミングモデル• GPUは多数のスレッドを並列実行できる外部演算デバイスとして扱われる• GPUで走るプログラムをカーネルと呼び、カーネルを実行すると、同一のカーネルを実行するスレッドが多数走る• 複数のスレッド群をまとめて、ブロックと呼ぶ。各ブロックは、
ブロック内のスレッドのみがアクセスできる共有メモリを持つ• カーネル実行時に、ブロック数及び各ブロックごとのスレッド数を指定する。概念的には、GPU上でブロック数×スレッド数の個数のスレッドが走ることになる– 実際には、Tesla C2050には448個しかCUDAコア（GPU演算コア）がな
いので、何千・何万のスレッドが同時に走るわけではない– 最適なブロック数・スレッド数というのはケースバイケース• 各スレッドには、ブロックIDとスレッドIDという固有のIDが振られる。これらは3次元配列のindex
3 Mar 2011

Copyright (C) 2011 X-Ability Co.,Ltd.
All rights reserved.

CUDAプログラミングとGPUの関係

CPUスレッドから起動するkernel(grid,
block)の中のblockの説明。GPUThreadはblockの中のWarpという単位
でスケジューリングされ、Warp内のthreadは同じ命令を実行する(SIMT)。
多くのWarpを使うことでメモリ遅延を隠蔽できる。32threads単位で動く。
3 Mar 2011

GlobalメモリはOff-chipなので、Global
memoryからchipに転送する(ホスト機GPU転送コストよりはダイブ低い）。
FermiはL1/L2 cacheが追加され、SMの数が16, SPの数が32になっている。

Copyright (C) 2011 X-Ability Co.,Ltd.
All rights reserved.

実習内容
下記、行列乗算のコードを書くことで、徐々に高速化を実感してもらいます。
1.
2.
3.
4.

普通に書くC++のコード
普通に書くCUDAのコード
チューニングしたCUDAのコード
CUDAのライブラリを使ったコード

時間があれば、、、
5. チューニングしたOpenCLのコードを3と比較
3 Mar 2011

Copyright (C) 2011 X-Ability Co.,Ltd.
All rights reserved.

マシンアクセス方法• Tesla C1060のマシン (kai)$ ssh gpuschoolXXX@133.30.112.247
GPUのメモリキャッシュが効かないマシン• Tesla C2050のマシン (ise)$ ssh gpuschoolXXX@133.30.112.247
GPUのメモリキャッシュが効くマシン※ gpuschoolXXX はアカウント名です

3 Mar 2011

Copyright (C) 2011 X-Ability Co.,Ltd.
All rights reserved.

メイン関数(main.cu)#include <sys/time.h>float* A = new float[M];#include <stdio.h>

float* B = new float[M];const int N = 512;float* C = new float[M];//512 x 512の行列乗算

const int M = N * N;for(int i=0; i<M; i++) A[i] = 1.0f;
for(int j=0; j<M; j++) B[j] = 2.0f;//時間計測用コード
double gettimeofday_sec()

double t1 = gettimeofday_sec();{
struct timeval tv;

matmul(A,B,C,N);

gettimeofday(&tv, NULL);

double t2 = gettimeofday_sec();// この行で関数を呼ぶ

return tv.tv_sec + (double)tv.tv_usec*1e-6;
float fAns= 1.0f * 2.0f * N;}float fDiff = 0.0f;
void matmul(float* A, float* B, float* C, int N);for(int k=0; k<M; k++) {int main(void)

float f = C[k] - fAns;{

fDiff += f * f;}// 以下3行はCUDAのコードのみで必要 ※CUDAランタイムライブラ
リの初期化する時間を節約する

printf("Time = %10.30f¥n", t2 - t1);

float* p;printf("Accuracy : %f¥n", sqrt( fDiff / M ) );cudaMalloc((void**)&p, sizeof(float));return 0;

cudaFree(p);3 Mar 2011}

Copyright (C) 2011 X-Ability Co.,Ltd.
All rights reserved.

普通に書くC++のコード(naive_cpu.cpp)#include <omp.h>
void matmul(float *A, float *B, float *C, int N){#pragma omp parallel for//※OpenMPによるスレッド並列

for(int i=0; i<N; i++){for(int j=0; j<N; j++){float sum = 0.0f;
for(int k=0; k<N; k++){sum += A[i*N+k]*B[k*N+j];}

コンパイルと実行

C[i*N+j] = sum;$ nvcc –O3 -Xcompiler -fopenmp main.cu naive_cpu.cpp}}$ export OMP_NUM_THREADS=4$ ./a.out}

3 Mar 2011

Copyright (C) 2011 X-Ability Co.,Ltd.
All rights reserved.

普通に書くCUDAのコード(naive_cuda.cu)__global__ void _matmul(float *A, float *B, float *C, int N){
int x = threadIdx.x + blockIdx.x * blockDim.x;

cudaMemcpy(devA, A, sizeof(float)*N*N,
cudaMemcpyHostToDevice);
cudaMemcpy(devB, B, sizeof(float)*N*N,
cudaMemcpyHostToDevice);

int y = threadIdx.y + blockIdx.y * blockDim.y;
float sum = 0.0f;//kernel execution

for(int k=0; k<N; k++){dim3 nThreads(16, 16);

sum += A[y*N+k] * B[k*N+x];

dim3 nBlocks(N/16, N/16);}
C[y*N+x] = sum;}_matmul <<< nBlocks, nThreads >>> (devA, devB, devC, N);
cudaMemcpy(C, devC, sizeof(float)*N*N,
cudaMemcpyDeviceToHost);//wrapper for _matmul kernel
void matmul(float *A, float *B, float *C, int N)

cudaFree(devA);{

cudaFree(devB);
float *devA, *devB, *devC;
cudaMalloc((void**)&devA, sizeof(float)*N*N);

cudaFree(devC);}

cudaMalloc((void**)&devB, sizeof(float)*N*N);
cudaMalloc((void**)&devC, sizeof(float)*N*N);

コンパイルと実行$ nvcc –O3 main.cu naive_cuda.cu$ ./a.out

3 Mar 2011

Copyright (C) 2011 X-Ability Co.,Ltd.
All rights reserved.

CUDAの最適化
今回関係あるのは2.と4.のみ（共有メモリ使う＆#pragma unroll）
グローバルメモリアクセスはcoalesce（複数スレッドからのメモリアクセ
スが１回のフェッチになるように）共有メモリ使うときはバンクコンフリクトをしないように

1.
2.

CUDA PROFILEのwarp_serializeで回数が見れる

条件分岐は減らす
loop unrollingは地味に有効__syncthreadsも減らす

3.
4.
5.

Block内のスレッド間を同期しないようにすれば必要なくなる

オフチップメモリのレイテンシの隠蔽

6.

warpを沢山使えば、特定のwarpが演算中に別のwarpが通信できる

etc

3 Mar 2011

Copyright (C) 2011 X-Ability Co.,Ltd.
All rights reserved.

チューニングしたCUDAのコード(cuda_opt.cu)#define blockSize 16//wrapper for _matmul kernel__global__ void _matmul(float *A, float *B, float *C, int N)

void matmul(float *A, float *B, float *C, int N){{
int bx = blockIdx.x;

float *devA, *devB, *devC;

int by = blockIdx.y;cudaMalloc((void**)&devA, sizeof(float)*N*N);

int tx = threadIdx.x;

cudaMalloc((void**)&devB, sizeof(float)*N*N);

int ty = threadIdx.y;cudaMalloc((void**)&devC, sizeof(float)*N*N);

int a = N*blockSize*by; //submatrix adress of Matrix A

cudaMemcpy(devA, A, sizeof(float)*N*N, cudaMemcpyHostToDevice);

int b = blockSize*bx;

cudaMemcpy(devB, B, sizeof(float)*N*N, cudaMemcpyHostToDevice);

float tmp = 0.0f;
for(int i=0; i<N; i+=blockSize){//kernel execution__shared__ float As[blockSize][blockSize];dim3 nThreads(blockSize, blockSize);__shared__ float Bs[blockSize][blockSize];dim3 nBlocks(N/blockSize, N/blockSize);As[ty][tx] = A[a + N*ty + tx];
Bs[ty][tx] = B[b + N*ty + tx];

_matmul <<< nBlocks, nThreads >>> (devA, devB, devC, N);__syncthreads();

cudaMemcpy(C, devC, sizeof(float)*N*N, cudaMemcpyDeviceToHost);#pragma unroll
for(int k=0; k<blockSize; k++){cudaFree(devA);

tmp += As[ty][k] * Bs[k][tx];cudaFree(devB);}__syncthreads();

cudaFree(devC);}

a += blockSize;b += blockSize*N;

コンパイルと実行

int c = N*blockSize*by + blockSize*bx;$ nvcc –O3 main.cu cuda_opt.cu}

C[c + N*ty + tx] = tmp;

$ ./a.out}

3 Mar 2011

Copyright (C) 2011 X-Ability Co.,Ltd.
All rights reserved.

CUDAのライブラリを使ったコード(cublas.cu)#include <cublas.h>//CALL SGEMM

void matmul(float *A, float *B, float *C, int N)

cublasSgemm('N', 'N', N, N, N, 1.0f, devA, N,
devB, N, 0.0f, devC, N);{
float *devA, *devB, *devC;

cublasGetMatrix(N, N, sizeof(*C), devC, N, C, N);//cublasInit();cublasFree(devA);//Allocate Memory

cublasFree(devB);

cublasAlloc(N*N, sizeof(*A), (void**)&devA);

cublasFree(devC);

cublasAlloc(N*N, sizeof(*B), (void**)&devB);}

cublasAlloc(N*N, sizeof(*C), (void**)&devC);//set matrix
cublasSetMatrix(N, N, sizeof(*A), A, N, devA, N);

コンパイルと実行

cublasSetMatrix(N, N, sizeof(*B), B, N, devB, N);$ nvcc –O3 main.cu cublas.cu -lcublas$ ./a.out

3 Mar 2011

Copyright (C) 2011 X-Ability Co.,Ltd.
All rights reserved.

事前計測タイム
512 x 512 1024 x 1024 2048 x 2048
naïve CPUnaïve CUDA
CUDA opt
CUBLAS

121.0
4.83
2.69
2.22

2086
30.4
12.3
7.75

CPU : Intel Core i7 860 @ 2.80GHz

2.37 x 105
232
76.8
34.8
unit : msec

GPU : NVIDIA Geforce GTX580 (Fermi Core)naïve CPU : OpenMP 4threads, CUDA opt : 16 x 16 blocked
結論：cublasのようなライブラリがあればそれを使った方がいいが、ない場合は頑張ってチューニングしましょう
3 Mar 2011

Copyright (C) 2011 X-Ability Co.,Ltd.
All rights reserved.

チューニングしたOpenCLのコード(ocl_MatMul.cl)• サンプル“oclMatrixMul”を改良– 時間の関係で、あらかじめ置いてあるコードを実行して、「チューニングしたCUDAのコード」と時間比較する• h,cpp,cl,Makefileが全て必要です• matrixMul_gold.cppはnaïve_cpu.ccと同様、
oclMatrixMulの比較用コードです

デモで説明します
3 Mar 2011

Copyright (C) 2011 X-Ability Co.,Ltd.
All rights reserved.

ここから座学メイン
1. CUDA_Occupancy_calculator
2. OpenCLのTips
3. Allinea DDT (debugger)
GPGPUコーディング入門マシン

4. 応用例：Amber11(MD)、XA‐CHEM‐SUITEXA‐CHEM‐SUITEの中のXA‐CUDA‐QMはCUDAで量子化学計算を加速するモジュール

3 Mar 2011

Copyright (C) 2011 X-Ability Co.,Ltd.
All rights reserved.

CUDA_Occupancy_calculator（１）3 Mar 2011

Copyright (C) 2011 X-Ability Co.,Ltd.
All rights reserved.

CUDA_Occupancy_calculator（２）• 以下の３つをいじると何がLimitになっているかを示してくれるエクセルマクロ– Threads Per Block– Registers Per Thread– Shared Memory Per Block (bytes)
３つのバランスが重要– ３つのパラメータはnvcc ‐‐ptxas‐option=‐vでも見れる• 100%だから最高の速度が出ているとは限らない– 対象となるアルゴリズムによる律速があるため（レジスタ使用量が異常に多い、共有メモリを多数必要とする、etc）
Occupancyだけで判断はできないが、参考にはなる。
3 Mar 2011

Copyright (C) 2011 X-Ability Co.,Ltd.
All rights reserved.

OpenCLの始め方• CUDA3.1以上のsdkを入れれば入ってる– NVIDIAサイトのOpenCL driver & sdkは不要• 動作Tips– まずliboclUtil.a(liboclUtil_x86_64.a)を作成• これがないとSDK内のサンプルコードがビルドできない
cd $OPENCL_SDK/OpenCL/common($OPENCL_SDK : デフォルトで/usr/local/cuda/sdk/OpenCL)
make– その後、例えばoclDeviceQueryを実行
cd $OPENCL_SDK/OpenCL/src/oclDeviceQuerymake
cd $OPENCL_SDK/OpenCL/bin/linux/release./oclDeviceQuery
3 Mar 2011

Copyright (C) 2011 X-Ability Co.,Ltd.
All rights reserved.

OpenCLの流れ（１）プラットフォーム取得 clGetPlatFormIDs()
デバイス取得 clGetDeviceIDs()
コンテキスト作成 clCreateContext()
コマンドキュー作成
clCreateCommandQueue()
5. プログラム作成
clCreateProgramWithBinary(), or
clCreateProgramWithSource()
6. カーネル作成 clCreateKernel()
1.
2.
3.
4.

3 Mar 2011

Copyright (C) 2011 X-Ability Co.,Ltd.
All rights reserved.

OpenCLの流れ（２）7. バッファオブジェクト作成 clCreateBuffer()
8. バッファ書込 clEnqueueWriteBuffer()
9. カーネル実行 clEnqueueNDRangeKernal()
10. バッファ読込 clEnqueueReadBuffer()
11. OpenCLオブジェクトリリース
clReleaseKernel(kernel), clReleaseProgram(program),
clReleaseMemObject(memobj),
clReleaseCommandQueue(command_queue),clReleaseContext(context)
3 Mar 2011

Copyright (C) 2011 X-Ability Co.,Ltd.
All rights reserved.

Allinea DDT 【デモ】nvcc –O3 cuda_opt.cu –g -G

直感的な
インター
フェース
‐g ‐Gの
コンパイ
ルオプシ
ョン必要
※代理店始
めました
3 Mar 2011

Copyright (C) 2011 X-Ability Co.,Ltd.
All rights reserved.

GPGPUコーディング入門マシン
1. Intel SSE/AVXとCUDAの併用による最高クラスのSIMD高速化コーディングが可能
2. OpenMP/MPIローカル並列マシン環境にてX-terminalによる
直感的かつシー ムレスなデバッグが可能
3. 次世代GPGPU言語であるOpenCLによる開発も可能

皆様に実習でお使い頂いたマシンとほぼ同じ（例）CPU : Intel Zeon 3.33GHz (X5680 6Core)GPGPU : NVIDIA Tesla C2050 x 2
Compiler : Intel Cluster Studio 2011
Intel Composer Xe, MKL, Intel MPI

CentOS 5.5 + GUI Debugger : Allinea DDT (X‐terminal)
3 Mar 2011

Copyright (C) 2011 X-Ability Co.,Ltd.
All rights reserved.

XA‐CHEM‐SUITE• XA‐CUDA‐QM
CUDAで量子化学計算(Quantum Mechanics : QM)を加速するモジュール• XA‐SSE‐QM
SSEで量子化学計算を高速化するモジュール• XA‐AVX‐QM
AVXで量子化学計算を高速化するモジュール

3 Mar 2011

Copyright (C) 2011 X-Ability Co.,Ltd.
All rights reserved.

ハートリー・フォック法のプロセス
2/ N
core
F  H   CaC*a 2      
a



1

 SCF
 H    P       
2



core
 (r)  (r )  (r' )  (r' ) H 
 G
( ab | cd )  
dr ' dr
core

a

b

c

d

| r  r' |Density Matrix : Initial Guessで初期値を作った後、非線形方程式を解いている間（SCFサイクル）アップデートされ続ける。 F C C  SC
Electron Replusion Integral J-matrix : Coulomb Potential J。理屈では一度計算してメインメモリにおけばいいのだが、O N 2  のため少し基底Nが大きくなると置けなくなる。ディスクI/Oは時間がかかり毎回演算するのがリーズナブルとなるが大変。
Electron Replusion Integral K-matrix : O N 4  、HF exchange K。J-matrixと同様の問題を抱えるが、使用レジスタの量が多く、メモリ少のSIMD型への実装が難しい。
3 Mar 2011

Copyright (C) 2011 X-Ability Co.,Ltd.
All rights reserved.

そもそもGPGPUを使って良い場合• 下記条件が揃っていて初めて意味がある– C/CUDAプログラミングに抵抗がなければハードを買うだ
けでよい– 理論、アルゴリズム、CPUソフトが成熟して高速化が見込
めない• 苦労して並列化したが最新理論で楽々抜かれた、、、– 高速で通信が少ないアルゴリズムが確立されている• 計算量が通信量より１～２桁大きい– 類似問題がGPUで加速できると分かっている

量子化学計算はこれらの条件を満たしている
3 Mar 2011

Copyright (C) 2011 X-Ability Co.,Ltd.
All rights reserved.

実アプリの高速化Tips
先の条件を満たした上で、凄く頑張るのは前提として、、、• CUDAだけだと対応部分の速度は出るが精度が足りな
い、メモリ制約があると言われる– 結局、倍精度＆メモリ沢山のCPU演算とのハイブリッド– ホスト側のプログラムの高速化が重要• SSE, AVXのintrinsicを使ってSIMD実装• CUDAで機能をFull実装するのは、ハードウェアの制約もあり新アルゴリズムを生み出す必要がありかなり大変– 結局、これもホスト側とのハイブリッド• アプリケーションによってはhost‐device通信が律速– 結局、これもCPU演算とのハイブリッドで、CPU‐GPUが同一ダイ上にあるプロセッサ（sandy‐bridgeなど）が有効
3 Mar 2011

Copyright (C) 2011 X-Ability Co.,Ltd.
All rights reserved.

Intel Sandy‐bridge• 2011年2月末現在、チップセットがリコール中（デスクトップは動いてますが、、、）のCPU‐GPUハイブリッ
ドプロセッサ
現時点でIntel Composer XEでGPU部分はコンパイルできな
いようである。言語は未定。

• AVXで256bitの要素をdouble x 4 / float x 8にSIMD的に計算できる
SSEの延長• Intel Intrinsics Guide
http://software.intel.com/en‐us/avx/
Intrinsicを使わないとSSEオプションつけてコンパイルしても十分に最適化されない
3 Mar 2011

Copyright (C) 2011 X-Ability Co.,Ltd.
All rights reserved.

AVX Document

3 Mar 2011

Copyright (C) 2011 X-Ability Co.,Ltd.
All rights reserved.

受講ありがとうございました

3 Mar 2011

Copyright (C) 2011 X-Ability Co.,Ltd.
All rights reserved.