Computing the inverse of a matrix using lapack in C

Question

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_?

Answer 1

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

#include <cstdio>  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];     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; } 

Answer 2

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      }