c++ - LAPACKE C++ 实矩阵求逆
问题描述
我试图在 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 的双精度矩阵。虽然程序似乎几乎就在那里,但逆还没有完全出现,它与 python linalg inverse 给我的有很大不同(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];
}
}
推荐阅读
- html - *ngIf else 一段时间后转到 else 模板
- javascript - 在 keyup 上设置范围滑块进度
- java - 序列化接收到的具有“LocalDateTime”数据类型的 JSON 对象
- python - 序列化数据类列表:str 总是调用 repr
- php - 我们如何通过限制相同列数MySql从表中获取数据
- fortran - Fortran 满足某些条件时如何跳过行?
- python - 基于单个字典中键的值的平均值和标准差
- android - 如果使用片段管理器和导航控制器,片段不可见
- java - 向/从此 Map 添加/删除时,如何扩展 java 的 HashMap 以添加钩子操作?
- python - 将 sigmoid 结果解释为神经网络中的概率