example - openmp c++

Descomposición de Cholesky con OpenMP

Question

Tengo un proyecto donde resolvemos el inverso de matrices densas definidas positivas grandes (más de 3000x3000) utilizando descomposición de Cholesky . El proyecto está en Java y utilizamos la biblioteca BLAS Colt BLAS . El perfil del código muestra que la descomposición de Cholesky es el cuello de botella.

Decidí intentar y paralelizar la descomposición de Cholesky usando OpenMP y usarlo como una DLL en Java (con JNA). Empecé con el código de descomposición de Cholesky en C del código de Rosetta .

Lo que noté es que los valores en una columna a excepción del elemento diagonal son independientes. Así que decidí calcular los elementos diagonales en serie y el resto de los valores de la columna en paralelo. También cambié el orden de los bucles para que el bucle interno se ejecute sobre las filas y el bucle externo sobre las columnas. La versión en serie es un poco más lenta que la de RosettaCode, pero la versión paralela es seis veces más rápida que la versión RosettaCode en mi sistema de 4 núcleos (8 HT). Usar el DLL en Java también acelera nuestros resultados seis veces. Aquí está el código:

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;
}

Puede encontrar el código completo para probar esto en http://coliru.stacked-crooked.com/a/6f5750c20d456da9

Inicialmente pensé que compartir falsamente sería un problema cuando los elementos restantes de una columna eran pequeños en comparación con el número de subprocesos, pero ese no parece ser el caso. Lo intenté

#pragma omp parallel for schedule(static, 8) // a cache line is 8 doubles

No he encontrado ejemplos claros de cómo paralelizar la descomposición de Choleskey. No sé si lo que he hecho es ideal. Por ejemplo, ¿funcionará bien en un sistema NUMA?

Tal vez un enfoque basado en la tarea es mejor en general? En las diapositivas 7-9 en http://courses.engr.illinois.edu/cs554/fa2013/notes/07_cholesky.pdf hay un ejemplo de descomposición paralela de cholesky utilizando "tareas de grano fino". Todavía no tengo claro cómo implementar esto.

Tengo dos preguntas, específicas y generales. ¿Tiene alguna sugerencia sobre cómo mejorar mi implementación de la descomposición de Cholesky con OpenMP? ¿Puede sugerir una implementación diferente de la descomposición de Cholesky con OpenMP, por ejemplo, con tareas?

Editar: como se solicita aquí está la función AVX que utilicé para calcular s . No ayudó

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;
}

Answer 1

Logré hacer funcionar SIMD con la descomposición de Cholesky. Hice esto utilizando mosaico de bucle como lo he usado antes en la multiplicación de matrices. La solución no fue trivial. Estos son los tiempos para una matriz de 5790x5790 en mi sistema 4 core / 8 HT Ivy Bridge (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%

El nuevo método es 25 veces más rápido para el doble y 40 veces más rápido para el simple. La eficiencia es de aproximadamente 35-40% de los FLOPS de máximo ahora. Con la multiplicación de matrices obtengo hasta 70% con AVX en mi propio código. No sé qué esperar de la descomposición de Cholesky. El algoritmo es parcialmente serial (al calcular el bloque diagonal, llamado triangle en mi código a continuación) a diferencia de la multiplicación de matrices.

Actualización: estoy dentro de un factor para 2 de MKL. No sé si debería estar orgulloso de eso o avergonzado por eso, pero aparentemente mi código aún se puede mejorar significativamente. Encontré una tesis de doctorado sobre esto que muestra que mi algoritmo de bloque es una solución común, así que pude reinventar la rueda.

Utilizo azulejos 32x32 para azulejos dobles y 64x64 para flotar. También reordené la memoria para que cada mosaico sea contiguo y sea su transposición. Definí una nueva función de producción de matriz. La multiplicación de la matriz se define como:

C_i,j = A_i,k * B_k,j //sum over k

Me di cuenta de que en el algoritmo de Cholesky hay algo muy similar

C_j,i = A_i,k * B_j,k //sum over k

Al escribir la transposición de las fichas, pude usar mi función optimizada para la multiplicación de matrices casi aquí (solo tuve que cambiar una línea de código). Aquí está la función principal:

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); 

Aquí están las otras funciones. Tengo seis funciones de producto para SSE2, AVX y FMA, cada una con versión doble y flotante. Solo muestro el de AVX y el doble aquí:

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);
}
}