%
% 列主元 LU 分解 
%    A(p,:) = L * U
%
clear all; close all;

n = 32;
A = randn(n);

[A_new,p] = PLU(A);
U = triu(A_new); 
L = A_new - U + eye(n);
fprintf('PLU: relerr=%.2e\n', norm(A(p,:)-L*U)/norm(A));

A_new = LU(A);
U = triu(A_new); 
L = A_new - U + eye(n);
fprintf(' LU: relerr=%.2e\n', norm(A-L*U)/norm(A));


%%%%% 列主元 LU 分解： A(p,:) = L * U
function [A,p] = PLU(A)
n = size(A,1);
p = 1 : n; 
for k = 1 : n-1
    [a_max,l] = max(abs(A(k:n,k)));
    if ( a_max == 0 )
        fprintf('Error: 第 %d 步的列主元为 0!\n', k); return;
    end
    l = l+k-1;
    if (l ~= k )
        t = A(k,:); A(k,:) = A(l,:); A(l,:) = t;
        tmp = p(k); p(k) = p(l); p(l) = tmp;
    end
    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

%%%%% LU 分解：A = L * U
function A = LU(A)
n = size(A,1); 
for k = 1 : n-1
    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