Apresentação de CUDA - Implementando uma convolução Bidimensional

Autor: Luke
Data:21/10/2010

Introdução

A evolução das Graphics Processor Units ( GPUs ) tem criado uma notável diferença de poder de execução em comparação com as convencionais CPUs. Esta evolução tem gerado uma grande demanda por desenvolvedores habilitados a utilizar o poder destas placas ao máximo. Porem esta nova tecnologia demanda um paradigma de programação em paralelo, um modelo de computação bem diferente dos comuns orientados a objetos e seriais. A linguagem utilizada pra programar estas placas é a CUDA, e esta será utilizada e exemplificada no decorrer deste artigo. A Imagem 1 é um levantamento do poder de processamento de 2001-2009 entre as GPUs e as CPUs.

/media/Attachments/courseIA366F2S2010/lc100373_6/fig1.png

Figure 1: Evolução entre GPUs e CPUs

O objetivo é demonstrar os passos envolvidos na criação de um programa em CUDA para usuários sem conhecimento de prévio de programação paralela ou CUDA. Alguns métodos exemplificados são uma compilação de conhecimento próprio adquirido durante o estudo de CUDA voltado para processamento de imagens, um simples exemplo de convolução bidimensional será incrementado desde sua implementação mais simples até uma implementação mais sofisticada com o intuito de acompanhar o leitor através das variadas etapas necessárias para extrair ótimos resultados do processamento em CUDA. O artigo tem um fundo acadêmico para iniciantes, portanto não levará o leitor ao ápice de aprimoramento do algoritmo, e nem tem como intuito definir um ápice para o algoritmo. A medida que os algoritmos são aprimorados sua legibilidade é perdida no processo, e a funcionalidade do código se perde entre as técnicas de aprimoramento.

Conceitos Básicos

Antes de entrar no CUDA é necessário entender a arquitetura que existente todo o paradigma é dado em função deste modelo e a programação é fortemente influenciada pela arquitetura, o conceito básico será mostrado, já que para entender o funcionamento dos códigos é necessário ter esta base. Inicialmente será exibido o funcionamento em serial do algoritmo de Convolução Bi-dimensional, este algoritmo será o único utilizado no artigo, portanto há necessidade de deixar claro seu funcionamento.

Convolução Bi-Dimensional

A convolução é um operador extremamente útil para o processamento de imagens, várias técnicas são implementadas utilizando apenas uma variação da máscara. O algoritmo define que um pixel é dado pela soma ponderada dele e de seus vizinhos, onde a ponderação é definida por uma máscara. O algoritmo 1 mostra como seu funcionamento é extremamente simples e por isso foi um candidato a ser utilizado neste estudo.

Algoritmo 1: Convolução Bi-Dimensional Serial

O algoritmo trabalha basicamente com quatro “FOR”s aninhados. Onde os dois primeiros ( Algorimo 1, linhas 9 e 10 ) percorrem todos os pixeis da imagem, e os dois últimos percorrem a vizinhança do pixel atual para obter a soma ponderada. Esse algoritmo dita uma certa seqüência de execução na qual os pixeis serão acessados, a figura 2 mostra a ordem em que uma suposta imagem de entrada seria percorrida, e os índices utilizados para uma suposta máscara 3x3 tendo o centro como origem.

/media/Attachments/courseIA366F2S2010/lc100373_6/fig3.png

Figure 2: Ordem do algoritmo serial

Analisando o algoritmo podemos calcular a quantidade de operações necessárias para obter a convolução. O algoritmo faz para cada pixel da imagem (L*K) multiplicações e (L*K-1) somas, olhando o conjunto como operações temos [(L*K)+(L*K-1)] operações por pixel. Adicionando o tamanho da imagem para encontrar a quantidade de pixeis temos, OPE=(M*N)[ (L*K)+(L*k-1)]. Esta equação irá gerar valores enormes dependendo de suas entradas, alguns valores foram gerados para mostrar o crescimento de operações variando a imagem de entrada e a mascara, a Tabela 1 e Figura 3 possuem os resultados destes cálculos.

Imagem Máscara Multiplicações Somas Operações
(100x100) 10000 9 90000 80000 170000
(500x500) 250000 9 2250000 2000000 4250000
(1000x1000) 1000000 9 9000000 8000000 17000000
(2000x2000) 4000000 9 36000000 32000000 68000000
(4000x4000) 16000000 9 144000000 128000000 272000000
(100x100) 10000 81 810000 800000 1610000
(500x500) 250000 81 20250000 20000000 40250000
(1000x1000) 1000000 81 81000000 80000000 161000000
(2000x2000) 4000000 81 324000000 320000000 644000000
(4000x4000) 16000000 81 1296000000 1280000000 2576000000
Table 1 : Calculo da quantidade de Operações

O algoritmo de convolução é uma boa escolha para implementação em GPU devido a esta alta quantidade de operações necessárias (Tabela 1 e Figura 3), e sua necessidade de utilizar dados na região na um pixel, no decorrer do artigo esta afirmação ficará mais clara.

/media/Attachments/courseIA366F2S2010/lc100373_6/fig4.png

Figure 3: Quantidade de Operações x Imagem

CUDA (Compute Unified Device Architecture)

CUDA é uma arquitetura de computação paralela desenvolvida pela NVIDIA, que permite melhorar drasticamente o desempenho computacional utilizando o poder das GPUs. Para possibilitar essa melhora é necessário entender bem esta arquitetura, e compreender os recursos que esta tecnologia disponibiliza. A grande diferença de desempenho é dada pela diferente filosofia dos processadores, uma CPU (Figura 4, a esquerda) tem como objetivo otimizar códigos seriais, já uma GPU (Figura 4, a direita) tem o objetivo de maximizar a quantidade de cálculos com ponto flutuantes, para isto os cálculos são feitos em massivas quantidades de threads.

Núcleos de Processamento/ Threads

Em uma CPU é utilizado uma sofisticada lógica de controle e um grande cache de memória para diminuir a latência de acesso a memória. Em uma CPU comum quatro núcleos de processamento estão organizados para ter um forte desempenho seqüencial.

/media/Attachments/courseIA366F2S2010/lc100373_6/fig5.png

Figure 4: Organização de uma CPU e uma GPU respectivamente

Voltando a GPU, o hardware toma vantagem da massiva quantidade de threads para encontrar trabalho quando a execução depende de acessos a memória de alta latência, isto diminui a lógica de controle necessária para cada thread (diminuindo o chip físico). Pequenos caches de memória (impacta o tamanho físico do chip) estão disponíveis para auxiliar no acesso a memória DRAM. Com isto temos uma grande área do chip disponível para os processadores de ponto flutuante. É importante tomar conhecimento destas grandes diferenças, já que isto determina se o desempenho será melhor para algumas tarefas em CPU caso sejam altamente serias, ou em GPU caso seja capaz de serem paralelizadas. Em sua organização a GPU possui um conjunto de processadores (Streaming Multiprocessor [SM], na Figura 6 isto é visível pelo agrupamento de cada linha composta por um cachê, uma unidade de controle e um conjunto de ALU) que compartilham uma unidade de controle e um pequeno cache de memória, o acesso a este cache esta restrito, para uso somente as threads dentro do mesmo SM, porem seu acesso é mais rápido do que a memória global (DRAM), este nome é dado pelo fato de todas as threads de todos os SM possuírem acesso, o custo deste múltiplo acesso é o alto tempo de acesso necessário.

Diferentes Memórias

A GPU trabalha com 5 tipos de memórias, sendo elas registradores, cacheada, global,constante e de textura. Todas possuem suas vantagens/desvantagens e limitações, um grande esforço durante a implementação de um código em CUDA é voltado ao acesso à memória. Uma etapa crucial para otimizar implementações está no uso adequado destas memórias. A solução mais simples e direta é utilizar os dados em memória global, enquanto esta implementação é bem fácil de ser desenvolvida, seu aprimoramento é extremamente restrito, os tempos de acesso aos dados serão altos devido a latência de acesso a esse tipo de memória e os ganhos na implementação serão baixos. Os registradores são extremamente limitados e só podem ser utilizados para variáveis unitárias, portanto acesso a esta memória e extremamente rápida. A memória cache também possui um tamanho restrito, ela aceita dados na forma de array, e seu acesso é extremamente rápido. Esta memória requer certa lógica de carregamento que é mais complexa que a global, isto é devido ao fato desta memória ser acessível somente dentro do próprio SM. A memória constante pode ser lida por qualquer SM porem não pode ser alterada, esta memória possui alta velocidade de acesso. Por ultimo temos a memória de textura, onde os dados estão armazenados na memória global e são cacheados transparentemente com eficiência.

Conceitos de Implementação

Para começar explorar implementações em CUDA há necessidade de definir como as estruturas definidas acima são acessíveis para os programadores. A estrutura de programação básica ira ser apresentada seguida de como acessar as threads e memórias. As estruturas serão exibidas de modo genérico para singularizar o uso de cada um destes itens.

Host/Device

Na programação seu código utiliza tanto a CPU quanto a GPU, cabe ao desenvolvedor escolher onde fazer a implementação de seu código para obter os melhores resultados. Para acessar o host ou o device as funções possuem especificação de onde deve ser executada, adicionalmente isto define de onde ela pode ser chamada. Existem, portanto três modificadores (Figura 5), “__device__” define uma função que será executada na GPU e possui funcionamento paralelo, esta função só pode ser invocada por uma função sendo executada dentro da GPU. O segundo modificador é “__global__”, uma função global é executada na GPU e pode ser invocada de fora da GPU, esta função é a interface de entrada/saída entre a CPU e GPU, esta função é executada em seriale tem acesso a memória do host e do device. Por ultimo o modificador “__host__” define uma função executada na CPU que só pode ser invocada na CPU, esta função definirá uma interface transparente para o restante do programa.

/media/Attachments/courseIA366F2S2010/lc100373_6/fig7.png

Figure 5: Modificadores de Funções e comunicação Host/Device respectivamente

Analisando as definições de chamadas/comunicação uma ordem de criação é definida, uma função host deve organizar os dados para uma função global e pode executar algum pré/pos processamento em serial à função global. Cabe a função host organizar os dados em memória da GPU, e ditar a configuração da GPU para o kernel (função global paralela). Nesta etapa é importante ressaltar que a transferência de dados da memória da CPU para a da GPU é demorada, portanto a quantidade de dados a serem transferidos deve ser limitada.

Configurando a GPU para processamento

O kernel é executado por um array de threads, onde todos os threads do array executam o mesmo kernel. Para que um thread consiga se diferenciar de outras ela possui um identificador único, com esse identificador a thread consegue acessar diferentes áreas da memória e trabalhar com dados diferentes enquanto executando o mesmo kernel, a figura 6 mostra um exemplo de como esta organização funcionaria e como o ID do thread poderia ser utilizado para acessar diferentes dados na memória.

/media/Attachments/courseIA366F2S2010/lc100373_6/fig8.png

Figure 6: Thread ID Array e memória

Adicionalmente para dar escalabilidade ao kernel é possível dividir o thread em diferentes grupos, estes grupos são chamados de Blocos (BLOCK). Um conjunto de blocos por sua vez compõe uma grade (GRID). Para facilitar o trabalho com dados multidimensionais uma grade pode ser dividida em uma matriz bidimensional de blocos, e um bloco pode ser dividido em uma matriz tridimensional de threads. Este arranjo é chamado de hierarquia de memória CUDA (Figura 7).

/media/Attachments/courseIA366F2S2010/lc100373_6/fig9.png

Figure 7: Hierarquia de Memória

Em código estas estruturas são facilmente definidas e acessadas. A configuração é feita por função que será executada na GPU, o tipo auxiliar “dim3” é um array de 3 inteiros, qualquer valor não definido é inicializado com 1. Abaixo o código 1 ilustra como configurar um kernel na GPU. Para acessar os índices atribuídos as threads/blocos, variáveis internas estão disponíveis às threads, a tabela 2 mostra estas variáveis e seus conteúdos.

Variável Modificadores Conteúdo
gridDim .x,.y Contém as dimensões da grade.
blockDim .x,.y,.z Contém as dimensões do bloco.
blockIdx .x,.y Índices de um bloco dentro da grade
threadIdx .x,.y,.z Índices de um thread dentro do bloco
Table 2 : Variáveis Auxiliares

Organizando dados em memória da GPU

A segunda obrigação da função host é copiar os dados da CPU para a GPU, como existem diferentes memórias cada memória possui uma maneira de inicialização, atribuição e definição. Explorar a eficiência de cada tipo de memória é de extrema importância para uma aplicação otimizada. A figura 8 mostra como as diferentes memórias são organizadas dentro da GPU, adicionalmente as setas indicam a direção de acesso aos dados.

/media/Attachments/courseIA366F2S2010/lc100373_6/fig10.png

Figure 8: Organização das Memórias

A criação de variáveis é igual a de um código em C, o que diferencia as variáveis são os qualificadores. Cada qualificador especifica o local onde a variável será criada. A tabela 3 esclarece as variações de qualificadores, em que memória a variável irá residir, quem possui acesso a variável e seu tempo de duração. O qualificador “__device__” é opcional caso outro qualificador também seja utilizado. Uma variável quando criada sem um qualificador é guardada em um registrador, com exceção de ponteiro, que serão criados na memória global. Conforme mostra a figura 8 a copia de dados do host para a memória global e constante é direta e bastante simples, as reais dificuldades são encontradas quando se utiliza as memórias compartilhadas. Como é necessário passar pelos threads, o carregamento dos dados deve ser feito durante a execução do kernel, o que se prova não ser tão trivial dependendo da implementação. Como o tamanho dos registradores e memória compartilhada são restritos cabe ao desenvolvedor conhecer o problema e segmentar apropriadamente os dados. A memória global possui amplo espaço porem é lenta, a memória compartilhada é rápida porem limitada e envolve um passo extra de carregar os dados durante a execução do kernel. A memória constante reside dentro da memória global, porem é cacheada, devido ao fato de não aceitar escrita isso a torna inviável para algumas variáveis que necessitam ser atualizadas.

Declaração Memória Escopo Duração
__device__ __local__ int LocalVar; Local Thread Thread
__device__ __shared__ int SharedVar; Compartilhada Bloco Bloco
__device__ int GlobalVar; Global Grade Aplicação
__device__ __constant__ int ConstantVar; Constante Grade Aplicação
Table 3: Qualificadores de Variáveis

Primeira Implementação

Agora que as estruturas do CUDA foram apresentadas, e como acessá-las um programa pode começar a ser concebido. O algoritmo de convolução apresentado deve agora ser re-pensado para utilizar um kernel. Não existe uma formula para transformar um código serial em paralelo, a seqüência de passos exibidas aqui é apenas um exemplo de como visualizar esta transformação. O primeiro passo após examinar o algoritmo serial é tentar encontrar a menor parte do algoritmo que executa com singularidade, esta parte deve conter o trecho de código que faz um calculo unitário ou de um conjunto relativamente pequeno de cálculos repetidas vezes. O segundo é encontrar como os laços fazem a distribuição destes vários cálculos, descobrir se estes necessitam de uma ordem específica ou se estão apenas sendo executados em ordem. Voltando rapidamente a figura 2 esta ordem já foi exposta durante a explicação do algoritmo, enquanto o algoritmo possui quatro laços serializados, os dois exteriores que percorrem a imagem é o que aparenta gerar o grande número de interações (já que geralmente a imagem será maior que a mascara). A idéia, portanto para otimizar o código é gerar um kernel que faça o calculo conjunto dos dois laços internos e que cada thread processe o equivalente a um pixel. Nesta primeira implementação simplicidade e legibilidade do código é priorizado, portanto os dados serão organizados em memória global a distribuição dos threads será bidimensional dentro dos blocos para facilitar a adaptação à imagem. A implementação no código 1 apresenta a função host para configurar a GPU, a tabela 4 explica os parâmetros e a tabela 5 as linhas do código.

Código da Função host

Code 1: Função host da primeira implementação
 1 #include <cuda.h>
 2 #define TILE_DIM_W 16
 3 #define TILE_DIM_H 32
 4 
 5 __host__ void cdConv_cpp( float *I, int wI, int hI,
 6                           float *M, int wM, int hM,
 7                           float *R ){
 8 
 9     float *dI;
10     float *dM;
11     float *dR;
12 
13     cudaMalloc((void **)&dI, wI*hI * sizeof(float));
14     cudaMalloc((void **)&dM, wM*hM * sizeof(float));
15     cudaMalloc((void **)&dR, wI*hI * sizeof(float));
16 
17     cudaMemcpy(dI, I, wI*hI * sizeof(float), cudaMemcpyHostToDevice);
18     cudaMemcpy(dM, M, wM*hM * sizeof(float), cudaMemcpyHostToDevice);
19     cudaMemset(dR, 0, wI*hI * sizeof(float));
20 
21     int fixedGridW = (wI % TILE_DIM_W == 0)? wI/TILE_DIM_W : wI/TILE_DIM_W+1;
22     int fixedGridH = (hI % TILE_DIM_H == 0)? hI/TILE_DIM_H : hI/TILE_DIM_H+1;
23     dim3 gridDim(fixedGridW,fixedGridH );
24     dim3 blockDim(TILE_DIM_W,TILE_DIM_H);
25 
26     cdConv_cd<<<gridDim, blockDim>>>(dI,wI,hI,dM,wM,hM,dR);
27 
28     cudaMemcpy(R, dR, wI*hI * sizeof(float), cudaMemcpyDeviceToHost);
29 
30     cudaFree(dR);
31     cudaFree(dI);
32     cudaFree(dM);
33 
34 }
Parâmetro Conteúdo
I Ponteiro para Imagem de Entrada na CPU
wI Largura da Imagem
hI Altura da Imagem
M Ponteiro para Máscara na CPU
wM Largura da Máscara
hM Altura da Máscara
R Ponteiro para Imagem de Saída na CPU
Table 4: Parâmetros da Função host
Linha(s)  
9-11 Cria ponteiros na memória global da GPU
13-15 Aloca os ponteiros na GPU com os tamanhos iguais a seus equivalentes na CPU
17-19 Copia os valores da CPU para a GPU
21-24 Calcula a quantidade de blocos necessários para cobrir a imagem. Tratando caso a imagem nao seja multipla do tamanho do bloco, exemplificado na figura 10
26 Invoca o kernel passando as imagens, e a mascara junto de seus tamanhos.
28 Copia a imagem calculada pelo kernel da GPU para a CPU
30-32 Libera a memória reservada na GPU.
Table 5: Explicação das linhas de código da função host
/media/Attachments/courseIA366F2S2010/lc100373_6/fig11.png

Figure 9: Desperdício de threads devido ao tamanho do bloco.

Código da Função global

Como foi definido que cada Thread cuidaria do calculo da mascara sobre um pixel, basta que o kernel calcule os dois loops internos do algoritmo 1. Neste kernel é importante notar como o acesso a memória global é direta e simples, enquanto já é possível notar como é importante entender bem o calculo dos índices a partir dos IDs.

Code 2: Função global da primeira implementação
 1 __global__ void cdConv_cd(float* dI, int wI, int hI,
 2                           float* dM, int wM, int hM,
 3                           float* dR){
 4 
 5 
 6     int Col = blockIdx.x * blockDim.x + threadIdx.x;
 7     int Row = (blockIdx.y * blockDim.y) + threadIdx.y ;
 8     int pX,pY,mY,mX;
 9 
10     float v = 0;
11 
12     if (Col >= wI || Row >= hI) return;
13 
14     for( int r = -(int)(hM/2) ; r <= (int)(hM/2) ; r++ ){
15 
16         pY = ((Row+r+hI)%hI)*wI;
17         mY = wM*(r+ (int)(hM/2));
18 
19         for( int c = -(int)(wM/2) ; c <= (int)(wM/2) ; c++ ){
20 
21             pX = (Col+c+wI)%wI;
22             mX = (c+(int)(wM/2));
23 
24             v = v+ dM[ mY+mX ]*dI[ pY+pX ];
25         }
26     }
27     dR[(Row*wI)+Col]= v;
28 }
Linha(s)  
6-7 Calculo dos índices da imagem para um grid com múltiplos blocos.
8-10 Variáveis auxiliares ( do tipo registrador )
12 Para a execução da Thread caso ela não esteja dentro da imagem ( problema exibido na figura 9 )
14-26 Implementa os dois loops que percorre a máscara fazendo as somas ponderadas. A figura 10 mostra como calcular o índice utilizando o ID da Thread.
27 Copia o valor calculado para a nova imagem.
Table 6: Explicação das linhas de código da função global
/media/Attachments/courseIA366F2S2010/lc100373_6/fig12.png

Figure 10: Calculo de Índices da imagem pelos IDs.

Aprimoramento para Memória Compartilhada

Uma importante porem complexa tarefa é otimizar o código para trabalhar com memória compartilhada para otimizar a quantidade de acessos a memória global feita.Devido a alta latência de acesso e concorrência de acesso das múltiplas Thread o desempenho será afetado. Para que a função seja adequada a utilizar memória compartilhada algumas alterações devem ser feitas às funções host e global.

Alterações no Host

A grande mudança no host será como ele configura a GPU, especificamente a grade. Nosso algoritmo requer que não só os dados pontuais estejam disponíveis para o calculo mas também uma região ao redor do pixel. Como as memórias compartilhadas dos blocos não são acessíveis entre si cabe a cada bloco carregar todos os dados necessários para sua execução, a figura 11 mostra como o bloco deve ser configurado. O antigo bloco deve aumentar por metade do tamanho da máscara, assim garantindo que todos os dados necessários serão carregados.

/media/Attachments/courseIA366F2S2010/lc100373_6/fig13.png

Figure 11: Dados a serem carregados por um bloco.

Alterações na global

Uma nova variável do tipo "__shared__" deve ser declarada agora, será nela que os dados serão carregados para executar os cálculos. A função agora se encarregará de copiar os dados da memória global para a memória compartilhada, é importante notar que um novo comando "__syncthreads()" é utilizado. Esta função garante que todas as threads que devem chegar a este ponto antes de seguir adiante, enquanto este comando possui um alto custo no desempenho ele é de extrema importância. Caso este comando não estivesse presente algumas threads poderiam seguir e tentariam utilizar um dado que ainda não foi carregado. Cada thread portanto deverá calcula qual pixel da imagem ele se encarregará de carregar, a figura 12 mostra como os pixeis serão carregados na memória compartilhada, após carregados as threads fora do conteúdo de execução( em vermelho na figura 12 ) devem ser paradas.

/media/Attachments/courseIA366F2S2010/lc100373_6/fig14.png

Figure 12: Organização na memória compartilhada e pixeis a serem executados.

Comparação de Desempenho

Para que o poder do CUDA fique claro as funções exibidas aqui serão comparadas com diferentes entradas para mostrar os ganhos. A tabela 7 mostra os resultado de comparação entre um algoritmo serial em CPU e o paralelo proposto acima.

Lado GPU CPU Speedup Compare
35 0.342 0.157 0.459 True
45 0.345 0.192 0.556 True
55 0.351 0.227 0.647 True
117 0.434 0.756 1.741 True
139 0.454 0.989 2.177 True
177 0.525 1.519 2.893 True
298 0.859 4.076 4.745 True
365 1.236 6.355 5.142 True
400 1.378 6.578 4.773 True
431 1.548 8.746 5.650 True
493 1.888 12.377 6.555 True
655 3.004 25.537 8.501 True
659 3.073 26.182 8.520 True
939 5.700 108.070 18.959 True
946 5.766 110.257 19.121 True
1046 6.840 142.819 20.880 True
1075 7.273 150.849 20.741 True
1110 7.653 159.756 20.875 True
1305 10.215 220.609 21.596 True
1499 12.884 290.674 22.561 True
1500 12.950 292.779 22.609 True
1577 14.193 322.031 22.689 True
1629 15.422 345.279 22.389 True
1747 17.280 398.863 23.083 True
1765 17.577 402.681 22.910 True
1964 21.728 499.397 22.984 True
2006 22.243 522.384 23.485 True
2021 22.628 530.112 23.427 True
2027 22.680 536.175 23.641 True
Table 7: Teste de desempenho CPU/GPU.

Conclusões

Fica claro o poder de processamento das GPU em comparação as CPU quando observamos os valores obtidos no teste de desempenho. Até mesmo a implementação mais simples atingiu um ganho na ordem de 16-17 vezes mais rápido do que a CPU, portanto não se mostra eficiente para pequenas quantidades de dados (pouco conteúdo de processamento), como já era de se esperar. Um bom ganho é possível se ser obtido utilizando o mínimo da GPU, mantendo a complexidade de implementação baixa.