Main Content

GPU における A\b のベンチマーク

この例では、GPU における線形計算の求解のベンチマークを実行する方法について考えます。A*x = bx を求めるための MATLAB® コードは非常にシンプルです。最も一般的な方法は、行列の左除算 (mldivide またはバックスラッシュ演算子 (\) とも呼ばれる) を使用して x (つまり、x = A\b) を計算するものです。

関連する例:

この例のコードは以下の関数に含まれています。

function results = paralleldemo_gpu_backslash(maxMemory)

計算には適切な行列サイズを選択することが非常に重要です。そのために、CPU および GPU で使用可能なシステム メモリの量 (単位: GB) を指定します。既定値は GPU で使用可能なメモリ量のみに基づくものであり、ここでは使用するシステムに適した値を指定することができます。

if nargin == 0
    g = gpuDevice;
    maxMemory = 0.4*g.AvailableMemory/1024^3;
end

ベンチマーク関数

ここで実行するのは、行列の左除算 (\) のベンチマークであり、CPU と GPU 間のデータ転送コストや行列作成時間などのパラメーターのベンチマークではありません。そのため、データ生成を線形計算の求解と切り離してから、求解のみにかかる時間を測定します。

function [A, b] = getData(n, clz)
    fprintf('Creating a matrix of size %d-by-%d.\n', n, n);
    A = rand(n, n, clz) + 100*eye(n, n, clz);
    b = rand(n, 1, clz);
end

function time = timeSolve(A, b, waitFcn)
    tic;
    x = A\b; %#ok<NASGU> We don't need the value of x.
    waitFcn(); % Wait for operation to complete.
    time = toc;
end

問題の規模の選択

他の多数の並列アルゴリズムと同様に、並列での線形計算の求解のパフォーマンスは行列のサイズによって大きく異なります。A\b のベンチマークなどの他の例と同様に、さまざまな行列サイズについてアルゴリズムのパフォーマンスを比較します。

% Declare the matrix sizes to be a multiple of 1024.
maxSizeSingle = floor(sqrt(maxMemory*1024^3/4));
maxSizeDouble = floor(sqrt(maxMemory*1024^3/8));
step = 1024;
if maxSizeDouble/step >= 10
    step = step*floor(maxSizeDouble/(5*step));
end
sizeSingle = 1024:step:maxSizeSingle;
sizeDouble = 1024:step:maxSizeDouble;

パフォーマンスの比較: ギガフロップス

パフォーマンスを測定するための値として 1 秒当たりの浮動小数点演算の回数 (フロップス) を使用します。そのようにすると、さまざまな行列サイズについてアルゴリズムのパフォーマンスを比較できるようになります。

行列サイズを指定すると、ベンチマーク関数によって行列 A と右辺 b が 1 回作成されてから、A\b の求解が複数回実行されて、所要時間の正確な測定値が取得されます。ここでは、HPC チャレンジの浮動小数点演算の回数を使用します。したがって、n 行 n 列の行列の場合、浮動小数点演算の回数は 2/3*n^3 + 3/2*n^2 になります。

この関数はハンドルで関数 "wait" に渡されます。CPU では、この関数を実行しても何も起こりません。GPU では、この関数は保留中のすべての演算が完了するまで待機します。このように待機することで、正確な時間測定値が得られます。

function gflops = benchFcn(A, b, waitFcn)
    numReps = 3;
    time = inf;
    % We solve the linear system a few times and calculate the Gigaflops
    % based on the best time.
    for itr = 1:numReps
        tcurr = timeSolve(A, b, waitFcn);
        time = min(tcurr, time);
    end

    % Measure the overhead introduced by calling the wait function.
    tover = inf;
    for itr = 1:numReps
        tic;
        waitFcn();
        tcurr = toc;
        tover = min(tcurr, tover);
    end
    % Remove the overhead from the measured time. Don't allow the time to
    % become negative.
    time = max(time - tover, 0);
    n = size(A, 1);
    flop = 2/3*n^3 + 3/2*n^2;
    gflops = flop/time/1e9;
end

% The CPU doesn't need to wait: this function handle is a placeholder.
function waitForCpu()
end

% On the GPU, to ensure accurate timing, we need to wait for the device
% to finish all pending operations.
function waitForGpu(theDevice)
    wait(theDevice);
end

ベンチマークの実行

これでセットアップがすべて完了したので、ベンチマークを容易に実行できるようになりました。ただし、計算が完了するまでに時間がかかることがあるため、それぞれの行列サイズのベンチマークが完了するたびに中間ステータス情報を出力するものとします。また、単精度計算と倍精度計算の両方のベンチマークを実行するため、すべての行列サイズに対するループ処理を関数でカプセル化します。

function [gflopsCPU, gflopsGPU] = executeBenchmarks(clz, sizes)
    fprintf(['Starting benchmarks with %d different %s-precision ' ...
         'matrices of sizes\nranging from %d-by-%d to %d-by-%d.\n'], ...
            length(sizes), clz, sizes(1), sizes(1), sizes(end), ...
            sizes(end));
    gflopsGPU = zeros(size(sizes));
    gflopsCPU = zeros(size(sizes));
    gd = gpuDevice;
    for i = 1:length(sizes)
        n = sizes(i);
        [A, b] = getData(n, clz);
        gflopsCPU(i) = benchFcn(A, b, @waitForCpu);
        fprintf('Gigaflops on CPU: %f\n', gflopsCPU(i));
        A = gpuArray(A);
        b = gpuArray(b);
        gflopsGPU(i) = benchFcn(A, b, @() waitForGpu(gd));
        fprintf('Gigaflops on GPU: %f\n', gflopsGPU(i));
    end
end

次に、単精度および倍精度でのベンチマークを実行します。

[cpu, gpu] = executeBenchmarks('single', sizeSingle);
results.sizeSingle = sizeSingle;
results.gflopsSingleCPU = cpu;
results.gflopsSingleGPU = gpu;
[cpu, gpu] = executeBenchmarks('double', sizeDouble);
results.sizeDouble = sizeDouble;
results.gflopsDoubleCPU = cpu;
results.gflopsDoubleGPU = gpu;
Starting benchmarks with 7 different single-precision matrices of sizes
ranging from 1024-by-1024 to 19456-by-19456.
Creating a matrix of size 1024-by-1024.
Gigaflops on CPU: 43.805496
Gigaflops on GPU: 78.474002
Creating a matrix of size 4096-by-4096.
Gigaflops on CPU: 96.459635
Gigaflops on GPU: 573.278854
Creating a matrix of size 7168-by-7168.
Gigaflops on CPU: 184.997657
Gigaflops on GPU: 862.755636
Creating a matrix of size 10240-by-10240.
Gigaflops on CPU: 204.404384
Gigaflops on GPU: 978.362901
Creating a matrix of size 13312-by-13312.
Gigaflops on CPU: 218.773070
Gigaflops on GPU: 1107.983667
Creating a matrix of size 16384-by-16384.
Gigaflops on CPU: 233.529176
Gigaflops on GPU: 1186.423754
Creating a matrix of size 19456-by-19456.
Gigaflops on CPU: 241.482550
Gigaflops on GPU: 1199.151846
Starting benchmarks with 5 different double-precision matrices of sizes
ranging from 1024-by-1024 to 13312-by-13312.
Creating a matrix of size 1024-by-1024.
Gigaflops on CPU: 34.902918
Gigaflops on GPU: 72.191488
Creating a matrix of size 4096-by-4096.
Gigaflops on CPU: 74.458136
Gigaflops on GPU: 365.339897
Creating a matrix of size 7168-by-7168.
Gigaflops on CPU: 93.313782
Gigaflops on GPU: 522.514165
Creating a matrix of size 10240-by-10240.
Gigaflops on CPU: 104.219804
Gigaflops on GPU: 628.301313
Creating a matrix of size 13312-by-13312.
Gigaflops on CPU: 108.826886
Gigaflops on GPU: 681.881032

パフォーマンスのプロット

ここで、結果をプロットし、単精度と倍精度の両方について CPU でのパフォーマンスと GPU でのパフォーマンスを比較します。

最初に、単精度でのバックスラッシュ演算子のパフォーマンスを確認します。

fig = figure;
ax = axes('parent', fig);
plot(ax, results.sizeSingle, results.gflopsSingleGPU, '-x', ...
     results.sizeSingle, results.gflopsSingleCPU, '-o')
grid on;
legend('GPU', 'CPU', 'Location', 'NorthWest');
title(ax, 'Single-precision performance')
ylabel(ax, 'Gigaflops');
xlabel(ax, 'Matrix size');
drawnow;

次に、倍精度でのバックスラッシュ演算子のパフォーマンスを確認します。

fig = figure;
ax = axes('parent', fig);
plot(ax, results.sizeDouble, results.gflopsDoubleGPU, '-x', ...
     results.sizeDouble, results.gflopsDoubleCPU, '-o')
legend('GPU', 'CPU', 'Location', 'NorthWest');
grid on;
title(ax, 'Double-precision performance')
ylabel(ax, 'Gigaflops');
xlabel(ax, 'Matrix size');
drawnow;

最後に、GPU と CPU を比較した場合のバックスラッシュ演算子の高速化を確認します。

speedupDouble = results.gflopsDoubleGPU./results.gflopsDoubleCPU;
speedupSingle = results.gflopsSingleGPU./results.gflopsSingleCPU;
fig = figure;
ax = axes('parent', fig);
plot(ax, results.sizeSingle, speedupSingle, '-v', ...
     results.sizeDouble, speedupDouble, '-*')
grid on;
legend('Single-precision', 'Double-precision', 'Location', 'SouthEast');
title(ax, 'Speedup of computations on GPU compared to CPU');
ylabel(ax, 'Speedup');
xlabel(ax, 'Matrix size');
drawnow;

end
ans = 

         sizeSingle: [1024 4096 7168 10240 13312 16384 19456]
    gflopsSingleCPU: [1x7 double]
    gflopsSingleGPU: [1x7 double]
         sizeDouble: [1024 4096 7168 10240 13312]
    gflopsDoubleCPU: [34.9029 74.4581 93.3138 104.2198 108.8269]
    gflopsDoubleGPU: [72.1915 365.3399 522.5142 628.3013 681.8810]