LAPACK 求矩阵的逆

#LAPACK 求矩阵的逆| 来源: 网络整理| 查看: 265

整理后转载！！！！

=======================================================================

I would like to be able to compute the inverse of a general NxN matrix in C/C++ using lapack.

My understanding is that the way to do an inversion in lapack is by using the dgetri function, however, I can't figure out what all of its arguments are supposed to be.

Here is the code I have:

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; }

How would you complete it to obtain the inverse of the 3x3 matrix M using dgetri_?

4 Answers

First, M has to be a two-dimensional array, like double M[3][3]. Your array is, mathematically speaking, a 1x9 vector, which is not invertible.

N is a pointer to an int for the order of the matrix - in this case, N=3.

A is a pointer to the LU factorization of the matrix, which you can get by running the LAPACK routine dgetrf.

LDA is an integer for the "leading element" of the matrix, which lets you pick out a subset of a bigger matrix if you want to just invert a little piece. If you want to invert the whole matrix, LDA should just be equal to N.

IPIV is the pivot indices of the matrix, in other words, it's a list of instructions of what rows to swap in order to invert the matrix. IPIV should be generated by the LAPACK routine dgetrf.

LWORK and WORK are the "workspaces" used by LAPACK. If you are inverting the whole matrix, LWORK should be an int equal to N^2, and WORK should be a double array with LWORK elements.

INFO is just a status variable to tell you whether the operation completed successfully. Since not all matrices are invertible, I would recommend that you send this to some sort of error-checking system. INFO=0 for successful operation, INFO=-i if the i'th argument had an incorrect input value, and INFO > 0 if the matrix is not invertible.

So, for your code, I would do something like this:

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 }

Regarding 1x9 or 3x3. There really isn't any difference in terms of memory layout. In fact the BLAS/LAPACK routines don't take 2d arrays. They take 1d arrays and make assumptions about how you've laid it out. Beware though that BLAS and LAPACK will assume FORTRAN ordering (rows change fastest) rather than C ordering. – MRocklin Oct 18 '12 at 18:28

Here is the working code for computing the inverse of a matrix using lapack in C/C++:

#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; }

Here is a working version of the above using OpenBlas interface to LAPACKE. Link with openblas library (LAPACKE is already contained)

#include #include "cblas.h" #include "lapacke.h" // inplace inverse n x n matrix A. // matrix A is Column Major (i.e. 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

【本文地址】

LAPACK 求矩阵的逆

LAPACK 求矩阵的逆

今日新闻

推荐新闻