我有一个项目，我们用Cholesky分解来求解大型（超过3000x3000)正定稠密矩阵的逆。该项目在Java，我们使用的是欧洲核子研究中心Colt BLAS图书馆。对代码的分析表明Cholesky分解是瓶颈。
我决定尝试使用OpenMP并行化Cholesky分解，并将其作为Java中的DLL（与JNA一起使用）。我从C语言中的Cholesky分解代码开始，代码来自Rosetta代码。
我注意到的是，一列中除了对角线元素以外的值都是独立的。所以我决定串行计算对角线元素，并行计算列的其余值。我还交换了循环的顺序，使内部循环在行上运行，外部循环在列上运行。串行版本比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.staked-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分解，你有什么建议吗？你能用OpenMP为Cholesky分解提供一种不同的实现方式吗？例如，用任务来实现Cholesky分解吗？
编辑:这里需要的是我用来计算
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;
}
我设法让SIMD处理Cholesky分解。我使用循环平铺来实现这一点，就像我以前在矩阵乘法中使用的那样。解决办法并不简单。以下是我的4 Core/8 HT Ivy Bridge系统上5790x5790矩阵的时间（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倍。效率约为目前峰值滤波器的35～40%。对于矩阵乘法，我在自己的代码中使用AVX可以达到70%。我不知道从Cholesky分解中能期待什么。与矩阵乘法不同，该算法是部分串行的（当计算对角线块时，在我下面的代码中称为
更新:我在MKL的2因子之内。我不知道我应该为此感到骄傲还是感到尴尬，但显然我的代码仍然可以得到显著的改进。我找到了一篇关于这个的博士论文，它表明我的块算法是一个通用的解决方案，所以我设法重新发明了轮子。
我使用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
通过编写转置的瓦片，我能够使用我的优化函数在这里几乎完全的矩阵乘法（我只需要改变一行代码）。主要功能如下:
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，每个与双和浮动版本。这里我只显示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);
}
}
