使用 BLAS 和 OpenMP 优化特征重构（矩阵 - 对角矩阵 - 矩阵）乘积 C++

Question

我编写了一个 C++ 代码来求解线性系统 A.x = b，其中 A 是一个对称矩阵，首先使用 LAPACK(E) 对矩阵 A = V.D.V^T 进行对角化（因为我稍后需要特征值），然后求解 @987654324 @ 当然V 是正交的。

现在我想尽可能优化最后一个操作，例如通过使用 (C)BLAS 例程和 OpenMP。

这是我的幼稚实现：

// Solve linear system A.X = B for X (V contains eigenvectors and D eigenvalues of A)
void solve(const double* V, const double* D, const double* B, double* X, const int& N)
{
    #ifdef _OPENMP
    #pragma omp parallel for
    #endif
    for (int i=0; i<N; i++)
    {
        for (int j=0; j<N; j++)
        {
            for (int k=0; k<N; k++)
            {
                X[i] += B[j] * V[i+k*N] * V[j+k*N] / D[k];
            }
        }
    }
}

所有数组都是 C 样式的数组，其中 V 的大小为 N^2，D 的大小为 N，B 的大小为 N，X 的大小为 @987654334 @（并用零初始化）。

目前，这种幼稚的实现非常缓慢，是代码的瓶颈。任何提示和帮助将不胜感激！

谢谢

编辑 感谢 Jérôme Richard 的回答和评论，我通过调用 BLAS 并使用 OpenMP 并行化中间循环进一步优化了他的解决方案。在 1000x1000 矩阵上，这个解决方案比他的命题快约 4 倍，它本身比我的幼稚实现快 1000 倍。

对于N=1000 和N=2000，我发现#pragma omp parallel for simd 子句在两台分别具有 4 核和 20 核的不同机器上比其他替代方案更快。

void solve(const double* V, const double* D, const double* B, double* X, const int& N)
{

    double* sum = new double[N]{0.};

    cblas_dgemv(CblasColMajor,CblasTrans,N,N,1.,V,N,B,1,0.,sum,1);

    #pragma omp parallel for simd
    for (int i=0; i<N; ++i)
    {
        sum[i] /= D[i];
    }

    cblas_dgemv(CblasColMajor,CblasNoTrans,N,N,1.,V,N,sum,1,0.,X,1);

    delete [] sum;
}

Answer 1

此代码目前高度内存限制。因此，生成的程序可能很难扩展（只要启用了编译器优化）。实际上，在大多数常见系统（例如 1 个插槽非 NUMA 处理器）上，RAM 吞吐量是核心之间的共享资源，也是一种稀缺资源。此外，内存访问模式效率低，可以提高代码的算法复杂度。

为了加快计算速度，可以交换 j 和 k 循环，以便连续读取 V。此外，除以V[i+k*N] 和D[k] 在最内层循环中成为常数。然后计算可以分解更快，因为B[j] 和V[j+k*N] 也不依赖于i。由于总和预计算，生成的算法在 O(n^2) 而不是 O(n^3) 中运行！

最后，omp simd 可用于帮助编译器向量化代码，使其更快！

请注意，_OPENMP 在这里似乎没用，因为当 OpenMP 被禁用或不支持时，编译器应该忽略 #pragma。

// Solve linear system A.X = B for X (V contains eigenvectors and D eigenvalues of A)
void solve(const double* V, const double* D, const double* B, double* X, const int& N)
{
    std::vector<double> kSum(N);

    #pragma omp parallel for
    for (int k=0; k<N; k++)
    {
        const double sum = 0.0;

        #pragma omp simd reduction(+:sum)
        for (int j=0; j<N; j++)
        {
            sum += B[j] * V[j+k*N];
        }

        kSum[k] = sum / D[k];
    }

    // Loop tiling can be used to speed up this section even more.
    // The idea is to swap i-based and j-based loops and work on thread-private copies
    // of X and finally sum the thread-private versions into a global X.
    // The resulting code should work on contiguous data and can even be vectorized.
    #pragma omp parallel for
    for (int i=0; i<N; i++)
    {
        double sum = X[i];

        for (int k=0; k<N; k++)
        {
            sum += kSum[k] * V[i+k*N];
        }

        X[i] = sum;
    }
}

新代码应该比原来的代码快几个数量级（但仍然受内存限制）。请注意，结果可能会有些不同（因为浮点运算并不是真正的关联），但我希望结果更准确。