1 5 9 3 6 7 2 4
+ + + +
+ +
+
第 1 ステージ
第 2 ステージ
第 3 ステージ
図4.34: リダクション計算のカスケード実行
図4.35: リダクション計算のカスケード実行
float sum = 0.0;
1
for(int j = 0; j < M; j++){
2
for(int i = 0; i < N; i++){
3
sum += z[i][j];
4
}
5
}
6
この計算も繰り返し順に値を足し込むから一見すると並列性がない。が、図4.35のように
7
隣接する4要素毎に加算を行えばカスケード計算が可能である。1次元のカスケード計算
8
が二分木構造であったのに対して2次元のそれは四分木構造になる。
9
図4.36を用いてい、この2次元カスケード計算についてさらに考察する。この図は、図
10
4.35の最初の計算：黒いメッシュの最下層のテクスチャから青いメッシュのひとつ上位の、
11
サイズが1/4のメッシュを求めるステップを抜き出したものである。今、図4.36の青いメッ
12
シュの画素点(m, n)に黒いメッシュの4点のデータの和を求めることを考える。(m, n)に
13
対応する4点の位置は(2m,2n)、(2m+ 1,2n)、(2m,2n+ 1)、(2m+ 1,2n+ 1)、である。
14
フラグメントシェーダーでは(2m,2n)、(2m+ 1,2n)、(2m,2n+ 1)、(2m+ 1,2n+ 1)の
15
位置の数値の和を求め、それを(m, n)の位置に書き込めばよい。ホストプログラムはこの
16
計算を図4.35のように階層的に繰り返せばよい。
17
この方法はデータサイズが大きい場合にはそれなりの高速化が期待できる。
18
m 2m 2m+1 n
2n 2n+1
0 0
図4.36: 2次元リダクション計算のカスケード実行
4.7.1 ホストプログラム
1
さて、またプログラム全体の見通しが悪くなってきた。処理内容も相当に複雑であるか
2
ら、以下ではこの章のまとめを兼ねて、この計算に関する主要なプログラムを全て載せる
3
こととする。
4
まず、宣言などを集めたヘッダーファイルAll.hは図4.37、図4.38の通りである。各種
5
システムヘッダーの読み込み、構造体の定義、頂点バッファを表すクラスArrayBuffer、読
6
み書き可能2次元テクスチャのクラスRWTexture2D、シェーダーを統括するクラスShader
7
を定義している。ArrayBufferクラスのメンバー関数 moveData()を新規に導入してい
8
る。これについては後述する。
9
main関数は図4.39の通りである。図4.2から変更ない。
10
initSystem()は図4.40の通りである。これは図4.3から変更ない。
11
2次元リダクション計算に関係するヘッダファイルの読み込みと大域変数の宣言は図4.41
12
の通りである。頂点バッファはinitData()とcompute()の両方からアクセスするため、
13
大域化した。
14
図4.42が、合計値計算のためのinitData()である。ここでは頂点データの座標値が−1
15
の場合のみ設定している。というもの、図4.35、図4.36に示したように、計算を行う矩形領
16
域の左辺、底辺は固定のままなのだが、右辺は徐々に左側へ移動し、上辺は徐々に下方へ移
17
✓ ✏
#include <iostream>
using namespace std;
#include <stdio.h>
#include <stdlib.h>
#if defined(WIN32)
# pragma comment(lib, "glew32.lib")
# include "glew.h"
# include "glut.h"
# include "glext.h"
#elif defined(__APPLE__) || defined(MACOSX)
# include <GLUT/glut.h>
#else
# define GL_GLEXT_PROTOTYPES
# include <GL/glut.h>
#endif
struct Position2D { float x;
float y;
};
struct RGBA { float r;
float g;
float b;
float a;
};
struct ArrayBuffer { GLuint bufID;
int size;
ArrayBuffer(float* data, int s, int n);
void moveData(float* data, int n); // 頂点データをGPUへ転送 };
struct RWTexture2D { GLuint texID;
GLint num;
RWTexture2D(int tnum, void* data, int w, int h);
void readData(void*, int w, int h);
✒}; ✑
図 4.37: 2 次 元 リ ダ ク ション 計 算 の た め の ヘッダ ファイ ル（ 図 4.38 へ 続 く、
ArrayBuffer::moveData()を新規で追加）
✓ ✏ struct Shader {
GLuint program;
Shader(const char* vsn, const char* fsn);
void use();
void bindArrayBuffer(const char* vname, ArrayBuffer* ap);
void setFloat(const char* vname, float val);
void bindTextureR(const char* vname, RWTexture2D* tp);
void bindTextureW(RWTexture2D* tp);
void run(GLenum mode, int n);
GLuint compileProgram(GLenum type, const GLchar *file);
void buildProgram(const GLchar *vsfile, const GLchar *fsfile);
✒}; ✑
図 4.38: 2次元リダクション計算のためのヘッダファイル（図4.37からの続き、リダク ション計算固有の部分はない）
✓ ✏
int main(int argc, char *argv[]) {
initSystem(argc,argv);
initData();
compute();
showResults();
return 0;
✒} ✑
図4.39: 2次元リダクション計算のためのmain()（リダクション計算固有の部分はない）
✓ ✏ void initSystem(int argc, char *argv[])
{
glutInit(&argc,argv);
glutInitDisplayMode(GLUT_RGB|GLUT_SINGLE);
glutInitWindowSize(width,height);
glutCreateWindow("Test Window");
glClearColor(0.0,0.0,0.0,0.0);
#if defined(WIN32) glewInit();
#endif
GLuint fb;
glGenFramebuffers(1, &fb);
glBindFramebuffer(GL_FRAMEBUFFER, fb);
sp = new Shader("shader.vert","shader.frag");
sp->use();
✒} ✑
図 4.40: 2次元リダクション計算のためのinitSystem()（リダクション計算固有の部分 はない）
✓ ✏
#include "All.h"
#include <math.h>
#include <time.h>
const int NUM_POINTS = 4;
const int width = 4096;
const int height = 4096;
RWTexture2D *texXp; // テクスチャ X へのポインタ
RGBA x[height][width]; // Xのテクスチャデータ
RWTexture2D *texZp; // テクスチャ Z へのポインタ
Shader *sp; // シェーダーへのポインタ
ArrayBuffer* abp; // 頂点バッファへのポインタ
Position2D pos[NUM_POINTS]; // 頂点バッファの配列 double start_clock, end_clock; // 実行時間計測のため
✒ ✑
✓ ✏ void initData()
{
pos[0].x = -1.0; pos[0].y = -1.0; // 頂点データの部分的な初期化 pos[1].y = -1.0; // 残りはcompute()内で計算 pos[3].x = -1.0;
//頂点バッファの生成とシェーダーへの結合
abp = new ArrayBuffer((float*)pos,2,NUM_POINTS);
sp->bindArrayBuffer("position",abp);
for(int h = height-1; h >= 0; h--){ // テクスチャ X の初期値の設定 for(int w = 0; w < width; w++){
x[h][w].r = w+h+0.1;
x[h][w].g = w+h+0.2;
x[h][w].b = w+h+0.3;
x[h][w].a = w+h+0.4;
} }
start_clock = clock();
texXp = new RWTexture2D(0,x,width,height); // Xを入力に設定 texZp = new RWTexture2D(1,NULL,width,height); // Zを出力に設定 sp->setFloat("width",width); // uniform変数 width の設定 sp->setFloat("height",height); // uniform変数 height の設定
✒} ✑
図4.42: 2次元リダクション計算のためのinitData()
✓ ✏ void compute(void)
{
float offset;
for(int size = width; size > 1; ){
offset = float(size)/float(width)-1.0;
pos[1].x = offset;
pos[2].x = offset; pos[2].y = offset;
pos[3].y = offset;
abp->moveData((float*)pos,NUM_POINTS);
sp->bindTextureR("tx",texXp);
sp->bindTextureW(texZp);
glClear(GL_COLOR_BUFFER_BIT);
sp->run(GL_POLYGON,NUM_POINTS);
size /= 2;
offset = float(size)/float(width)-1.0;
pos[1].x = offset;
pos[2].x = offset;
pos[2].y = offset;
pos[3].y = offset;
abp->moveData((float*)pos,NUM_POINTS);
sp->bindTextureR("tx",texZp);
sp->bindTextureW(texXp);
glClear(GL_COLOR_BUFFER_BIT);
sp->run(GL_POLYGON,NUM_POINTS);
size /= 2;
}
✒} ✑
図4.43: 2次元リダクション計算のためのcompute()
動しながらリダクション演算を行うため、初期設定では左辺、底辺に関する座標値のみを設
1
定している。右辺、上辺は関数compute()内で設定する。また、テクスチャとシェーダー
2
プログラムの関連付けを行う関数呼び出しsp->bindTextureR()、sp->bindTextureW()
3
も全て compute() へ移動させた。
4
リダクション計算の中核は、図4.43の関数compute()である。
5
図4.36を思い出そう。まず、画像全体を表す矩形領域は、
左下の座標点：(−1, −1), 右上の座標点：(1, 1)
である。それに対して、青いメッシュの矩形領域は、縦横の辺の長さが半分になるから
✓ ✏ void ArrayBuffer::moveData(float* data, int n){
glBindBuffer(GL_ARRAY_BUFFER, bufID);
glBufferSubData(GL_ARRAY_BUFFER, 0, sizeof(float)*size*n, data);
✒} ✑
図 4.44: 頂点データをGPUへ転送する関数 である。さらに緑のメッシュの矩形領域は、
左下の座標点：(−1, −1), 右上の座標点：(0.5, 0.5)
である。これを繰り返し、その矩形領域に含まれる画素点が1個になったところで計算を
1
終了させる。関数 compute()のforループは、矩形のサイズを1／2に減じながら実行
2
する繰り返しである。
3
関数 compute()の中を少し詳しく述べると、図4.43のプログラム辺：
4
offset = float(size)/float(width)-1.0;
5
pos[1].x = offset;
6
pos[2].x = offset;
7
pos[2].y = offset;
8
pos[3].y = offset;
9
はその矩形領域の右辺と上辺の座標値の設定を行なっている部分である。
10
設定された頂点座標値は
11
abp->moveData((float*)pos,NUM_POINTS);
12
によってGPUへ転送する。関数moveData()の実装は図4.44の通りである。ここにOpenGL
13
の関数呼び出し：
14
glBufferSubData(G ARRAY BUFFER, 0, sizeof(float)*size*n, data)17
15
は、すでにGPU内に領域確保されている頂点バッファに頂点データを一斉転送する18。
16
以上が矩形領域の設定である。
17
一旦、矩形が設定されたならば、関数呼び出し：
18
sp->bindTextureR("tx",texXp);
19
sp->bindTextureW(texZp);
20
glClear(GL_COLOR_BUFFER_BIT);
21
sp->run(GL_POLYGON,NUM_POINTS);
22
17void glBufferSubData(GLenum target, GLintptr offset, GLsizeiptr size, const GLvoid * data);
18一斉転送ではなく、glMapBuffer()19、glUnmapBuffer()20を用いることも可能であるが、それについ ての解説は省略する。
✓ ✏ void showResults()
{
RGBA result;
texXp->readData(&result,1,1);
end_clock = clock();
double total = (end_clock-start_clock)/CLOCKS_PER_SEC;
double gflops = (4.0*(4096*4096-1)) / (total * 1e9);
start_clock = clock();
RGBA result2;
result2.r = result2.g = result2.b = result2.a = 0.0;
for(int h = height-1; h >= 0; h--){
for(int w = 0; w < width; w++){
result2.r += x[h][w].r;
result2.g += x[h][w].g;
result2.b += x[h][w].b;
result2.a += x[h][w].a;
} }
end_clock = clock();
double total2 = (end_clock-start_clock)/CLOCKS_PER_SEC;
double gflops2 = (4.0*(4096*4096-1)) / (total2 * 1e9);
printf("GPU: %f %f %f %f\n",result.r,result.g,result.b,result.a);
printf("CPU: %f %f %f %f\n",result2.r,result2.g,result2.b,result2.a);
printf("%le %le %le %le\n",total,gflops,total2,gflops2);
✒} ✑
図 4.45: リダクション計算のためのshowResults()
は、テクスチャの設定、フレームバッファの初期化、描画の実行を行う。この一連の実行
1
で図4.35の黒いメッシュの4画素点の値の和は、対応する青いメッシュの画素点に格納さ
2
れる。
3
OpenGLによるGPGPUの基本はピンポン計算であるから、図4.43の後半では、2枚
4
のテクスチャの役割を入れ替え、矩形の一片のサイズを半分に減じて、同様の計算を繰り
5
返す。
6
このピンポン計算を矩形に含まれる画素点が1個になるまで繰り返す。
7
図4.45は、計算結果、実行時間を出力するプログラムである。
8
GPUの計算結果はGPUのメモリ内のテクスチャの左下の1画素分の求められている
9
から、
10
texXp->readData(&result,1,1);
1
によってCPUの1画素分の変数 resultに転送している。
2
CPUによる合計計算は、単純な二重ループで実装した。
3
出力は、まずGPUとCPUの計算結果を
4
printf("GPU: %f %f %f %f\n",result.r,result.g,result.b,result.a);
5
printf("CPU: %f %f %f %f\n",result2.r,result2.g,result2.b,result2.a);
6
によって比較出力する。それに続きて実行時間、演算性能を計算し、出力する。この計算
7
に必要な総演算数は、4.0*(4096*4096-1)である21。GPU内での実際の演算数も同じで
8
ある。
9
4.7.2 シェーダープログラム
10
バーテックスシェーダーのソースプログラムは図4.46の通りである。これは図4.11と
11
同じものである。
12
フラグメントシェーダーのソースプログラムは図4.47の通りである。図4.36で見たよ
13
うに、フラグメントシェーダーの画素点が(m, n)であるとき、テクスチャ上の(2m,2n)、
14
(2m+ 1,2n)、(2m,2n+ 1)、(2m+ 1,2n+ 1)の位置の値の和を求めねばならない。そこで、
15
座標値を保持する変数texCoord、texCoord00、texCoord10、texCoord01、texCoord00
16
を用いて、テクスチャ上の座標を計算している。プログラム中の定数 0.5 は座標値をよ
17
り正確に計算するための調整用である。
18
4.7.3 実行結果
19
プログラムを実行するとshowResults()が以下のような計算結果を出力する。
20
GPU: 68704378880.000 68706054144.000 68707733504.000 68709408768.000
21
CPU: 68711055360.000 68711055360.000 68711055360.000 68711055360.000
22
GPUとCPUの計算結果が異なる。これはGPUのリダクション計算とCPUの単純な二重
23
ループの計算の計算順序の違いに依る。float型は有効桁数がせいぜい7桁程度しかない
24
ため、7桁以上大きさの異なる数の加算では小さい方の数が無視される。単純な二重ルー
25
プによる積算ではそれが起きやすいが、階層的なリダクション計算ではそれが起きにくい。
26
実際、テクスチャサイズを小さくすると、GPUとCPUの計算結果は完全に一致するこ
27
とを確認している。
28
21正確に言えば、CPUによる単純な二重ループでの演算数は4.0*4096*4096であって、4個だけ無駄な 演算があるが、もちろん演算性能の数値ではほとんど無視できる。
✓ ✏
#version 120
attribute vec2 position;
void main(void) {
gl_Position = vec4(position,0.0,1.0);
✒} ✑
図4.46: リダクション計算のためのバーテックシェーダーソースプログラム
✓ ✏
#version 120
uniform sampler2D tx;
uniform float width;
uniform float height;
void main(void) {
vec2 texCoord = 2*(gl_FragCoord.xy-vec2(0.5));
vec2 delta = vec2(0.5/width,0.5/height);
vec2 texCoord00
= vec2(texCoord.x/width,texCoord.y/height)+delta;
vec2 texCoord10
= vec2((texCoord.x+1)/width,texCoord.y/height)+delta;
vec2 texCoord01
= vec2(texCoord.x/width,(texCoord.y+1)/height)+delta;
vec2 texCoord11
= vec2((texCoord.x+1)/width,(texCoord.y+1)/height)+delta;
vec4 x00 = texture2D(tx,texCoord00);
vec4 x10 = texture2D(tx,texCoord10);
vec4 x01 = texture2D(tx,texCoord01);
vec4 x11 = texture2D(tx,texCoord11);
gl_FragColor = x00+x10+x01+x11;
✒} ✑
図 4.47: リダクション計算のためのフラグメントシェーダーソースプログラム
表4.4: リダクション計算の実行時間と実行速度
GPU CPU
時間 速度 時間 速度
1.69E-01 4.00E-01 2.13E-02 3.15E+00 時間の単位は秒、速度の単位はGFLOPS
4.7.4 実行速度
1
図4.4は、GPUによる計算とCPUによる計算の実行時間、演算性能である。テクスチャ
2
サイズは、図4.41のプログラム中に埋め込んでいるように、40962とした。前節、前々節
3
の計算では繰り返しを多数回行うことで演算性能を稼ぐことができたが、この節のリダク
4
ション演算では合計値を1回計算しただけの時間を測定した。
5
4.7.5 演算性能の評価
6
図4.4によれば、GPUの演算性能は約0.4GFLOPSであるが、CPUのそれは3.15GFLOPS
7
であった。
8
総演算数が少ない上に単純な計算ではないことから、GPUの性能がCPUの約1/10で
9
あることは仕方がない。むしろ、1/10で済んだことを良しとすべきである。特殊な目的
10
でない限り、リダクション計算はそれほど頻度の高い計算ではなく、多用されることはな
11
い。リダクション計算が全体の計算の中でそれほどネックにならないことを確認できたと
12
考えるべきである。
13
この章では、OpenGLを用いたGPGPUを見てきた。
14
現在はCUDAやOpenCLといったGPGPU専用の環境が利用できるから、この章の
15
内容がGPGPUの主流になることはなく、むしろ古い技術とみなすべきである。しかし、
16
「なぜグラフィックス専用ハードウェアが高速計算に利用できるのか」という問いへのひと
17
つの答えになっていると考える。
18
最近の3D CGレンダリングは大規模計算を必要としている。ここで述べた内容はその
19
用途では依然として重要な技術でありうる。
20