Exponenta Banner

Сравнение A\b на графическом процессоре

Открыть скрипт

В этом примере рассматривается, как мы можем сравнить решение линейной системы на графическом процессоре. Код MATLAB ®, для решения x в A*x = b очень просто. Чаще всего мы используем матричное левое деление, также известное как mldivide или оператор обратной косой черты (\), чтобы вычислить x (то есть x = A\b).

Похожие примеры:

Код, показанный в этом примере, можно найти в этой функции:

function results = paralleldemo_gpu_backslash(maxMemory)

Важно выбрать подходящий размер матрицы для расчетов. Мы можем сделать это, указав объем системной памяти в ГБ, доступный для CPU и GPU. Значение по умолчанию основано только на объеме памяти, доступном на графическом процессоре, и можно задать значение, подходящее для вашей системы.

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

Функция бенчмаркинга

Мы хотим бенчмарковать матричное левое деление (\), а не стоимость передачи данных между центральным процессором и графическим процессором, время создания матрицы или другие параметры. Поэтому мы отделяем генерацию данных от решения линейной системы и измеряем только время, необходимое для выполнения последней.

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

Выбор размера задачи

Как и в случае с большим количеством других параллельных алгоритмов, эффективность решения линейной системы параллельно сильно зависит от размера матрицы. Как видно из других примеров, таких как Benchmarking 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;

Сравнение эффективности: Gigaflops

Мы используем количество операций с плавающей точкой в секунду в качестве меры эффективности, потому что это позволяет сравнивать эффективность алгоритма для различных размеров матрицы.

Учитывая размер матрицы, функция бенчмаркинга создает матрицу A и правую сторону b один раз, а затем решает A\b несколько раз, чтобы получить точную меру времени, которое требуется. Мы используем счетчик операций с плавающей точкой HPC Challenge, так что для матрицы n на n мы считаем операции с плавающей точкой как 2/3*n^3 + 3/2*n^2.

Функция передается в указателе в функцию 'wait'. На центральном процессоре эта функция ничего не делает. На графическом процессоре эта функция ожидает завершения всех ожидающих операций. Ожидание таким образом обеспечивает точность синхронизации.

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;

Наконец, рассмотрим ускорение работы оператора backslash при сравнении графический процессор с центральный процессор.

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]

Документация по Parallel Computing Toolbox

Поддержка