В этом примере используется «Игра жизни» Конвея для демонстрации того, как можно выполнять трафаретные операции с помощью графического процессора.
Многие операции массива могут быть выражены как «операция трафарета», где каждый элемент выходного массива зависит от небольшой области входного массива. Примеры включают конечные различия, свертку, медианную фильтрацию и методы конечных элементов. Этот пример использует «Игру жизни» Конвея для демонстрации двух способов выполнения трафаретной операции на GPU, начиная с кода в электронной книге Клеве Молера «Эксперименты в MATLAB».
«Игра жизни» следует нескольким простым правилам:
Ячейки расположены в 2D сетке
На каждом шаге судьба каждой клетки определяется жизненной силой её восьми ближайших соседей
Любая клетка с ровно тремя живыми соседями оживает на следующем шаге
Живая клетка с ровно двумя живыми соседями остается живой на следующем шаге
Все другие клетки (включая клетки с более чем тремя соседями) умирают на следующем шаге или остаются пустыми
Таким образом, «трафаретом» в данном случае является область 3x3 вокруг каждого элемента. Вот несколько примеров обновления ячейки:

Этот пример является функцией, разрешающей использование вложенных функций:
function paralleldemo_gpu_stencil()
Первоначальная популяция клеток создается на 2D сетке, где примерно 25% мест живы.
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');

Электронная книга Эксперименты в MATLAB предоставляет начальную реализацию, которая может быть использована для сравнения. Эта версия полностью векторизирована, обновляя все ячейки в сетке за один проход за поколение.
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. Алгоритм остается неизменным. Обратите внимание, что 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).
Просмотр расчетов в updateGrid очевидно, что одни и те же операции применяются в каждом местоположении сетки независимо. Это говорит о том, что arrayfun может использоваться для проведения оценки. Однако каждая ячейка должна знать о своих восьми соседях, нарушая элементарную независимость. Каждый элемент должен иметь доступ к полной сетке, а также работать независимо.
Решение состоит в использовании вложенной функции. Вложенные функции, даже используемые с 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 здесь: расширение размера. Нам нужно было пройти только векторы строк и столбцов, и они были автоматически расширены в полную сетку. Эффект такой, как мы назвали:
[cols,rows] = meshgrid(cols,rows);
в рамках arrayfun звоните. Это экономит как некоторые вычисления, так и некоторую передачу данных между CPU-памятью и GPU-памятью.
В этом примере простая трафаретная операция, «Игра жизни» Конвея, была реализована на GPU с использованием arrayfun и переменные, объявленные в родительской функции. Этот метод может быть использован для реализации ряда трафаретных операций, включая алгоритмы конечных элементов, свертки и фильтры. Он также может использоваться для доступа к элементам в справочной таблице, определенной в родительской функции.
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