clear
clc
A=input(\'输入对称正定矩阵A=\')
B=input(\'输入自由项B=\')
n=length(A(:,1));
for k=1:n
   if (det(A(1:k,1:k))<=0)
       input(\'矩阵不是正定矩阵,请重新运行程序\')
   end
end
%分解A=L*L\'
for i=1:n
   t=0;
   for s=1:i-1
       t=t+L(i,s)^2;
   end
   L(i,i)=sqrt(A(i,i)-t);
   for k=i+1:n
       u=0;
       for s=1:i-1
           u=u+L(i,s)*L(k,s);
       end
       L(k,i)=(A(k,i)-u)/L(i,i);
   end
end
%分解AX=B为Ly=B  Lx=y

%求y
for i=1:n
   r=0;
   for k=1:i-1
       r=r+L(i,k)*y(k);
   end
   y(i)=(B(i)-r)/L(i,i);
end
%求x
for i=n:-1:1
   q=0;
   for k=i+1:n
       q=q+L(k,i)*x(k);
   end
   x(i)=(y(i)-q)/L(i,i);
end
x\' %若要输出行向量，此处改为x即可

function pngfanggen(n)
A=[4 2 -4 0 2 4 0 0;
   2 2 -1 -2 1 3 2 0;
  -4 -1 14 1 -8 -3 5 6;
   0 -2 1 6 -1 -4 -3 3;
   2 1 -8 -1 22 4 -10 -3;
   4 3 -3 -4 4 11 1 -4;
   0 2 5 -3 -10 1 14 2;
   0 0 6 3 -3 -4 2 19];
B = [0 -6 20 23 9 -22 -15 45];
L=zeros(8,1);
for j=1:n
    L(j,j)=sqrt(A(j,j)-sum(L(j,:)*L(j,:)\'));
    for i=j+1:n
        L(i,j)=(A(i,j)-sum(L(i,:)*L(j,:)\'))/L(j,j);
    end
end


Y=zeros(1,n);
for i=1:n
    Y(i)=(B(i)-sum(L(i,:)*Y\'))/L(i,i);
end

X(n)=Y(n)/L(n, n);
for i=n-1:-1:1
    X(i)=(Y(i)-sum(L(:,i)\'*X\'))/L(i,i);
end

disp(\'平方根法求得L=\');
disp(L);
disp(\'该线性方程组的解为X=\');disp(X\');
end