我试图在C++ LAPACKE中反转一个真实矩阵。对于复杂矩阵,我具有相同的功能,并且可以正常工作。但是实际情况给出了错误的答案。这是我的功能:
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的 double 矩阵反转。虽然程序似乎快要存在了,但是逆函数还不完全存在,它与python linalg逆函数给我的有很大不同(python就在这里,因为我将矩阵乘以逆函数,结果非常接近于恒等式)。在LAPACKE输出中,我将矩阵乘以它的逆,得到对角线为1,但非对角线的值最高为0.17,与0相比,这是巨大的。有没有办法使LAPACKE程序提供更好的结果?谢谢!
最佳答案
对于行列式较大的矩阵,可以重新缩放输入,计算逆,然后重新缩放输出。这是一个非常简单的Python示例,您的比例因子应为1/25左右,以便获得(1/25)24 = 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];
}
}
关于c++ - LAPACKE C++实矩阵求逆,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/59901747/