用OpenMP进行Cholesky分解
我有一个项目，我们可以用Cholesky分解来求解大的（超过3000×3000）正定密集matrix的逆。 该项目是在Java中，我们使用的是使用CERN 马驹BLAS库 。 分析代码表明Cholesky分解是瓶颈。
我决定尝试使用OpenMP并行化Cholesky分解，并将其用作Java中的DLL（使用JNA）。 我从Rosetta Code的 C中的Cholesky分解代码开始。
我注意到，除了对angular元素之外，列中的值是独立的。 所以我决定并行计算串行对angular元素和列的其余值。 我也交换了循环的顺序，以便内部循环遍历行和遍历列的外部循环。 串行版本比RosettaCode稍慢， 但并行版本比我的4核（8 HT）系统上的RosettaCode版本快6倍。 在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/6f5750c20d456da9find完整的testing代码
我最初认为，当一列的其余元素与线程的数量相比较小时，假分享会成为一个问题，但似乎并非如此。 我试过了
#pragma omp parallel for schedule(static, 8) // a cache line is 8 doubles
我还没有find明确的例子来说明如何将Choleskey分解并行化。 我不知道我做了什么是理想的。 例如，它会在NUMA系统上运行吗？
一般而言，基于任务的方法可能更好？ 在http://courses.engr.illinois.edu/cs554/fa2013/notes/07_cholesky.pdf的幻灯片7-9中，有一个使用“细粒度任务”的平行胆固醇分解的例子。 目前还不清楚如何实现这一点。
我有两个问题，具体和一般。 对于如何提高OpenMP Cholesky分解的实现性，您有什么build议吗？ 你可以build议一个不同的实施Cholesky分解与OpenMP例如与任务？
编辑：这里要求的是我用来计算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; }
我pipe理SIMD与Cholesky分解。 我使用了循环平铺，就像之前在matrix乘法中使用过的一样。 解决scheme不是微不足道的。 以下是我的4核/ 8 HT Ivy Bridge系统上5790x5790matrix的时间（eff = GFLOPS /（峰值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％左右。 通过matrix乘法，AVX在我自己的代码中可以达到70％。 我不知道从乔列斯基分解会有什么期望。 与matrix乘法不同，该algorithm是部分串行的（当计算对angular线块时，我的代码中称为triangle ）。
更新：我在MKL的2个因素之内。 我不知道我应该为此感到骄傲还是为此感到尴尬，但显然我的代码仍然可以大大改善。 我在这发现了一篇博士论文 ，这说明我的分块algorithm是一个通用的解决scheme，所以我设法重新发明了方向盘。
我使用双32×32瓷砖和64×64瓷砖浮法。 我也重新sorting每个瓷砖的内存是连续的，是它的转置。 我定义了一个新的matrix生产函数。 matrix乘法定义如下：
C_i,j = A_i,k * B_k,j //sum over k
我意识到在Choleskyalgorithm中有一些非常相似的东西
C_j,i = A_i,k * B_j,k //sum over k
通过编写瓷砖的转置，我几乎可以完全使用我的优化函数进行matrix乘法（我只需要更改一行代码）。 这里是主要function：
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);
以下是其他function。 我有六个SSE2，AVX和FMA的产品function，每个都有双重版本和浮动版本。 我只展示了一个AVX和双倍在这里：
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); } }
