Ich habe die komprimierte Sparse Column (csc) -Darstellung der nxn-Dreiecksmatrix A mit Nullen auf der Hauptdiagonale und möchte nach b in auflösen

(A + I)' * x = b

Dies ist die Routine, die ich habe, um dies zu berechnen:

void backsolve(const int*__restrict__ Lp,
               const int*__restrict__ Li,
               const double*__restrict__ Lx,
               const int n,
               double*__restrict__ x) {
  for (int i=n-1; i>=0; --i) {
      for (int j=Lp[i]; j<Lp[i+1]; ++j) {
          x[i] -= Lx[j] * x[Li[j]];
      }
  }
}

Somit bwird über das Argument übergeben xund von der Lösung überschrieben. Lp, Li, LxSind jeweils die Zeile, Indizes und die Datenzeiger in dem Standard csc Darstellung schwach besetzte Matrizen. Diese Funktion ist der oberste Hotspot im Programm mit der Zeile

x[i] -= Lx[j] * x[Li[j]];

der größte Teil der Zeit verbracht. Kompilieren mit gcc-8.3 -O3 -mfma -mavx -mavx512fgibt

backsolve(int const*, int const*, double const*, int, double*):
        lea     eax, [rcx-1]
        movsx   r11, eax
        lea     r9, [r8+r11*8]
        test    eax, eax
        js      .L9
.L5:
        movsx   rax, DWORD PTR [rdi+r11*4]
        mov     r10d, DWORD PTR [rdi+4+r11*4]
        cmp     eax, r10d
        jge     .L6
        vmovsd  xmm0, QWORD PTR [r9]
.L7:
        movsx   rcx, DWORD PTR [rsi+rax*4]
        vmovsd  xmm1, QWORD PTR [rdx+rax*8]
        add     rax, 1
        vfnmadd231sd    xmm0, xmm1, QWORD PTR [r8+rcx*8]
        vmovsd  QWORD PTR [r9], xmm0
        cmp     r10d, eax
        jg      .L7
.L6:
        sub     r11, 1
        sub     r9, 8
        test    r11d, r11d
        jns     .L5
        ret
.L9:
        ret

Laut vtune,

vmovsd  QWORD PTR [r9], xmm0

ist der langsamste Teil. Ich habe fast keine Erfahrung mit der Montage und weiß nicht, wie ich diesen Vorgang weiter diagnostizieren oder optimieren kann. Ich habe versucht, mit verschiedenen Flags zu kompilieren, um SSE, FMA usw. zu aktivieren / deaktivieren, aber nichts hat funktioniert.

Prozessor: Xeon Skylake

Frage Was kann ich tun, um diese Funktion zu optimieren?

Dies sollte ziemlich stark vom genauen Sparsity-Muster der Matrix und der verwendeten Plattform abhängen. Ich habe einige Dinge mit gcc 8.3.0und Compiler-Flags -O3 -march=native(die sich -march=skylakeauf meiner CPU befinden) im unteren Dreieck dieser Matrix der Dimension 3006 mit Einträgen ungleich Null von 19554 getestet. Hoffentlich liegt dies etwas in der Nähe Ihres Setups, aber ich hoffe auf jeden Fall, dass diese Ihnen eine Vorstellung davon geben können, wo Sie anfangen sollen.

Für das Timing habe ich Google / Benchmark mit dieser Quelldatei verwendet . Es definiert, benchBacksolveBaselinewelche Benchmarks die in der Frage angegebene Implementierung und benchBacksolveOptimizedwelche Benchmarks die vorgeschlagenen "optimierten" Implementierungen sind. Es gibt auch benchFillRhsseparate Benchmarks für die Funktion, die in beiden verwendet wird, um einige nicht vollständig triviale Werte für die rechte Seite zu generieren. Um die Zeit der "reinen" Backsolves zu erhalten, sollte die benötigte Zeit benchFillRhsabgezogen werden.

1. Streng rückwärts iterieren

Die äußere Schleife in Ihrer Implementierung durchläuft die Spalten rückwärts, während die innere Schleife die aktuelle Spalte vorwärts durchläuft. Es scheint konsistenter zu sein, jede Spalte auch rückwärts zu durchlaufen:

for (int i=n-1; i>=0; --i) {
    for (int j=Lp[i+1]-1; j>=Lp[i]; --j) {
        x[i] -= Lx[j] * x[Li[j]];
    }
}

Dies ändert kaum die Assembly ( https://godbolt.org/z/CBZAT5 ), aber die Benchmark-Timings zeigen eine messbare Verbesserung:

------------------------------------------------------------------
Benchmark                        Time             CPU   Iterations
------------------------------------------------------------------
benchFillRhs                  2737 ns         2734 ns      5120000
benchBacksolveBaseline       17412 ns        17421 ns       829630
benchBacksolveOptimized      16046 ns        16040 ns       853333

Ich gehe davon aus, dass dies durch einen irgendwie vorhersehbareren Cache-Zugriff verursacht wird, aber ich habe mich nicht weiter damit befasst.

2. Weniger Lasten / Speicher in der inneren Schleife

Da A ein unteres Dreieck ist, haben wir i < Li[j]. Daher wissen wir, dass x[Li[j]]sich dies aufgrund der Änderungen x[i]in der inneren Schleife nicht ändern wird. Wir können dieses Wissen mithilfe einer temporären Variablen in unsere Implementierung einbringen:

for (int i=n-1; i>=0; --i) {
    double xi_temp = x[i];
    for (int j=Lp[i+1]-1; j>=Lp[i]; --j) {
        xi_temp -= Lx[j] * x[Li[j]];
    }
    x[i] = xi_temp;
}

Dadurch wird gcc 8.3.0der Speicher von innerhalb der inneren Schleife direkt nach seinem Ende in den Speicher verschoben ( https://godbolt.org/z/vM4gPD ). Der Benchmark für die Testmatrix auf meinem System zeigt eine kleine Verbesserung:

------------------------------------------------------------------
Benchmark                        Time             CPU   Iterations
------------------------------------------------------------------
benchFillRhs                  2737 ns         2740 ns      5120000
benchBacksolveBaseline       17410 ns        17418 ns       814545
benchBacksolveOptimized      15155 ns        15147 ns       887129

3. Rollen Sie die Schleife ab

Während clangbereits nach der ersten vorgeschlagenen Codeänderung mit dem Abrollen der Schleife begonnen wird, ist dies gcc 8.3.0immer noch nicht der Fall . Versuchen wir es also, indem wir zusätzlich bestehen -funroll-loops.

------------------------------------------------------------------
Benchmark                        Time             CPU   Iterations
------------------------------------------------------------------
benchFillRhs                  2733 ns         2734 ns      5120000
benchBacksolveBaseline       15079 ns        15081 ns       953191
benchBacksolveOptimized      14392 ns        14385 ns       963441

Beachten Sie, dass sich auch die Baseline verbessert, da die Schleife in dieser Implementierung ebenfalls abgewickelt wird. Unsere optimierte Version profitiert auch ein wenig vom Abrollen der Schleife, aber vielleicht nicht so sehr, wie es uns gefallen hat. Wenn man sich die generierte Assembly ( https://godbolt.org/z/_LJC5f ) ansieht , scheint es, als gccwäre sie mit 8 Abrollvorgängen ein wenig weit gegangen. Für mein Setup kann ich mit nur einem einfachen manuellen Abrollen tatsächlich etwas besser abschneiden. Lassen Sie also die Flagge -funroll-loopswieder fallen und implementieren Sie das Abrollen mit so etwas:

for (int i=n-1; i>=0; --i) {
    const int col_begin = Lp[i];
    const int col_end = Lp[i+1];
    const bool is_col_nnz_odd = (col_end - col_begin) & 1;
    double xi_temp = x[i];
    int j = col_end - 1;
    if (is_col_nnz_odd) {
        xi_temp -= Lx[j] * x[Li[j]];
        --j;
    }
    for (; j >= col_begin; j -= 2) {
        xi_temp -= Lx[j - 0] * x[Li[j - 0]] +
                   Lx[j - 1] * x[Li[j - 1]];
    }
    x[i] = xi_temp;
}

Damit messe ich:

------------------------------------------------------------------
Benchmark                        Time             CPU   Iterations
------------------------------------------------------------------
benchFillRhs                  2728 ns         2729 ns      5090909
benchBacksolveBaseline       17451 ns        17449 ns       822018
benchBacksolveOptimized      13440 ns        13443 ns      1018182

Andere Algorithmen

Alle diese Versionen verwenden immer noch dieselbe einfache Implementierung der Rückwärtslösung für die Struktur mit geringer Dichte. Inhärent kann das Arbeiten mit solchen Strukturen mit geringer Dichte erhebliche Probleme mit dem Speicherverkehr haben. Zumindest für Matrixfaktorisierungen gibt es ausgefeiltere Methoden, die mit dichten Submatrizen arbeiten, die aus der spärlichen Struktur zusammengesetzt sind. Beispiele sind übernodale und multifrontale Methoden. Ich bin ein bisschen verschwommen, aber ich denke, dass solche Methoden diese Idee auch auf das Layout anwenden und dichte Matrixoperationen für niedrigere dreieckige Rückwärtslösungen verwenden werden (zum Beispiel für Faktorisierungen vom Cholesky-Typ). Es könnte sich also lohnen, diese Art von Methoden zu untersuchen, wenn Sie nicht gezwungen sind, sich an die einfache Methode zu halten, die direkt auf der spärlichen Struktur funktioniert. Siehe zum Beispiel diese Umfrage von Davis.

Sie können einige Zyklen rasieren, indem Sie unsignedanstelle intder Indextypen verwenden, was >= 0ohnehin sein muss:

void backsolve(const unsigned * __restrict__ Lp,
               const unsigned * __restrict__ Li,
               const double * __restrict__ Lx,
               const unsigned n,
               double * __restrict__ x) {
    for (unsigned i = n; i-- > 0; ) {
        for (unsigned j = Lp[i]; j < Lp[i + 1]; ++j) {
            x[i] -= Lx[j] * x[Li[j]];
        }
    }
}

Beim Kompilieren mit dem Compiler-Explorer von Godbolt wird etwas anderer Code für den Innerloop angezeigt, wodurch die CPU-Pipeline möglicherweise besser genutzt wird. Ich kann nicht testen, aber Sie könnten es versuchen.

Hier ist der generierte Code für die innere Schleife:

.L8:
        mov     rax, rcx
.L5:
        mov     ecx, DWORD PTR [r10+rax*4]
        vmovsd  xmm1, QWORD PTR [r11+rax*8]
        vfnmadd231sd    xmm0, xmm1, QWORD PTR [r8+rcx*8]
        lea     rcx, [rax+1]
        vmovsd  QWORD PTR [r9], xmm0
        cmp     rdi, rax
        jne     .L8