-9

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

在 C 中使用 lapack 计算矩阵的逆

谢谢你。

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

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

我究竟做错了什么？

4

2 回答 2

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

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

于 2016-06-24T15:47:48.740 回答
0

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

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

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

于 2019-06-07T09:58:55.697 回答