function P2()
clc
randn('seed',0);
n = 30;
A = randn(n,n); 
b = randn(n,1);
[L,U,rowpiv,colpiv] = QueenLU(A);
[p,q] = Convert(rowpiv,colpiv);
disp(['p = ' sprintf('%3d',p)])
disp(['q = ' sprintf('%3d',q)])
disp(' ')
I = eye(n,n); P = I(p,:); Q = I(q,:);
errLU = norm(P*A*Q - L*U,1)
y = L\b(p);
z = U\y;
x = z(q);
errSolve = norm(A\b - x)

function [p,q] = ConvertCVL(rowpiv,colpiv)
% rowpiv and colpiv are row vectors of integers with length n-1
% p and q are row vectors of integers of length n with the
% property that 
%    P = I(p,:) = Pnm1*..*P1  Pk = I with rows k and rowpiv(k) swapped
%    Q = I(q,:) = Q1*..*Qnm1  Qk = I with rows k and colpiv(k) swapped


n = length(rowpiv)+1;
p = 1:n;
q = 1:n;
for k=1:n-1
    p([rowpiv(k) k]) = p([k rowpiv(k)]);
end

% Note the reverse order...
for k=n-1:-1:1
    q([colpiv(k) k]) = q([k colpiv(k)]);
end

%-2 if q is wrong, i.e., if you had k=1:n-1 in the second loop


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);

