% JULIA(C,K,V) draws the Julia set with the following parameters:
% c is a complex number used in the map f(z) = z.^2+ c.
% k gives the number of iterations
% v determines the number of points on the x-axis.
% JULIA uses c = 0.2+0.65i, k = 14, v = 500.
% default settings
if nargin < 3
c = 0.2+0.65i;
k = 200;
v = 500;
end
% c=gsingle(c);
% radius of the circle beyond which every point diverges
r = max(abs(c),2);
% divide the x-axis
d = linspace(-r,r,v);
% create the matrix A containing complex numbers
A = ones(v,1)*d+i*(ones(v,1)*d)';
% create the point matrix
B = zeros(v,v);
gA=gsingle(A);
gB=gsingle(B);
%%%CPU
tic
% iteration
for s = 1:k
B = B+(abs(A)<=r);
% the map
% A=sin(A)+exp(A)+gones(v,v).*c;
A =power(A,5)+power(A,2)+ones(v,v).*c;
end;
toc
% plot settings
imagesc(B);
%%%GPU
tic;
gsync;
% iteration
for s = 1:k
gB = gB+(abs(gA)<=r);
% the map
% A=sin(A)+exp(A)+gones(v,v).*c;
gA =power(gA,5)+power(gA,2)+gones(v,v).*c;
end;
gsync;
toc
isequal(B-gB)
imagesc(gB);


Elapsed time is 69.214393 seconds.
Elapsed time is 1.130295 seconds.
ans =
1