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_?
Here is the working code for computing the inverse of a matrix using lapack in C/C++:
Here is a working version of the above using OpenBlas interface to LAPACKE. Link with openblas library (LAPACKE is already contained)
Example:
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:
Here is a working version of Spencer Nelson's example above. One mystery about it is that the input matrix is in row-major order, even though it appears to call the underlying fortran routine
dgetri
. I am led to believe that all the underlying fortran routines require column-major order, but I am no expert on LAPACK, in fact, I'm using this example to help me learn it. But, that one mystery aside:The input matrix in the example is singular. LAPACK tries to tell you that by returning a
3
in theerrorHandler
. I changed the9
in that matrix to a19
, getting anerrorHandler
of0
signalling success, and compared the result to that fromMathematica
. The comparison was also successful and confirmed that the matrix in the example should be in row-major order, as presented.Here is the working code:
I built and ran it as follows on a Mac:
I did an
nm
on the LAPACKE library and found the following:and it looks like there is a LAPACKE [sic] wrapper that would presumably relieve us of having to take addresses everywhere for fortran's convenience, but I am probably not going to get around to trying it because I have a way forward.
EDIT
Here is a working version that bypasses LAPACKE [sic], using LAPACK fortran routines directly. I do not understand why a row-major input produces correct results, but I confirmed it again in Mathematica.
built and run like this:
EDIT
The mystery no longer appears to be a mystery. I think the computations are being done in column-major order, as they must, but I am both inputting and printing the matrices as if they were in row-major order. I have two bugs that cancel each other out so things look row-ish even though they're column-ish.