function LUTool_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
%% LUTool.m
function [Aout,pout,qout,Lout,Uout] = LUTool(name,n,pivs,paws,wantgif)
    % LUTool Gaussian elimination demonstration.
    %
    %   LUTool(A) computes a unit lower triangular matrix L, an upper
    %   triangular matrix U, and permutation vectors p and q so that
    %
    %       A(p,q) = L*U
    %
    %   The LUTool interface lets you select the matrix family and order,
    %   choose the Gaussian elimination pivoting strategy,
    %   increase or decrease the matrix order, and speed up or
    %   slow down the animation.
    %
    %   LUTool, with no input arguments, uses the defaults in LUTool.m .
    %
    %   LUTool(A), with one input argument, uses the matrix A .
    %
    %   LUTool('name',n,'pivot',paws), with four input arguments, uses
    %       name    matrix family
    %       n       order
    %       pivot   pivoting strategy, 'none', 'partial' or 'complete', 
    %       paws    unit pause
    %   Names of matrix families are in the string array mats in LUTool.m .
    %
    %   [A,p,q,L,U] = LUTool(..) returns five outputs with A(p,q) = L*U.
    % 
    %   Example:
    %       [A,p,q,L,U] = LUTool('magic',5,'partial',0.05)
    %
    %   Blog: https://blogs.mathworks.com/cleve/2025/02/08.
    %
    %   Feb. 11, 2025
    
    %   Copyright 2013-2025 Cleve Moler
    
    mats = {"    ","magic","pascal","golub","wilkinson","frank", ...
        "rando","fiedler","redheff","redheff_3","binomial","gcdmat"};

    lufigure
    if nargin == 0
        name = 'frank';
        n = 2;
        A = lugallery(name,n);
        pivs = 'none';
        paws = 0.05;
        wantgif = false;
    elseif nargin == 1
        A = name;
        name = 'mymat';
        [n,~] = size(A);
        pivs = 'partial';
        paws = .02;
        wantgif = false;
    else
        A = lugallery(name,n);
        if nargin < 5
            wantgif = false;
        end
    end
    set(gcf,'userdata',{name,n,pivs,paws})
    lucontrols(mats,pivs,wantgif);
    T = luinit(A);
    lupawser('start')

    [L,U,p,q] = LUTooler(A,T,pivs,name);
    
    if nargout > 0
        Aout = A; pout = p; qout = q; Lout = L; Uout = U;
    end
    lupawser(80)
    lupawser('wrap')
end
%% LUTooler.m
function [L,U,pv,qv] = LUTooler(A,T,pivs,name)
    Lcolor = nytcolors(3);
    Ucolor = nytcolors(4);
    Pivotcolor = nytcolors(2);
    ElimColor = nytcolors(5);
    [n,~] = size(A);
    titl = text(.40,.65+.031*n,[name '(' int2str(n) ')'], ...
        'fontsize',16,'fontweight','bold','horiz','center', ...
        'interpreter','none');
    axis off
    stop = findobj('string','X');
    lupawser(40)

    % Elimination
    
    n = size(A,1);
    pv = 1:n;
    qv = 1:n;
    dsig = 1;
    for k = 1:n
    
    %   If column is zero, abort.

        if all(all(A(k:n,k:n)==0))
            for l = k:n
                for i = l+1:n
                    set(T(i,l),'back',Lcolor)
                end
                for j = l:n
                    set(T(l,j),'back',Ucolor)
                end
            end
            continue
        end
    
        p = 0;
        q = 0;
        while p < k || q < k || p > n || q > n

            switch pivs

                case "none"
                    p = k;
                    q = k;
                
                case "partial"
                    p = find(abs(A(k:n,k)) == max(abs(A(k:n,k))),1)+k-1;
                    q = k;
                
                case "complete"
                    [p,q] = find(abs(A(k:n,k:n)) == max(max(abs(A(k:n,k:n)))),1);
                    p = p(1)+k-1;
                    q = q(1)+k-1;
                    
            end
            set(T(p,q),'back',Pivotcolor)
        end
        lupawser(40)
          
        u = 15;
        w = .9;
        x = (u-4-n)/2;  
        y = (u-n)/2; 

    %   Swap columns

        if q > k
            A(:,[q,k]) = A(:,[k,q]);
            qv([q,k]) = qv([k,q]);
            dsig = -dsig;
            for s = .10:.10:1
                for i = 1:n
                    set(T(i,k),'Pos',[(x+k+s*(q-k)) (y+n-i) w w]/u)
                    set(T(i,q),'Pos',[(x+q+s*(k-q)) (y+n-i) w w]/u) 
                end
                lupawser(0.33)
            end
            T(:,[q,k]) = T(:,[k,q]);
            for i = 1:n
                set(T(i,k),'string',luprint(A(i,k)),'userdata',[i k])
                set(T(i,q),'string',luprint(A(i,q)),'userdata',[i q])
            end
            lupawser(1)
        end
    
    %   Swap rows
    
        if p > k
             A([p,k],:) = A([k,p],:);
             pv([p,k]) = pv([k,p]);
             dsig = -dsig;
             for s = .10:.10:1
                 for j = 1:n
                     set(T(k,j),'Pos',[(x+j) (y+n-(k+s*(p-k))) w w]/u)
                     set(T(p,j),'Pos',[(x+j) (y+n-(p+s*(k-p))) w w]/u)
                 end
                 lupawser(0.33)
             end
                
             T([p,k],:) = T([k,p],:);
             lupawser(1)
             for j = k:n
                 set(T(k,j),'string',luprint(A(k,j)),'userdata',[k j])
                 set(T(p,j),'string',luprint(A(p,j)),'userdata',[p j])
             end
             lupawser(10)
         end

         if all(A(k,k) ~= 0)

    %        Compute multipliers in L
    
             for i = k+1:n
                 A(i,k) = A(i,k)/A(k,k);
                 set(T(i,k),'string',luprint(A(i,k)),'back',Lcolor)
                 lupawser(1)
             end

             for j = k:n
                 set(T(k,j),'string',luprint(A(k,j)),'back',Ucolor)
                 drawnow
             end
             lupawser(4)
  
    %        Elimination    
    
             for j = k+1:n
                 for i = k+1:n
                     A(i,j) = A(i,j) - A(i,k)*A(k,j);
                     set(T(i,j),'string',luprint(A(i,j)),'back',ElimColor)
                     if A(i,j) == 0
                         lupawser(0.01)
                     else
                         lupawser(1)
                     end
                 end
             end
             i = k+1:n;
             set(T(i,i),'back','w')
             lupawser(40)
         end
         lupawser(4)
         if stop.Value
             break
         end
     end
    
    % Seperate L and U into two matrices
    
    for s = 0.10:0.10:1.0
        for j = 1:n
            for i = 1:n
                if i <= j
                    set(T(i,j),'Pos',[(x+j+1.0*s) (y+n+1.0*s-i) w w]/u)
                else
                    set(T(i,j),'Pos',[(x+j-1.0*s) (y+n-1.0*s-i) w w]/u)
                end
            end
        end
        set(titl,'pos',get(titl,'pos')+[0 s/30 0])
        lupawser(1)
    end

    % Insert ones on diagonal of L
    
    for j = 1:n
        c = lupushbut('1',[(x+j-1) (y+n-j-1) 1 1]/u, '');
        set(c,'BackgroundColor',Lcolor)
    end
    lupawser(1)
    
    % Pivot Indices

    for i = 1:n
        pvi = lupushbut(int2str(pv(i)),[x+i-1 y-2.5 0.5 .75]/u, '' );
        set(pvi,'BackgroundColor',Lcolor)
        qvi = lupushbut(int2str(qv(i)),[x+i+1 y+n+2 0.5 .75]/u, '');
        set(qvi,'BackgroundColor',Ucolor)  
    end
    lupawser(1)

    % Determinant

    dt = dsig*prod(diag(A));
    if min(abs(diag(A))) < 10*n*eps(max(abs(diag(A))))
        dtpr = '0.0';
    else
        dtpr = luprint(dt);
    end
    dbut = lupushbut(dtpr,[x-0.5 y+n+2 1.5 0.75]/u, '');
    set(dbut,'BackgroundColor',Ucolor)
    lupawser(1)

    L  = tril(A,-1) + eye(size(A));
    U  = triu(A);
end
%% lib\golub.m
function A = golub(n)
    %GOLUB  Badly conditioned integer test matrices.
    %   GOLUB(n) is the product of two random integer n-by-n matrices,
    %   one of them unit lower triangular and one unit upper triangular.
    %   LU factorization without pivoting fails to reveal that such
    %   matrices are badly conditioned.
    s = 10;
    L = tril(round(s*randn(n)),-1)+eye(n);
    U = triu(round(s*randn(n)),1)+eye(n);
    A = L*U;
end
%% lib\lucontrols.m
function lucontrols(mats,pivs,wantgif)
    % Control buttons
    u = 15;
    x = u-3;  
    y = u-2.5;
    B(1) = lupopup(        mats, [ x y-1/4 3/2 1]/u, @lupopupcb);
    B(2) = lupushbut(  "repeat", [ x y-1 3/2 1]/u, @lucontrolscb);
    B(3) = lutoggbut(    "none", [ x y-2 3/2 1]/u, @lucontrolscb);
    B(4) = lutoggbut( "partial", [ x y-3 3/2 1]/u, @lucontrolscb);
    B(5) = lutoggbut("complete", [ x y-4 3/2 1]/u, @lucontrolscb);
    B(6) = lupushbut(  "slower", [ x y-5 3/2 1]/u, @lucontrolscb);
    B(7) = lupushbut(  "faster", [ x y-6 3/2 1]/u, @lucontrolscb);
    B(8) = lupushbut(      "n-", [ x y-7 3/2 1]/u, @lucontrolscb);
    B(9) = lupushbut(      "n+", [ x y-8 3/2 1]/u, @lucontrolscb);
    B(10) = lupushbut(   "help", [ x y-9 3/2 1]/u, @luinfo);
    B(11) = lupushbut(   "info", [ x y-10 3/2 1]/u, @luinfo);
    B(12) = lutoggbut(    "gif", [ x y-11 3/2 1]/u, @lugifcb);
    B(13) = lutoggbut("X", [u-3/4 u-3/4 3/4 3/4]/u, @lustopcb);
    fud = get(gcf,'userdata');
    [name,~,~,~] = deal(fud{:});
    nyt2 = nytcolors(2);
    if isequal(name,'mymat')
        B(1).Value = 1;
    else
        B(1).Value = find([mats{:}]==name);
    end
    B(1).BackgroundColor = nyt2;
    switch pivs
        case "none", B(3).BackgroundColor = nyt2;
        case "partial", B(4).BackgroundColor = nyt2;
        case "complete", B(5).BackgroundColor = nyt2;
    end
    if wantgif
        B(12).Value = 1;
        B(12).BackgroundColor = nyt2;
    end
end
%% lib\lucontrolscb.m
function lucontrolscb(arg,~)
    fud = get(gcf,'userdata');
    [name,n,pivs,paws] = deal(fud{:});
    switch arg.String
        case 'n-',       n = max(1,n-1);
        case 'n+',       n = min(15,n+1);
        case 'slower',   paws = 1.5*paws;
        case 'faster',   paws = paws/1.5;
        case 'repeat'  
        otherwise,       pivs = arg.String;
    end
    LUTool(name,n,pivs,paws)
end         
%% lib\lufigure.m
function lufigure
    close
    pos = get(0,'screensize')-[0 0 0 40];
    set(gcf,'name','LUTool', ...
        'numbertitle','off', ...
        'toolbar','none', ...
        'menubar','none', ...
        'inverthardcopy','off', ...
        'paperpositionmode','auto', ...
        'position',pos)
end
%% lib\luframe.m
function luframe(arg1,arg2)
    persistent filename
    im = frame2im(getframe(gcf));
    gif_shrink = 5/16;
    im = imresize(im,gif_shrink);
    [imind,cm] = rgb2ind(im,256);
    if isequal(arg1,'wrap')
        delay = 5;
        imwrite(imind,cm,filename,'gif','WriteMode','append','DelayTime',delay);
    elseif ischar(arg1)
        filename = arg1;
        delay = 2;
        imwrite(imind,cm,filename,'gif','Loopcount',arg2,'DelayTime',delay)
    else
        imwrite(imind,cm,filename,'gif','WriteMode','append','DelayTime',arg1)
    end
end
%% lib\lugallery.m
function A = lugallery(name,n)
    try
        A = feval(name,n);
    catch
        A = gallery(name,n);
    end
end
%% lib\lugifcb.m
function lugifcb(arg,~)
    if arg.Value == 0
        arg.BackgroundColor = 'w';
    else
        arg.BackgroundColor = nytcolors(2);
    end
end
%% lib\luinfo.m
function luinfo(arg,~)
    switch arg.String
        case "help"
            doc("LUTool")
        case "info"
            gcfud = get(gcf,'userdata');
            name = gcfud{1};
            w = which(name);
            if isempty(w)
                helpwin(['private/' name])
            else
                helpwin(name)
            end
    end
end         
%% lib\luinit.m
function T = luinit(A)
    delete(findobj('type','text'))
    ch = flipud(get(gcf,'children'));
    delete(ch(15:end))
    [m,n] = size(A);
    u = 15;
    x = (u-4-n)/2; 
    y = (u-n)/2; 
    T = repmat(uicontrol('style','text'),m,n);
    for k = 1:m
        for j = 1:n
            T(k,j) = lutext(luprint(A(k,j)),[x+j y+m-k 1 1]/u);
            set(T(k,j),'userdata',A(k,j),'backgroundcolor','w')       
         end
    end
end    
%% lib\lupawser.m
function lupawser(paws)
    gifobj = findobj('String','gif');
    gif = gifobj.Value == 1;
    if isequal(paws,'start')
        if gif
            gcfud = get(gcf,'userdata');
            n = gcfud{2};
            fname = ['html/' gcfud{1} '_' num2str(n) '_' gcfud{3} '.gif'];
            luframe(fname,inf)
        end
    elseif isequal(paws,'wrap')
        if gif
            luframe('wrap')
        end
    else
        gcfud = get(gcf,'userdata');
        pause(paws*gcfud{4})
        if gif
            luframe(paws*gcfud{4})
        end
    end
end
%% lib\lupopup.m
 function hout = lupopup(string,position,callback)
    position(3:4) = 0.9*position(3:4);
    h = uicontrol(Style = "popupmenu", ...
        String = string, ...
        Units = "normalized", ...
        Position = position, ...
        Callback = callback, ...
        Fontname = "LucidaConsole", ...
        Fontsize = 16, ...
        Fontweight = "bold", ...
        BackgroundColor = 'w');
    if nargout > 0
        hout = h;
    end
end
%% lib\lupopupcb.m
function lupopupcb(arg,~)
    fud = get(gcf,'userdata');
    [~,n,pivs,paws] = deal(fud{:});
    name = arg.String{arg.Value};
    LUTool(name,n,pivs,paws)
end         
%% lib\luprint.m
function s = luprint(x)
    % LUTool sprintf
    if abs(x-round(x)) < 100*eps(x)
        f = '%8.0f';
    elseif abs(x) >= 100 && abs(x) < 1000
        f = '%8.1f';
    elseif abs(x) < 0.01 || abs(x) >= 1000
        f = '%8.0e';
    else
        f = '%8.2f';
    end 
    s = sprintf(f,x);
end
%% lib\lupushbut.m
function hout = lupushbut(string,position,callback)
    position(3:4) = 0.9*position(3:4);
    h = uicontrol(Style = "pushbutton", ...
        String = string, ...
        Units = "normalized", ...
        Position = position, ...
        Callback = callback, ...
        Fontname = "LucidaConsole", ...
        Fontsize = 16, ...
        Fontweight = "bold", ...
        BackgroundColor = 'w');
    if nargout > 0
        hout = h;
    end
end
%% lib\lustopcb.m
function lustopcb(arg,~)
    if arg.Value
        arg.BackgroundColor = nytcolors(2);
        pause(0.5)
        close(gcf)
    else
        arg.BackgroundColor = 'w';
    end    
end         
%% lib\lutext.m
function h = lutext(string,position)
    position(3:4) = 0.9*position(3:4);
    h = uicontrol(Style = "text", ...
        String = string, ...
        Units = "normalized", ...
        Position = position, ...
        Fontname = "LucidaConsole", ...
        Fontsize = 16, ...
        Fontweight = "bold", ...
        BackgroundColor = 'w');
end
%% lib\lutoggbut.m
function togg = lutoggbut(str,position,callback)
    position(3:4) = 0.9*position(3:4);
    togg = uicontrol(Style = "toggle", ...
        String = str, ...
        Units = "normalized", ...
        Position = position, ...
        Callback = callback, ... 
        Fontname = "LucidaConsole", ...
        Fontsize = 16, ...
        Fontweight = "bold", ...
        BackgroundColor = 'w');
end
           
%% lib\nytcolors.m
function nyt = nytcolors(k)
   % Colors from New York Times Games.
   nytco = uint8( ...
      [184   228    68
       248   211     0
       121   160   238
       184   113   198
       192   192   192]);
   if nargin == 0
       nyt = nytco;
   else
       nyt = nytco(k,:);
   end
end
%% lib\redheff_3.m
function R = redheff_3(n)
    % R = redheffs(n) is redheffer(n) with first and last columns
    % swapped to improve sparsity.  See "Computing Mertens #3" in
    % https://blogs.mathworks.com/cleve/2024/10/22/mobius-mertens-and-redheffer/
    R = gallery("redheff",n);
    R(:,[n 1]) = R(:,[1 n]);
end
%% EOF