ldm.m

MATLAB科学计算与工程分析源代码源程序4

function [L,D,M]=ldm(A)
% 此函数用来求解矩阵A的LDM'分解
% 其中L,M均为单位下三角矩阵，D为对角矩阵
% 编写日期：2007-5-13

[m,n]=size(A);
if m~=n
    error('输入矩阵不是方阵，请正确输入矩阵！');
    return;
end
D(1,1)=A(1,1);
for i=1:n
    L(i,i)=1;
    M(i,i)=1;
end
L(2:n,1)=A(2:n,1)/D(1,1);
M(2:n,1)=A(1,2:n)'/D(1,1);

for j=2:n
    v(1)=A(1,j);
    for i=2:j
        v(i)=A(i,j)-L(i,1:i-1)*v(1:i-1)';
    end
    for i=1:j-1
        M(j,i)=v(i)/D(i,i);
    end
    D(j,j)=v(j);
    L(j+1:n,j)=(A(j+1:n,j)-L(j+1:n,1:j-1)*v(1:j-1)')/v(j);
end