为什么矩阵的逆在C中的LAPACKE中如此缓慢++

首先,让我抱怨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 可能。我还添加了一些计时器读数用于基准测试。除此之外,代码的逻辑是不变的。

现在在一个不错的工作站上进行了基准测试。它列出了转换所花费的时间相对于所花费的总时间的百分比 矩阵_反转_ lapack 换言之,我测量转换开销:

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

拉帕克所花费的时间可以很好地缩放为N ^3. ,如预期(数据未显示)。对于N=3000,反转矩阵的时间约为16秒,约为5秒 ^6. s(5微秒),其中N=10。

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

那么,什么OP(或任何有类似问题的人)可能做错了,从而获得了“不可接受的开销转换值”?这是我的短名单

1. 汇编不当
2. 矩阵初始化不正确(用于测试)
3. 基准测试不当

1.编制不当

如果忘记在发布模式下编译,最终会出现优化的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很快返回,状态与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似乎认为他处理2D阵列的方法是好的,而Lapack将2D阵列打包为1D阵列是一种奇怪而过时的方法。不,拉帕克是对的。

如果将二维阵列定义为 vector<vector<double>> ,获得一个优点:代码简单。这是有代价的。这样一个矩阵的每一行都是与其他行分开分配的。因此,矩阵100乘100可以存储在100个完全不同的存储器块中。这会对缓存(和预取器)利用率产生不良影响。Lack(和其他线性代数包)在单个连续数组中强制数据的紧凑化。这样可以最大限度地减少缓存和预取器未命中。如果OP从一开始就使用这种方法,他可能会获得他们现在为转换支付的1-3%以上的收益。

这种紧凑化可以通过至少三种方式来实现。

• 为2D矩阵编写一个自定义类,内部数据存储在1D数组中,并方便访问成员函数(例如: operator () ),或者找一个做这件事的图书馆
• 为编写自定义分配器 std::vector (或者找一个图书馆)。这个分配器应该从预先分配的1D向量中分配内存,该向量与Lapack使用的数据存储模式完全匹配
• 使用 std::vector<double*> 并且初始化具有指向预分配的1D阵列的适当元素的地址的指针。

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

EDIT:替代库

另一种方法是使用 已知的 因为它是高度优化的。Lapack本身可以被重新归类为具有许多实现的标准接口,并且可能会发生OP使用未优化接口的情况。选择哪个库可能取决于OP感兴趣的硬件/软件平台,并且可能随时间变化。

就目前(2021年年中)而言,一个不错的建议是:

• 拉帕克 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可能值得一试。

一种更容易(安装、运行)的方式可能是使用并行CPU库,例如Plasma。Plsama是Lapack兼容的,它是由包括Jack Dongarra在内的一个大型团队开发的,它也应该很容易在本地编译,因为它提供了CMake脚本。

一个基于并行CPU的多核实现在多大程度上优于LU分解的单线程实现的例子可以在这里找到,例如: https://cse.buffalo.edu/faculty/miller/Courses/CSE633/Tummala-Spring-2014-CSE633.pdf (简短回答:对于大小为1000的矩阵,为5到15次)。