%
%  LU 分解的紧凑方式
%
clearvars; 
close all;

n = 32;
A = randn(n);

[L,U] = LU(A);

fprintf('||A-L*U||=%.2e\n', norm(L*U-A));
subplot(1,3,1), spy(A), title('A');
subplot(1,3,2), spy(L), title('L');
subplot(1,3,3), spy(U), title('U');
shg


%%%%% LU 分解
function [L,U] = LU(A)
[m,n] = size(A);
if ( m ~= n )
    fprintf('\n Error: A 不是方阵!\n'); return; 
end

L = eye(n); 
U = zeros(n);
for k = 1 : n
    U(k,k) = A(k,k) - sum(L(k,1:k-1)'.*U(1:k-1,k));
    if (U(k,k) == 0)
        fprintf('\n Error: 第 %d 个主元＝0 !\n',k); return; 
    end
    for j = k+1 : n
        U(k,j) = A(k,j) - sum(L(k,1:k-1)'.*U(1:k-1,j));
        L(j,k) = (A(j,k) - sum( L(j,1:k-1)'.*U(1:k-1,k) ))/U(k,k);
    end
end
end

