% Ali Dorostkar

clearvars
N = 100;
nc = feature('numcores');
parfortimes = zeros(nc,1);
drange_times = zeros(nc,1);
spmd_times = zeros(nc,1);
y = zeros(N,1);
ys = zeros(N,1);

%% Serial computation
tic;
for n = 1:N
    ys(n) = max(svd(randn(n)));
end
serialTime = toc;

a = gcp;
delete(a);
%% Initialize cluster with different number of workers
for i=1:nc
    a = parpool('local',i);
    nw = a.NumWorkers;
    y = zeros(N,1);
    tic;
    parfor n = 1:N
        y(n) = max(svd(randn(n)));
    end
    parfortimes(i) = toc;
%     if (norm(ys - y) > 1e-10)
%         error('The computed answers are wrong');
%     end
    
    y = zeros(N,1);
    tic;
    for n = drange(1:N)
        y(n) = max(svd(randn(n)));
    end
    drange_times(i) = toc;
%     if (norm(ys - y) > 1e-10)
%         error('The computed answers are wrong');
%     end
    
    y = zeros(N,1);
    tic;
    spmd
        NN = floor(N/nw);
        start = (labindex-1)*NN+1;
        fin   = 0;
        if labindex == a.NumWorkers
            fin   = N;
        else
            fin   = labindex*NN;
        end
        for n=start:fin
            y(n) = max(svd(randn(n)));
        end
    end
    spmd_times(i) = toc;
%     if (norm(ys - y) > 1e-10)
%         error('The computed answers are wrong');
%     end
    
    delete(a);
end

%% plot results
semilogy(parfortimes, '-*');
xlabel('Number of workers')
ylabel('Time to compute')
figure;
plot(parfortimes/serialTime, '-*');
hold on
plot(drange_times/serialTime, '-ro');
plot(spmd_times/serialTime, '-cd');
hold off
xlabel('Number of workers')
ylabel('Speedup relative to serial run')
legend('parfor', 'drange','spmd')
