C算法反转1024x1024矩阵？答案

【问题标题】：C algorithm to invert 1024x1024 matrix?C算法反转1024x1024矩阵？
【发布时间】：2016-10-26 07:53:47
【问题描述】：

这是我在这个网站上的第一个问题。我（拼命地）试图在我的程序中反转一个大矩阵。我想使用 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

我做错了什么？

【问题讨论】：

你对另一个答案有什么问题？
该帖子包含有效的 C 代码。
请考虑包含一些您已经编写（或从该链接使用）的代码。或者考虑更详细地描述运行代码所必须的配置/设置。读者需要做的事情越少（例如点击并阅读另一个链接）越好；你更有可能得到有意义的答案。
Computing the inverse of a matrix using lapack in C的可能重复
stackoverflow.com/questions/3519959/… 中的答案由@Olaf 评论：There is no language C/C++. The code you show is C++. But the question is about C.。确实，#include <cstdio> 是 C++，而 #include <stdio.h> 是 C。这就是编译器抱怨 cstdio 既不是文件也不是目录的原因。

标签： c algorithm matrix lapack inversion

【解决方案1】：

就个人而言，我尝试使用以下两种方法实现矩阵求逆：

使用伴随方法的矩阵求逆。
使用 Gauss-Jordan 方法的矩阵求逆。

我发现，在这两个实现中，Gauss-Jordan 表现出色。我尝试了 100x100 矩阵，并在我的机器上不到 2 秒就得到了结果。虽然没有尝试过 1000x1000。不知道其他更好的求逆算法。 Gauss-Jordan 实现起来并不复杂。

【讨论】：

【解决方案2】：

你可以验证这个 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

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

【讨论】：