LDLT 分解法求解对称线性方程组
朱松盛 041002045 南京师范大学
一、 LDLT 分解法原理
利用矩阵 A 的
分解来解Ax = b，的方法称为
分解法。当求解方程组的系数矩阵是对称矩阵时，则用
分解法可以简化程序设计并减少计算量。
当矩阵 A 的各阶顺序主子式不为零时，A 有唯一的 Doolittle 分解 A = LU，其中
，
。
此时，当然有
，所以矩阵U 的对角线元素
，（
），将矩阵U 每行依次提出
，则有
，其中
，
因而
，显然这种A 的分解也是唯一的，当A 为对称矩阵时，由
，得
，由分解的唯一性，有
，即
。
由此可得，若对称矩阵A 的各阶顺序主子式不为零时，A 可唯一分解为
，其中
，
，
为L 的转置矩阵。
当A 有
分解时，利用矩阵运算法则及相等原理易得计算
及
的
公式
小学单位换算公式大全免费下载公式下载行测公式大全下载excel公式下载逻辑回归公式下载
为
，
为减少乘法次数，引入辅助量
，则上面公式可写成
，
（1）
当设计程序时，为了减少存储空间，又由于矩阵A 在经过计算后就不需要再使用了，因此可以采用原位存储的方法。也就是可以将
存于矩阵A 的对角线元素即
中，
存于矩阵A 的上三角部分即
中，而
则可以存于矩阵A 的下三角部分即
中。
当需要求解对称线性方程组
时，由
分解，有
，这可以分解成三个简单的方程组
，
和
由此易得求解公式为
（2）
同样，向量b 在经过计算后也不再需要，因此也可以采用相同的方法处理。可以在计算
时将
存于
中，然后计算
时也将
存于
中，最后计算
时还是存于
中，由
输出计算结果即可。
公式（1）和（2）就是
分解的求解公式。
用
分解法求解对称方程组的计算量约为
，这比用Gauss 消元法求解的计算量要少。对一般的对称方程组，会出现因某个
而不能进行
分解的情况，不过当矩阵A 的各阶顺序主子式都不为零时，则
分解总能进行到底，从而可用
分解法来求解。
二、 LDLT 分解法程序主要部分（C 语言）
程序中所有矩阵的下标都是从0～n－1，与公式中的1～n不同，这只是习惯上的不同，其余都是一致的。
程序主要部分有两个模块LDLTDCMP 和LDLTBKSB。第一个是
分解，第二个是解方程组，A 经LDLTDCMP计算后输出给LDLTBKSB 使用。两个模块分开，可用于解不同b，相同A 的多个对称线性方程组。
void LDLTDCMP (unsigned int n, double * *a)
{
int k;
int m;
int i;
for (k = 0; k < n; k++)
{
/* Calculate d[k], and saved in a[k][k] */
for (m = 0; m < k; m++)
a[k][k] = a[k][k] - a[m][k] * a[k][m];
if (a[k][k] == 0)
{
printf ("\n\nERROR: LDL\' decompose failed !!\n");
printf (" Main submatrix of every rank of A cannot be zero.\n\n");
return;
}
for (i = k + 1; i < n; i++)
{
/* temporary varible u[i][k] is saved in a[k][i]*/
for (m = 0; m < k; m++)
a[k][i] = a[k][i] - a[m][i] * a[k][m];
/* Calculate l[i][k], and saved in a[i][k]*/
a[i][k] = a[k][i] / a[k][k];
}
}
}
void LDLTBKSB (unsigned int n, double * *a, double *b)
{
int i;
int k;
for (i = 0; i < n; i++)
{
/* Calculate y[i], and saved in b[i] */
for (k = 0; k < i; k++)
b[i] = b[i] - a[i][k] * b[k];
}
for (i = n - 1; i >= 0; i--)
{
/* Calculate z[i], and saved in b[i] */
b[i] = b[i] / a[i][i];
/* Calculate x[i], and saved in b[i] */
for (k = i + 1; k < n; k++)
b[i] = b[i] - a[k][i] * b[k];
}
}
最终的测试程序还加入了参数的输入输出（都是文件读写），和一些信息提示。
三、 LDLT 分解法程序的检验
用
，
检验程序，得结果为
，与实际计算结果一致。
四、 LDLT 分解法程序的使用说明
程序目录结构：
|-c_source C 源代码文件，能编译测试相应的算法函数
|-data 输入数据文本文件
|-func C 的函数代码（可被引用于其他程序中）
|-output 输出文件，包括执行文件、目标文件和输出结果的文本文件
|-spec 说明文挡
1． 输入数据文件
格式
pdf格式笔记格式下载页码格式下载公文格式下载简报格式下载
输入数据文件名为LDLTPARA.DAT（在./data/下），为文本文件，文件内容格式如下：
# 以'#'开头的行将被忽略掉, 是注释行. 不要有空行.
# 第一个有用行必须以'n ='开头, 接着是系数矩阵的维数.
n = 3
# 接着输入增广矩阵（也就是系数矩阵A 和b. A 在左边, b 在右边）.
# 每一行应该有 n+1 个数, n 个是矩阵A 的, 一个b 的.
# A b
1 0.5 0.5 1
0.5 1 0.5 -2
0.5 0.5 1 3
# 按维数输完A 和b 后, 接下来的所有行都会被忽略.
2． 输出数据文件
输出数据文件名为LDLT_KEY.DAT（在./output/下），为文本文件，文件内容格式如下：
Original parameters:
n = 3
1 0.5 0.5 1
0.5 1 0.5 -2
0.5 0.5 1 3
After calculated:
1 0.5 0.5 1
0.5 0.75 0.25 -5
0.5 0.333 0.667 5
Answers x[]:
1 -5 5
其中，Original parameters 输出的是方程组原来的参数（维数和增广矩阵），After calculated 输出的是经过计算后矩阵A 和b 的内容，最后输出计算结果即方程组的解。
浙公网安备 33021202002483