%%
%MATLAB on GPU basics

gpuDevice %find properties of installed GPU card

%Over 200 operations (see http://www.mathworks.com/help/distcomp/run-built-in-functions-on-a-gpu.html) 
%are available as built-in MATLAB functions that can be executed directly on the GPU by 
%providing an input argument of the type gpuArray. These GPU-enabled functions are overloaded - 
%in other words, they operate differently depending on the data type of the arguments passed to them.

A=rand(5000);
tic, eig(A); time_cpu=toc %how fast is diagonalization on multicore CPU
B=gpuArray(A); %push random matrix A into fast memory on GPU card
class(A), class(B)
tic, eig(B); time_gpu=toc, speedup=time_cpu/time_gpu %find speedup after diagonalization on GPU

%%
%Proper time measurement should also include transfer between fast memory
%on GPU and slow memory on CPU
tic;
A1=rand(5000);
B1 = fft(A1);
time_cpu = toc
tic;
A1=rand(5000);
A1 = gpuArray(A1);
B1 = fft(A1);
B1 = gather(B1);
time_gpu = toc
speedUp = time_cpu/time_gpu;
disp(speedUp) %speedup below 1 means that it is inefficient to use GPU

%%
%Creating Mandelbrot Set (an example of fractal) using MATLAB
%See also http://www.mathworks.com/moler/exm/chapters/mandelbrot.pdf 
%in C. Moler, "Experiments with MATLAB"

%CPU version:
maxIterations = 500;
gridSize = 1000;
xlim = [-0.748766713922161, -0.748766707771757]; % define zoomed part of Mandelbrot Set
ylim = [ 0.123640844894862,  0.123640851045266]; % define zoomed part of Mandelbrot Set

% 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 = zeros( size(z0) );

% Calculate
z = z0;
for n = 0:maxIterations
    z = z.*z + z0;
    inside = abs( z )<=2;
    count = count + inside;
end
count = log( count+1 );

% Plot
cpuTime = toc( t );
set( gcf, 'Position', [200 200 600 600] );
figure(1);
imagesc( x, y, count );
axis image
colormap( [jet();flipud( jet() );0 0 0] );
title( sprintf( '%1.2fsecs (without GPU)', cpuTime ) );

%%
%GPU version I: 
%gpuArray provides GPU versions of many functions that can be used to 
%create data arrays, including the linspace, logspace, and meshgrid functions needed here. 
%Similarly, the count array is initialized directly on the GPU using the
%function gpu.GPUArray.zeros; other than these simple changes to the data initialization, 
%the algorithm is unchanged. 

maxIterations = 500;
gridSize = 1000;
xlim = [-0.748766713922161, -0.748766707771757]; % define zoomed part of Mandelbrot Set
ylim = [ 0.123640844894862,  0.123640851045266]; % define zoomed part of Mandelbrot Set

% 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 = zeros( 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+1 );

% Plot
naiveGPUTime = toc( t );
figure(2);
imagesc( x, y, count )
axis image
title( sprintf( '%1.2fsecs (naive GPU) = %1.1fx faster', ...
    naiveGPUTime, cpuTime/naiveGPUTime ) )


%%
%GPU version II:
%The algorithm is operating equally on every element of the input, so we can 
%place the code in a helper function and call it using arrayfun. For GPU array 
%inputs, the function used with arrayfun gets compiled into native GPU code. 
%In this case we placed the loop in processMandelbrotElement.m which looks as follows:

%function count = processMandelbrotElement(x0,y0,maxIterations)
%z0 = complex(x0,y0);
%z = z0;
%count = 1;
%while (count <= maxIterations) ...
%      && ((real(z)*real(z) + imag(z)*imag(z)) <= 4)
%    count = count + 1;
%    z = z*z + z0;
%end
%count = log(count);

%Note that an early abort has been introduced since this function only processes a single element. 
%For most views of the Mandelbrot Set a significant number of elements stop very early and this can 
%save a lot of processing. The for loop has also been replaced by a while loop because they are usually 
%more efficient. This function makes no mention of the GPU and uses no GPU-specific features - it is standard 
%MATLAB code.

%Using arrayfun means that instead of many thousands of calls to separate GPU-optimized operations (at 
%least 6 per iteration), we make one call to a parallelized GPU operation that performs the whole 
%calculation:

% 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( @processMandelbrotElement, xGrid, yGrid, maxIterations );

% Show
gpuArrayfunTime = toc( t );
figure(3);
imagesc( x, y, count )
axis image
title( sprintf( '%1.2fsecs (GPU arrayfun) = %1.1fx faster', ...
    gpuArrayfunTime, cpuTime/gpuArrayfunTime ) );