使用C中的lapack计算矩阵的逆矩阵

我希望能够使用lapack计算C / C ++中一般NxN矩阵的逆。

我的理解是在lapack中进行反转的方法是使用dgetri函数,但是,我无法弄清楚它的所有参数应该是什么。

这是我的代码:

 void dgetri_(int* N, double* A, int* lda, int* IPIV, double* WORK, int* lwork, int* INFO); int main(){ double M [9] = { 1,2,3, 4,5,6, 7,8,9 }; return 0; } 

你将如何完成它以使用dgetri_获得3x3矩阵M的逆?

首先,M必须是一个二维数组,如double M[3][3] 。 从数学上讲,你的数组是1×9向量,它是不可逆的。

  • N是指向矩阵顺序的int的指针 – 在这种情况下,N = 3。

  • A是指向矩阵的LU分解的指针,您可以通过运行LAPACK例程dgetrf获得该dgetrf

  • LDA是矩阵“前导元素”的整数,如果您只想反转一小块,它可以让您选择更大矩阵的子集。 如果要反转整个矩阵,LDA应该等于N.

  • IPIV是矩阵的枢轴索引,换句话说,它是要交换哪些行以反转矩阵的指令列表。 IPIV应该由LAPACK例程dgetrf生成。

  • LWORK和WORK是LAPACK使用的“工作空间”。 如果要反转整个矩阵,LWORK应该是一个等于N ^ 2的int,而WORK应该是一个带有LWORK元素的双数组。

  • INFO只是一个状态变量,用于告诉您操作是否成功完成。 由于并非所有矩阵都是可逆的,我建议您将其发送到某种错误检查系统。 对于成功操作,INFO = 0,如果第i个参数具有不正确的输入值,则INFO = -i,如果矩阵不可逆,则INFO> 0。

所以,对于你的代码,我会做这样的事情:

 int main(){ double M[3][3] = { {1 , 2 , 3}, {4 , 5 , 6}, {7 , 8 , 9}} double pivotArray[3]; //since our matrix has three rows int errorHandler; double lapackWorkspace[9]; // dgetrf(M,N,A,LDA,IPIV,INFO) means invert LDA columns of an M by N matrix // called A, sending the pivot indices to IPIV, and spitting error // information to INFO. // also don't forget (like I did) that when you pass a two-dimensional array // to a function you need to specify the number of "rows" dgetrf_(3,3,M[3][],3,pivotArray[3],&errorHandler); //some sort of error check dgetri_(3,M[3][],3,pivotArray[3],9,lapackWorkspace,&errorHandler); //another error check } 

这是使用C / C ++中的lapack计算矩阵的逆矩阵的工作代码:

 #include  extern "C" { // LU decomoposition of a general matrix void dgetrf_(int* M, int *N, double* A, int* lda, int* IPIV, int* INFO); // generate inverse of a matrix given its LU decomposition void dgetri_(int* N, double* A, int* lda, int* IPIV, double* WORK, int* lwork, int* INFO); } void inverse(double* A, int N) { int *IPIV = new int[N+1]; int LWORK = N*N; double *WORK = new double[LWORK]; int INFO; dgetrf_(&N,&N,A,&N,IPIV,&INFO); dgetri_(&N,A,&N,IPIV,WORK,&LWORK,&INFO); delete IPIV; delete WORK; } int main(){ double A [2*2] = { 1,2, 3,4 }; inverse(A, 2); printf("%f %f\n", A[0], A[1]); printf("%f %f\n", A[2], A[3]); return 0; } 

以下是使用OpenBlas接口到LAPACKE的上述工作版本。 与openblas库链接(已包含LAPACKE)

 #include  #include "cblas.h" #include "lapacke.h" // inplace inverse nxn matrix A. // matrix A is Column Major (ie firts line, second line ... *not* C[][] order) // returns: // ret = 0 on success // ret < 0 illegal argument value // ret > 0 singular matrix lapack_int matInv(double *A, unsigned n) { int ipiv[n+1]; lapack_int ret; ret = LAPACKE_dgetrf(LAPACK_COL_MAJOR, n, n, A, n, ipiv); if (ret !=0) return ret; ret = LAPACKE_dgetri(LAPACK_COL_MAJOR, n, A, n, ipiv); return ret; } int main() { double A[] = { 0.378589, 0.971711, 0.016087, 0.037668, 0.312398, 0.756377, 0.345708, 0.922947, 0.846671, 0.856103, 0.732510, 0.108942, 0.476969, 0.398254, 0.507045, 0.162608, 0.227770, 0.533074, 0.807075, 0.180335, 0.517006, 0.315992, 0.914848, 0.460825, 0.731980 }; for (int i=0; i<25; i++) { if ((i%5) == 0) putchar('\n'); printf("%+12.8f ",A[i]); } putchar('\n'); matInv(A,5); for (int i=0; i<25; i++) { if ((i%5) == 0) putchar('\n'); printf("%+12.8f ",A[i]); } putchar('\n'); } 

例:

 % g++ -I [OpenBlas path]/include/ example.cpp [OpenBlas path]/lib/libopenblas.a % a.out +0.37858900 +0.97171100 +0.01608700 +0.03766800 +0.31239800 +0.75637700 +0.34570800 +0.92294700 +0.84667100 +0.85610300 +0.73251000 +0.10894200 +0.47696900 +0.39825400 +0.50704500 +0.16260800 +0.22777000 +0.53307400 +0.80707500 +0.18033500 +0.51700600 +0.31599200 +0.91484800 +0.46082500 +0.73198000 +0.24335255 -2.67946180 +3.57538817 +0.83711880 +0.34704217 +1.02790497 -1.05086895 -0.07468137 +0.71041070 +0.66708313 -0.21087237 -4.47765165 +1.73958308 +1.73999641 +3.69324020 -0.14100897 +2.34977565 -0.93725915 +0.47383541 -2.15554470 -0.26329660 +6.46315378 -4.07721533 -3.37094863 -2.42580445 

这是上面Spencer Nelson的一个工作版本。 关于它的一个谜团是输入矩阵是行主要顺序,即使它似乎调用了基础fortran例程dgetri 。 我被引导相信所有潜在的fortran例程都需要列主要顺序,但我不是LAPACK的专家,事实上,我正在使用这个例子来帮助我学习它。 但是,除了一个谜团:

示例中的输入矩阵是单数。 LAPACK尝试通过在errorHandler返回3来告诉您。 我将该矩阵中的9更改为19 ,得到0信号成功的errorHandler ,并将结果与Mathematica的结果进行比较。 比较也是成功的,并确认示例中的矩阵应按行显示主要顺序。

这是工作代码:

 #include  #include  #include  int main() { int N = 3; int NN = 9; double M[3][3] = { {1 , 2 , 3}, {4 , 5 , 6}, {7 , 8 , 9} }; int pivotArray[3]; //since our matrix has three rows int errorHandler; double lapackWorkspace[9]; // dgetrf(M,N,A,LDA,IPIV,INFO) means invert LDA columns of an M by N matrix // called A, sending the pivot indices to IPIV, and spitting error information // to INFO. also don't forget (like I did) that when you pass a two-dimensional // array to a function you need to specify the number of "rows" dgetrf_(&N, &N, M[0], &N, pivotArray, &errorHandler); printf ("dgetrf eh, %d, should be zero\n", errorHandler); dgetri_(&N, M[0], &N, pivotArray, lapackWorkspace, &NN, &errorHandler); printf ("dgetri eh, %d, should be zero\n", errorHandler); for (size_t row = 0; row < N; ++row) { for (size_t col = 0; col < N; ++col) { printf ("%g", M[row][col]); if (N-1 != col) { printf (", "); } } if (N-1 != row) { printf ("\n"); } } return 0; } 

我在Mac上构建并运行如下:

 gcc main.c -llapacke -llapack ./a.out 

我在LAPACKE库上做了一个nm ,发现了以下内容:

 liblapacke.a(lapacke_dgetri.o): U _LAPACKE_dge_nancheck 0000000000000000 T _LAPACKE_dgetri U _LAPACKE_dgetri_work U _LAPACKE_xerbla U _free U _malloc liblapacke.a(lapacke_dgetri_work.o): U _LAPACKE_dge_trans 0000000000000000 T _LAPACKE_dgetri_work U _LAPACKE_xerbla U _dgetri_ U _free U _malloc 

而且看起来有一个LAPACKE [sic]包装器可能会让我们不得不为了fortran的方便而到处寻找地址,但我可能不会去尝试它,因为我有前进的方向。

编辑

这是一个绕过LAPACKE [原文如此]的工作版本,直接使用LAPACK fortran例程。 我不明白为什么行主要输入会产生正确的结果,但我在Mathematica中再次证实了这一点。

 #include  #include  int main() { int N = 3; int NN = 9; double M[3][3] = { {1 , 2 , 3}, {4 , 5 , 6}, {7 , 8 , 19} }; int pivotArray[3]; //since our matrix has three rows int errorHandler; double lapackWorkspace[9]; /* from http://www.netlib.no/netlib/lapack/double/dgetrf.f SUBROUTINE DGETRF( M, N, A, LDA, IPIV, INFO ) * * -- LAPACK routine (version 3.1) -- * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * November 2006 * * .. Scalar Arguments .. INTEGER INFO, LDA, M, N * .. * .. Array Arguments .. INTEGER IPIV( * ) DOUBLE PRECISION A( LDA, * ) */ extern void dgetrf_ (int * m, int * n, double * A, int * LDA, int * IPIV, int * INFO); /* from http://www.netlib.no/netlib/lapack/double/dgetri.f SUBROUTINE DGETRI( N, A, LDA, IPIV, WORK, LWORK, INFO ) * * -- LAPACK routine (version 3.1) -- * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * November 2006 * * .. Scalar Arguments .. INTEGER INFO, LDA, LWORK, N * .. * .. Array Arguments .. INTEGER IPIV( * ) DOUBLE PRECISION A( LDA, * ), WORK( * ) */ extern void dgetri_ (int * n, double * A, int * LDA, int * IPIV, double * WORK, int * LWORK, int * INFO); // dgetrf(M,N,A,LDA,IPIV,INFO) means invert LDA columns of an M by N matrix // called A, sending the pivot indices to IPIV, and spitting error information // to INFO. also don't forget (like I did) that when you pass a two-dimensional // array to a function you need to specify the number of "rows" dgetrf_(&N, &N, M[0], &N, pivotArray, &errorHandler); printf ("dgetrf eh, %d, should be zero\n", errorHandler); dgetri_(&N, M[0], &N, pivotArray, lapackWorkspace, &NN, &errorHandler); printf ("dgetri eh, %d, should be zero\n", errorHandler); for (size_t row = 0; row < N; ++row) { for (size_t col = 0; col < N; ++col) { printf ("%g", M[row][col]); if (N-1 != col) { printf (", "); } } if (N-1 != row) { printf ("\n"); } } return 0; } 

建立并运行如下:

 $ gcc foo.c -llapack $ ./a.out dgetrf eh, 0, should be zero dgetri eh, 0, should be zero -1.56667, 0.466667, 0.1 1.13333, 0.0666667, -0.2 0.1, -0.2, 0.1 

编辑

这个谜团似乎不再是一个谜。 我认为计算是按照列主要顺序进行的,因为它们必须如此,但我输入和打印矩阵就像它们按行主顺序一样。 我有两个相互抵消的错误,所以即使它们是列列的,它们看起来仍然是行列式的。