Ich habe ein Projekt, in dem wir die Inverse von großen (über 3000x3000) positiv definitiven dichten Matrizen mit lösenCholesky-Zersetzung. Das Projekt ist in Java und wir verwenden das CERNColt BLAS Bibliothek. Das Profilieren des Codes zeigt, dass die Cholesky-Zerlegung der Engpass ist.

Ich beschloss, die Cholesky-Zerlegung mit OpenMP zu parallelisieren und als DLL in Java (mit JNA) zu verwenden. Ich habe mit dem Cholesky-Dekompositionscode in C von begonnenRosetta Code.

Was mir aufgefallen ist, dass die Werte in einer Spalte mit Ausnahme des diagonalen Elements unabhängig sind. Also habe ich beschlossen, die diagonalen Elemente seriell und die restlichen Werte der Spalte parallel zu berechnen. Ich habe auch die Reihenfolge der Schleifen vertauscht, sodass die innere Schleife über die Zeilen und die äußere Schleife über die Spalten verläuft. Die Serienversion ist etwas langsamer als die von RosettaCodeAber die parallele Version ist sechsmal schneller als die RosettaCode-Version auf meinem 4-Kern-System (8 HT). Die Verwendung der DLL in Java beschleunigt unsere Ergebnisse ebenfalls um das Sechsfache. Hier ist der Code:

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

Den vollständigen Code zum Testen finden Sie unterhttp://coliru.stacked-crooked.com/a/6f5750c20d456da9

Anfangs dachte ich, dass falsches Teilen ein Problem sein würde, wenn die verbleibenden Elemente einer Spalte im Vergleich zur Anzahl der Threads klein wären, aber das scheint nicht der Fall zu sein. Ich habe es versucht

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

Ich habe keine klaren Beispiele für die Parallelisierung der Choleskey-Zerlegung gefunden. Ich weiß nicht, ob das, was ich getan habe, ideal ist. Funktioniert es beispielsweise auf einem NUMA-System?

Vielleicht ist ein aufgabenbasierter Ansatz im Allgemeinen besser? In den Folien 7-9 umhttp://courses.engr.illinois.edu/cs554/fa2013/notes/07_cholesky.pdf Es gibt ein Beispiel für eine parallele Cholesky-Zerlegung mit "feinkörnigen Aufgaben". Mir ist noch nicht klar, wie ich das umsetzen soll.

Ich habe zwei spezifische und allgemeine Fragen. Haben Sie Vorschläge zur Verbesserung meiner Implementierung von Cholesky Decomposition mit OpenMP? Können Sie eine andere Implementierung von Cholesky Decomposition mit OpenMP vorschlagen, z. mit aufgaben?

Bearbeiten: wie hier angefordert ist die AVX-Funktion, die ich zum Berechnen verwendet habes. Es hat nicht geholfen

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