%
% Cholesky 分解
%
clear all; close all;

n = 32;
A = randn(n);
A = A'*A;

L = Cholesky(A);

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

%%%%% Cholesky Decomposition 
function L = Cholesky(A)
n = size(A,1);
L = zeros(n);
for j = 1 : n
    L(j,j) = A(j,j) - dot(L(j,1:j-1),L(j,1:j-1));
    L(j,j) = sqrt(L(j,j));
    for i = j+1 : n
        L(i,j) = (A(i,j) - dot(L(i,1:j-1),L(j,1:j-1)) )/L(j,j);
    end
end
end
