cholesky.asv

用matlab编写的cholesky分解程序

I=eye(M);
for j=1:M
    for i=1:M
        if i==j+1
        H(i,j)=1;
        else H(i,j)=0;
        end
    end
end
J=H';
A=[I,H;J,I] ; 
A=[1,0,0,0,0,0,0,0;0,1,0,0,1,0,0,0;0,0,1,0,0,1,0,0;0,0,0,
n=length(A)
lambda=eig(A);
flagdc=0;
flagzd=0;
while(flagdc~=n*(n-1)/2)|(flagzd~=n)
    for i=1:n
        for j=i+1:n
            if(A(i,j)==A(j,i))
                flagdc=flagdc+1;
            end
        end
    end
    if(flagdc~=n*(n-1)/2)
        disp('矩阵A不是对称矩阵，请重新输入矩阵A:')
        A=input('A=\n');
        disp('矩阵的阶数为:')
        n=length(A)
        flagdc=0;
    end
    if flagdc==n*(n-1)/2
        disp('矩阵A的特征值是:')
        lambda=eig(A)
        for i=1:n
            if(lambda(i)>0)
                flagzd=flagzd+1;
            end
        end
        if flagzd~=n
            disp('矩阵A不是正定矩阵。请重新输入矩阵A:')
            A=input('A=\n');
            disp('矩阵的阶数为：')
            n=length(A)
            flagzd=0;
            flagdc=0;
        end
    end
end
l(1,1)=sqrt(A(1,1));
for j=2:n
    l(j,1)=A(j,1)/l(1,1);
end
for i=2:n
    for j=i:n
        p(j,i)=0;
        for k=1:i-1
            p(j,i)=p(j,i)+l(i,k)*l(j,k);
            t(j,i)=A(j,i)-p(j,i);
            if (j==i)
                l(j,i)=A(j,i)-p(j,i);
                if (j==i)
                    l(j,i)=sqrt(t(j,i));
                else
                    l(j,i)=t(j,i)/l(i,i);
                end
            end
        end
    end
end