C++中的并行稠密矩阵对角化

Question

我已经使用OPENMP编写了一个并行的C++代码。我一直在使用EIGEN来处理线性代数。然而，对于稠密矩阵的特征求解器，Eigen并不能并行化。我已经用openmp构建了OPENBLAS，但是当使用更多线程(使用LAPACKE_zheevd函数(时，我没有看到任何性能提升。

你有什么建议？

这是我一直在使用的代码：

#include <stdlib.h>
#include <stdio.h>
#include <omp.h>
// C-wrapper to the Fortran Lapack lib.
#include <lapacke.h>
#include <iostream>
#include <chrono>
int main()
{
const int N = 1024;

printf("Matrix size=%in", N);
printf("Number of threads=%in", omp_get_max_threads());

int LDA=N, info,i,j;
double w[N];

lapack_complex_double * a;

a = (lapack_complex_double *) malloc( N*LDA*sizeof(lapack_complex_double) );

srand(999); /* Dense random matrix */
for (i = 0; i < N; i++ )
for(j = 0; j < N; j++ )
a[i*LDA + j] = lapack_make_complex_double(rand() , rand());
// LAPACKE_zheevd: computes all
// eigenvalues and eigenvectors of a
// complex Hermitian matrix A using divide
// and conquer algorithm
auto start = std::chrono::system_clock::now();
info = LAPACKE_zheevd( LAPACK_ROW_MAJOR, 'V', 'L', N, a, LDA, w );
auto duration = std::chrono::duration_cast< std::chrono::milliseconds >( std::chrono::system_clock::now() - start);
std::cout << duration.count() << std::endl;
/* Check for convergence */

if( info > 0 ) {
printf( "The algorithm failed to compute eigenvalues.n" );
exit( 1 );
}
/* Print the extreme eigenvalues */
printf("Smallest eigen value=%6.2fn",w[0]);
printf("Largest value=%6.2fnn",w[N-1]);
}

Answer 1

LAPACKE_zheevd是LAPACK基元，而不是BLAS基元。前者提供高级线性代数函数(如线性求解器、奇异值分解等(。后者提供非常基本的线性代数函数，如(密集(矩阵乘法、向量矩阵乘积、点积等。LAPACK库在内部使用BLAS库来尽可能地加快速度(通常通过重写算法来使用高度优化的矩阵乘法(。

OpenBLAS主要提供BLAS原语，而不是真正的LAPACK原语。它只实现了几个LAPACK函数，因此使用了它们的内部并行矩阵乘法。关于目标输入，所提供的LAPACK函数的实现可能远不是最佳的(只有当矩阵乘法足够大以便并行运行时，它才是快速的(。大多数函数都来自NetLib实现，它是一个标准的开源LAPACK包。关于ZPut-it，OpenBLAS是一个非常好的BLAS实现，但不是一个很好的LAPACK实现(尽管它可以比NetLib实现提供实质性的改进(。

事实上，大多数LAPACK函数的并行化非常困难。有几个高技能的研究&工程团队几十年来一直致力于该主题，以便编写快速并行LAPACK原语。AFAIK，结果喜忧参半，尽管到目前为止已经做出了努力，但对于一些原语来说，可伸缩性相当令人失望。这是一项艰巨的任务，尤其是因为他们需要关心算法的数值稳定性和多核上的平衡计算。

除此之外，还有一个巨大的问题：内存限制算法无法扩展。事实上，很少有内核能够饱和计算节点的大部分内存带宽。问题是许多(部分(LAPACK基元都是内存绑定的。更不用说，由于记忆墙，这在未来可能不会更好。

---

特别是关于LAPACKE_zheevd，OpenBLAS的一个开放问题表明ZHEEV的实现实际上效率很低，使用多个线程往往会降低性能。"英特尔MKL"可更有效地实现这一点(尽管可扩展性仍然不太好(。尝试PLASMA(专门为多核架构设计(这样的替代库可能是个好主意。另一种解决方案是基于关于该主题的最新研究论文来实现这一点(通常通过实现基于任务的平铺奇异值分解(，但这远不是一件容易的事。