Main Content

GPU でのステンシル演算

スクリプトを開く

この例では Conway の "ライフ ゲーム" を使用し、GPU を使ってステンシル演算をどのように実行するのかを示します。

多くの配列演算は "ステンシル演算" として表現でき、そこでは出力配列の各要素が入力配列の小さな領域によって決まります。例としては、有限差分、畳み込み、メディアン フィルター処理、有限要素法などが挙げられます。この例では Conway の "ライフ ゲーム" を使用して、GPU でステンシル演算を実行する 2 つの方法を示します。まずは、Cleve Moler の電子書籍 "Experiments with MATLAB" にあるコードから開始します。

"ライフ ゲーム" ではいくつかの簡単なルールに従います。

  • セルは 2 次元グリッドに並べられる

  • 各ステップで、それぞれのセルの運命は隣接する 8 つのセルの生命力によって決定される

  • あるセルが、生命をもつ 3 つのセルと隣接していると次のステップで生命を得る

  • 生命をもつセルが、生命をもつ 2 つのセルと隣接していると、次のステップで生命が維持される

  • 他のセルはすべて (生命をもつ隣接セルが 3 つを超える場合を含め) 次のステップで生命を失うか、空のままとなる

したがって、この場合の "ステンシル" は、各要素を囲む 3 行 3 列の領域となります。以下に、セルがどのように更新されるかについていくつかの例を挙げます。

この例は、入れ子関数の使用が可能な関数です。

function paralleldemo_gpu_stencil()

ランダムな初期集団の生成

セルの初期集団は、およそ 25% の位置が生命をもつように 2 次元グリッド上に作成されます。

gridSize = 500;
numGenerations = 100;
initialGrid = (rand(gridSize,gridSize) > .75);
gpu = gpuDevice();

% Draw the initial grid
hold off
imagesc(initialGrid);
colormap([1 1 1;0 0.5 0]);
title('Initial Grid');

ライフ ゲームの実行

電子書籍 "Experiments with MATLAB" には、比較に使用できる初期実装が提示されています。このバージョンは完全にベクトル化されており、グリッドのすべてのセルを世代ごとに 1 回の受け渡しで更新します。

    function X = updateGrid(X, N)
        p = [1 1:N-1];
        q = [2:N N];
        % Count how many of the eight neighbors are alive.
        neighbors = X(:,p) + X(:,q) + X(p,:) + X(q,:) + ...
            X(p,p) + X(q,q) + X(p,q) + X(q,p);
        % A live cell with two live neighbors, or any cell with
        % three live neighbors, is alive at the next step.
        X = (X & (neighbors == 2)) | (neighbors == 3);
    end

grid = initialGrid;
% Loop through each generation updating the grid and displaying it
for generation = 1:numGenerations
    grid = updateGrid(grid, gridSize);

    imagesc(grid);
    title(num2str(generation));
    drawnow;
end

次に、ゲームを再度実行し、各世代にどのくらい時間を要するかを測定します。

grid = initialGrid;
timer = tic();

for generation = 1:numGenerations
    grid = updateGrid(grid, gridSize);
end

cpuTime = toc(timer);
fprintf('Average time on the CPU: %2.3fms per generation.\n', ...
    1000*cpuTime/numGenerations);

Average time on the CPU: 11.323ms per generation.

下記の各バージョンの妥当性を確認するために、この結果を保持します。

expectedResult = grid;

ライフ ゲームを GPU で実行するための変換

ライフ ゲームを GPU で実行するために、gpuArrayを使用して初期集団を GPU に送信します。アルゴリズムは変更されません。確実に GPU での計算が完了してからタイマーを停止する目的で、wait (GPUDevice)が使用されることに注意してください。これは、時間を正確に測定する際にのみ必要になります。

grid = gpuArray(initialGrid);
timer = tic();

for generation = 1:numGenerations
    grid = updateGrid(grid, gridSize);
end

wait(gpu); % Only needed to ensure accurate timing
gpuSimpleTime = toc(timer);

% Print out the average computation time and check the result is unchanged.
fprintf(['Average time on the GPU: %2.3fms per generation ', ...
    '(%1.1fx faster).\n'], ...
    1000*gpuSimpleTime/numGenerations, cpuTime/gpuSimpleTime);
assert(isequal(grid, expectedResult));

Average time on the GPU: 1.655ms per generation (6.8x faster).

GPU 用の要素単位バージョンの作成

関数 updateGrid の計算を見ると、明らかに各グリッド位置で同じ演算が独立に適用されています。このことから、評価の実行にarrayfunを使用できる可能性が示されます。ただし、各セルには隣接する 8 つのセルの情報が必要であるため、要素単位の独立性はありません。各要素は独立して機能しつつ、同時にフル グリッドにアクセスできることが必要です。

解決策は、入れ子関数を使用することです。入れ子関数は、arrayfunで使用される場合でも、親関数で宣言された変数にアクセスできます。つまり、各セルは先行するタイムステップのグリッド全体を読み取って、これにインデックスを付けることができます。

grid = gpuArray(initialGrid);

    function X = updateParentGrid(row, col, N)
        % Take account of boundary effects
        rowU = max(1,row-1);  rowD = min(N,row+1);
        colL = max(1,col-1);  colR = min(N,col+1);
        % Count neighbors
        neighbors ...
            = grid(rowU,colL) + grid(row,colL) + grid(rowD,colL) ...
            + grid(rowU,col)                   + grid(rowD,col) ...
            + grid(rowU,colR) + grid(row,colR) + grid(rowD,colR);
        % A live cell with two live neighbors, or any cell with
        % three live neighbors, is alive at the next step.
        X = (grid(row,col) & (neighbors == 2)) | (neighbors == 3);
    end


timer = tic();

rows = gpuArray.colon(1, gridSize)';
cols = gpuArray.colon(1, gridSize);
for generation = 1:numGenerations
    grid = arrayfun(@updateParentGrid, rows, cols, gridSize);
end

wait(gpu); % Only needed to ensure accurate timing
gpuArrayfunTime = toc(timer);

% Print out the average computation time and check the result is unchanged.
fprintf(['Average time using GPU arrayfun: %2.3fms per generation ', ...
    '(%1.1fx faster).\n'], ...
    1000*gpuArrayfunTime/numGenerations, cpuTime/gpuArrayfunTime);
assert(isequal(grid, expectedResult));

Average time using GPU arrayfun: 0.795ms per generation (14.2x faster).

arrayfunのもう 1 つの新しい機能もここで使用されていることに注意してください。すなわち、次元拡張です。行ベクトルと列ベクトルを渡すだけで、それが自動的にフル グリッドに拡張されます。その効果は、

[cols,rows] = meshgrid(cols,rows);

をarrayfun呼び出しの一部として呼び出した場合に似ています。これによって計算と、CPU メモリおよび GPU メモリ間のデータ転送が両方とも一部軽減されます。

まとめ

この例では、簡単なステンシル演算である Conway の "ライフ ゲーム" を、arrayfunおよび親関数で宣言された変数を使用して GPU 上に実装しました。この手法は、有限要素アルゴリズム、畳み込み、フィルターなど、さまざまなステンシル演算の実行に使用できます。また、親関数で定義されたルックアップ テーブル内の要素へのアクセスにも使用できます。

fprintf('CPU:          %2.3fms per generation.\n', ...
    1000*cpuTime/numGenerations);
fprintf('Simple GPU:   %2.3fms per generation (%1.1fx faster).\n', ...
    1000*gpuSimpleTime/numGenerations, cpuTime/gpuSimpleTime);
fprintf('Arrayfun GPU: %2.3fms per generation (%1.1fx faster).\n', ...
    1000*gpuArrayfunTime/numGenerations, cpuTime/gpuArrayfunTime);

CPU:          11.323ms per generation.
Simple GPU:   1.655ms per generation (6.8x faster).
Arrayfun GPU: 0.795ms per generation (14.2x faster).

end

参考

| |

関連する例

詳細