Sinopse do artigo: BEMxGPU

Achei um artigo muito interessante,

Takahashi, T.; Hamada, T. (2009): “GPU-accelerated boundary element method for Helmhotz’ equation in three dimensions“.

Mesmo para quem conhece Elementos de Contorno, vale a pena ler o artigo. Aprendi várias coisas sobre GPU com ele, mesmo sem ter entendido muito a parte de Elementos de Contorno (apesar de ser minha área de pesquisa):

Peak performance

O artigo fala que “um processador escalar é capaz de executar 2 operações aritméticas por ciclo de clock. Nossa GPU tem 1.17 GHz e 12 multiprocessadores, num total de 96 processadores escalares. Portanto, a peak performance é de 228.1 Gflop/s”.

O “portanto“ é novo pra mim. Eu não sabia que era possível calcular a peak performance da máquina a partir desses dados… Como dá pra ver o clock da minha GPU no manual (1,296 GHz), e eu sei que ela tem 30 multiprocessadores, agora sei que a peak performance dela é 622.1 Gflop/s.

Fonte de citação

O artigo revê vários conceitos da programação CUDA e da tecnologia GPGPU. Isso é bom não só pra reforçar nosso know-how, mas também é sempre bom ter fontes adicionais pra citar nos nossos artigos. Frases de efeito que o artigo usa que talvez você queira citar:

• Um warp tem mesmo 32, e não 16 threads;
• O manual CUDA recomenda blocos com múltiplos de 64 threads para máxima performance;
• A velocidade de floating-point operations throughput atual das GPUs é até 20 vezes maior que da CPU;
• Programas de engenharia têm sido experimentados na GPU antes de CUDA, com OpenGL, DirectX, Cg e dispositivos dedicados.

Custo computacional

O artigo revisa o custo computacional de importantes métodos matemáticos. Por exemplo, ele fala que o custo computacional da solução direta de sistemas lineares (como escalonamentos, pivoteamentos, eliminação de Gauss) é da ordem de O(N³), e o consumo de memória é O(N²). Ele diz isso pra justificar a escolha de métodos iterativos no trabalho, como o gradiente conjugado, cujo custo computacional cai pra O(N²) e o consumo de memória cai pra O(N). Boa justificativa pro pessoal que mexe com gradiente conjugado & cia! O artigo ressalva, contudo, que o custo nesses métodos iterativos é proporcional ao número de iterações necessários para o método convergir, e que portanto um pré-condicionamento da matriz é sempre necessário para minimizar esse custo.

Em outro ponto no final do artigo, eles confirmam essa afirmação comparando um pré-condicionador (point-Jacobi), que permitiu convergência em 543 passos e 7498 s, com outro (block-diagonal) que convergiu em 446 passos em 9095 s (não, o resultado não está invertido). Essa diferença de 21% na eficiência justifica investir no estudo de pré-condicionadores.

E já que estamos falando de sistemas lineares: o artigo também alerta que sistemas lineares podem apresentar autovalores fictícios, que podem não só desacelerar a convergência do solver como degenerar os resultados. Há métodos para remover esses autovalores, como o de Burton-Miller, citado no trabalho.

Benchmark justo

Quando comparamos um código em GPU com CPU, os códigos têm que ser o mais parecidos possível para que a comparação seja justa. Para garantir que a comparação seja justa, o artigo faz uma coisa que eu nunca tinha visto: ele conta o número de operações que o algoritmo executa: quantas adições/subtrações, multiplicações, e uso de funções transcendentais. Aí é só garantir que CPU e GPU executem o mesmo número de operações, e a comparação terá sido justa. De quebra, aprendi que há uma forma de executar divisões implementada em hardware! É a função __fdividef(alfa,beta), que equivale a alfa/beta (eu só conhecia __sinf e __cosf).

Precisão mista

O artigo traz uma explicação bem detalhada sobre como fazer precisão mista, que serve para usar precisão dupla só quando for necessário ou para emular precisão dupla em GPUs que não têm esse recurso. Tem até um trecho de código explicando como fazer isso:

static __device__ __inline__ void dsacc0(float *ah, float *al, float b)

{

float t1 = *ah + b;

float e = t1 – *ah;

float t2 = ((b-e) + (*ah – (t1 – e))) + *al;

*ah = t1 + t2;

*al = t2 – (*ah – t1);

}

Todos sabemos que precisão dupla ainda é um ponto fraco das GPGPUs. A maioria só tem precisão simples, e nas que têm precisão dupla, o desempenho é dezenas de vezes inferior do que o da precisão simples. No entanto, precisão dupla, que é inútil em cálculos gráficos, é essencial em problemas de engenharia como o método dos elementos de contorno (assunto do artigo).

Threads por bloco

Na página 18, os autores contam em detalhes como fizeram para otimizar o número de threads por bloco.

• Nthreads deve ser múltiplo de 64, de acordo com o manual do CUDA;
• Cada bloco tem somente 16kb de memória compartilhada;
• Os múltiplos de 64 que não excedem 16kb (neste problema) são: 64, 128, 192, 256 e 320;
• Blocos com esses números de threads usariam respectivamente 1408, 2816, 4416, 5632 e 7360 registradores (neste problema). Como o número máximo de registradores é 8192, esse não é um parâmetro a se preocupar (neste problema);
• Dois ou mais blocos são recomendados por multiprocessador, para compensar a carga de latência;
• Como os blocos de 192, 256 e 320 threads resultariam em menos que isso, essas opções são descartadas;
• 64 threads resultaria em 5 blocos por multiprocessador, mas o artigo afirma que a eficiência de caching de dados diminui com o tamanho do bloco;
• Conclusão: o tamanho ótimo do bloco é 128 threads (de largura).

É interessante ler a discussão inteira. Dá várias idéias sobre como escolher o tamanho de bloco nos nossos próprios problemas.

Padrões de floats

Devido à imprecisão dos resultados da GPU, não é uma boa idéia usá-la na integração numérica de singularidades. Como essa integração é necessária no método dos elementos de contorno, os autores fizeram com que a parte singular das integrais fosse calculada pela CPU e a parte não-singular pela GPU. O resultado da integral era a soma dos resultados dos dois dispositivos. Aí a surpresa: GPUs não seguem a mesma especificação IEEE-754 sobre aritmética de ponto flutuante que as CPUs seguem! Assim, a partir de um dado número de elementos, a precisão do resultado da integral “quebrava”, isto é, a integral não melhorava se o integrador fosse refinado.

Observe como a GPU em precisão simples perde precisão mesmo com o aumento do refinamento.

Isso parece ter sido tão frustrante que os autores desistiram de usar precisão simples a partir daí (apesar dela ser incrivelmente mais rápida na GPU).

Limite da GPU

Meus programas sempre têm um limite de tamanho de dados que eu posso usar, mas nunca investiguei isso a fundo. No fim, não sei a razão desse limite. Nesse artigo, um gráfico pouco usual é plotado. No eixo y, em vez de milissegundos, são mostrados Gflop/s.

Figura 5: Performance da GPU em tempo total versus tempo de kernel.

Como a porcentagem de capacidade de computação do kernel vai sempre se aproximando da do programa completo e o total tende a uma assíntota (claro, a GPU tem uma peak performance), os autores concluíram que o programa deles era limitado pela aritmética, e não pela largura de banda de memória.

Boa leitura!