% function Checkerboard_mzip
%   MATLAB zip file, a self-extracting MATLAB archive.
%   Usage: Run this file to recreate the original directory.

    fname = mfilename;
    fin = fopen([fname '.m'],'r');
    dname = fname(1:find(fname=='_',1,'last')-1);
    mkdir(dname);
    mkdir([dname filesep 'lib'])
    addpath(dname)
    addpath([dname filesep 'lib'])    
    
    L = fgetl(fin);
    while length(L) < 2 || ~isequal(L(1:2),'%%')
        L = fgetl(fin);
    end
    while ~isequal(L,'%% EOF')
        F = [dname filesep L(4:end)];
        disp(F)
        fout = fopen(F,'w');
        L = fgetl(fin);
        while length(L) < 2 || ~isequal(L(1:2),'%%')
            fprintf(fout,'%s\n',L);
            L = fgetl(fin);
        end
        fclose(fout);
    end
    fclose(fin);
% end
%% imagem.m
function imagem(A)
    p = parula;
    p(end,:) = [.95 .95 .95];
    colormap(p)
    dpos = get(0,'defaultFigurePosition');
    pos = 3/4*dpos.*[1 1 .5 .5];
    set(gcf,'pos',pos)
    imagesc(A)
    axis image
end
%% plotem.m
function plotem(s1,s2,logylim)
    co = get(gca,'colororder');
    dpos = get(0,'defaultFigurePosition');
    pos = 3/4*dpos;
    set(gcf,'pos',pos)
    clf
    if nargin == 1
        plot(s1,'.-','markersize',10,'color',co(1,:));
        set(gca,'ytick',[0 s1(1)/2 s1(1)])
        axis padded
    else
        if nargin == 2
            h1 = semilogy(s1,'.-','markersize',10,'color',co(1,:));
            axis padded
            ym = s1(1);
        else
            h1 = semilogy(s1,'.-','markersize',10,'color',co(1,:));
            hold on
            h2 = semilogy(s2,'o-','markersize',3,'color',co(2,:));
            hold off
            axis padded
            set(gca,'ylim',10.^logylim)
            ym = max(s1(1),s2(1));
        end
        em = eps(ym);
        ax = xlim;
        line(ax,[ym ym],'color','k')
        line(ax,[em em],'color','k')
        if nargin < 3
            legend(h1, 's1', ...
                'location','best')
        else
            legend([h1 h2], {'s1','s2'}, ...
                'location','best')
        end
    end
end
%% timefit.m
t = [ ...
    0.0053    0.0097
    0.0164    0.0336
    0.0339    0.0737
    0.0845    0.1824
    0.1537    0.3203
    0.2693    0.4800
    0.5271    0.9695
    0.8001    1.2523
    1.4335    2.1614
    2.1234    3.1476
    2.9524    4.2736
    4.0253    5.4319
    5.4193    7.2862
    6.9529    9.1811
    8.7620   11.5007
   10.8354   14.1059
   13.0821   16.8855
   15.6420   20.2465
   18.8214   24.0954
   22.1447   28.0566];
k = (1:20)';
A = [ones(20,1) k k.^2 k.^3];
c = A\t;
ratio = c(end,2)/c(end,1)
f  = A*c;
m = 250;
n = m*k;
plot(n,t,'o',n,f,'-','linewidth',1.5)
axis padded
xlabel('n','fontweight','bold')
ylabel('t','fontweight','bold')
title('time, (seconds)')
legend({'s1','s2','f1','f2'},'location','northwest')
text(3000,25,sprintf('ratio = %5.2f',ratio), ...
    'fontweight','bold')

%% timem.m
m = 250;
co = get(gca,'colororder');
clf
axis([0 21*m 0 6])
hold on
shg
for j = 1:6
    t = zeros(0,2);
    for k = 1:20
        n = m*k;
        A = gallery('minij',n);
        rnk = rank(A);
        title(rnk)
        tic , s1 = svd(A); t(k,1) = toc;
        tic, [~,s2,~] = svd(A,'vector'); t(k,2) = toc;
        c = (1:k)';
        y = 1000*t(c,:)./(c.^3);
        line(m*c,y, ...
            'marker','o', ...
            'linewidth',1.5, ...
            'linestyle','-', ...
            'color',co(j,:))
        drawnow
    end
    print -dpng times.png
    t
end

%% EOF