clc,clear

% A = [4 -1 1; -1 4.25 2.75; 1 2.75 3.5];
A = gallery('lehmer',6);


L = zeros(size(A));
for n = 1:size(A,1)
    
    if n==1
       L(n,n) = sqrt( A(n,n) ); 
       L(n+1:end,n) = A(n,n+1:end)/L(n,n);
       continue;
    end
    
    L(n,n) = sqrt( A(n,n)-sum(L(n,1:n-1).^2) );
    
    if n~=size(A,1) % 避免矩阵溢出
        M = L(n+1:end,1:n-1).*(ones(size(A,1)-n,1)*L(n,1:n-1)); %为下一步计算l_ij中的求和项
        L(n+1:end,n) = ( A(n+1:end,n)- sum(M,2) )/L(n,n);
    end
    
end