C algorithm to invert 1024x1024 matrix?

Question

this is my firs question on this website. I am (desperately) trying to invert a big matrix in my program. I would like to use lapack in order to do this, and I found this thread which looks quite promising, but I think it's in C++ language. Could you help me out?

Computing the inverse of a matrix using lapack in C

Thank you.

Update: you are right, the answer is a bit unclear. After incorporating the program i posted into mine, I get the following error:

mymatrixmalloc_2.c:15:18: fatal error: cstdio: non existing File or directory 
#include <cstdio>
              ^
compilation terminated.

I guess the prolem is that I havent's properly installed the llapack library, or that I am including it while compiling.

This is how I installed the library (from terminal, I have Ubuntu):

sudo apt-get install build-essential
sudo apt-get install liblapack*
sudo apt-get install libblas*

And this is how I compile:

davide@John:~$ gcc -Wall -lm -llapack -lblas mymatrixmalloc_2.c -o mymatrixmalloc_2.exe

What am I doing wrong?

Please consider including some code you have written (or used from that link). Or consider describing in more detail the configuration/setup you have to run the code. The less readers have to do (e.g. click and read another link), the better; more likely you'll get a meaningful answer. — FishStix, Jun 23 '16 at 20:29
Possible duplicate of [Computing the inverse of a matrix using lapack in C](http://stackoverflow.com/questions/3519959/computing-the-inverse-of-a-matrix-using-lapack-in-c) — m69's been on strike for years, Jun 24 '16 at 04:13
The answer in https://stackoverflow.com/questions/3519959/computing-the-inverse-of-a-matrix-using-lapack-in-c is commented by @Olaf : `There is no language C/C++. The code you show is C++. But the question is about C.`. Indeed, `#include ` is C++, while `#include ` is C. This is the reason why the compiler complains that cstdio is neither a file nor a directory. — francis, Jun 24 '16 at 13:30

Niklas Rosencrantz · Answer 1 · 2016-06-24T17:56:42.993

You can verify that this C algorithm performs a correct inversion of a small matrix:

gcc main.c -llapacke -llapack
dac@dac-Latitude-E7450 ~/C/gnu> ./a.out 
dgetrf eh, 0, should be zero
dgetri eh, 0, should be zero
0.6, -0.7
-0.2, 0.4⏎

The above number are verified for the 2*2 matrix in this sample program:

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

int main() {
    int N = 2;
    int NN = 4;
    double M[2][2] = {{4, 7},
                      {2, 6}};
    int pivotArray[2];
    int errorHandler;
    double lapackWorkspace[4];
    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;
}

Now all we need is to define a larger matrix 1024*1024 and invert it same way.

#include <stdio.h>
#include <lapacke.h>
int main() {
    int N = 1024;
    int NN = N*N;
    double M[N][N];
    for(int i=0;i<N;i++) {
        for(int j=0;j<N;j++) {
            M[i][j] =  0;
            if(i==j)
                M[i][j] =  1;
        }
    }
    int pivotArray[N];
    int errorHandler;
    double lapackWorkspace[N*N];
    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);
    return 0;
}

To run the above code, I also had to increase my stack size on Linux:

ulimit -s 65532

The matrix the above code uses is the unit matrix which is its own inverse. You can also use any other matrix that has an inverse, and invert it twice to check that you get the original matrix.

score 0 · Answer 2 · answered Jun 07 '19 at 09:58

Personally, I tried implementing matrix inverse using following two methods:

Matrix inverse using adjoint method.
Matrix inverse using Gauss-Jordan method.

And I found that, among these two implementations, Gauss-Jordan excelled. I tried it for 100x100 matrices, and got results in less than 2 sec on my machine. Not tried 1000x1000 though. Don't know other better algorithms for finding inverse. Gauss-Jordan is not that complicated to implement.

C algorithm to invert 1024x1024 matrix?

2 Answers2