我有 n 个任意 px 1 向量 x_i 和 pxk 矩阵 A_i 和 npxp 半正定矩阵 S_i，其中一些（通常大多数）*S_i* 是相同的（例如只有两个不同的 S 矩阵，一个正定矩阵适用到 i=1,..., n-1 和半定 S for i=n)。我想对所有 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'_i=S_i 其中 D_i 是对角矩阵。）
n 的范围可以从几到几千甚至更多，p 通常小于 10，k 通常小于 100。S可以包含带零的行和列，它最初形成为 _BB' 其中 B 是下三对角矩阵具有非负对角元素。但这些 B 不可用。
我想知道在速度和准确性方面实现这一目标的最佳方法是什么，并且更加重视准确性？
目前我正在为不同的 S 使用我自己编写的 LDL 分解函数，反转 L 并计算转换，因为我一直在处理具有少量不同 S 和大 n 的情况。如果我理解正确，那么只求解线性方程组而不在精度方面显式反转矩阵会更明智，但就速度而言，它似乎是更好的选择？
我正在使用 Fortran，我宁愿使用一些 LAPACK 函数进行分解，而不是我自己的 LDL 分解，我不能保证它是稳定的等，但我不确定哪个是好方法。