Tengo la representación de columna dispersa comprimida (csc) de la matriz triangular inferior nxn A con ceros en la diagonal principal, y me gustaría resolver b en

(A + I)' * x = b

Esta es la rutina que tengo para calcular esto:

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

Por lo tanto, bse pasa a través del argumento xy se sobrescribe con la solución. Lp, Li, LxSon, respectivamente, las hileras, índices, y los punteros de datos en la representación csc estándar de matrices dispersas. Esta función es el punto de acceso superior en el programa, con la línea

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

siendo la mayor parte del tiempo dedicado. Compilando con gcc-8.3 -O3 -mfma -mavx -mavx512fda

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

Según vtune,

vmovsd  QWORD PTR [r9], xmm0

Es la parte más lenta. Casi no tengo experiencia con el ensamblaje, y no sé cómo diagnosticar u optimizar aún más esta operación. He intentado compilar con diferentes indicadores para habilitar / deshabilitar SSE, FMA, etc., pero nada ha funcionado.

Procesador: Xeon Skylake

Pregunta ¿Qué puedo hacer para optimizar esta función?

Esto debería depender bastante del patrón de dispersión exacto de la matriz y la plataforma que se utiliza. Probé algunas cosas con gcc 8.3.0los indicadores del compilador -O3 -march=native(que está -march=skylakeen mi CPU) en el triángulo inferior de esta matriz de dimensión 3006 con 19554 entradas distintas de cero. Espero que esto esté algo cerca de su configuración, pero en cualquier caso espero que esto pueda darle una idea de dónde comenzar.

Para el tiempo utilicé google / benchmark con este archivo fuente . Define benchBacksolveBaselinequé puntos de referencia la implementación dada en la pregunta y benchBacksolveOptimizedqué puntos de referencia las implementaciones "optimizadas" propuestas. También hay benchFillRhsque compara por separado la función que se utiliza en ambos para generar algunos valores no completamente triviales para el lado derecho. Para obtener el tiempo de las soluciones "puras", benchFillRhsse debe restar el tiempo que lleva.

1. Iterando estrictamente al revés

El bucle externo en su implementación itera a través de las columnas hacia atrás, mientras que el bucle interno itera a través de la columna actual hacia adelante. Parece que sería más consistente iterar a través de cada columna también hacia atrás:

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

Esto apenas cambia el ensamblaje ( https://godbolt.org/z/CBZAT5 ), pero los tiempos de referencia muestran una mejora considerable:

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

Supongo que esto es causado por un acceso a la caché más predecible, pero no lo examiné mucho más.

2. Menos cargas / tiendas en bucle interno

Como A es triangular inferior, tenemos i < Li[j]. Por lo tanto, sabemos que x[Li[j]]no cambiará debido a los cambios x[i]en el bucle interno. Podemos poner este conocimiento en nuestra implementación mediante el uso de una variable temporal:

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

Esto hace que gcc 8.3.0mover la tienda a la memoria desde el interior del bucle interno directamente después de su finalización ( https://godbolt.org/z/vM4gPD ). El punto de referencia para la matriz de prueba en mi sistema muestra una pequeña mejora:

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

3. Desenrollar el bucle

Si bien clangya comienza a desenrollar el bucle después del primer cambio de código sugerido, gcc 8.3.0aún no lo ha hecho. Así que vamos a intentarlo pasando adicionalmente -funroll-loops.

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

Tenga en cuenta que la línea de base también mejora, ya que el ciclo en esa implementación también se desenrolla. Nuestra versión optimizada también se beneficia un poco del desenrollado del bucle, pero quizás no tanto como nos hubiera gustado. Mirando el ensamblaje generado ( https://godbolt.org/z/_LJC5f ), parece que gccpodría haber ido un poco lejos con 8 desenrollamientos. Para mi configuración, de hecho, puedo hacerlo un poco mejor con solo un desenrollado manual simple. Así que suelte la bandera -funroll-loopsnuevamente e implemente el desenrollado con algo como esto:

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

Con eso mido:

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

Otros algoritmos

Todas estas versiones todavía usan la misma implementación simple de la resolución hacia atrás en la estructura de matriz dispersa. Inherentemente, operar en estructuras de matriz dispersas como estas puede tener problemas significativos con el tráfico de memoria. Al menos para las factorizaciones matriciales, existen métodos más sofisticados que operan en submatrices densas que se ensamblan a partir de la estructura dispersa. Los ejemplos son métodos supernodales y multifrontales. Estoy un poco confuso en esto, pero creo que tales métodos también aplicarán esta idea al diseño y usarán operaciones de matriz densas para soluciones triangulares inferiores hacia atrás (por ejemplo, para factorizaciones de tipo Cholesky). Por lo tanto, podría valer la pena analizar ese tipo de métodos, si no se ve obligado a seguir el método simple que funciona directamente en la estructura dispersa. Ver por ejemplo esta encuestapor Davis.

Puede reducir algunos ciclos utilizando en unsignedlugar de intlos tipos de índice, que deben ser de >= 0todos modos:

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

Compilando con compilador explorador de Godbolt de código se muestra ligeramente diferente para el Innerloop, lo que podría hacer un mejor uso de la tubería de la CPU. No puedo probar, pero puedes intentarlo.

Aquí está el código generado para el bucle interno:

.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