LAPACK 求矩阵的逆

original source: http://stackoverflow.com/questions/3519959/computing-the-inverse-of-a-matrix-using-lapack-in-c

Function:

inv(A)

C++ Code (Edited):

#include <cstdio>
#include <cblas.h>

// pre install :
// sudo apt-get install build-essential
// sudo apt-get install liblapack*
// sudo apt-get install libblas*
// compile command : 
// g++ test.cpp -llapacke -llapack
// ./a.out
// May need to adjust stack memory upper limit:
// ulimit -s 1024000 # 1GB

extern "C" {
    // LU decomoposition of a general matrix
    void dgetrf_(int* M, int *N, double* A, int* lda, int* IPIV, int* INFO);
    // void sgetrf_(int* M, int *N, float* 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 sgetri_(int* N, float* A, int* lda, int* IPIV, float* 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(){
    
    int N = 1000;
    double M[N*N];
    
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            M[i*N+j] = 0;
            if (i == j)
                M[i*N+j] = 1;
        }
    }
    inverse(M, N);

    printf("%f %f\n", M[0], M[N+1]); // should be 1 1 
    printf("%f %f\n", M[1], M[N]); // 0 0

    return 0;
}