Настройка

Значения ниже задают высоко масштабируемую часть Множества Мандельброта в овраге между основной кардиоидой и p/q лампой с ее левой стороны от него.

1000x1000 сетка действительных частей (X) и мнимых частей (Y) создается между этими пределами, и алгоритм Мандельброта выполнен с помощью итераций в каждом местоположении сетки. Для этого конкретного местоположения 500 итераций будет достаточно, чтобы полностью представить изображение.

maxIterations = 500;
gridSize = 1000;
xlim = [-0.748766713922161, -0.748766707771757];
ylim = [ 0.123640844894862,  0.123640851045266];

Множество Мандельброта в MATLAB

Ниже реализация Множества Мандельброта с помощью стандартных команд 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 обеспечивает версии графического процессора многих функций, которые можно использовать, чтобы создать массивы данных, включая linspacelogspace, и 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 ) );

Работа с CUDA

В Экспериментах в 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 может быть адаптирован, чтобы использовать оборудование графического процессора:

title('The Mandelbrot Set on a GPU')