BLOG C++

“I want it all and I want it now”

Diversas vezes temos de programar um algoritmo para realizar um processamento e queremos uma resposta rápida. Nesta situação, muitos programadores saem a escovar bits trocando multiplicações por 2 por shifts, enchem o programa de ponteiros e se esquecem das três soluções iniciais que devem ser pensadas:

1) Usar um algoritmo eficiente: por melhor que seja a sua implementação de bubble-sort a complexidade do algoritmo é O(n²) ao passo que existem algoritmos para ordenação com complexidade O(n log n).

2) Esquecem de ligar o otimizador

3) Usam construções que impedem a otimização do código.

Eu diria que hoje em dia, escovar bits é o último recurso, para ser usado somente quando você vai tentar salvar o mundo…

No lugar disso, devemos seguir as recomendações anteriores, sempre escrevendo o código de forma que o otimizador consiga dar jeito. Não é difícil: em geral, quanto mais legível é o programa mais fácil é para o otimizador melhorar o código.

Vamos ver um exemplo de uma rotina simples para multiplica matrizes n x n com n grande, digamos 1000 – ou seja, matrizes de um milhão de elementos…

Multiplicação de Matrizes

O primeiro passo é a escolha do algoritmo. O algoritmo básico que aprendemos em Álgebra possui complexidade O(n³). Embora existam algoritmos com complexidade menor, tais como os algoritmos de Strassen (complexidade O(n^2.807)) e o de Coppersmith–Winograd (complexidade O(n^2.376)), eles só se tornam mais rápidos que o algoritmo comum a partir de n realmente grande, como apontado neste artigo. Portanto, fiquemos com o algoritmo comum, mostrado no código abaixo. As matrizes são inicializadas pseudo-randomicamente, mas nesse caso isso não é importante.

#include <cstdlib>

int main(int argc, char *argv[])
{
  const int N = 1000;
  typedef double Linha[N];
  Linha *a = new Linha[N],
        *b = new Linha[N],
        *c = new Linha[N];

  srand(1);
  for( int i = 0; i < N; i++ )
    for( int j = 0; j < N; j++ ) {
      a[i][j] = rand()/(rand()+0.1);
      b[i][j] = rand()/(rand()+0.1);
      c[i][j] = 0;
  }  

  for( int i = 0; i < N; i++ )
    for( int j = 0; j < N; j++ )
      for( int k = 0; k < N; k++ )
        c[i][j] += a[i][k] * b[k][j];

  delete [] a;
  delete [] b;
  delete [] c;
}

O programa acima, compilado no Dev-C++ 4.9.9.2 usando o g++ 3.4.2 (um pouco ultrapassado…) e sem nenhuma modificação nos parâmetros defaults do compilador, roda em 18.50 segundos (média de 4 execuções sem nenhum outro programa executando e descartando a primeira execução – por causa do cache de disco). O problema aqui é que por default o otimizador do g++ fica desligado. Então, vamos ligá-lo passando alguns parâmetros:

-O3 -funroll-loops -fprefetch-loop-arrays -march=pentium4 -msse2

O resultado melhora um pouco: 15.95 segundos – um ganho de 13,7%. Pouco? Eu também acho, mas o interessante é que veio sem trabalho nenhum… Na verdade, ligando apenas o otimizador (-O3) já há um ganho significativo: 16.16 segundos. Para entender o que essas outras opções significam recomendo olhar a documentação do g++.

Para mostrar que escovar bits pura e simplesmente não resolve, fiz duas alterações no programa acima que não resultaram em nenhum ganho de desempenho, na verdade até piorou um pouco…

1) Alocação estática ds matrizes

static double a[N][N],
              b[N][N],
              c[N][N];

// O restante do código não muda, exceto pelo delete do final

Essa opção tem o inconveniente de alocar a memória o tempo todo da execução. Essa versão levou 16.07 segundos…

2) Uso de ponteiros duplos

int main(int argc, char *argv[])
{
  const int N = 1000;
  double **a = new double*[N],
        **b = new double*[N],
        **c = new double*[N];

  for( int i = 0; i < N; i++ ) {
    a[i] = new double[N];
    b[i] = new double[N];
    c[i] = new double[N];
  }

  srand(1);
  for( int i = 0; i < N; i++ )
    for( int j = 0; j < N; j++ ) {
      a[i][j] = rand()/(rand()+0.1);
      b[i][j] = rand()/(rand()+0.1);
      c[i][j] = 0;
    }  

  for( int i = 0; i < N; i++ )
    for( int j = 0; j < N; j++ )
      for( int k = 0; k < N; k++ )
        c[i][j] += a[i][k] * b[k][j]; 

  for( int i = 0; i < N; i++ ) {
    delete [] a[i];
    delete [] b[i];
    delete [] c[i];
  }

  delete [] a;
  delete [] b;
  delete [] c;
}

Note que é necessário alocar (e desalocar) espaço para cada uma das linhas da matriz. Essa versão levou 16.95 segundos… E se passarmos para aritmética de ponteiros o resultado só piora: 17.54 segundos. Porque isso ocorre? Por que o otimizador não consegue otimizar muito bem ponteiros…

E agora? Bem, agora é hora de olhar o código e pensar.  De cara, vemos que na linha

c[i][j] += a[i][k] * b[k][j]

apenas k varia no loop mais interno. Porém, isso o otimizador também nota. Olhando para a[i][k], e já  que no loop interno apenas k varia, estamos lendo posições sequenciais na memória. Porém, em b[k][j] isso não ocorre. Acessos sequenciais na memória são muito bons, porque o processador pode utilizar bem o cache L1 e o prefetching. Olhando o comando acima, vemos que se transpusermos a matriz b (b^t[j][k] = b[k][j]) , o acesso seria sequencial. Será que isso faz diferença? Faz e muita… não vou colocar o código para transpor a matriz, é simplesmente  dois loops aninhados trocando b[i][j] com b[j][i], e mais para frente veremos que não vamos precisar disso. Assim, a linha mais interna do loop fica:

c[i][j] += a[i][k] * b[j][k]

E o tempo de execução? Pasmem… 2.61 segundos… Inacredtável, não? 7x mais rápido!

Mas transpor a matriz em outros casos pode não ser tão simples. Melhor procurar outra opção. Olhando os índices da linha original, o melhor seria se no loop interno quem variasse fosse o j.  a[i][k] se tornaria constante e tanto b quanto c seriam acessados sequencialmente. Mas é possível trocar a ordem dos loops? Sim, e isso é uma otimização conhecida de muitos otimizadores, porém não está presente no g++3.4.2 usado pelo Dev-C++. A boa notícia é que está presente nas versões mais recentes do g++, embora não ainda para Windows.

Enquanto essa otimização não vem, trocamos a ordem do loop para ver o resultado:

#include <cstdlib>

int main(int argc, char *argv[])
{
  const int N = 1000;
  typedef double Linha[N];
  Linha *a = new Linha[N],
        *b = new Linha[N],
        *c = new Linha[N];

  srand(1);
  for( int i = 0; i < N; i++ )
    for( int j = 0; j < N; j++ ) {
      a[i][j] = rand()/(rand()+1);
      b[i][j] = rand()/(rand()+1);
      c[i][j] = 0;
  }  

  for( int i = 0; i < N; i++ )
    for( int k = 0; k < N; k++ )
      for( int j = 0; j < N; j++ )
        c[i][j] += a[i][k] * b[k][j];

  delete [] a;
  delete [] b;
  delete [] c;
}

Tempo: 2.79 segundos, bem semelhante à transposição da matriz. E de qualquer forma, essa modificação não afeta em nada a clareza do código. E o melhor, há compiladores que já a fazem.

Por último, podemos ainda usar outro compilador. Para este exemplo , testei com o compilador Intel C++, reconhecidamente um dos mais rápidos do mercado. No caso do primeiro e do último programa não houve diferença nesses casos pois esse compilador faz essa otimização de trocar a ordem dos loops. Opções:

/O3 /QxHost

Resultado: primeiro programa, 0.62 segundos!!! Quase 30x mais rápido! E dessa vez, sem termos de mudar uma linha no código! Em parte, essa diferença pode ser explicada também pelo fato de que o compilador Intel faz um bom uso das instruções vetoriais presentes nos processadores Intel. Já a versão com ponteiros para ponteiros levou 0.80 segundos… parece pouco, mas ela foi quase 30% mais lenta…

Resumindo:

Sempre é bom olhar as alternativas antes de sair otimizando código na mão para usar ponteiros (que trazem a falsa sensação de velocidade), e principalmente, entender o que o compilador é capaz de otimizar e o que ele não otimiza, o que afeta a velocidade do processador  onde o programa irá executar e quais recursos ele oferece. Uma simples troca de posição dos loops deu um ganho enorme! Uma simples troca de compilador deu um ganho maior ainda…

Até o próximo post!