function [L,X]= Cholesky(A,b)
N=size(A);%记录矩阵的大小
m=N(1);%矩阵的行数
n=N(2);%矩阵的列数
%判断输入矩阵是否为正定矩阵
if m~=n
disp('输入的矩阵不是方阵,请重新输入')
return;
end
%判断输入的矩阵是否为正定矩阵
d=eig(A);
for i=1:n
if d(i)<0
disp('输入的矩阵不是正定矩阵,请重新输入')
return;
else
break;
end
end
L=zeros(m,n);%初始化L矩阵
L(1,1)=sqrt(A(1,1));
for k=2:n
L(k,1)=A(k,1)/L(1,1);
end
for k=2:n
L(k,k)=sqrt(A(k,k)-L(k,1:k-1)*L(k,1:k-1)');
L(k+1:m,k)=(A(k+1:m,k)-L(k+1:m,1:k-1)*L(k,1:k-1)')/L(k,k);
if L(k,k)==0
break
end
end
% 前代法
y=L\b;
X=L'\y;
## 测试代码部分
M=diag(rand(100,1));
U=orth(rand(100,100));
A=U'*M*U;
b=rand(100,1);
[L,X]=Cholesky(A,b)
eps=norm(A*X-b);
display(eps)