%
%  LU 分解的紧凑方式
%
clear all; 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

