我在C中编写一个需要矩阵和向量乘法的算法。 我有一个矩阵Q (W x W),它是通过将矢量J (1 x W)的转置与其自身相乘并添加使用标量a缩放的单位矩阵I而 创建的 。
Q = [(J ^ T)* J + aI]。
然后我必须将Q的倒数乘以向量G以得到向量M.
M =(Q ^( – 1))* G.
我正在使用cblas和clapack来开发我的算法。 当使用随机数(类型float)填充矩阵Q并使用例程sgetrf_和sgetri_反转时 ,计算的逆是正确的 。
但是当矩阵Q是对称的时,你乘以(J ^ T)x J就是这种情况,计算出的逆是错误的! 。
我知道从C调用lapack例程时数组的行主要(在C中)和列主要(在FORTRAN中)格式,但对于对称矩阵,这应该不是问题,因为A ^ T = A.
我已经在下面附加了我的C函数代码用于矩阵求逆。
我相信有更好的方法来解决这个问题。 谁能帮我这个?
使用cblas的解决方案会很棒……
谢谢。
void InverseMatrix_R(float *Matrix, int W) { int LDA = W; int IPIV[W]; int ERR_INFO; int LWORK = W * W; float Workspace[LWORK]; // - Compute the LU factorization of a M by N matrix A sgetrf_(&W, &W, Matrix, &LDA, IPIV, &ERR_INFO); // - Generate inverse of the matrix given its LU decompsotion sgetri_(&W, Matrix, &LDA, IPIV, Workspace, &LWORK, &ERR_INFO); // - Display the Inverted matrix PrintMatrix(Matrix, W, W); } void PrintMatrix(float* Matrix, int row, int colm) { int i,k; for (i =0; i < row; i++) { for (k = 0; k < colm; k++) { printf("%g, ",Matrix[i*colm + k]); } printf("n"); } }
我不知道BLAS或LAPACK,所以我不知道是什么原因导致这种行为。
但是,对于给定forms的矩阵,计算逆是非常容易的。 对此的重要事实是
(J^T*J)^2 = (J^T*J)*(J^T*J) = J^T*(J*J^T)*J = * (J^T*J)
其中表示内部产品(如果组件是真实的 – 复杂组件的规范双线性forms,但是你可能不会考虑转置而是共轭转置,你会回到内部产品)。
要概括,
(J^T*J)^n = ()^(n-1) * (J^T*J), for n >= 1.
让我们用S
表示对称方阵(J^T*J)
,用q
表示标量
。 那么,对于一般来说,绝对值足够大a != 0
( |a| > q
):
(a*I + S)^(-1) = 1/a * (I + a^(-1)*S)^(-1) = 1/a * (I + ∑ (-1)^k * a^(-k) * S^k) k>0 = 1/a * (I + (∑ (-1)^k * a^(-k) * q^(k-1)) * S) k>0 = 1/a * (I - 1/(a+q)*S) = 1/a*I - 1/(a*(a+q))*S
除了a = 0
和a = -q
,该公式(通过解析性)保持所有a
,这可以通过计算来validation
(a*I + S) * (1/a*I - 1/(a*(a+q))*S) = I + 1/a*S - 1/(a+q)*S - 1/(a*(a+q))*S^2 = I + 1/a*S - 1/(a+q)*S - q/(a*(a+q))*S = I + ((a+q) - a - q)/(a*(a+q))*S = I
使用S^2 = q*S
该计算也比首先找到LU分解更简单,更有效。
您可能想尝试Armadillo ,这是一个易于使用的LAPACK C ++包装器。 它提供了几个与逆相关的函数:
3×3矩阵求逆的示例,请访问sgetri.f了解更多信息
//__CLPK_integer is typedef of int //__CLPK_real is typedef of float __CLPK_integer ipiv[3]; { //Compute LU lower upper factorization of matrix __CLPK_integer m=3; __CLPK_integer n=3; __CLPK_real *a=(float *)this->m1; __CLPK_integer lda=3; __CLPK_integer info; sgetrf_(&m, &n, a, &lda, ipiv, &info); } { //compute inverse of a matrix __CLPK_integer n=3; __CLPK_real *a=(float *)this->m1; __CLPK_integer lda=3; __CLPK_real work[3]; __CLPK_integer lwork=3; __CLPK_integer info; sgetri_(&n, a, &lda, ipiv, work, &lwork, &info); }
以上就是c/c++开发分享用CBLAS / LAPACK进行C中的对称矩阵求逆相关内容,想了解更多C/C++开发(异常处理)及C/C++游戏开发关注计算机技术网(www.ctvol.com)!)。
本文来自网络收集,不代表计算机技术网立场,如涉及侵权请联系管理员删除。
ctvol管理联系方式QQ:251552304
本文章地址:https://www.ctvol.com/c-cdevelopment/562340.html