为什么我在 C++ 中使用 LAPACKE 进行矩阵求逆如此缓慢：MAGMA Alternative and set up

Question

我正在使用 LAPACK 对矩阵进行求逆：我进行了引用传递，即处理地址。下面是带有输入矩阵和由其地址引用的输出矩阵的函数。

问题是我不得不将 F_matrix 转换为 1D array，我认为这在运行时级别上浪费了性能：我可以找到什么方法来摆脱它补充任务很耗时我想如果我打电话很多次函数 matrix_inverse_lapack.

更新

有人建议我在我的 MacOS Big Sur 11.3 上安装 MAGMA，但我在设置它时遇到了很多困难。

我有一张 AMD Radeon Pro 5600M 显卡。我已经默认安装了 Big Sur 版本的所有 Framework OpenCL（也许我这么说是错误的）。任何人都可以告诉安装 MAGMA 的过程。我看到 http://magma.maths.usyd.edu.au/magma/ 上存在一个 MAGMA 软件，但它真的很贵而且不符合我的要求：我只需要所有的 SDK（头文件和库），如果可能的话用我的 GPU 卡构建。我已经在我的 MacOS 上安装了所有英特尔 OpenAPI SDK。也许，我可以 link 把它放到 MAGMA 安装中。

我看到另一个 link https://icl.utk.edu/magma/software/index.html 哪里 MAGMA 似乎是 public : 有 none link 上面的非免费版本，是吗没有吗？

Answer 1

首先让我抱怨OP没有提供所有必要的数据。该程序几乎完成，但不是 minimal, reproducible example。这很重要，因为 (a) 它浪费时间，并且 (b) 它隐藏了潜在的相关信息，例如。关于矩阵初始化。其次，OP 没有提供任何关于编译的细节，这又可能是相关的。最后但同样重要的是，OP 没有检查 Lapack 函数可能出现的错误的状态代码，这对于正确解释结果也很重要。

让我们从一个最小的可重现示例开始：

#include <lapacke.h>
#include <vector>
#include <chrono>
#include <iostream>

using Matrix = std::vector<std::vector<double>>;

std::ostream &operator<<(std::ostream &out, Matrix const &v)
{
    const auto size = std::min<int>(10, v.size());
    for (int i = 0; i < size; i++)
    {
        for (int j = 0; j < size; j++)
        {
            out << v[i][j] << "\t";
        }
        if (size < std::ssize(v)) out << "...";
        out << "\n";
    }
    return out;
}

void matrix_inverse_lapack(Matrix const &F_matrix, Matrix &F_output, std::vector<int> &IPIV_buffer,
                           std::vector<double> &matrix_buffer)
{
    //  std::cout << F_matrix << "\n";
    auto t0 = std::chrono::steady_clock::now();

    const int N = F_matrix.size();

    for (int i = 0; i < N; i++)
    {
        for (int j = 0; j < N; j++)
        {
            auto idx = i * N + j;
            matrix_buffer[idx] = F_matrix[i][j];
        }
    }

    auto t1 = std::chrono::steady_clock::now();
    // LAPACKE routines
    int info1 = LAPACKE_dgetrf(LAPACK_ROW_MAJOR, N, N, matrix_buffer.data(), N, IPIV_buffer.data());
    int info2 = LAPACKE_dgetri(LAPACK_ROW_MAJOR, N, matrix_buffer.data(), N, IPIV_buffer.data());
    auto t2 = std::chrono::steady_clock::now();

    for (int i = 0; i < N; i++)
    {
        for (int j = 0; j < N; j++)
        {
            auto idx = i * N + j;
            F_output[i][j] = matrix_buffer[idx];
        }
    }
    auto t3 = std::chrono::steady_clock::now();

    auto whole_fun_time = std::chrono::duration<double>(t3 - t0).count();
    auto lapack_time = std::chrono::duration<double>(t2 - t1).count();
    //   std::cout << F_output << "\n";
    std::cout << "status: " << info1 << "\t" << info2 << "\t" << (info1 == 0 && info2 == 0 ? "Success" : "Failure")
              << "\n";
    std::cout << "whole function:            " << whole_fun_time << "\n";
    std::cout << "LAPACKE matrix operations: " << lapack_time << "\n";
    std::cout << "conversion:                " << (whole_fun_time - lapack_time) / whole_fun_time * 100.0 << "%\n";
}

int main(int argc, const char *argv[])
{
    const int M = 5;  // numer of test repetitions

    const int N = (argc > 1) ? std::stoi(argv[1]) : 10;
    std::cout << "Matrix size = " << N << "\n";

    std::vector<int> IPIV_buffer(N);
    std::vector<double> matrix_buffer(N * N);

    // Test matrix_inverse_lapack M times
    for (int i = 0; i < M; i++)
    {
        Matrix CO_CL(N);
        for (auto &v : CO_CL) v.resize(N);

        int idx = 1;
        for (auto &v : CO_CL)
        {
            for (auto &x : v)
            {
                x = idx + 1.0 / idx;
                idx++;
            }
        }
        matrix_inverse_lapack(CO_CL, CO_CL, IPIV_buffer, matrix_buffer);
    }
}

在这里，operator<< 有点矫枉过正，但对于任何想要半手动验证代码是否有效（通过取消注释第 26 和 58 行）并确保代码正确的人来说可能有用衡量其性能很重要。

代码可以用

编译

g++ -std=c++20 -O3  main.cpp -llapacke

该程序依赖于一个外部库，lapacke，需要安装它，头文件 + 二进制文件，用于编译代码和运行。

我的代码与 OP 的有点不同：它更接近于“现代 C++”，因为它避免使用裸指针；我还在 matrix_inverse_lapack 中添加了外部缓冲区，以抑制内存分配器和释放器的持续启动，这是一个以可衡量的方式减少 2D-1D-2D 转换开销的小改进。我还必须初始化矩阵并找到一种方法来在 OP 的脑海中读取 N 的值可能是什么。我还添加了一些计时器读数用于基准测试。除此之外，代码的逻辑没有改变。

现在在像样的工作站上进行基准测试。它列出了转换所用时间相对于 matrix_inverse_lapack 所用总时间的百分比。换句话说，我测量转换开销：

 N =   10, 3.5%   
 N =   30, 1.5%   
 N =  100, 1%   
 N =  300, 0.5%   
 N = 1000, 0.35%  
 N = 3000, 0.1%

Lapack 花费的时间很好地缩放为 N³，正如预期的那样（数据未显示）。 N = 3000 时矩阵求逆的时间约为 16 秒，N = 10 时约为 5^-6 秒（5 微秒）。

我假设即使是 3% 的开销也是完全可以接受的。我相信 OP 使用大小大于 100 的矩阵，在这种情况下，1% 或以下的开销当然是可以接受的。

那么 OP（或有类似问题的任何人）可能做错了什么以获得“不可接受的开销转换值”？这是我的候选名单

编译不当
不正确的矩阵初始化（用于测试）
基准测试不当

1.编译不当

如果忘记在 Release 模式下编译，优化的 Lapacke 将与未优化的转换竞争。在我的机器上，当 N = 20 时，开销达到 33% 的峰值。

2。不正确的矩阵初始化（用于测试）

如果像这样初始化矩阵：

        for (auto &v : CO_CL)
        {
            for (auto &x : v)
            {
                x = idx; // rather than, eg., idx + 1.0/idx
                idx++;
            }
        }

那么矩阵是奇异的，lapack returns 相当快，状态不为0。这增加了转换部分的相对重要性。但是奇异矩阵不是人们想要反转的（这是不可能的）。

3.基准测试不当

这是 N = 10 的程序输出示例：

 ./a.out 10 
 Matrix size = 10
 status: 0  0   Success
 whole function:            0.000127658
 LAPACKE matrix operations: 0.000126783
 conversion:                0.685425%
 status: 0  0   Success
 whole function:            1.2497e-05
 LAPACKE matrix operations: 1.2095e-05
 conversion:                3.21677%
 status: 0  0   Success
 whole function:            1.0535e-05
 LAPACKE matrix operations: 1.0197e-05
 conversion:                3.20835%
 status: 0  0   Success
 whole function:            9.741e-06
 LAPACKE matrix operations: 9.422e-06
 conversion:                3.27482%
 status: 0  0   Success
 whole function:            9.939e-06
 LAPACKE matrix operations: 9.618e-06
 conversion:                3.2297%

可以看到第一次调用lapack函数比后续调用要多10倍的时间。这是一个相当稳定的模式，好像 Lapack 需要一些时间来进行自我初始化。它会严重影响小 N 的测量。

4.还能做什么？

OP apperas 相信他的 2D 数组方法很好，Lapack 将 2D 数组打包成 1D 数组的方法很奇怪而且过时。不对，Lapack才是对的

如果将二维数组定义为 vector<vector<double>>，则可以获得一个优势：代码简单。这是有代价的。这种矩阵的每一行都与其他行分开分配。因此，一个 100×100 的矩阵可以存储在 100 个完全不同的存储块中。这对缓存（和预取器）利用率有不良影响。 Lapck（和其他线性代数包）强制将数据压缩到一个连续的数组中。这是为了最大限度地减少缓存和预取器未命中。如果 OP 从一开始就使用这种方法，他可能会获得比他们现在为转换支付的 1-3% 以上的收益。

这种压缩至少可以通过三种方式实现。

为二维矩阵编写自定义 class，将内部数据存储在一维数组中并方便地访问成员函数（例如：operator ()），或者找一个库来做这件事
为 std::vector 编写自定义分配器（或查找库）。这个分配器应该从一个预先分配的一维向量中分配内存，该向量与 Lapack
使用 std::vector<double*> 并使用指向预分配一维数组的适当元素的地址初始化指针。

上述每个解决方案都会强制对周围代码进行一些更改，OP 可能不想这样做。一切都取决于代码的复杂性和预期的性能提升。

编辑：替代库

另一种方法是使用已知高度优化的库。 Lapack 本身可以被视为具有许多实现的标准接口，OP 可能会使用未优化的接口。选择哪个库可能取决于 hardware/software OP 感兴趣的平台，并且可能会随时间变化。

就目前（2021 年年中）而言，一个不错的建议是：

Lapack https://www.netlib.org/lapack/
图集https://en.wikipedia.org/wiki/Automatically_Tuned_Linear_Algebra_Software http://math-atlas.sourceforge.net/
OpenBlas https://www.openblas.net/
岩浆https://developer.nvidia.com/magma
等离子https://bitbucket.org/icl/plasma/src/main/

如果 OP 使用大小至少为 100 的矩阵，那么面向 GPU 的 MAGMA 可能值得尝试。

一个更简单的（安装，运行ning）方法可能是使用并行 CPU 库，例如等离子体。 Plsama 是 Lapack 兼容的，它是由包括 Jack Dongarra 在内的一大批人开发的，它也应该很容易在本地编译它，因为它提供了 CMake 脚本。

基于并行 CPU 的多核实现可以在多大程度上优于 LU 分解的单线程实现的示例可以在此处找到：https://cse.buffalo.edu/faculty/miller/Courses/CSE633/Tummala-Spring-2014-CSE633.pdf（简短回答：5对于大小为 1000 的矩阵是 15 倍）。

为什么我在 C++ 中使用 LAPACKE 进行矩阵求逆如此缓慢：MAGMA Alternative and set up

Why my inversions of matrices are such slow with LAPACKE in C++ : MAGMA Alternative and set up

c++

matrix

inversion

lapack

magma

更新