ldltduicheng.m

来自「这是线性议程组的直接解法」· M 代码 · 共 33 行

function [L,D,LT,x]=LDLTduicheng(A,b)
% A为对称正定矩阵
n=length(A);
L=eye(n);
for k=1:n
    if k==1
        D(k,k)=A(k,k);
        L(k+1:n,k)=A(k+1:n,k)/D(k,k);
    else
        D(k,k)=A(k,k)-(L(k,1:k-1).^2*diag(D(1:k-1,1:k-1)));
        if k~=n
            L(k+1:n,k)=(A(k+1:n,k)-L(k+1:n,1:k-1)*(L(k,1:k-1)*diag(D(1:k-1,1:k-1))))/D(k,k);
        end
    end
end
LT=L';
for k=1:n
    if k==1
        y(k,1)=b(k,1)/L(k,k);
    else
        y(k,1)=(b(k,1)-L(k,1:k-1)*y(1:k-1,1))/L(k,k);
    end
end
LTT=D*LT;
for k=n:-1:1
    if k==n
        x(k,1)=y(k,1)/LTT(k,k);
    else
        x(k,1)=(y(k,1)-LTT(k,k+1:n)*x(k+1:n,1))/LTT(k,k);
    end
end
end