C算法反转1024x1024矩阵？

Question

这是我在这个网站上的第一个问题。我(拼命)试图在我的程序中反转一个大矩阵。我想使用 lapack 来做到这一点，我发现这个线程看起来很有前途，但我认为它是用 C++ 语言编写的。你能帮我一下吗？

Computing the inverse of a matrix using lapack in C

谢谢。

更新:你是对的，答案有点不清楚。将我发布的程序合并到我的程序中后，出现以下错误:

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

我想问题是我没有正确安装llapack库，或者我在编译时包含它。

这就是我安装库的方式(从终端，我有 Ubuntu):

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

这就是我的编译方式:

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

我做错了什么？

Answer 1

您可以验证此 C 算法是否能正确执行小矩阵的求逆:

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⏎   

以上数字在本示例程序中针对 2*2 矩阵进行了验证:

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

现在我们需要定义一个更大的矩阵 1024*1024 并以同样的方式对其求逆。

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

要运行上面的代码，我还必须增加 Linux 上的堆栈大小:

ulimit -s 65532

上面代码使用的矩阵是单位矩阵，它是它自己的逆矩阵。您还可以使用具有逆矩阵的任何其他矩阵，并将其反转两次以检查是否获得原始矩阵。