В этом примере показано, как простая, известная математическая проблема, Множество Мандельброта, может быть выражена в коде MATLAB®. Используя Parallel Computing Toolbox™ этот код затем адаптируется, чтобы использовать оборудование графического процессора тремя способами:
Используя существующий алгоритм, но с данными графического процессора, как введено
Используя arrayfun
выполнять алгоритм на каждом элементе независимо
Используя интерфейс MATLAB/CUDA, чтобы запустить некоторый существующий CUDA/C ++ код
Значения ниже задают высоко масштабируемую часть Множества Мандельброта в овраге между основной кардиоидой и p/q лампой с ее левой стороны от него.
1000x1000 сетка действительных частей (X) и мнимых частей (Y) создается между этими пределами, и алгоритм Мандельброта выполнен с помощью итераций в каждом местоположении сетки. Для этого конкретного местоположения 500 итераций будет достаточно, чтобы полностью представить изображение.
maxIterations = 500; gridSize = 1000; xlim = [-0.748766713922161, -0.748766707771757]; ylim = [ 0.123640844894862, 0.123640851045266];
Ниже реализация Множества Мандельброта с помощью стандартных команд MATLAB, работающих на центральном процессоре. Это основано на коде, предоставленном в "Экспериментах Клева Молера с MATLAB" электронная книга.
Это вычисление векторизовано таким образом, что каждое местоположение обновляется целиком.
% Setup t = tic(); x = linspace( xlim(1), xlim(2), gridSize ); y = linspace( ylim(1), ylim(2), gridSize ); [xGrid,yGrid] = meshgrid( x, y ); z0 = xGrid + 1i*yGrid; count = ones( size(z0) ); % Calculate z = z0; for n = 0:maxIterations z = z.*z + z0; inside = abs( z )<=2; count = count + inside; end count = log( count ); % Show cpuTime = toc( t ); fig = gcf; fig.Position = [200 200 600 600]; imagesc( x, y, count ); colormap( [jet();flipud( jet() );0 0 0] ); axis off title( sprintf( '%1.2fsecs (without GPU)', cpuTime ) );
gpuArray
Когда данные об обнаружении MATLAB по графическому процессору, вычисления с теми данными выполняются на графическом процессоре. Класс gpuArray
обеспечивает версии графического процессора многих функций, которые можно использовать, чтобы создать массивы данных, включая linspace
logspace
, и meshgrid
функции необходимы здесь. Точно так же count
массив инициализируется непосредственно на графическом процессоре с помощью функционального ones
.
С этими изменениями в инициализации данных вычисления будут теперь выполняться на графическом процессоре:
% Setup t = tic(); x = gpuArray.linspace( xlim(1), xlim(2), gridSize ); y = gpuArray.linspace( ylim(1), ylim(2), gridSize ); [xGrid,yGrid] = meshgrid( x, y ); z0 = complex( xGrid, yGrid ); count = ones( size(z0), 'gpuArray' ); % Calculate z = z0; for n = 0:maxIterations z = z.*z + z0; inside = abs( z )<=2; count = count + inside; end count = log( count ); % Show count = gather( count ); % Fetch the data back from the GPU naiveGPUTime = toc( t ); imagesc( x, y, count ) axis off title( sprintf( '%1.3fsecs (naive GPU) = %1.1fx faster', ... naiveGPUTime, cpuTime/naiveGPUTime ) )
Отмечая, что алгоритм работает одинаково с каждым элементом входа, мы можем поместить код в функцию помощника и вызвать его с помощью arrayfun
. Для входных параметров графического процессора массивов функция используется с arrayfun
скомпилирован в нативный код графического процессора. В этом случае мы поместили цикл в pctdemo_processMandelbrotElement.m
:
function count = pctdemo_processMandelbrotElement(x0,y0,maxIterations) z0 = complex(x0,y0); z = z0; count = 1; while (count <= maxIterations) && (abs(z) <= 2) count = count + 1; z = z*z + z0; end count = log(count);
Обратите внимание на то, что раннее аварийное прекращение работы было введено потому что этот функциональные процессы только один элемент. Для большинства представлений Множества Мандельброта значительное количество элементов останавливается очень рано, и это может сохранить большую обработку. for
цикл был также заменен while
цикл, потому что они обычно более эффективны. Эта функция не упоминает о графическом процессоре и не использует специфичных для графического процессора функций - это - стандартный код MATLAB.
Используя arrayfun
средние значения, что вместо многих тысяч вызовов, чтобы разделить оптимизированные графическим процессором операции (по крайней мере 6 на итерацию), мы выполняем один вызов к параллелизированной операции GPU, которая выполняет целое вычисление. Это значительно уменьшает наверху.
% Setup t = tic(); x = gpuArray.linspace( xlim(1), xlim(2), gridSize ); y = gpuArray.linspace( ylim(1), ylim(2), gridSize ); [xGrid,yGrid] = meshgrid( x, y ); % Calculate count = arrayfun( @pctdemo_processMandelbrotElement, ... xGrid, yGrid, maxIterations ); % Show count = gather( count ); % Fetch the data back from the GPU gpuArrayfunTime = toc( t ); imagesc( x, y, count ) axis off title( sprintf( '%1.3fsecs (GPU arrayfun) = %1.1fx faster', ... gpuArrayfunTime, cpuTime/gpuArrayfunTime ) );
В Экспериментах в MATLAB улучшенная производительность достигается путем преобразования основного алгоритма для C-MEX-функции. Если вы готовы сделать, некоторые работают на C/C++, то можно использовать Parallel Computing Toolbox, чтобы вызвать предзаписанное использование ядер CUDA данные MATLAB. Вы делаете это с parallel.gpu.CUDAKernel
функция.
CUDA/C ++ реализация алгоритма обработки элемента был написан от руки в pctdemo_processMandelbrotElement.cu
: Это должно затем быть вручную скомпилировано с помощью компилятора Nvidia NVCC, чтобы произвести уровень ассемблера pctdemo_processMandelbrotElement.ptx
(.ptx
обозначает "Параллельный язык выполнения Потока").
CUDA/C ++ код немного более включен, чем версии MATLAB, которые мы видели до сих пор, из-за отсутствия комплексных чисел на C++. Однако сущность алгоритма неизменна:
__device__ unsigned int doIterations( double const realPart0, double const imagPart0, unsigned int const maxIters ) { // Initialize: z = z0 double realPart = realPart0; double imagPart = imagPart0; unsigned int count = 0; // Loop until escape while ( ( count <= maxIters ) && ((realPart*realPart + imagPart*imagPart) <= 4.0) ) { ++count; // Update: z = z*z + z0; double const oldRealPart = realPart; realPart = realPart*realPart - imagPart*imagPart + realPart0; imagPart = 2.0*oldRealPart*imagPart + imagPart0; } return count; }
Один поток графического процессора требуется для местоположения во Множестве Мандельброта с потоками, сгруппированными в блоки. Ядро указывает, насколько большой блок потока, и в коде ниже мы используем это, чтобы вычислить количество требуемых блоков потока. Это затем становится GridSize
.
% Load the kernel cudaFilename = 'pctdemo_processMandelbrotElement.cu'; ptxFilename = ['pctdemo_processMandelbrotElement.',parallel.gpu.ptxext]; kernel = parallel.gpu.CUDAKernel( ptxFilename, cudaFilename ); % Setup t = tic(); x = gpuArray.linspace( xlim(1), xlim(2), gridSize ); y = gpuArray.linspace( ylim(1), ylim(2), gridSize ); [xGrid,yGrid] = meshgrid( x, y ); % Make sure we have sufficient blocks to cover all of the locations numElements = numel( xGrid ); kernel.ThreadBlockSize = [kernel.MaxThreadsPerBlock,1,1]; kernel.GridSize = [ceil(numElements/kernel.MaxThreadsPerBlock),1]; % Call the kernel count = zeros( size(xGrid), 'gpuArray' ); count = feval( kernel, count, xGrid, yGrid, maxIterations, numElements ); % Show count = gather( count ); % Fetch the data back from the GPU gpuCUDAKernelTime = toc( t ); imagesc( x, y, count ) axis off title( sprintf( '%1.3fsecs (GPU CUDAKernel) = %1.1fx faster', ... gpuCUDAKernelTime, cpuTime/gpuCUDAKernelTime ) );
Этот пример показал три пути, которыми алгоритм MATLAB может быть адаптирован, чтобы использовать оборудование графического процессора:
Преобразуйте входные данные, чтобы быть на графическом процессоре с помощью gpuArray
, оставление без изменений алгоритма
Используйте arrayfun
на gpuArray
введите, чтобы выполнить алгоритм на каждом элементе входа независимо
Используйте parallel.gpu.CUDAKernel
запускать некоторый существующий CUDA/C ++ использование кода данные MATLAB
title('The Mandelbrot Set on a GPU')