使用 OpenMP 进行 Cholesky 分解
Posted
技术标签:
【中文标题】使用 OpenMP 进行 Cholesky 分解【英文标题】：Cholesky decomposition with OpenMP 【发布时间】：2014-04-24 03:03:44 【问题描述】：我有一个项目，我们使用Cholesky Decomposition 求解大型（超过 3000x3000）正定稠密矩阵的逆矩阵。该项目使用 Java，我们使用的是 CERN Colt BLAS library。分析代码表明 Cholesky 分解是瓶颈。
我决定尝试使用 OpenMP 并行化 Cholesky 分解，并将其用作 Java 中的 DLL（使用 JNA）。我从Rosetta Code 开始使用 C 中的 Cholesky 分解代码。
我注意到，除了对角元素之外，列中的值是独立的。所以我决定串行计算对角元素，并行计算列的其余值。我还交换了循环的顺序，以便内循环在行上运行，外循环在列上运行。在我的 4 核 (8 HT) 系统上，串行版本比 RosettaCode 的版本稍慢，但并行版本比 RosettaCode 版本快六倍。在 Java 中使用 DLL 可以加快我们的结果六次也是如此。代码如下：
double *cholesky(double *A, int n)
double *L = (double*)calloc(n * n, sizeof(double));
if (L == NULL)
exit(EXIT_FAILURE);
for (int j = 0; j <n; j++)
double s = 0;
for (int k = 0; k < j; k++)
s += L[j * n + k] * L[j * n + k];
L[j * n + j] = sqrt(A[j * n + j] - s);
#pragma omp parallel for
for (int i = j+1; i <n; i++)
double s = 0;
for (int k = 0; k < j; k++)
s += L[i * n + k] * L[j * n + k];
L[i * n + j] = (1.0 / L[j * n + j] * (A[i * n + j] - s));
return L;
你可以在http://coliru.stacked-crooked.com/a/6f5750c20d456da9找到完整的测试代码
我最初认为当列的剩余元素与线程数相比较小时，错误共享将是一个问题，但事实似乎并非如此。我试过了
#pragma omp parallel for schedule(static, 8) // a cache line is 8 doubles
我还没有找到关于如何并行化 Choleskey 分解的明确示例。我不知道我所做的是否理想。例如，它能否在 NUMA 系统上正常工作？
也许一般来说，基于任务的方法更好？在http://courses.engr.illinois.edu/cs554/fa2013/notes/07_cholesky.pdf 的幻灯片 7-9 中，有一个使用“细粒度任务”的并行 Cholesky 分解示例。我还不清楚如何实现这一点。
我有两个问题，具体的和一般的。您对如何使用 OpenMP 改进 Cholesky 分解的实现有什么建议吗？您能否建议的不同实现，例如有任务吗？
编辑：这里要求的是我用来计算 s 的 AVX 函数。它没有帮助
double inner_sum_AVX(double *li, double *lj, int n)
__m256d s4;
int i;
double s;
s4 = _mm256_set1_pd(0.0);
for (i = 0; i < (n & (-4)); i+=4)
__m256d li4, lj4;
li4 = _mm256_loadu_pd(&li[i]);
lj4 = _mm256_loadu_pd(&lj[i]);
s4 = _mm256_add_pd(_mm256_mul_pd(li4, lj4), s4);
double out[4];
_mm256_storeu_pd(out, s4);
s = out[0] + out[1] + out[2] + out[3];
for(;i<n; i++)
s += li[i]*lj[i];
return s;
【问题讨论】：
您的加速很好，我不认为仅使用 OpenMP 可以获得其他性能。您可以尝试使用 AVX/SSE 来计算s。也许可以做一些改进，但那将是数学上的..
@user3018144，我同意 6x 已经相当不错了。我想主要问题是我是否会在 NUMA 系统上获得相同的加速，或者是否可以改进单线程代码（超线程有很大帮助的事实告诉我它可以）。关于 s 上的 AVX/SSE 的好点。我已经考虑了几天，但还没有尝试过。最好使用 SIMD 同时在多行上执行此操作，但对角线使其变得困难。
如果我错了，请纠正我，但您似乎正在使用 omp 并行化内部循环。如果你想让多个线程并行计算，你不想启动很多短命的线程，而是保持一些与 CPU 数量相似的线程连续忙碌。我会尝试并行化 外循环，这样线程开销（创建、调度、运行、终止）会更低。
@EOF，要是这么简单就好了……每一列都取决于它之前所有列的值。它们必须按顺序计算。但是除了第一个元素之外，列中的值可以并行完成。
@EOF，现在没有 OpenMP 运行时会在并行区域结束时杀死工作线程。相反，所有线程都保存在一个池中，并在进入新的并行区域时（廉价地）被调用。 MSVC 的 OpenMP 运行时使用 Windows 本机线程池实现，因此以最小的开销获得最大的性能。
【参考方案1】：
我设法让 SIMD 与 Cholesky 分解一起工作。我使用循环平铺来做到这一点，就像我之前在矩阵乘法中使用的那样。解决方案并非易事。以下是我的 4 核/8 HT Ivy Bridge 系统上 5790x5790 矩阵的时间（eff = GFLOPS/(peak GFLOPS)）：
double floating point peak GFLOPS 118.1
1 thread time 36.32 s, GFLOPS 1.78, eff 1.5%
8 threads time 7.99 s, GFLOPS 8.10, eff 6.9%
4 threads+AVX time 1.36 s, GFLOPS 47.64, eff 40.3%
4 threads MKL time 0.68 s, GFLOPS 95.14, eff 80.6% // from LAPACKE_dpotrf
single floating point peak GFLOPS 236.2
1 thread time 33.88 s, GFLOPS 1.91, eff 0.8%
8 threads time 4.74 s, GFLOPS 13.64, eff 5.8%
4 threads+AVX time 0.78 s, GFLOPS 82.61, eff 35.0%
新方法双倍速度提高 25 倍，单倍速度提高 40 倍。现在的效率大约是峰值 FLOPS 的 35-40%。使用矩阵乘法，我在自己的代码中使用 AVX 获得了 70%。我不知道 Cholesky 分解会发生什么。与矩阵乘法不同，该算法是部分串行的（在计算对角块时，在下面的代码中称为 triangle）。
更新：我在 MKL 的 2 分之内。我不知道我应该为此感到自豪还是感到尴尬，但显然我的代码仍然可以显着改进。我在上面找到了PhD thesis，这表明我的块算法是一种常见的解决方案，所以我设法重新发明了***。
我使用 32x32 瓷砖作为双精度瓷砖，使用 64x64 瓷砖作为浮动瓷砖。我还将每个图块的内存重新排序为连续的并成为它的转置。我定义了一个新的矩阵生产函数。矩阵乘法定义为：
C_i,j = A_i,k * B_k,j //sum over k
我意识到在 Cholesky 算法中有一些非常相似的东西
C_j,i = A_i,k * B_j,k //sum over k
通过编写图块的转置，我几乎可以将优化后的函数用于矩阵乘法here（我只需要更改一行代码）。这是主要功能：
reorder(tmp,B,n2,bs);
for(int j=0; j<nb; j++)
#pragma omp parallel for schedule(static) num_threads(ncores)
for(int i=j; i<nb; i++)
for(int k=0; k<j; k++)
product(&B[stride*(nb*j+k)],&B[stride*(nb*i+k)],&B[stride*(nb*i+j)],bs);
triangle(&B[stride*(nb*j+j)], bs);
#pragma omp parallel for schedule(static)
for(int i=j+1; i<nb; i++)
block(&B[stride*(nb*i+j)],&B[stride*(nb*j+j)],bs);
reorder_inverse(B,tmp,n2,bs);
这是其他功能。我有 SSE2、AVX 和 FMA 的六个产品函数，每个函数都有 double 和 float 版本。我这里只展示一个用于 AVX 和 double 的：
template <typename Type>
void triangle(Type *A, int n)
for (int j = 0; j < n; j++)
Type s = 0;
for(int k=0; k<j; k++) s+= A[k*n+j]*A[k*n+j];
//if((A[j * n + j] - s)<0) printf("asdf3 j %d, %f %f\n", j, A[j * n + j] - s, sqrt(A[j * n + j] - s));
A[j*n+j] = sqrt(A[j*n+j] - s);
Type fact = 1.0/A[j*n+j];
for (int i = j+1; i<n; i++)
Type s = 0;
for(int k=0; k<j; k++) s+=A[k*n+i]*A[k*n+j];
A[j*n+i] = fact * (A[j*n+i] - s);
template <typename Type>
void block(Type *A, Type *B, int n)
for (int j = 0; j <n; j++)
Type fact = 1.0/B[j*n+j];
for (int i = 0; i<n; i++)
Type s = 0;
for(int k=0; k<j; k++)
s += A[k*n+i]*B[k*n+j];
A[j*n+i] = fact * (A[j*n+i] - s);
template <typename Type>
void reorder(Type *A, Type *B, int n, int bs)
int nb = n/bs;
int stride = bs*bs;
//printf("%d %d %d\n", bs, nb, stride);
#pragma omp parallel for schedule(static)
for(int i=0; i<nb; i++)
for(int j=0; j<nb; j++)
for(int i2=0; i2<bs; i2++)
for(int j2=0; j2<bs; j2++)
B[stride*(nb*i+j) + bs*j2+i2] = A[n*bs*i + j*bs + n*i2 + j2];
template <typename Type>
void reorder_inverse(Type *A, Type *B, int n, int bs)
int nb = n/bs;
int stride = bs*bs;
//printf("%d %d %d\n", bs, nb, stride);
#pragma omp parallel for schedule(static)
for(int i=0; i<nb; i++)
for(int j=0; j<nb; j++)
for(int i2=0; i2<bs; i2++)
for(int j2=0; j2<bs; j2++)
B[n*bs*i + j*bs + n*i2 + j2] = A[stride*(nb*i+j) + bs*j2+i2];
extern "C" void product32x32_avx(double *a, double *b, double *c, int n)
for(int i=0; i<n; i++)
__m256d t1 = _mm256_loadu_pd(&c[i*n + 0]);
__m256d t2 = _mm256_loadu_pd(&c[i*n + 4]);
__m256d t3 = _mm256_loadu_pd(&c[i*n + 8]);
__m256d t4 = _mm256_loadu_pd(&c[i*n + 12]);
__m256d t5 = _mm256_loadu_pd(&c[i*n + 16]);
__m256d t6 = _mm256_loadu_pd(&c[i*n + 20]);
__m256d t7 = _mm256_loadu_pd(&c[i*n + 24]);
__m256d t8 = _mm256_loadu_pd(&c[i*n + 28]);
for(int k=0; k<n; k++)
__m256d a1 = _mm256_set1_pd(a[k*n+i]);
__m256d b1 = _mm256_loadu_pd(&b[k*n+0]);
t1 = _mm256_sub_pd(t1,_mm256_mul_pd(a1,b1));
__m256d b2 = _mm256_loadu_pd(&b[k*n+4]);
t2 = _mm256_sub_pd(t2,_mm256_mul_pd(a1,b2));
__m256d b3 = _mm256_loadu_pd(&b[k*n+8]);
t3 = _mm256_sub_pd(t3,_mm256_mul_pd(a1,b3));
__m256d b4 = _mm256_loadu_pd(&b[k*n+12]);
t4 = _mm256_sub_pd(t4,_mm256_mul_pd(a1,b4));
__m256d b5 = _mm256_loadu_pd(&b[k*n+16]);
t5 = _mm256_sub_pd(t5,_mm256_mul_pd(a1,b5));
__m256d b6 = _mm256_loadu_pd(&b[k*n+20]);
t6 = _mm256_sub_pd(t6,_mm256_mul_pd(a1,b6));
__m256d b7 = _mm256_loadu_pd(&b[k*n+24]);
t7 = _mm256_sub_pd(t7,_mm256_mul_pd(a1,b7));
__m256d b8 = _mm256_loadu_pd(&b[k*n+28]);
t8 = _mm256_sub_pd(t8,_mm256_mul_pd(a1,b8));
_mm256_storeu_pd(&c[i*n + 0], t1);
_mm256_storeu_pd(&c[i*n + 4], t2);
_mm256_storeu_pd(&c[i*n + 8], t3);
_mm256_storeu_pd(&c[i*n + 12], t4);
_mm256_storeu_pd(&c[i*n + 16], t5);
_mm256_storeu_pd(&c[i*n + 20], t6);
_mm256_storeu_pd(&c[i*n + 24], t7);
_mm256_storeu_pd(&c[i*n + 28], t8);
【讨论】：
重新发明***并不丢人。它只是表明你的想法与其他在你之前做过的有成就的人一样。你仍然必须弄清楚。 你能不能写一个使用这段代码的例子？我想我想通了，但我不确定要调用什么参数。 bs = 块大小，nb = 块数，对吧？ @ТимофейЛомоносов，我的部分代码我还不能发布，但这里是主要功能coliru.stacked-crooked.com/a/9c00d5ac7332e1c8 @ТимофейЛомоносов，这里是AVX的产品功能coliru.stacked-crooked.com/a/4c934a4775dcd2f1 @ТимофейЛомоносов，如果您想要 SSE2 和 FMA 的产品功能，请告诉我，但这应该足以让您弄清楚。如果我有时间，我会清理我无法发布的代码并将整个事情公开。以上是关于使用 OpenMP 进行 Cholesky 分解的主要内容，如果未能解决你的问题，请参考以下文章
三十分钟理解：矩阵Cholesky分解，及其在求解线性方程组矩阵逆的应用