我有n个任意的px1向量x_I、pxk矩阵A_I和n个pxp半正定矩阵S_I，其中一些(通常大多数)*S_I*是相同的(例如只有两个不同的S矩阵，一个正定适用于I=1，…，n-1和I=n的半定S)。我想对所有x_I和A_I进行以下线性变换：
x_i*=inv(L_i)x_i和
A_i*=inv(L_i)A_i，其中
L_i是下三角矩阵，使得L_i L’_i＝S_i(或者甚至更好的L_i D_i L’_id＝S_
n可以在几到几千甚至更多的范围内，p通常小于10并且k通常小于100S可以包含具有零的行和列，并且它最初形成为_BB'，其中B是具有非负对角元素的下三对角矩阵。不过这些B没有。
我想知道，在速度和准确性方面，在更注重准确性的情况下，实现这一目标的最佳方式是什么？
目前，我正在对不同的S使用我自己编写的LDL分解函数，反转L并计算变换，因为我一直在处理具有少量不同S和大量n的情况。如果我理解正确，那么就精度而言，只求解线性方程组而不显式反转矩阵会更明智，但就速度而言，这似乎是更好的选择？
我使用的是Fortran，我宁愿使用一些LAPACK函数进行因子分解，而不是我自己的LDL分解，我不能保证它是稳定的等等，但我不确定哪种方法是好的。
我想你需要的是Cholesky变换和使用Cholesky转换形式的解算器。使用LAPACK绝对是个好主意。寻找?potrf(Cholesky分解)和?potrs(求解Cholesky形式的线性方程组)函数。请参阅《LAPACK用户指南》的线性方程部分。