Chcę, aby zoptymalizować tę krótką pętlę

Question

Chciałbym, aby zoptymalizować tę prostą pętlę:

unsigned int i; 
while(j-- != 0){ //j is an unsigned int with a start value of about N = 36.000.000 
    float sub = 0; 
    i=1; 
    unsigned int c = j+s[1]; 
    while(c < N) { 
     sub += d[i][j]*x[c];//d[][] and x[] are arrays of float 
     i++; 
     c = j+s[i];// s[] is an array of unsigned int with 6 entries. 
    } 
    x[j] -= sub;      // only one memory-write per j 
} 

Pętla ma czas realizacji wynosi około jedną sekundę z 4000 MHz AMD Bulldozer. Myślałem o SIMD i OpenMP (których zwykle używam, aby uzyskać większą prędkość), ale ta pętla jest rekurencyjna.

Wszelkie sugestie?

Answer 1

zgadzam się z transpozycji do lepszego buforowania (ale zobacz moje komentarze na tym na końcu), a tam jest więcej do zrobienia, więc zobaczymy, co możemy zrobić z pełną funkcją ...

funkcja Original, odsyłające (z pewnymi sprzątania dla mojego zdrowia psychicznego):

void MultiDiagonalSymmetricMatrix::CholeskyBackSolve(float *x, float *b){ 
    //We want to solve L D Lt x = b where D is a diagonal matrix described by Diagonals[0] and L is a unit lower triagular matrix described by the rest of the diagonals. 
    //Let D Lt x = y. Then, first solve L y = b. 

    float *y = new float[n]; 
    float **d = IncompleteCholeskyFactorization->Diagonals; 
    unsigned int *s = IncompleteCholeskyFactorization->StartRows; 
    unsigned int M = IncompleteCholeskyFactorization->m; 
    unsigned int N = IncompleteCholeskyFactorization->n; 
    unsigned int i, j; 
    for(j = 0; j != N; j++){ 
     float sub = 0; 
     for(i = 1; i != M; i++){ 
      int c = (int)j - (int)s[i]; 
      if(c < 0) break; 
      if(c==j) { 
       sub += d[i][c]*b[c]; 
      } else { 
       sub += d[i][c]*y[c]; 
      } 
     } 
     y[j] = b[j] - sub; 
    } 

    //Now, solve x from D Lt x = y -> Lt x = D^-1 y 
    // Took this one out of the while, so it can be parallelized now, which speeds up, because division is expensive 
#pragma omp parallel for 
    for(j = 0; j < N; j++){ 
     x[j] = y[j]/d[0][j]; 
    } 

    while(j-- != 0){ 
     float sub = 0; 
     for(i = 1; i != M; i++){ 
      if(j + s[i] >= N) break; 
      sub += d[i][j]*x[j + s[i]]; 
     } 
     x[j] -= sub; 
    } 
    delete[] y; 
} 

ze względu na komentarz na temat podziału równolegle dając impuls prędkości (mimo że tylko o (N)), jestem przy założeniu, że sama funkcja jest wywoływana wiele . Dlaczego więc przydzielać pamięć? Po prostu zaznacz x jako __restrict__ i zmień y na x wszędzie (__restrict__ jest rozszerzeniem GCC, pobranym z C99.Możesz użyć do tego celu define.Może biblioteka ma już jeden).

Podobnie, chociaż domyślam się, że nie można zmienić podpisu, można sprawić, że funkcja będzie przyjmować tylko jeden parametr i modyfikować go. b nigdy nie jest używany, gdy ustawiono x lub y. Oznaczałoby to również, że możesz pozbyć się gałęzi w pierwszej pętli, która działa ~ N * M razy. Użyj memcpy na początku, jeśli musisz mieć 2 parametry.

A dlaczego jest d tablica wskaźników? Czy to musi być? Wydaje się to zbyt głębokie w oryginalnym kodzie, więc nie będę go dotykał, ale jeśli istnieje możliwość spłaszczenia zapisanej tablicy, będzie to zwiększenie prędkości, nawet jeśli nie możesz jej transponować (pomnóż, dodaj, dereferencja jest szybsza niż dereferencja, dodawanie, dereferencja).

więc nowy kod:

void MultiDiagonalSymmetricMatrix::CholeskyBackSolve(float *__restrict__ x){ 
    // comments removed so that suggestions are more visible. Don't remove them in the real code! 
    // these definitions got long. Feel free to remove const; it does nothing for the optimiser 
    const float *const __restrict__ *const __restrict__ d = IncompleteCholeskyFactorization->Diagonals; 
    const unsigned int *const __restrict__ s = IncompleteCholeskyFactorization->StartRows; 
    const unsigned int M = IncompleteCholeskyFactorization->m; 
    const unsigned int N = IncompleteCholeskyFactorization->n; 
    unsigned int i; 
    unsigned int j; 
    for(j = 0; j < N; j++){ // don't use != as an optimisation; compilers can do more with < 
     float sub = 0; 
     for(i = 1; i < M && j >= s[i]; i++){ 
      const unsigned int c = j - s[i]; 
      sub += d[i][c]*x[c]; 
     } 
     x[j] -= sub; 
    } 

    // Consider using processor-specific optimisations for this 
#pragma omp parallel for 
    for(j = 0; j < N; j++){ 
     x[j] /= d[0][j]; 
    } 

    for(j = N; (j --) > 0;){ // changed for clarity 
     float sub = 0; 
     for(i = 1; i < M && j + s[i] < N; i++){ 
      sub += d[i][j]*x[j + s[i]]; 
     } 
     x[j] -= sub; 
    } 
} 

Dobrze to wygląda porządniej, a brak alokacji pamięci i zmniejszona rozgałęzienia, jeśli nic innego, to impuls. Jeśli możesz zmienić s na dodatkową wartość UINT_MAX na końcu, możesz usunąć więcej gałęzi (zarówno sprawdzanie i<M, które ponownie uruchamiają ~ N * M razy).

Teraz nie możemy utworzyć więcej pętli równolegle i nie możemy łączyć pętli. Teraz podbicie będzie, jak zasugerowano w drugiej odpowiedzi, przestawić d. Poza tym ... praca konieczna do przestawienia d ma dokładnie takie same problemy z pamięcią podręczną jak praca do wykonania pętli. I potrzebowałaby pamięci przydzielonej. Niedobrze. Jedynymi opcjami do dalszej optymalizacji są: zmiana samej struktury IncompleteCholeskyFactorization->Diagonals, co prawdopodobnie będzie oznaczać wiele zmian, lub znalezienie innego algorytmu, który działa lepiej z danymi w tej kolejności.

Jeśli chcesz pójść dalej, twoje optymalizacje będą musiały wpłynąć na znaczną część kodu (nie jest to złe, chyba że istnieje dobry powód, dla którego Diagonals jest tablicą wskaźników, wydaje się, że może to zrobić przy pomocy refaktor).

Answer 2

że może chcesz transponować macierz D - oznacza przechowywać go w taki sposób, że można wymieniać indeksy - make I indeks zewnętrzna:

sub += d[j][i]*x[c]; 

zamiast

sub += d[i][j]*x[c]; 

Powinno to spowodować lepszą wydajność pamięci podręcznej.

Answer 3

Chcę dać odpowiedź na moje własne pytanie: Zła wydajność została spowodowana przez brakujący konflikt pamięci podręcznej ze względu na fakt, że (przynajmniej) Win7 wyrównuje duże bloki pamięci do tej samej granicy. W moim przypadku, dla wszystkich buforów adresy miały to samo wyrównanie (bufferadress% 4096 był taki sam dla wszystkich buforów), więc przypadają one na ten sam bufor pamięci podręcznej L1. Zmieniłem alokację pamięci, aby wyrównać bufory do różnych granic, aby uniknąć konfliktów w pamięci podręcznej i uzyskać przyspieszenie czynnika 2. Dziękuję za wszystkie odpowiedzi, szczególnie za odpowiedzi od Dave'a!