Lapack 很可能没有任何计算行列式的例程。但是，我们可以使用 LU、QR 或 SVD 分解来计算它。我更喜欢使用 LU 分解。现在 lapack 使用一些 dgetrf 子程序将矩阵 A 分解为带有一些 IPIV 数组的 PLU 格式。我不太清楚如何处理这些信息。为了计算行列式，我只需将 U 矩阵的对角元素相乘。但是 PLU 格式中的 L 和 U 是什么以及如何提取它们。我正在用 C 语言编程。
【问题讨论】：
标签： lapack
Lapack 的 dgetrf() 计算一般 M×N 矩阵 A 的 A=P*L*U 分解。假设可逆方阵 A，其行列式可以计算为乘积：
U 是一个上三角矩阵。因此，它的行列式是对角元素的乘积，恰好是输出A 的对角元素。确实，见how the output A is defined：
退出时，因式分解
A = P*L*U中的因子L和U； L的单位对角线元素不存储。
L 是一个下三角矩阵，具有未存储的单位对角元素。因此，其行列式始终为 1。
P 是一个置换矩阵，编码为转置的乘积（即 2 循环或交换）。确实，请参阅dgetri() 以了解它的使用方式。因此，它的行列式是 1 或 -1，这取决于转置的数量是偶数还是奇数。因此，P 的行列式可以计算为：
int j;
double detp=1.;
for( j=0;j<n;j++){
if(j+1!=ipiv[j]){
// j+1 : following feedback of ead : ipiv is from Fortran, hence starts at 1.
// hey ! This is a transpose !
detp=-detp;
}
}
这种方法的复杂性主要由使用部分旋转的高斯消除的成本决定，即 O(2/3n^3)。
您可能会使用 dgetc2() 或 QR 分解转向完全旋转以提高准确性。正如潘等人在Algebraic and Numerical Techniques for the Computation of Matrix Determinants 中发出的信号。等，结合方程 4.8、4.9 和命题 4.1，行列式上的最终误差可能像 ed=(a+eps*a*n^4)^{n-1}*eps*an^5=a^n*(1+eps*n^4)^{n-1}*n^5*eps 这样缩放，其中 eps 是双精度（大约 1e-13），a 是矩阵中所有元素的最大幅度A 和 n 是矩阵的大小。这意味着计算出的行列式对于“大”矩阵不是很重要：见表格，特别是使用 PLU 分解时的相对误差！本文还提供了一种算法来跟踪误差的传播并产生更好的误差估计。
你也可以试试Faddeev–Le Verrier algorithm...
【讨论】：
1...n 它应该可能是 j+1!=ipiv[j] 另请参阅 stackoverflow.com/a/50285530/5769463
L 是一个单位对角矩阵，所以它的行列式总是统一的。对于P，它是1 或-1，具体取决于排列的奇数/偶数，但ipiv 是交换数组而不是排列数组，因此您可以将以下简单的Python 循环转换为C
det = 1
for i in range(len(ipiv)):
det *= u(i,i)
if ipiv[i] != i:
det *= -1
但请注意，如果矩阵的条件在条件方面存在问题，则乘法可能会溢出或下溢。相反，只需使用 SVD。
【讨论】：