function q=ldlt(A,NA,b)
% LDLT分解
% A 一维存储的系数矩阵
% NA 对角线元素在一维矩阵中的位置
% b 右端项
% A=[a d g
% d e h
% g h i];
% 则 A=[a d e g h i]; NA=[1 3 6]; 是以下三角的行存储的
n = size(NA,2);
for k=2:n
u = zeros(k-1,1);
Lk = k - (NA(k)-NA(k-1)) + 1;
u(Lk:(k-1)) = A((NA(k-1)+1):(NA(k)-1));
for i=(Lk+1):(k-1)
Li = i - (NA(i)-NA(i-1)) + 1;
u(i) = u(i) - A((NA(i-1)+1):(NA(i)-1))*u(Li:(i-1));
end
for i=Lk:(k-1)
A(NA(k)-k+i) = u(i)/A(NA(i));
end
A(NA(k)) = A(NA(k)) - A((NA(k-1)+1):(NA(k)-1))*u(Lk:(k-1));
end
D=A(NA);
L=A;L(NA)=ones(1,n);
%求解方程L*y=b
y=b; s=2;
for k=2:n
Lk = k - (NA(k)-NA(k-1)) + 1;
for j=Lk:k-1
y(k)=y(k)-L(s)*y(j);
s=s+1;
end
s=s+1;
end
%求解D*z=y,D是对角阵
z=zeros(n,1);
for k=1:n
z(k)=y(k)/D(k);
end
%求解方程L'*q=z
q=z; s=NA(n)-1;
for k=n-1:-1:1
for j=k:-1:k-(NA(k+1)-NA(k)-1)+1
q(j)=q(j)-q(k+1)*L(s);
s=s-1;
end
s=s-1;
end
end