Shader OSL Multithread

Julho 24, 2011

Essas últimas semanas eu estava tentando implementar suporte a múltiplas threads do shader OSL. O BRL-CAD usa pthreads para paralelizar o trabalho.

A primeira vez que eu pus o shader OSL pra rodar, tive problemas, provavelmente devido à leitura e escritas simultâneas a dados na memória. Eu tinha adiado esse problema e estava trabalhando apenas com uma thread.

O que a minha função osl_render faz é preencher uma estrutura chamada RenderInfo com dados como ponto de interseção, a normal, um ponteiro para o shader dessa superfície, etc. Preenchida essa estrutura, ela faz uma chamada ao método QueryColor da classe OSLRender. O objeto dessa classe é global, o que significa que está em uma posição de memória compartilhada pelas threads. Como a chamada do método QueryColor modifica o objeto, temos o problema do acesso simultâneo a esse método.

Primeira tentativa

A primeira coisa que eu tentei, foi proteger a chamada dessa função com semáforos. O BRL-CAD tem uma biblioteca interna de utilitários, entre eles a implementação de um semáforo.

Basicamente, o que eu fiz foi o seguinte:

OSLRender oslr;
...
bu_semaphore_acquire(BU_SEM_SYSCALL);
oslr.QueryColor();
bu_semaphore_release(BU_SEM_SYSCALL);

BU_SEM_SYSCALL na verdade é uma macro representando um inteiro. Quando uma thread chama a função bu semaphore acquire(BU_SEM_SYSCALL), ela bloqueia o acesso a outras threads para o mesmo inteiro.

Com isso o código passou a funcionar para mais de uma thread! Porém, fui conservador demais ao bloquear a chamada do QueryColor. Creio que 90% do trabalho da função osl_render seja feito lá. Logo, colocando um semáforo nessa função, a paralelização vai por água abaixo. De fato, fiz testes com 1 e 2 processadores e os tempos para renderizar uma imagem de teste foram exatamente os mesmos.

Segunda tentativa

Fui estudar um pouco mais o código do OSL e descobri que ele tem suporte a aplicativos multithread. Basicamente, cada thread deve manter uma estrutura local chamada ThreadInfo. Ao chamar o sistema de shaders do OSL, passamos essa informação, que será modificada.

Gambiarra!

Para manter uma ThreadInfo para cada thread, precisamos identificá-las de alguma forma. O problema é que não temos acsso a essa informação diretamente. O único jeito que encontrei sem ter que mudar a estrutura interna do raytracer do BRL-CAD foi a seguinte: toda vez que a função osl_render é chamada, uma estrutura chamada struct application é passada. Essa estrutura tem um campo chamado resource, que é um ponteiro para um pedaço de memória exclusivo para cada thread.

Esse pedaço é alocado na criação das threads e não muda ao longo de suas vidas. Assim, usar o endereço desse pedaço de memória é um jeito (feio) de identificá-las. Outro problema é que, embora haja uma função de inicialização do shader (chamado osl_setup), as threads são criadas após isso, o que me obriga a inicializá-las na função osl_render mesmo.

A solução na qual cheguei foi a seguinte:

/* Every time a thread reaches osl_render for the 
   first time, we save the address of their own
   buffers, which is an ugly way to identify them */
std::vector<struct resource *> visited_addrs;
/* Holds information about the context necessary to
   correctly execute a shader */
std::vector<void *> thread_infos;

int osl_render(struct application *ap, ...){
    ...
    bu_semaphore_acquire(BU_SEM_SYSCALL);
    
    /* Check if it is the first time this thread is
       calling this function */
    bool visited = false;
    for(size_t i = 0; i < visited_addrs.size(); i++){
        if(ap->a_resource == visited_addrs[i]){
            visited = true;
            thread_info = thread_infos[i];
            break;
        }
    } 
    if(!visited){
        visited_addrs.push_back(ap->a_resource);
        /* Get thread specific information from 
           OSLRender system */
        thread_info = oslr->CreateThreadInfo();
        thread_infos.push_back(thread_info);
    }
    bu_semaphore_release(BU_SEM_SYSCALL);
    ...
}

Classe thread-safe

Com as ThreadInfo's, podemos nos despreocupar com o acesso concorrente ao objeto da classe OSLRender. Para ter uma garantia de que não será feita escrita nesse objeto, podemos declarar os métodos como const, como no exemplo abaixo:

struct A {
    int x;
    // OK
    void set_x(int _x) { x = _x; };
    // ERRO, tentando modificar membro do objeto 
    // num método const
    void const_set_x(int _x) const { x = _x; };
};

Testes

Com essa nova implementação, os tempos melhoraram bastante com mais processadores. A tabela abaixo mostra os tempos de execução para renderizar uma cena de exemplo, usando até 4 processadores.

Gráfico 1: Tempo x no. de processadores

Os ganhos foram praticamente lineares, sendo que para 4 processadores a medida de paralelização foi de 1.25 (se fosse perfeitamente linear seria 1.0).

Próximos passsos

Falei com meu mentor sobre o desenvolvimento daquele modo de visualização incremental, mas parece que ele não gostou muito da ideia. Isso exigiria uma mudança bastante grande no sistema, pois o renderizador do BRL-CAD foi desenvolvido para ray-tracing e não path-tracing.

Por enquanto estou estudando maneiras de adaptar o código OSL para suportar ray-tracing, mas não sei se isso é viável. A propósito, eu ainda confundo bastante os termos ray-tracing, path-tracing, photon mapping e todos esses algoritmos de iluminação e pretendo em breve escrever um post cobrindo esses tópicos bem por cima.

Anúncios

GDB Gnu Debugger

Julho 17, 2011

Acho que fui apresentado ao GDB logo que aprendi a programar, mas nunca me acostumei a usá-lo. Talvez por causa da maratona, sempre fui adepto do debug com printf, exceto quando há falha de segmentação. Mas mesmo nesse caso eu acho o valgrind mais eficaz.

Entretanto, dizem que essa é uma má prática e que o GDB pode ser uma ferramenta essencial na hora de encontrar causas de bugs. Pessoalmente, não me lembro de nenhuma vez em que o GDB me ajudou a achar um bug que eu não teria achado de outra forma, tão ou mais rapidamente quanto.

Mas isso deve ser porque eu não aprendi a usá-lo direito :P Em uma nova tentativa de melhorar meu uso com essa ferramenta, decidi organizar informações em um post de forma que eu possa vir a aprender novas funcionalidades.

Exemplo Básico

Para poder rodar código com GDB, devemos passar a flag -g ou -ggdb.

Dado um exemplo teste.c

#include <stdlib.h>
int main (int argc, char **argv){

    int i, v[5], n = atoi(argv[1]);

    for(i=n; i>=0; i++)
        v[i] = i;

    return 0;
}

Podemos compilá-lo com:

$ gcc teste.c -g

Gerando o executável ./a.out. Aí invocamos o gdb com o comando:

$ gdb ./a.out

Sem os possíveis argumentos do programa. Esses argumentos serão passados quando formos rodar o programa dentro do gdb, através do comando run.

(gdb) run 4

O programa acima tem um bug que causará falha de segmentação (ou algum erro relacionado ao acesso indevido de memória). No meu caso a saída foi:

Starting program: /home/kunigami/workspace/c/a.out 

Program received signal SIGSEGV, Segmentation fault.
0x00000000004004d9 in main () at teste.c:6
6	        v[i] = i;

Ele acusa o erro na linha 6 e mostra o trecho de código correspondente. Nesse caso, podemos analisar o código com a função print (ou simplesmente ‘p’) do gdb:

(gdb) p v[i]
Cannot access memory at address 0x7ffffffff000
(gdb) p i
$2 = 832

Vemos, ao tentar imprimir v[i], que tal posição da memória é inválida e isso porque o índice i está além das 5 posições alocadas.

Breakpoints

breakpoint (b) — coloca um marcador em uma linha do código para que o programa pare. A partir daí podemos executar o código passo-a-passo para uma análise mais detalhada.

Há varias maneiras de se colocar esse marcador. Uma é você dizer diretamente o número da linha, como por exemplo:

(gdb) b 123 

Colocará um breakpoint na linha 123. Se há vários arquivos fontes que formam o executável, dá pra especificar em qual arquivo quer se colocar o breakpoint. Se por exemplo quisermos colocar no código foo.c, fazemos:

(gdb) b foo.c:123

Também dá para passar o nome da função.

(gdb) b foo.c:bar

Para listar os breakpoints setados,

info b

Associado a cada breakpoint está um identificador, que pode ser usado se quisermos remover um dado breakpoint através de

delete <id do breakpoint>

Breakpoints condicionais

Se uma dada função é chamada muitas vezes e sabemos em quais condições um erro acontecerá, poderíamos fazer a seguinte gambiarra no código:

...
if(<condição>){
/* linha L do código */
}
...

E colocar o breakpoint na linha L. Ao invés disso, podemos colocar um breakpoint condicional. Dado o id do breakpoint, a sintaxe é a seguinte:

condition <id do breakpoint> <condição>

A condição é exatamente a mesma que você colocaria dentro do if.

Andando pelo código

Uma vez que o breakpoint fez a execução do programa parar, temos vários comandos para analisar detalhadamente as próximas linhas que serão executadas.

* step (s) — vai para a próxima instrução. Se for uma função, entra nela (step into).

* next (n) — vai para a próxima instrução. Se for uma função, não entra nela (step over).

* continue (c) — continua a execução normal do código (irá para a próxima vez que encontrar um breakpoint)

* finish (sem atalho) — sai da função atual (step out)

GDB e multi-threading

Quando lidamos com código multi-thread, o GDB começa a fazer mais falta. Vamos analisar o seguinte exemplo, adaptado daqui, usando pthreads:

threads.c

#include <stdio.h>
#include <stdlib.h>
#include <pthread.h>

void *foo(void *id)
{
    int i, s = 0;
    for(i = 0; i < 10000000; i++)
        s = (s + i) % 1234567;
    printf("[%d] s = %d\n", *((int*)id), s);
}

main()
{
     pthread_t t1, t2;
     int  id1 = 1, id2 = 2;

     /* Create independent threads each of which
        will execute function */
     (void) pthread_create(&t1, NULL, foo, (void*) &id1);
     (void) pthread_create(&t2, NULL, foo, (void*) &id2);

     /* Wait till threads are complete before main
        continues. Unless we wait we run the risk
        of executing an exit which will terminate   
        the process and all threads before the 
        threads have completed. */
     pthread_join(t1, NULL);
     pthread_join(t2, NULL); 

     exit(0);
}

Nesse código a thread inicial cria duas novas threads e bloqueia. Essas duas threads chamam a função foo, onde ficam em um loop. Quando ambas saem da função, a thread inicial desbloqueia e termina o código. Para compilá-lo,

gcc threads.c -lpthread -g

Vamos rodar o gdb e colocar um breakpoint na função foo, mas apenas para uma das threads. Podemos usar o breakpoint condicional.

(gdb) b foo
Breakpoint 1 at 0x400680: file teste.c, line 7.
(gdb) info b
Num     Type           Disp Enb Address            What
1       breakpoint     keep y   0x0000000000400680 in foo at teste.c:7
(gdb) condition 1 (*(int*)id) == 1

Aí rodamos o código

(gdb) run
Starting program: /home/kunigami/workspace/c/threads/a.out 
[Thread debugging using libthread_db enabled]
[New Thread 0x7ffff783c700 (LWP 30676)]
[Switching to Thread 0x7ffff783c700 (LWP 30676)]

Breakpoint 1, foo (id=0x7fffffffe2fc) at teste.c:7
7	    int i, s = 0;

Aqui ele avisa que trocou para outra thread (“Switching to Thread….“). Podemos visualizar as threads através do comando info threads.

(gdb) info threads
[New Thread 0x7ffff703b700 (LWP 30677)]
  3 Thread 0x7ffff703b700 (LWP 30677)  clone () at ../sysdeps/unix/sysv/linux/x86_64/clone.S:84
* 2 Thread 0x7ffff783c700 (LWP 30676)  foo (id=0x7fffffffe2fc) at teste.c:7
  1 Thread 0x7ffff7fcf700 (LWP 30673)  clone () at ../sysdeps/unix/sysv/linux/x86_64/clone.S:84

O que a tabela acima mostra é que a thread com id 1 é a thread inicial que está sendo clonada para criar a thread 3. A thread 2 está na função foo, onde ativou o breakpoint e é ela quem o GDB está analisando, indicado pelo asterisco na frente do id. Note que essa tabela provavelmente vai ser diferente a cada execução. Depende de como o sistema operacional está escalonando a execução das tarefas, por exemplo.

Vamos dar uns steps até que a variável i do loop valha, por exemplo, 3.

(gdb) p i
$2 = 3

Agora, vamos trocar para a thread 3 através do comando thread 3 e imprimir o valor de i.

(gdb) thread 3 
[Switching to thread 3 (Thread 0x7ffff703b700 (LWP 30779))]
#0  foo (id=0x7fffffffe2f8) at teste.c:8
8	    for(i = 0; i < 100000000; i++)
(gdb) p i
$7 = 2

Novamente, note que o estágio em que a thread 3 se encontra dependerá muito. Se ela ainda não tiver entrado no loop, dê alguns steps.

Conclusão

A possibilidade de se analisar várias threads ao mesmo tempo torna o GDB uma ferramenta muito mais atraente para se debugar código multi-thread do que o uso de printf‘s.

Referências

[1] http://www.delorie.com/gnu/docs/gdb/gdb_25.html


Biblioteca OpenMP

Janeiro 7, 2011

Depois de estudar um pouco a biblioteca MPI e aprender suas principais funções, decidi investir um tempo na biblioteca OpenMP. Como eu já tinha comentado em outro post, a diferença principal entre as duas é que o MPI usa o paradigma de memória distribuída enquanto no OpenMP a memória é compartilhada.

Dispor de memória compartilhada pode facilitar bastante o trabalho. Um exemplo é a divisão balanceada das iterações entre os processadores e com trabalho totalmente distribuído (no referido post a distribuição era feita através de um processador mestre). Com o MPI chegamos a uma solução com comunicação não-bloqueante, mas que executava um número de iterações maior do que o previsto.

Com o OpenMP podemos manter um pool central de iterações. Conforme cada processador termine de executar seu bloco de iterações, ele adquire exclusividade sobre esse pool, atualiza o número de iterações restantes e libera o recurso.

Fiz a simulação de uma iteração usando a função sleep com tempos aleatórios.

Geração de números aleatórios

Em geral, um requerimento para reproducibilidade de experimentos em que se usam números aleatórios é que toda vez que executamos o código, a sequência gerada seja a mesma. Para isso a biblioteca padrão stdlib oferece a função srand, que recebe como parâmetro uma semente. Se a semente for fixada, a sequência terá essa propriedade.

Quando estamos trabalhando com paralelização de código, uma propriedade desejável é que a união das sequências geradas em paralelo correspondam à sequência que seria gerada sequencialmente. Por exemplo, se com um processo geramos 6 números aleatórios entre 1 e 10, {6, 6, 8, 9, 7, 9}, então essa mesma sequência deveria ser gerada para 2 processos. Poderia ser algo como {6, 8, 7} e {6, 9, 9} (na verdade depende da ordem em que cada processador invocou o gerador – o importante é que a união das sequências seja igual à sequência original).

Com isso, os resultados obtidos serão essencialmente os mesmos entre a versão sequencial e a versão paralela, restando comparar apenas o tempo de execução.

O problema do gerador da biblioteca padrão C é que não temos acesso ao estado do gerador. Assim, não dá para garantir que quando cada processador gerar um número aleatório, o gerador esteja com o mesmo estado em todos os processadores. A biblioteca Boost possui uma implementação de geradores de números aleatórios, com a vantagem de que o gerador é um objeto ao qual temos acesso.

Assim, podemos colocá-lo como variável compartilhada entre os processadores e colocar a geração do número aleatório dentro de uma região crítica do Open MP, que é uma região onde apenas um processador pode executar por vez (quem tentar acessar fica bloqueado até que o processador na região crítica termine). Com isso garantimos que o estado do gerador está sempre atualizado para todos os processadores. Essa abordagem pode ser conferida na seção de código abaixo.

Resultados computacionais

Rodei o código abaixo em uma máquina quad-core e obtive os seguintes tempos de execução (medidos com a função time do linux), em função do número de threads utilizadas:

Tempos de execução

Observamos que a aceleração é quase linear. A terceira coluna da tabela mostra um tempo “normalizado” crescente com o aumento do número de processadores, mas isso se deve ao trecho não paralelizado do código. De forma geral, a economia de tempo foi bastante satisfatória. Porém, trata-se de um exemplo muito simples. Pretendo ainda aplicar essa estratégia em uma implementação mais complexa para ver os ganhos obtidos.

Código

#include <stdio.h>
#include <iostream>
#include <algorithm>

#include <boost/random/mersenne_twister.hpp>
#include <boost/random/uniform_int.hpp>
#include <boost/random/variate_generator.hpp>
#include <omp.h>

/* A otenção do número aleatório deve ser protegido */
int random(boost::mt19937 &gen, int lb, int ub){
    int val;
    #pragma omp critical
    {
        boost::uniform_int< > dist(lb, ub);
        boost::variate_generator<boost::mt19937&, 
                                 boost::uniform_int<> > r(gen, dist);
        val = r();
    }
    return val;
}
int iteration(boost::mt19937 &gen){
    int rtime = random(gen, 1, 10);
    sleep(rtime);
}
/* A atualização das iterações restantes deve ser 
 * uma area protegida */
bool get_iteration_block(int *iter){
    bool stat;
    #pragma omp critical
    {
        if (*iter <= 0){
            stat = false;
        }
        else {
            *iter = *iter - 1;
            stat = true;
        }
    }
    return stat;
}
int main (){

    boost::mt19937 gen;
    gen.seed(0);

    int iter = 20;  // numero inicial de iteracoes
    #pragma omp parallel shared(iter, gen) num_threads(2) 
    {
        while (get_iteration_block(&iter))
            iteration(gen);         
    }
    return 0;
}

Para compilar código OpenMP no linux podemos usar o gcc com a flag -fopenmp

g++ -fopenmp <arquivo.cpp>

No código acima, o gerador mt19937 é baseado no algoritmo gerador de números pseudo-aleatórios Mersenne Twister.


Ordenação Bitônica Paralela

Dezembro 24, 2010

No último post apresentei o algoritmo de ordenação bitônica. Aqui apresentarei uma versão paralela da ordenação bitônica, que na verdade é um híbrido da ordenação bitônica com algoritmo de ordenação usual como o quicksort. Para simplificar, vamos considerar novamente que o tamanho do vetor é uma potência de 2, assim como o número de processadores. Rotulamos os processadores de 0 a p-1. Dividimos o vetor de entrada em partes iguais de tamanho n/p. Vamos definir m = n/p para facilitar a notação. Assim, o processador 0 contém os elementos 0 a m-1, o processador 1 os elementos de m a 2m-1 e assim por diante.

Divisão Bitônica

Primeiramente, vamos relembrar como era feita a divisão bitônica. Dada uma sequência bitônica, consideramos os pares (0, n/2), (1, n/2+1), … (n/2-1, n-1) para gerar as sequências A1 e A2.

Os elementos 0 a m-1, sob posse do processador 0, devem formar pares com os elementos de n/2 a m+n/2-1, que estão sob posse do processador p/2.

Vamos considerar um exemplo com p=8. Podemos abstrair o tamanho do vetor n, considerando que cada processador tem um bloco de tamanho m. Numeremos os blocos de acordo com o processador. Assim, para gerar as sequências A1 e A2 na divisão bitônica, devemos parear o bloco 0 com o bloco 4, o bloco 1 com o 5, o 2 com o 6 e o bloco 3 com o 7, conforme a Figura abaixo.

Divisão bitônica de um vetor com 8 blocos.

Na recursão, dividiremos o vetor em duas metades, uma do bloco 0 ao bloco 3 e a outra do bloco 4 ao bloco 7. Para fazer a divisão bitônica da primeira metade, pareamos os blocos 0 e 2 e os blocos 1 e 3. Na segunda metade o pareamento é análogo. No último nível da recursão formamos os pares com os adjacentes, ou seja, (0,1), (2,3), (4,5) e (6,7).

Observe que para gerar a subsequência correspondente ao processador 0, decorrente da primeira divisão bitônica, só precisamos da subsequência do processador 4, ficando independente do restante das outras subsequências. Na segunda divisão, o processador 0 só precisou da informação em 2 e na terceira apenas da informação em 1. Fica fácil ver que com trocas de mensagens entre pares a cada divisão bitônica, podemos obter a sequência ordenada no final.

Para determinar qual o par correto de cada processador em cada divisão bitônica, podemos usar a representação binária do número de cada processador. É possível mostrar que na primeira divisão bitônica, pareamos os processadores cujo rótulo difere apenas no bit mais significativo. Por exemplo, (000 e 100), (001 e 101), (010 e 110) e (011 e 111). Na segunda divisão, aqueles que diferem apenas no segundo bit mais significativo, ou seja, (000 e 010), (001 e 011), (100 e 110) e (101 e 111). O procedimento é análogo parra os níveis subsequentes.

Gerando uma sequência bitônica

Como no caso sequencial, partimos de uma sequência bitônica. Portanto, precisamos transformar uma sequência arbitrária em bitônica antes de aplicar o método acima. A ideia geral era, em um dado nível, quebrar a sequência em 2, ordenar a primeira metade em ordem crescente e a outra metade em descrescente, juntar e aplicar a divisão bitônica.

Podemos transformar a recursão em um processo iterativo. Inicialmente ordenamos conjuntos de blocos consecutivos de tamanho 2. Assim, por exemplo, para p = 8, devemos ordenar os pares de blocos (0,1) e (4, 5) em ordem crescente e os pares (2, 3) e (6, 7) em ordem descrescente. Na próxima iteração ordenamos o bloco (0,1,2,3) em ordem crescente e (4,5,6,7) em ordem descrescente. Finalmente, na terceira iteração ordenamos (0,1,2,3,4,5,6,7) em ordem crescente. A Tabela abaixo sumariza esse procedimento, onde l é o tamanho dos conjuntos de blocos.


Estratégia bottom-up para obter sequências bitônicas.

Podemos determinar quando ordernar uma sequência em ordem crescente ou decrescente, novamente de acordo com a representação binária. No primeiro nível, os processadores com bit menos significativo 0 ordenam em ordem crescente e os com bit menos significativo 1, em ordem decrescente. No segundo nível, o segundo bit mais significativo é usado para determinar a ordem.

Bitonic split vs. Merge split

Observe que em um certo estágio das divisões bitônicas, a sequência em questão vai estar toda contida em um só processador. A princípio poderíamos continuar as divisões para ordenar esse trecho, mas como não podemos mais explorar o paralelismo, é mais viável usar um algoritmo O(n log n) ao invés do O(n log^2 n).

Dá para fazer algo melhor ainda! Se mantivermos as sequências em cada processador sempre ordenadas, quando chegarmos ao ponto em que a sequência está contida em um só processador, não precisaremos fazer nada!

Para trabalharmos com as sequências ordenadas dentro de cada processador, precisamos modificar o algoritmo que gera a sequência A1 e A2. Ao gerar essas sequências entre um par de processadores, ao invés de receber uma sequência crescente e outra decrescente, por hipótese recebemos duas sequências crescentes. Aí então executamos a etapa de intercalação do merge sort e distribuímos a metade menor para A1 e a metade maior para A2. Observe que os elementos de A1 e A2 são os mesmos se tivéssemos executado a divisão bitônica, só mudando a ordem deles. Para podermos usar a hipótese das sequências ordenadas, basta que cada processador ordene sua sequência inicialmente.

Por exemplo, considere a sequência bitônica {2, 4, 6, 8, 7, 5, 3, 1} e 2 processadores A e B, sendo que {2, 4, 6, 8} está no processador A e {7, 5, 3, 1} no processador B. A divisão usual resultaria em {2, 4, 3, 1} e {7, 5, 6, 8}. A proposta é que A tivesse {2, 4, 6, 8}, B tivesse {1, 3, 5, 7} (ordenado!) e a divisão resultasse em {1, 2, 3, 4} e {5, 6, 7, 8}.

Complexidade

Como foi descrito, cada processador deve ordenar sua sequência localmente. Usando um algoritmo eficiente como o quicksort, obtemos a complexidade O(m log m). Já a ordenação bitônica consiste de O(log p) estágios, cada qual usado para obter uma sequência bitônica. Para tanto usamos O(log p) trocas de mensagens entre um dado processador e seus pares. Tanto a troca de mensagens quanto a operação de intercalação é O(m). Assim, a complexidade total do algoritmo pode ser dada por .

Código

Decidi implementar a versão paralela usando C++ e MPI. Sacrifiquei um pouco de legibilidade para o código ficar menor. Esse template que estou usando tem pouco espaço horizontal, o que é bem ruim para postar códigos :(

#include <algorithm>
#include <cmath>
#include <cstdio>
#include "mpi.h"

/* Faz a intercalacao entre v e u e guarda em r */
void Merge(int *v, int *u, int n, int *t){
    int i=0, j=0, k=0;
    while(i < n && j < n)
        t[k++] = (v[i] < u[j]) ? v[i++] : u[j++];
    while(i < n) t[k++] = v[i++];
    while(j < n) t[k++] = u[j++];
}
/* Ordena conjuntos de blocos de tamanho set_size */
void BitonicSplit(int myrank, int *v, int n, 
                  int set_size, bool inc){
    
    int u[n], t[2*n];
    int set_dim = log(set_size)/M_LN2 + 1e-8;
    int eor_bit = 1 << (set_dim - 1);
    for(int i = 0; i < set_dim; i++){
        int pair = myrank ^ eor_bit;
        MPI_Status st;
        MPI_Send(v, n, MPI_INT, pair, 0,
                 MPI_COMM_WORLD);
        MPI_Recv(u, n, MPI_INT, pair, 0, 
                 MPI_COMM_WORLD, &st);

        Merge(v, u, n, t);

        int offset = ((inc && myrank < pair) || 
                      (!inc && myrank > pair)) ? 0 : n;
        for(int j=0; j<n; j++)
            v[j] = t[j+offset];
        eor_bit >>= 1;
    }   
}
void BitonicSort(int myrank, int np, int *v, 
                 int n){
    for(int ssize = 2; ssize <= np; ssize *= 2)
        BitonicSplit(myrank, v, n, ssize, 
                     !(myrank & ssize));
}
int main (int argc, char *argv[]){

    int myrank, np;
 
    MPI_Init(&argc, &argv);
    MPI_Comm_size(MPI_COMM_WORLD, &np);
    MPI_Comm_rank(MPI_COMM_WORLD, &myrank);

    /* Para simplificar, considere que a distribuição 
       dos dados já foi feita p/ cada processador */
    int t[ ] = {9, 12, 16, 23, 26, 39, 42, 61,
                43, 17, 14, 13, 12, 7, 6, 5};

    int n = 16/np;
    int *v = t + myrank*n;

    /* Ordena localmente */
    std::sort(v, v+n);

    BitonicSort(myrank, np, v, n);

    /* Cada processador tem a parcela correspondente
       do vetor ordenado total */
    printf("[%d]", myrank);
    for(int i=0; i<n; i++)
        printf(" %d", v[i]);
    printf("\n");

    MPI_Finalize();

    return 0;
}

Ordenação bitônica

Dezembro 17, 2010

Aprendi vários algoritmos de ordenação nos cursos da faculdade e na maratona. Entre eles estão o Bubble Sort, o Insertion Sort, o Selection Sort, o Heap Sort, Merge Sort e o Quick Sort. Porém, um no qual eu não tinha ouvido falar é o Bitonic Sort ou Ordenação Bitônica (em tradução livre).

Uma sequência monotonicamente crescente é equivalente a uma sequência não-descrescente, ou seja, seus elementos nunca diminuem, mas podem ficar iguais. A definição de monotonicamente decrescente é análoga. Já uma sequência bitônica, é uma sequência monotonicamente crescente seguida de uma monotonicamente decrescente (ou uma rotação cíclica dessa sequência).

Assim, por exemplo, {1, 2, 2, 3, 7, 6, 4} é uma sequência bitônica, formada pela sequência monotonicamente crescente {1, 2, 2, 3} e pela monotonicamente decrescente {7, 6, 4}. Como {3, 7, 6, 4, 1, 2, 2} é uma rotação cíclica da sequência original, então ela é também bitônica.

 

Figura 1: Exemplo de sequência bitônica

Para facilitar a descrição do algoritmo, vamos supor que o tamanho da sequência é uma pontência de 2.

Bitonic Split

Bitonic Split ou Divisão Bitônica (novamente, uma tradução livre), consiste em dividir a sequência bitônica A de tamanho n em duas novas sequências bitônicas, A1 e A2, de tamanho n/2, com a propriedade adicional de que todos os elementos de A1 são menores ou iguais a quaisquer elementos de A2.

Para isso, criamos a sequência A1 e A2 da seguinte forma:

 

Figura 2: Ilustração da divisão bitônica da Figura 1. A1 está destacado em azul, A2 em verde.

Por exemplo, se o vetor for {1, 5, 6, 9, 8, 7, 3, 0}, teremos A1 = {1, 5, 3, 0} e A2 = {8, 7, 6, 9}. Podemos fazer essa divisão recursivamente, até que restem sequências de apenas um elemento. Essa estratégia está esquematizada na figura a seguir:

Observe que pela propriedade dos elementos da subárvore da esquerda serem sempre menores do que os da direita, no nível das folhas o vetor estará ordenado. Se estivéssemos fazendo essas divisões no próprio vetor teríamos algo do tipo:

Sequência em cada etapa do algoritmo.

Não podemos aplicar esse algoritmo a qualquer subsequência porque usamos a hipótese de que a sequência de entrada era bitônica. Portanto, devemos transformar a sequência de entrada em uma sequência bitônica. Para tanto usaremos um algoritmo baseado no seguinte argumento indutivo: Suponha que sabemos transformar em bitônica cada uma das metades da sequência original. Como, a partir disso, transformar a sequência original em bitônica?

Como cada metade é bitônica, usamos a divisão bitônica para ordená-las. A primeira metade ordenamos em ordem não-decrescente e a segunda metade em ordem não-crescente (é só trocar A1 por A2 no algoritmo). Pronto! A sequência original passou a ser bitônica.

Só resta mostrar a base da indução, que é quando temos 4 elementos, por exemplo {a1, a2, a3, a4}. Nesse caso, quebramos em {a1, a2} e {a3, a4}. Se a1 > a2, trocamos esses elementos de posição. Se a3 < a4, também trocamos eles de posição. Temos então uma sequência bitônica.

Complexidade

Vamos analisar primeiro a complexidade de ordenar uma sequência bitônica com a divisão bitônica recursiva. A cada nível quebramos o array em 2, com um processamento linear para construir A1 e A2. Não há processamento após chamar a recursão. Temos então O(log n) níveis gastando O(n) cada. Portanto, temos O(n log n).

Agora a complexidade do algoritmo para transformar uma sequência qualquer em bitônica. Em cada nó quebramos o array em 2 e após retornar da recursão, é preciso usar o algoritmo anterior para ordenar cada metade, gastando O(n/2 log n/2 + n/2 log n/2) = O(n log n). É fácil ver que em cada nível gastamos O(n log n) e temos O(log n) níveis. Logo, a complexidade total do algoritmo é .

Esse algoritmo poderia ser melhorado, observando que é possível ordenar uma sequência bitônica em O(n) (exercício), o que resultaria na complexidade ótima de ordenação que é O(n log n). Porém, esse algoritmo na forma original é interessante para ser paralelizado. Porém, esse é um assunto que discutiremos em um futuro post.

Código

Eu ia implementar esse algoritmo em python, mas como na Wikipedia já existe essa implementação, decidi fazer uma em C++. Para simplificar, o código supõe que o tamanho do vetor é uma potência de 2. Se não for o caso, basta encontrar o maior elemento do vetor e inserir várias cópias dele no final até completar um tamanho que satisfaça essa propriedade. Depois de ordenado, é só jogar fora os últimos elementos.

#include <stdio.h>
#include <algorithm>

// Ordena sequência bitonica
void BitonicSplit(int *v, int n, bool inc){
    if (n <= 1) return;

    for(int i=0; i<n/2; i++)
        if ((inc && v[i] > v[i+n/2])||
            (!inc && v[i] < v[i+n/2]))
            std::swap(v[i], v[i+n/2]);

    BitonicSplit(v, n/2, inc);
    BitonicSplit(v + n/2, n/2, inc);
}
// Ordena sequência arbitrária
void BitonicSort(int *v, int n, bool inc){
    if (n <= 1) return;

    BitonicSort(v, n/2, inc);        // Ordena cresc.
    BitonicSort(v + n/2, n/2, !inc); // Ordena decresc.
    BitonicSplit(v, n, inc);
}

int main (){

    int v[ ] = {5, 1, 6, 9, 8, 7, 0, 3};
    int n = 8;

    BitonicSort(v, n, true);

    for(int i=0; i<n; i++)
        printf("%d ", v[i]);
    printf("\n");
}

Referência

[1] Parallel Programming with MPI – Peter S. Pacheco
[2] Wikipedia – Bitonic Sorter


Comunicação Não Bloqueante

Dezembro 10, 2010

A forma mais comum de comunicação entre dois processos no MPI é através de Sends e Receives, implementados pelas funções MPI_Send e MPI_Recv, respectivamente. Em teoria essas funções são bloqueantes, no sentido de que ao enviar uma mensagem com MPI_Send, o programa para sua execução até que o processador para o qual ele enviou a mensagem receba sua mensagem com um MPI_Recv.

Imagine um cenário simples com 2 processos que queiram trocar mensagens entre si. Localmente, ambos executam um MPI_Send seguido de um MPI_Recv. O diagrama abaixo ilustra esse caso:


Diagrama ilustrando potencial deadlock

No caso em que o send e o receive são bloqueantes, teremos o chamado deadlock (no diagrama de dependência, ele é caracterizado por um ciclo dirigido), ou seja, os processadores ficam esperando infinitamente.


Exemplo de deadlock no dia-a-dia (foto recebida por email).

Felizmente, a maior parte das implementações do padrão MPI, provê um buffer para que o send não seja bloqueante. Ao chamar MPI_Send, a mensagem é copiada em um buffer e a execução do programa continua. Aí quando o destinatário chama MPI_Recv, ele lê a mensagem do buffer. Em implementações onde não há o esse buffer, deve-se usar MPI_Bsend, onde o buffer deve ser explicitamente alocado.

Porém, a chamada MPI_Recv continua bloqueante. Dependendo da aplicação, essa característica pode ser custosa, pois enquanto espera uma mensagem, o processador fica ocioso. Recentemente me deparei com uma situação que sofria desse problema.

Suponha três processos p1, p2, p3, sendo que p1 é o processo que gerencia a execução do algoritmo (leitura de entrada, repassar tarefa para p2 e p3, além de processar os resultados de p2 e p3). Ao processo p1 é dado um certo número de iterações que deverão executadas por p1, p2 e p3. Como cada iteração executada leva um tempo diferente, usamos a estratégia de distribuição por demanda.

Assim, p1 inicialmente envia um bloco de iterações a p2 e p3 para que estes executem. Enquanto isso, ele trata de executar também algumas iterações. Ao final dessa execução, ele comunica com p2 e p3 para verificar se eles já terminaram suas iterações. Se sim, envia mais um bloco de iterações para que eles executem. Caso contrário, volta para executar mais iterações, fazendo a verificação novamente, até que o número especificado de iterações seja executado.

Para fazer a comunicação com p2 (ou p3), o processo p1 pode chamar MPI_Recv, esperando que p2 retorne seus resultados. Depois, ele decide se ainda há iterações para serem executadas. Então, envia com o MPI_Send um ‘ok’, em caso positivo, ou um ‘stop’, em caso negativo. Por outro lado, o processo p2, ao terminar de executar seu bloco de iterações, envia seus resultados com o MPI_Send e espera a resposta de p1 para saber se deve parar ou continuar, através do MPI_Recv.

É muito provável que algum dos processos fique esperando os outros terminarem suas tarefas, enquanto poderia estar executando mais iterações. Nesse caso, operações de send e receive não bloqueantes seriam mais adequadas. Felizmente, a biblioteca MPI oferece essa alternativa.

Sends e Receives não bloqueantes

No MPI, as funções que implementam send e receive não bloqueantes são denominadas MPI_Isend e MPI_Irecv, respectivamente. Elas não bloqueiam, mas também não efetuam de fato a operação de send/receive, apenas iniciam o processo.

Ao chamar uma dessas funções, é retornado um objeto do tipo ‘request’, que indica o status do send ou do receive. Esse ‘request’ pode ser verificado basicamente de duas maneiras: através da função MPI_Wait ou da função MPI_Test. A primeira bloqueia a execução até que o send ou receive se complete. A segunda retorna uma flag que é true se a operação completou, ou false caso contrário.

O exemplo apresentado anteriormente pode ser adaptado para usar comunicação não bloqueantes. Agora, o processador p1 distribui o bloco de iterações e começa a executar as suas. Ao final, ele chama o MPI_Irecv para indicar que está esperando uma mensagem de p2 e p3. Aí então ele verifica através do MPI_Test se a sua mensagem já chegou. Em caso positivo, ele executa o MPI_Isend dizendo se é para parar ou continuar a execução. Caso não tenha chegado a mensagem, ele simplesmente executa mais um bloco de iterações e volta para verificar novamente.

No outro lado, o processo p2 (ou p3), ao terminar suas iterações, envia uma mensagem para p1 com MPI_Isend e então invoca MPI_Irecv para indicar que está esperando a resposta de p1. Em seguida, chama MPI_Test para determinar se a mensagem já chegou. Em caso negativo, volta a executar mais um bloco de iterações. Em caso positivo, verifica se a mensagem é para que ele pare de executar.

Experimentos computacionais

Implementei as duas estratégias mencionadas acima, simulando a execução de iteração com um sleep de um tempo aleatório entre 1 e 5 segundos. Fixei 100 iterações e variei os tamanhos de bloco em 1, 5, 10 e 20. A tabela abaixo sumariza os tempos de execução obtidos e os tempos de espera.

O tempo de espera é a soma do tempo de chamada de cada função de send ou receive. Mais especificamente, eu meço o tempo com MPI_Wtime antes e depois de chamar cada uma dessas funções, e somo a diferença desses tempos.

Para medir o tempo total, foi medido apenas no processador mestre, mas para garantir que a medida corresponde ao tempo do processador que terminou por último, usei a função MPI_Barrier, que bloqueia todo mundo que a chama enquanto houver que ainda não a chamou.

Tabela: Comparação de tempo de execução total e tempo ocioso total entre a estratégia bloqueante e a não bloqueante para diversos tamanhos de blocos.

Verificamos, como era esperado, que o uso de chamadas não bloqueantes resulta em tempo ocioso nulo. O tempo de execução total tende a aumentar com o tamanho dos blocos.

Observe que para blocos pequenos, a estratégia não-bloqueante se deu melhor, devido aos tempos de processador ocioso da estratégia bloqueante. Porém, para blocos maiores, a estratégia bloqueante é mais rápida. Uma explicação para isso é que na estratégia não-bloqueante, o número de iterações realizadas pode superar a quantidade estabelecida de iterações, pois quando o MPI_Test retorna falso, o processador vai executar outro bloco de iterações, mesmo que nesse momento o número de iterações já tenha sido atingido.

Da forma como o tempo de processador ocioso é medido, está incluindo o tempo de envio e recebimento das mensagens. Note então que esse tempo é irrelevante perto do tempo de execução do algoritmo. Assim, usar blocos de tamanho pequeno, em detrimento do aumento de trocas de mensagem, e comunicação não-bloqueante parece uma boa ideia.


MPI – Trocas de mensagens

Outubro 8, 2010

Nesse post vou comentar sobre alguns aspectos de comunicação entre nós, através da biblioteca MPI, a qual introduzi neste outro post.

Comunicação ponto a ponto

Como diz o próprio nome, a biblioteca MPI (interface de passagem de mensagens), tem como característica principal a comunicação entre nós através de mensagens. O modo mais simples de trocar mensagens, é através das funções de envio e recebimento de mensagens, respectivamente MPI_Send e MPI_Recv.

A ideia é muito próxima do envio de cartas do mundo real: o remetente deve especificar para a função MPI_Send, a mensagem e o destinatário. O destinatário, por sua vez, especifica para a função MPI_Recv o remetente de quem ele está esperando a carta. Um detalhe é que a função MPI_Recv é bloqueante, ou seja, o destinatário fica bloqueado até receber a mensagem que está esperando.

Broadcast

Se um processo p quiser mandar uma mensagem para todos os outros, pode utilizar a função MPI_Send para cada outro processo. Se o envio de cada uma dessas mensagens levar um tempo t e houver n processos, o tempo necessário para que todos os processos tenham a informação passada pelo processo p, será da ordem de t(n – 1). Além de não ser muito eficiente, ainda há o empecilho de o processo p estar concentrando todo o trabalho em si, o que não é interessante, uma vez que nosso objetivo é paralelizar o código.

Uma ideia é usar uma topologia em árvore para a distribuição da mensagem. Suponha que n = 8 e o processo 0 quer enviar a mensagem. A troca de mensagens segue o esquema da figura abaixo.


Topogia de árvore

Inicialmente o nó 0 envia a mensagem para o nó 1. Em um segundo momento, o nó 0 envia a mensagem para o nó 2, ao mesmo tempo em que o nó 1 envia a mensagem para o nó 3. Agora os nós 0, 1, 2 e 3 possuem a informação que o nó 0 queria transmitir. No terceiro estágio, o nó 0 envia para 0 4, o nó 1 para o 5, o 2 para o 6 e o 3 para o 7. Assim, todos os nós receberam a mensagem em 3 estágios. Não é difícil ver que o tempo gasto será proporcional à altura da árvore, que por sua vez é O(log n). Assim, o tempo total gasto é O(t log n). Na MPI, essa estratégia é implementada pela função MPI_Bcast.

Reduce

Imagine que paralelizamos um laço que fazia a soma de um vetor, distribuindo um pedaço para cada um dos processos. Depois de calculada a soma, é preciso juntar os resultados. No caso da soma, basta que cada nó envie seu resultado para o processo principal e que este os some para obter a resposta final. Assim como a versão simples do broadcast, esse esquema é proporcional ao número de nós.

Aqui, também é possível usarmos a topologia de árvore da seguinte forma: cada nó envia seu resultado para o seu predecessor na árvore. Pelo exemplo da figura anterior, o nó 0 e o nó 4 vão enviar para o nó 0. Na verdade não é preciso que o nó 0 envie o resultado para si mesmo, basta que ele receba a resposta do nó 4 e some com o valor que ele já tinha calculado. Assim, cada nó intermediário fará a soma dos seus dois filhos, de forma que no final o nó 0 conterá a soma total. Novamente conseguimos reduzir o tempo para O(t log n). Na MPI, essa abordagem é utilizada pela função MPI_Reduce.

Reduce + Broadcast

Ao realizar o reduce, o nó 0 conterá a soma total, mas os outros nós terão apenas somas parciais. Se quiséssemos que todos os nós tivessem a soma total, uma possibilidade seria o nó 0 fazer um broadcast para todos os outros nós, ou seja, um reduce seguido de um broadcast. Há um meio mais eficiente de fazer isso, através da topologia borboleta.

Topologia Borboleta

Para entender melhor esse esquema, consideremos a informação que o nó 0 envia (em verde). No primeiro estágio ele envia sua informação para o nó 4. Logo, os nós 0 e 4 têm a informação inicial do nó 0. No segundo estágio, o nó 0 envia sua informação para o nó 2, enquanto o nó 4 envia sua informação para o nó 6. Agora os nós 0, 2, 4 e 6 têm a informação do nó 0. Finalmente no estágio 3, o nó 0 envia para o nó 1, o 2 para o 3, o 4 para o 5 e o 6 para o 7. Assim, todos têm a informação inicial do nó 0. Seguindo a mesma lógica para os outros nós, concluímos que ao final do terceiro estágio, todos os nós terão informação dos outros nós. Se cada nó somar o seu conteúdo com o que receber, no final todos os nós terão as somas totais. A função MPI_Allreduce implementa essa abordagem.