c++ - 使用BLAS和OpenMP优化本征重组(矩阵-对角矩阵-矩阵)产品C++

Question

我编写了一个C++代码来解决线性系统A.x = b，其中A是一个对称矩阵，方法是首先使用LAPACK(E)对角矩阵A = V.D.V^T(因为以后需要特征值)，然后求解x = A^-1.b = V^T.D^-1.V.b，当然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的大小为 N(并用零初始化)。
目前，这种幼稚的实现非常慢，并且是代码的瓶颈。任何提示和帮助将不胜感激!
谢谢
编辑
感谢JérômeRichard的回答和评论，我通过调用BLAS并将中间循环与OpenMP并行化来进一步优化了他的解决方案。在1000x1000的矩阵上，此解决方案的速度是他的提议的4倍左右，而这个提议本身比我的幼稚实现快1000倍。
我发现 #pragma omp parallel for simd子句比 N=1000和 N=2000分别在具有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;
    }
}

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