function [L,D,M]=ldm(A)
% 此函數(shù)用來求解矩陣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