基于OpenMP的Cholesky分解
我有一个项目，我们用Cholesky分解求解大型正定稠密矩阵的逆(超过3000x3000)。该项目是用Java编写的，我们使用的是CERN 柯尔特布拉斯图书馆。对代码进行分析表明，Cholesky分解是瓶颈。
我决定使用OpenMP将Cholesky分解并行化，并将其作为Java (与JNA一起使用)使用。我从Rosetta码 C中的Cholesky分解代码开始。
我注意到的是，除了对角线元素之外，列中的值是独立的。因此，我决定以串行方式计算对角线元素，并并行计算列的其余值。我还交换了循环的顺序，以便内部循环在行上运行，外部循环在列上运行。串行版本比RosettaCode 的版本稍慢，但并行版本比我的4核(8 HT)系统上的RosettaCode版本快6倍，在Java中使用DLL的也使我们的结果提高了6倍。以下是代码：
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系统上运行良好吗？
也许以任务为基础的方法总体上更好？在cholesky.pdf的幻灯片7-9中，有一个使用“细粒度任务”并行cholesky分解的例子。我还不清楚如何实现这一点。
我有两个问题，一个是具体的，一个是一般性的。对于如何用OpenMP改进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;
}回答 1
Stack Overflow用户
发布于 2014-04-14 15:08:24
我想让SIMD和Cholesky分解一起工作。我用循环平铺来做这个，就像我以前在矩阵乘法中所用的一样。解决办法并非微不足道。下面是我的4核心/8 HT常春藤桥系统上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分解会带来什么。该算法是部分串行的(在计算对角线块时，在下面的代码中称为triangle )，与矩阵乘法不同。
更新：I在2 MKL的因子范围内。我不知道我应该为此感到骄傲还是感到尴尬，但显然我的代码仍然可以得到显著的改进。我在上面找到了一个PhD论文，这表明我的块算法是一个常见的解决方案，所以我成功地重新发明了轮子。
我使用32x32块表示double，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);
}
}https://stackoverflow.com/questions/22479258
复制相似问题
腾讯云开发者
