Eu tenho a representação da coluna esparsa compactada (csc) da matriz triangular inferior nxn A com zeros na diagonal principal e gostaria de resolver b em

(A + I)' * x = b

Esta é a rotina que tenho para calcular isso:

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

Assim, bé passado através do argumento xe é substituído pela solução. Lp, Li, LxSão, respectivamente, a linha, índices e ponteiros de dados na representação csc padrão de matrizes esparsas. Esta função é o principal ponto de acesso do programa, com a linha

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

sendo a maior parte do tempo gasto. Compilar com gcc-8.3 -O3 -mfma -mavx -mavx512fdá

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

De acordo com vtune,

vmovsd  QWORD PTR [r9], xmm0

é a parte mais lenta. Quase não tenho experiência com montagem e não sei como diagnosticar ou otimizar ainda mais esta operação. Eu tentei compilar com diferentes sinalizadores para ativar / desativar SSE, FMA, etc, mas nada funcionou.

Processador: Xeon Skylake

Pergunta O que posso fazer para otimizar esta função?

Isso deve depender um pouco do padrão exato de dispersão da matriz e da plataforma que está sendo usada. Testei algumas coisas com gcc 8.3.0sinalizadores de compilador -O3 -march=native(que estão -march=skylakena minha CPU) no triângulo inferior dessa matriz de dimensão 3006 com entradas 19554 diferentes de zero. Espero que isso esteja um pouco próximo da sua configuração, mas, em qualquer caso, espero que eles possam lhe dar uma idéia de por onde começar.

Por tempo, usei google / benchmark com este arquivo de origem . Ele define benchBacksolveBaselinequais benchmarks a implementação dada na pergunta e benchBacksolveOptimizedquais benchmarks as implementações "otimizadas" propostas. Também há benchFillRhsquais benchmarks separadamente da função usada em ambos para gerar alguns valores não completamente triviais para o lado direito. Para obter o tempo dos backsolves "puros", o tempo que benchFillRhsleva deve ser subtraído.

1. Iterando estritamente para trás

O loop externo em sua implementação percorre as colunas para trás, enquanto o loop interno percorre a coluna atual para a frente. Parece que seria mais consistente iterar através de cada coluna também:

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

Isso mal muda a montagem ( https://godbolt.org/z/CBZAT5 ), mas os tempos de referência mostram uma melhoria mensurável:

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

Suponho que isso seja causado por um acesso de cache de alguma forma mais previsível, mas não o examinei muito mais.

2. Menos cargas / lojas no loop interno

Como A é triangular inferior, temos i < Li[j]. Portanto, sabemos que x[Li[j]]isso não será alterado devido às alterações x[i]no loop interno. Podemos colocar esse conhecimento em nossa implementação usando uma variável temporária:

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

Isso faz com que gcc 8.3.0o armazenamento seja transferido para a memória de dentro do loop interno para diretamente após seu término ( https://godbolt.org/z/vM4gPD ). A referência para a matriz de teste no meu sistema mostra uma pequena melhoria:

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

3. Desenrole o loop

Embora clangjá comece a desenrolar o loop após a primeira alteração de código sugerida, gcc 8.3.0ainda não o fez. Então, vamos tentar, passando 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

Observe que a linha de base também melhora, pois o loop nessa implementação também é desenrolado. Nossa versão otimizada também se beneficia um pouco do desenrolar do loop, mas talvez não tanto quanto gostaríamos. Olhando para o assembly gerado ( https://godbolt.org/z/_LJC5f ), parece que gccpode ter ido um pouco longe com 8 desenrolamentos. Para minha configuração, na verdade, posso melhorar um pouco com apenas um desenrolar manual simples. Então solte a bandeira -funroll-loopsnovamente e implemente o desenrolamento com algo como isto:

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

Com isso eu meço:

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

Outros algoritmos

Todas essas versões ainda usam a mesma implementação simples da solução reversa na estrutura da matriz esparsa. Inerentemente, operar em estruturas de matriz esparsas como estas pode ter problemas significativos com o tráfego de memória. Pelo menos para fatorações matriciais, existem métodos mais sofisticados, que operam em submatrizes densas que são montadas a partir da estrutura esparsa. Exemplos são métodos supernodais e multifrontais. Estou um pouco confuso com isso, mas acho que esses métodos também aplicarão essa idéia ao layout e usarão operações de matriz densa para soluções triangulares inferiores para trás (por exemplo, para fatorações do tipo Cholesky). Portanto, pode valer a pena examinar esse tipo de método, se você não for forçado a seguir o método simples que funciona diretamente na estrutura esparsa. Veja, por exemplo, esta pesquisapor Davis.

Você pode cortar alguns ciclos usando, em unsignedvez de, intos tipos de índice, que devem ser >= 0assim:

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

A compilação com o explorador de compilador da Godbolt mostra um código ligeiramente diferente para o loop interno, potencialmente fazendo melhor uso do pipeline da CPU. Não posso testar, mas você pode tentar.

Aqui está o código gerado para o loop 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