function P1()
% Tests QueenLU
clc, clear, randn('seed',0), format short

% Example 1.
disp('Example 1.')
n = 5;
A = randn(n,n); A(1,1) = .001;
[L,U,rowpiv,colpiv] = QueenLU(A)

% Example 2.
disp('Example 2')
disp(' ')
n = 20;
randn('seed',0)
A = randn(n,n);
[L,U,rowpiv,colpiv] = QueenLU(A);
disp(['rowpiv = ' sprintf('%3d',rowpiv)])
disp(['colpiv = ' sprintf('%3d',colpiv)])


function [L,U,rowpiv,colpiv] = QueenLU(A)
% A is an n-by-n matrix
% L is unit lower triangular
% U is upper triangular
% rowpiv and colpiv are integer (n--vectors with the property that
% PAQ = LU where
%   P = P_{n-1}*...*P1     Pk = I with rows rowpiv(k) and k swapped
%   Q = Q1*...*Q_{n-1}     Qk = I with rows colpiv(k) and k swapped
%
[n,n] = size(A);
rowpiv = 1:n-1;
colpiv = 1:n-1;
for k=1:n-1
    [diagM,tau] = max(abs(diag(A(k:n,k:n)))); tau = tau+k-1;
    if diagM >= norm(A(k:n,k),inf) && diagM >= norm(A(k,k:n),inf)
        rowpiv(k) = tau; colpiv(k) = tau;
        A([k tau],:) = A([tau k],:);
        A(:,[k tau]) = A(:,[tau k]);
    else
        [colM,mu]     = max(abs(A(k:n,k))); mu = mu+k-1;
        [rowM,lambda] = max(abs(A(k,k:n))); lambda = lambda+k-1;
        if colM>=rowM 
            % Interchange rows...
            rowpiv(k) = mu; colpiv(k) = k; 
            A([mu k],:) = A([k mu],:);
        else
            % Interchange columns...
            rowpiv(k) = k; colpiv(k) = lambda; 
            A(:,[lambda k]) = A(:,[k lambda]);
        end
    end
    if A(k,k) ~=0
       A(k+1:n,k) = A(k+1:n,k)/A(k,k);
       A(k+1:n,k+1:n) = A(k+1:n,k+1:n) - A(k+1:n,k)*A(k,k+1:n);
    end
end
L = eye(n,n) + tril(A,-1);
U = triu(A);

% Bound Comment
% If in step k we choose pivot from diagonal
% then norm(L(k+1:n,k),inf) <= 1 and norm(U(k,k+1:n)/U(k,k),inf) <= 1
% If in step k we choose pivot from A(k:n,k)
% then norm(L(k+1:n,k),inf) <= 1 and norm(U(k,k+1:n),inf) > 1 maybe
% If in step k we choose pivot from A(k,k:n)
% then norm(L(k+1:n,k),inf) > 1 maybe and norm(U(k,k+1:n)/U(k,k),inf) <= 1








    