0

LAPACK用来逆矩阵:我做了一个引用传递,即通过处理地址。下面的函数带有一个输入矩阵和一个由它们的地址引用的输出矩阵。

问题是我必须将其转换F_matrix1D array,我认为这是在运行时级别的性能浪费:我可以找到哪种方法来摆脱这个耗时的补充任务我想如果我打电话很多次功能matrix_inverse_lapack

在相关功能下方:

// Passing Matrixes by Reference
void matrix_inverse_lapack(vector<vector<double>> const &F_matrix, vector<vector<double>> &F_output) {

  // Index for loop and arrays
  int i, j, ip, idx;

  // Size of F_matrix
  int N = F_matrix.size();

  int *IPIV = new int[N];

  // Statement of main array to inverse
  double *arr = new double[N*N];

  // Output Diagonal block
  double *diag = new double[N];

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

  // LAPACKE routines
  int info1 = LAPACKE_dgetrf(LAPACK_ROW_MAJOR, N, N, arr, N, IPIV);
  int info2 = LAPACKE_dgetri(LAPACK_ROW_MAJOR, N, arr, N, IPIV);

 for (i = 0; i<N; i++){
    for (j = 0; j<N; j++){
      idx = i*N + j;
      F_output[i][j] = arr[idx];
    }
  }

  delete[] IPIV;
  delete[] arr;
}

例如,我这样称呼它:

vector<vector<double>> CO_CL(lsize*(2*Dim_x+Dim_y), vector<double>(lsize*(2*Dim_x+Dim_y), 0));

... some code

matrix_inverse_lapack(CO_CL, CO_CL);

反转的表现不是预期的,我认为这是由于2D -> 1D我在函数中描述的这种转换matrix_inverse_lapack

更新

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

我有AMD Radeon Pro 5600M显卡。我已经默认安装了 Big Sur 版本的所有框架 OpenCL(也许我这样说是错的)。任何人都可以告诉安装 MAGMA 的程序。我在http://magma.maths.usyd.edu.au/magma/上看到了一个 MAGMA 软件,但它真的很贵而且不符合我想要的:我只需要所有的 SDK(头文件和库) ,如果可能的话,用我的 GPU 卡构建。我已经在我的 MacOS 上安装了所有英特尔 OpenAPI SDK。也许,我可以将它链接到 MAGMA 安装。

我看到另一个链接https://icl.utk.edu/magma/software/index.html MAGMA 似乎是公开的:上面的非免费版本没有链接,不是吗?

4

1 回答 1

1

首先让我抱怨 OP 没有提供所有必要的数据。该程序几乎完成,但它不是一个最小的、可重现的示例。这很重要,因为(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 3(数据未显示)。对于 N = 3000,反转矩阵的时间约为 16 秒,对于 N = 10,约为 5 -6 s(5 微秒)。

我认为即使是 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 数组时很奇怪且过时。不,Lapack 是对的。

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

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

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

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

编辑:替代库

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

至于现在(2021 年年中),一个不错的建议是:

如果 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 次)。

于 2021-07-09T12:32:24.880 回答