LDL_equ.m

% 函数功能：实现Ax=b的改进平方根法，其中A为n阶对称正定矩阵，将其分解为LDL'的形式
% 输入：矩阵A，b
% 输出：得到的解矩阵x
% 使用范围：改进平方根法比平方根法（GG'）少了开方运算，比高斯消去法和Doolittle分解法（LU）少了一半的乘除法
% 是求解系数矩阵为对称正定矩阵的线性方程组最有效的方法之一

function x = LDL_equ(A, b)

[~, n] = size(A);
%% 方阵的LDL分解(改进平方根分解)
y = b; % 把b放入增广矩阵，分解后得到y
A(2: n, 1) = A(2: n, 1) / A(1, 1); % 先求第一列的L矩阵，第一行不变不用求
A(2, 2: n) = A(2, 2: n) - A(2, 1) * A(1, 2: n); % 求第二行的上三角矩阵的值
y(2, :) = y(2, :) - A(2, 1) * y(1, :); % 用同样的方法对第二行的b分解
for k = 3: n % 第k行循环
    A(k, 2: k - 1) = (A(2: k - 1, k) ./ diag(A(2: k - 1, 2: k - 1)))'; % 此时计算L矩阵
    A(k, k: n) = A(k, k: n) - A(k, 1: k - 1) * A(1: k - 1, k: n); % 计算U矩阵，此时k<=j
    y(k, :) = y(k, :) - A(k, 1: k - 1) * y(1: k - 1, :); % 用同样的方法对b分解
end

%% 以下转化为求LDL'x=b方程组，拆分求解
% 解对角线矩阵方程组Dz=y
z = y ./ diag(A);
% 解上三角矩阵方程组L'x=z
A(eye(n) == 1) = 1; % 把A的对角线元素变为1，冒充L矩阵使用
x = UTri_equ(A', z);

end