LAPACK 求矩阵的逆

整理后转载！！！！

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

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 1×9 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 1×9 or 3×3. 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 <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+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 <stdio.h>
#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<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');
}

Example:

% 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 

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 the errorHandler. I changed the 9 in that matrix to a 19, getting an errorHandler of 0 signalling success, and compared the result to that from Mathematica. 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:

#include <stdio.h>
#include <stddef.h>
#include <lapacke.h>

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

I built and ran it as follows on a Mac:

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

I did an nm on the LAPACKE library and found the following:

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

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.

#include <stdio.h>
#include <stddef.h>

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

built and run like this:

$ 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

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.

今天的文章LAPACK 求矩阵的逆分享到此就结束了，感谢您的阅读。