Como o BLAS consegue um desempenho tão extremo?

108

Por curiosidade, decidi comparar a minha própria função de multiplicação de matrizes versus a implementação do BLAS ... Fiquei, para dizer, o menos surpreendido com o resultado:

Implementação personalizada, 10 testes de multiplicação de matriz 1000x1000:

Took: 15.76542 seconds.

Implementação do BLAS, 10 tentativas de multiplicação de matriz 1000x1000:

Took: 1.32432 seconds.

Isso está usando números de ponto flutuante de precisão única.

Minha implementação:

template<class ValT>
void mmult(const ValT* A, int ADim1, int ADim2, const ValT* B, int BDim1, int BDim2, ValT* C)
{
    if ( ADim2!=BDim1 )
        throw std::runtime_error("Error sizes off");

    memset((void*)C,0,sizeof(ValT)*ADim1*BDim2);
    int cc2,cc1,cr1;
    for ( cc2=0 ; cc2<BDim2 ; ++cc2 )
        for ( cc1=0 ; cc1<ADim2 ; ++cc1 )
            for ( cr1=0 ; cr1<ADim1 ; ++cr1 )
                C[cc2*ADim2+cr1] += A[cc1*ADim1+cr1]*B[cc2*BDim1+cc1];
}

Eu tenho duas perguntas:

  1. Dado que uma multiplicação matriz-matriz diz: nxm * mxn requer n * n * m multiplicações, portanto, no caso acima de 1000 ^ 3 ou 1e9 operações. Como é possível no meu processador de 2.6 GHz o BLAS fazer 10 * 1e9 operações em 1,32 segundos? Mesmo que as multiplicações fossem uma única operação e não houvesse mais nada sendo feito, isso levaria cerca de 4 segundos.
  2. Por que minha implementação é muito mais lenta?
DeusAduro
fonte
17
O BLAS foi otimizado de um lado e do outro por um especialista na área. Suponho que ele esteja aproveitando a unidade de ponto flutuante SIMD em seu chip e fazendo vários truques para melhorar o comportamento do cache também ...
dmckee --- ex-moderador gatinho
3
Ainda assim, como você faz operações 1E10 em um processador de 2,63E9 ciclos / segundo em 1,3 segundos?
DeusAduro
9
Várias unidades de execução, pipe-lining e Single Instruction Multiple Data ((SIMD), o que significa fazer a mesma operação em mais de um par de operandos ao mesmo tempo). Alguns compiladores podem direcionar as unidades SIMD em chips comuns, mas você quase sempre precisa ativá-los explicitamente, e isso ajuda a saber como tudo funciona ( en.wikipedia.org/wiki/SIMD ). O seguro contra falhas de cache é quase certamente a parte difícil.
dmckee --- ex-moderador gatinho de
13
A suposição está errada. Existem algoritmos melhores conhecidos, consulte Wikipedia.
MSalters de
2
@DeusAduro: Em minha resposta para Como escrever um produto de matriz de matriz que pode competir com Eigen? Publiquei um pequeno exemplo sobre como implementar um produto matriz-matriz eficiente de cache.
Michael Lehn

Respostas:

141

Um bom ponto de partida é o grande livro The Science of Programming Matrix Computations, de Robert A. van de Geijn e Enrique S. Quintana-Ortí. Eles fornecem uma versão de download grátis.

O BLAS está dividido em três níveis:

  • O nível 1 define um conjunto de funções de álgebra linear que operam apenas em vetores. Essas funções se beneficiam da vetorização (por exemplo, do uso de SSE).

  • As funções de nível 2 são operações de vetor de matriz, por exemplo, algum produto de vetor de matriz. Essas funções podem ser implementadas em termos de funções de Nível 1. No entanto, você pode aumentar o desempenho dessas funções se puder fornecer uma implementação dedicada que faça uso de alguma arquitetura de multiprocessador com memória compartilhada.

  • As funções de nível 3 são operações como o produto matriz-matriz. Novamente, você pode implementá-los em termos de funções de Nível2. Mas as funções de Nível 3 realizam operações O (N ^ 3) em dados O (N ^ 2). Portanto, se sua plataforma tem uma hierarquia de cache, você pode aumentar o desempenho se fornecer uma implementação dedicada otimizada para cache / amigável para cache . Isso é muito bem descrito no livro. O principal impulso das funções do Nível 3 vem da otimização do cache. Esse aumento excede significativamente o segundo aumento do paralelismo e outras otimizações de hardware.

A propósito, a maioria (ou mesmo todas) as implementações de BLAS de alto desempenho NÃO são implementadas em Fortran. ATLAS é implementado em C. GotoBLAS / OpenBLAS é implementado em C e suas partes críticas de desempenho em Assembler. Apenas a implementação de referência do BLAS é implementada em Fortran. No entanto, todas essas implementações BLAS fornecem uma interface Fortran de forma que ele pode ser vinculado ao LAPACK (LAPACK ganha todo o seu desempenho do BLAS).

Compiladores otimizados desempenham um papel menor a este respeito (e para GotoBLAS / OpenBLAS o compilador não importa de todo).

A implementação de IMHO no BLAS usa algoritmos como o algoritmo Coppersmith – Winograd ou o algoritmo Strassen. Não tenho certeza sobre o motivo, mas este é o meu palpite:

  • Talvez não seja possível fornecer uma implementação otimizada de cache desses algoritmos (ou seja, você perderia mais do que ganharia)
  • Esses algoritmos não são numericamente estáveis. Como o BLAS é o kernel computacional do LAPACK, isso é impossível.

Editar / Atualizar:

O artigo novo e inovador para este tópico são os artigos BLIS . Eles são excepcionalmente bem escritos. Para minha palestra "Noções básicas de software para computação de alto desempenho", implementei o produto matriz-matriz seguindo seu artigo. Na verdade, implementei várias variantes do produto matriz-matriz. As variantes mais simples são inteiramente escritas em C simples e têm menos de 450 linhas de código. Todas as outras variantes apenas otimizam os loops

    for (l=0; l<MR*NR; ++l) {
        AB[l] = 0;
    }
    for (l=0; l<kc; ++l) {
        for (j=0; j<NR; ++j) {
            for (i=0; i<MR; ++i) {
                AB[i+j*MR] += A[i]*B[j];
            }
        }
        A += MR;
        B += NR;
    }

O desempenho geral do produto matriz-matriz depende apenas desses loops. Cerca de 99,9% do tempo é gasto aqui. Nas outras variantes, usei intrínsecos e código assembler para melhorar o desempenho. Você pode ver o tutorial passando por todas as variantes aqui:

ulmBLAS: Tutorial sobre GEMM (produto Matrix-Matrix)

Junto com os papéis do BLIS, torna-se bastante fácil entender como bibliotecas como Intel MKL podem obter tal desempenho. E por que não importa se você usa armazenamento principal de linha ou coluna!

Os benchmarks finais estão aqui (chamamos nosso projeto ulmBLAS):

Benchmarks para ulmBLAS, BLIS, MKL, openBLAS e Eigen

Outra edição / atualização:

Eu também escrevi um tutorial sobre como o BLAS é usado para problemas de álgebra linear numérica, como resolver um sistema de equações lineares:

Fatoração LU de alto desempenho

(Esta fatoração LU é, por exemplo, usada pelo Matlab para resolver um sistema de equações lineares.)

Espero encontrar tempo para estender o tutorial para descrever e demonstrar como realizar uma implementação paralela altamente escalonável da fatoração LU como no PLASMA .

Ok, aqui está: Codificando uma Fatoração de LU paralela otimizada de cache

PS: Eu também fiz alguns experimentos para melhorar o desempenho do uBLAS. Na verdade, é muito simples aumentar (sim, brincar com as palavras :)) o desempenho do uBLAS:

Experiências em uBLAS .

Aqui está um projeto semelhante com BLAZE :

Experiências em BLAZE .

Michael Lehn
fonte
3
Novo link para “Benchmarks para ulmBLAS, BLIS, MKL, openBLAS e Eigen”: apfel.mathematik.uni-ulm.de/~lehn/ulmBLAS/#toc3
Ahmed Fasih
Acontece que o ESSL da IBM usa uma variação do algoritmo Strassen - ibm.com/support/knowledgecenter/en/SSFHY8/essl_welcome.html
ben-albrecht
2
a maioria dos links estão mortos
Aurélien Pierre
Um PDF do TSoPMC pode ser encontrado na página do autor, em cs.utexas.edu/users/rvdg/tmp/TSoPMC.pdf
Alex Shpilkin
Embora o algoritmo Coppersmith-Winograd tenha uma boa complexidade de tempo no papel, a notação Big O esconde uma constante muito grande, então ela só começa a se tornar viável para matrizes ridiculamente grandes.
DiehardThe Tryhard
26

Portanto, em primeiro lugar, o BLAS é apenas uma interface com cerca de 50 funções. Existem muitas implementações concorrentes da interface.

Em primeiro lugar, mencionarei coisas que não estão relacionadas:

  • Fortran vs C, não faz diferença
  • Algoritmos de matriz avançados, como Strassen, implementações não os usam porque não ajudam na prática

A maioria das implementações divide cada operação em matriz de pequena dimensão ou operações de vetor de maneira mais ou menos óbvia. Por exemplo, uma grande multiplicação de matriz de 1000x1000 pode ser quebrada em uma sequência de multiplicações de matriz de 50x50.

Essas operações de pequena dimensão de tamanho fixo (chamadas de kernels) são codificadas no código de montagem específico da CPU usando vários recursos da CPU de seu destino:

  • Instruções de estilo SIMD
  • Paralelismo de nível de instrução
  • Percepção de cache

Além disso, esses kernels podem ser executados em paralelo entre si, usando vários threads (núcleos de CPU), no padrão de design de redução de mapa típico.

Dê uma olhada no ATLAS, que é a implementação do BLAS de código aberto mais comumente usada. Ele tem muitos kernels concorrentes diferentes e, durante o processo de construção da biblioteca ATLAS, ele faz uma competição entre eles (alguns até são parametrizados, portanto, o mesmo kernel pode ter configurações diferentes). Ele tenta diferentes configurações e, a seguir, seleciona a melhor para o sistema de destino específico.

(Dica: é por isso que se você estiver usando ATLAS, é melhor construir e ajustar a biblioteca manualmente para sua máquina específica do que usar uma pré-construída.)

Andrew Tomazos
fonte
ATLAS não é mais a implementação BLAS de código aberto mais comumente usada. Foi superado por OpenBLAS (um fork do GotoBLAS) e BLIS (uma refatoração do GotoBLAS).
Robert van de Geijn
1
@ ulaff.net: Talvez. Isso foi escrito há 6 anos. Eu acho que a implementação BLAS mais rápida atualmente (na Intel, é claro) é Intel MKL, mas não é open source.
Andrew Tomazos
14

Primeiro, existem algoritmos mais eficientes para multiplicação de matrizes do que aquele que você está usando.

Em segundo lugar, sua CPU pode fazer muito mais do que uma instrução por vez.

Sua CPU executa 3-4 instruções por ciclo e, se as unidades SIMD forem usadas, cada instrução processa 4 flutuadores ou 2 duplos. (é claro que este número também não é preciso, já que a CPU normalmente pode processar apenas uma instrução SIMD por ciclo)

Terceiro, seu código está longe de ser ideal:

  • Você está usando ponteiros brutos, o que significa que o compilador deve assumir que eles podem ser apelidos. Existem palavras-chave específicas do compilador ou sinalizadores que você pode especificar para dizer ao compilador que eles não têm apelidos. Como alternativa, você deve usar outros tipos de ponteiros não processados, que cuidam do problema.
  • Você está alterando o cache realizando uma travessia ingênua de cada linha / coluna das matrizes de entrada. Você pode usar o bloqueio para realizar o máximo de trabalho possível em um bloco menor da matriz, que cabe no cache da CPU, antes de passar para o próximo bloco.
  • Para tarefas puramente numéricas, o Fortran é praticamente imbatível, e C ++ exige muita persuasão para atingir uma velocidade semelhante. Isso pode ser feito e há algumas bibliotecas que o demonstram (normalmente usando modelos de expressão), mas não é trivial e não acontece simplesmente .
Jalf
fonte
Obrigado, adicionei o código correto de restrição de acordo com a sugestão de Justicle, não vi muitas melhorias, gostei da ideia em bloco. Por curiosidade, sem saber o tamanho do cache da CPU, como seria um código ideal correto?
DeusAduro
2
Você não. Para obter o código ideal, você precisa saber o tamanho do cache da CPU. Obviamente, a desvantagem disso é que você está efetivamente codificando seu código para melhor desempenho em uma família de CPUs.
jalf de
2
Pelo menos o loop interno evita cargas com passos largos. Parece que isso foi escrito para uma matriz que já está sendo transposta. É por isso que é "apenas" uma ordem de magnitude mais lento que o BLAS! Mas sim, ainda está se debatendo por causa da falta de bloqueio de cache. Tem certeza de que o Fortran ajudaria muito? Acho que tudo o que você ganharia aqui é que restrict(sem aliasing) é o padrão, ao contrário de C / C ++. (E, infelizmente, ISO C ++ não tem uma restrictpalavra - chave, então você deve usar __restrict__em compiladores que a fornecem como uma extensão).
Peter Cordes,
11

Não sei especificamente sobre a implementação do BLAS, mas existem alogoritmos mais eficientes para a Multiplicação de Matriz que tem complexidade melhor do que O (n3). Um bem conhecido é Strassen Algorithm

softveda
fonte
8
O Algoritmo de Strassen não é usado em números por dois motivos: 1) Não é estável. 2) Você economiza alguns cálculos, mas isso vem com o preço que você pode explorar hierarquias de cache. Na prática, você até perde desempenho.
Michael Lehn
4
Para a implementação prática do Algoritmo Strassen fortemente construído sobre o código-fonte da biblioteca BLAS, há uma publicação recente: " Strassen Algorithm Reloaded " em SC16, que atinge um desempenho superior ao BLAS, mesmo para o tamanho do problema 1000x1000.
Jianyu Huang
4

A maioria dos argumentos para a segunda questão - montador, divisão em blocos etc. (mas não menos do que N ^ 3 algoritmos, eles são realmente superdesenvolvidos) - desempenham um papel. Mas a baixa velocidade do seu algoritmo é causada essencialmente pelo tamanho da matriz e pelo infeliz arranjo dos três loops aninhados. Suas matrizes são tão grandes que não cabem imediatamente na memória cache. Você pode reorganizar os loops de forma que o máximo possível seja feito em uma linha no cache, reduzindo drasticamente as atualizações do cache (BTW, a divisão em pequenos blocos tem um efeito analógico, melhor se os loops sobre os blocos forem organizados de forma semelhante). Segue-se uma implementação de modelo para matrizes quadradas. No meu computador, seu consumo de tempo foi de cerca de 1:10 em comparação com a implementação padrão (como a sua). Em outras palavras: nunca programe uma multiplicação de matriz ao longo do "

    void vector(int m, double ** a, double ** b, double ** c) {
      int i, j, k;
      for (i=0; i<m; i++) {
        double * ci = c[i];
        for (k=0; k<m; k++) ci[k] = 0.;
        for (j=0; j<m; j++) {
          double aij = a[i][j];
          double * bj = b[j];
          for (k=0; k<m; k++)  ci[k] += aij*bj[k];
        }
      }
    }

Mais uma observação: Esta implementação é ainda melhor no meu computador do que substituir tudo pela rotina BLAS cblas_dgemm (experimente no seu computador!). Mas muito mais rápido (1: 4) está chamando dgemm_ da biblioteca Fortran diretamente. Acho que essa rotina não é de fato Fortran, mas código assembler (não sei o que está na biblioteca, não tenho os fontes). Totalmente incerto para mim é por que cblas_dgemm não é tão rápido, uma vez que, que eu saiba, é apenas um invólucro para dgemm_.

Wolfgang Jansen
fonte
3

Esta é uma aceleração realista. Para obter um exemplo do que pode ser feito com o assembler SIMD sobre o código C ++, veja alguns exemplos de funções de matriz do iPhone - elas eram mais de 8x mais rápidas do que a versão C e nem mesmo são montagens "otimizadas" - ainda não há revestimento interno são operações de pilha desnecessárias.

Além disso, seu código não está " restrito correto " - como o compilador sabe que quando ele modifica C, não está modificando A e B?

Justicle
fonte
Claro, se você chamou a função como mmult (A ..., A ..., A); você certamente não obteria o resultado esperado. Mais uma vez, embora eu não estivesse tentando superar / reimplementar o BLAS, apenas vendo o quão rápido ele realmente é, então a verificação de erros não estava em mente, apenas a funcionalidade básica.
DeusAduro de
3
Desculpe, para ser claro, o que estou dizendo é que se você colocar "restringir" em seus ponteiros, obterá um código muito mais rápido. Isso ocorre porque toda vez que você modifica C, o compilador não precisa recarregar A e B - acelerando drasticamente o loop interno. Se você não acredita em mim, verifique a desmontagem.
Artigo de
@DeusAduro: Isso não é verificação de erro - é possível que o compilador não consiga otimizar os acessos ao array B [] no loop interno porque ele pode não ser capaz de descobrir que os ponteiros A e C nunca se chamam de B array. Se houvesse aliasing, seria possível que o valor na matriz B mudasse enquanto o loop interno estivesse em execução. Elevar o acesso ao valor B [] para fora do loop interno e colocá-lo em uma variável local pode permitir que o compilador evite acessos contínuos a B [].
Michael Burr de
1
Hmmm, então tentei primeiro usar a palavra-chave '__restrict' no VS 2008, aplicada a A, B e C. Isso não mostrou nenhuma alteração no resultado. No entanto, mover o acesso a B, do loop mais interno para o loop externo, melhorou o tempo em ~ 10%.
DeusAduro de
1
Não tenho certeza sobre o VC, mas com o GCC você precisa habilitar -fstrict-aliasing. Também há uma explicação melhor para "restringir" aqui: cellperformance.beyond3d.com/articles/2006/05/…
Justicle de
2

Com relação ao código original em multiplicação MM, a referência de memória para a maioria das operações é a principal causa do mau desempenho. A memória está funcionando de 100 a 1000 vezes mais lenta que o cache.

A maior parte da velocidade vem do emprego de técnicas de otimização de loop para esta função de loop triplo na multiplicação MM. Duas técnicas de otimização de loop principais são usadas; desenrolando e bloqueando. Com relação ao desenrolamento, desenrolamos os dois loops mais externos e os bloqueamos para reutilização de dados no cache. O desenrolamento do loop externo ajuda a otimizar o acesso aos dados temporariamente, reduzindo o número de referências à memória para os mesmos dados em momentos diferentes durante toda a operação. Bloquear o índice de loop em um número específico ajuda a reter os dados no cache. Você pode escolher otimizar para cache L2 ou cache L3.

https://en.wikipedia.org/wiki/Loop_nest_optimization

Pari Rajaram
fonte
-24

Por muitas razões.

Em primeiro lugar, os compiladores Fortran são altamente otimizados e a linguagem permite que sejam assim. C e C ++ são muito flexíveis em termos de manipulação de array (por exemplo, o caso de ponteiros que se referem à mesma área de memória). Isso significa que o compilador não pode saber com antecedência o que fazer e é forçado a criar um código genérico. No Fortran, seus casos são mais simplificados e o compilador tem melhor controle do que acontece, permitindo que ele otimize mais (por exemplo, usando registradores).

Outra coisa é que o Fortran armazena coisas em colunas, enquanto C armazena dados em linhas. Eu não verifiquei seu código, mas tome cuidado com a forma como você executa o produto. Em C, você deve fazer a varredura de linha: desta forma, você faz a varredura de seu array ao longo da memória contígua, reduzindo as perdas de cache. A falta de cache é a primeira fonte de ineficiência.

Terceiro, depende da implementação do blas que você está usando. Algumas implementações podem ser escritas em assembler e otimizadas para o processador específico que você está usando. A versão netlib foi escrita em fortran 77.

Além disso, você está fazendo muitas operações, a maioria delas repetidas e redundantes. Todas essas multiplicações para obter o índice são prejudiciais para o desempenho. Eu realmente não sei como isso é feito no BLAS, mas existem muitos truques para evitar operações caras.

Por exemplo, você pode retrabalhar seu código desta forma

template<class ValT>
void mmult(const ValT* A, int ADim1, int ADim2, const ValT* B, int BDim1, int BDim2, ValT* C)
{
if ( ADim2!=BDim1 ) throw std::runtime_error("Error sizes off");

memset((void*)C,0,sizeof(ValT)*ADim1*BDim2);
int cc2,cc1,cr1, a1,a2,a3;
for ( cc2=0 ; cc2<BDim2 ; ++cc2 ) {
    a1 = cc2*ADim2;
    a3 = cc2*BDim1
    for ( cc1=0 ; cc1<ADim2 ; ++cc1 ) {
          a2=cc1*ADim1;
          ValT b = B[a3+cc1];
          for ( cr1=0 ; cr1<ADim1 ; ++cr1 ) {
                    C[a1+cr1] += A[a2+cr1]*b;
           }
     }
  }
} 

Experimente, tenho certeza de que você salvará algo.

Em sua pergunta nº 1, o motivo é que a multiplicação de matrizes escala como O (n ^ 3) se você usar um algoritmo trivial. Existem algoritmos que escalam muito melhor .

Stefano Borini
fonte
36
Esta resposta está completamente errada, desculpe. As implementações do BLAS não são escritas em fortran. O código crítico de desempenho é escrito em assembly, e os mais comuns hoje em dia são escritos em C acima disso. Além disso, o BLAS especifica a ordem das linhas / colunas como parte da interface e as implementações podem lidar com qualquer combinação.
Andrew Tomazos
10
Sim, esta resposta está completamente errada. Infelizmente, está cheio de absurdos comuns, por exemplo, a alegação de que o BLAS foi mais rápido por causa do Fortran. Ter 20 (!) Avaliações positivas é uma coisa ruim. Agora, esse absurdo se espalha ainda mais por causa da popularidade do Stackoverflow!
Michael Lehn
12
Acho que você está confundindo a implementação de referência não otimizada com as implementações de produção. A implementação de referência serve apenas para especificar a interface e o comportamento da biblioteca e foi escrita em Fortran por motivos históricos. Não é para uso em produção. Na produção, as pessoas usam implementações otimizadas que exibem o mesmo comportamento da implementação de referência. Estudei o interior do ATLAS (que suporta o Octave - Linux "MATLAB") que posso confirmar em primeira mão está escrito em C / ASM internamente. As implementações comerciais também são quase certas.
Andrew Tomazos,
5
@KyleKanos: Sim, aqui está a fonte do ATLAS: sourceforge.net/projects/math-atlas/files/Stable/3.10.1 Pelo que sei, é a implementação BLAS portátil de código aberto mais comumente usada. Está escrito em C / ASM. Fabricantes de CPU de alto desempenho, como a Intel, também fornecem implementações BLAS especialmente otimizadas para seus chips. Eu garanto que as partes de baixo nível da biblioteca Intels são escritas em (duuh) x86 assembly, e tenho certeza que as partes de nível médio seriam escritas em C ou C ++.
Andrew Tomazos,
9
@KyleKanos: Você está confuso. O Netlib BLAS é a implementação de referência. A implementação de referência é muito mais lenta do que as implementações otimizadas (consulte a comparação de desempenho ). Quando alguém diz que está usando o netlib BLAS em um cluster, não significa que está realmente usando a implementação de referência do netlib. Isso seria simplesmente bobo. Significa apenas que eles estão usando um lib com a mesma interface do netlib blas.
Andrew Tomazos,