Tengo un proyecto donde resolvemos el inverso de matrices densas definidas positivas grandes (más de 3000x3000) usandoDescomposición Cholesky. El proyecto está en Java y usamos el CERNBiblioteca 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). Comencé con el código de descomposición Cholesky en C deCódigo Rosetta.

Lo que noté es que los valores en una columna, excepto el elemento diagonal, son independientes. Entonces 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 ligeramente más lenta que la de RosettaCodepero 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). El uso de la 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 enhttp://coliru.stacked-crooked.com/a/6f5750c20d456da9

Inicialmente pensé que el intercambio falso 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?

¿Quizás un enfoque basado en tareas es mejor en general? En las diapositivas 7-9 ahttp://courses.engr.illinois.edu/cs554/fa2013/notes/07_cholesky.pdf Hay un ejemplo de descomposición paralela de Cholesky usando "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, p. con tareas?

Editar: según lo solicitado aquí es la función AVX que solía calculars. 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;
}