拉帕克C++实矩阵反演

Question

我试图在 LAPACKE C++反转一个实矩阵。我对复杂矩阵具有相同的函数，并且它有效。但真实案例给出了错误的答案。这是我的函数：

void inv(std::vector<std::vector<double>> &ans, std::vector<std::vector<double>> MAT){
int N = MAT.size();
int *IPIV = new int[N];
double * arr = new double[N*N];
for (int i = 0; i<N; i++){
for (int j = 0; j<N; j++){
int idx = i*N + j;
arr[idx] = MAT[i][j];
}
}
LAPACKE_dgetrf(LAPACK_ROW_MAJOR, N, N, arr, N, IPIV);
LAPACKE_dgetri(LAPACK_ROW_MAJOR, N, arr, N, IPIV);
for (int i = 0; i<N; i++){
for (int j = 0; j<N; j++){
int idx = i*N + j;
ans[i][j] = arr[idx];
}
}
delete[] IPIV;
delete[] arr;
}

我尝试反转 24 x 24 的双精度矩阵。虽然程序似乎几乎就在那里，但逆向还没有完全存在，它与 python linalg 逆向给我的有很大不同(python 就在这里，因为我将矩阵与逆相乘，结果非常接近 indentity(。在 LAPACKE 输出中，我将矩阵乘以其逆矩阵，我得到对角线为 1，但非对角线的值高达 0.17，与 0 相比，这是巨大的。有没有办法使LAPACKE程序产生更好的结果？谢谢！

Answer 1

对于具有大行列式的矩阵，您可以重新缩放输入，计算逆运算，然后重新缩放输出。这是一个非常简单的 Python 示例，您的比例因子应该是 1/25 左右，以获得 (1/25(²⁴=2.8e-34 的总乘数，从而使输入矩阵行列式约为 1000。

import numpy as np
scale = 0.5
i = np.array([[1,2],[3,4]])
print(i)
print(np.linalg.det(i))
print("-----------------------------------")
x = np.multiply(i, [scale]) # rescale matrix
print(x)
print(np.linalg.det(x))     # determinant should be less
print("-----------------------------------")
y = np.linalg.inv(x)
print(y)
print(np.linalg.det(y))
print("-----------------------------------")
o = np.multiply(y, [scale]) # rescale matrix
print(o)
print(np.linalg.det(o))
print(np.dot(i, o))

您可以使用scale并验证任何值是代码返回正确的倒置输入矩阵。

所以你的代码将是

double scale = 1.0/25.0;
for (int i = 0; i<N; i++){
for (int j = 0; j<N; j++){
int idx = i*N + j;
arr[idx] = scale*MAT[i][j];
}
}
....
for (int i = 0; i<N; i++){
for (int j = 0; j<N; j++){
int idx = i*N + j;
ans[i][j] = scale * arr[idx];
}
}