Alcides Carneiro de Araújo Neto. Comparação Sequência-Família em GPU - PDF
Download "Alcides Carneiro de Araújo Neto. Comparação Sequência-Família em GPU"
Ana do Carmo da Mota Belo
1 Alcides Carneiro de Araújo Neto Comparação Sequência-Família em GPU Campo Grande MS 2014
2 Alcides Carneiro de Araújo Neto Comparação Sequência-Família em GPU Dissertação de Mestrado apresentada à Faculdade de Computação da Universidade Federal de Mato Grosso do Sul UFMS, como parte dos requisitos para obtenção do título de Mestre em Ciência da Computação. Universidade Federal de Mato Grosso do Sul UFMS Faculdade de Computação FACOM Orientadora: Prof. a Dr. a Nahri Balesdent Moreano Campo Grande MS 2014
3 Resumo Durante as últimas décadas o volume de informações biológicas em algumas bases de dados cresceu em um ritmo quase exponencial. Ferramentas como o HMMER podem encontrar sequências biológicas homólogas a uma família de sequências modelada estatisticamente por um profile HMM utilizando o algoritmo de Viterbi. Dada a complexidade quadrática desse algoritmo, esse procedimento pode consumir longos tempos de execução dependendo da quantidade de sequências, do tamanho do profile HMM e do hardware utilizado. Esse trabalho descreve o desenvolvimento de uma solução em GPU, de alto desempenho, para o problema de determinar se uma nova sequência biológica é homóloga a uma família de sequências conhecida. A solução implementada alcançou desempenho compatível ou superior ao HMMER. Palavras-chaves: Alinhamento sequência-família, algoritmo de Viterbi, HMMER, GPU, CUDA, acelerador. i
4 Abstract Over the past few decades the amount of biological data in some databases grew up in an almost exponential rate. Tools such as HMMER use the Viterbi algorithm to find biological sequences that are homologue to a family of sequences represented by a statistical model called profile HMM. Due to the quadratic time complexity of the Viterbi algorithm, this search procedure can demand long execution times depending on database size, profile HMM length, and hardware used. The purpose of this project is to design a high performance GPU solution for the problem of finding out if a new biological sequence is homologue to a known family of sequences. The implemented solution reached a performance compatible or superior to HMMER. Keywords: Sequence-profile alignment, Viterbi algorithm, HMMER, GPU, CUDA, accelerator. ii
5 Sumário 1 Introdução Motivação Objetivos do Trabalho Organização do Texto Plataforma de Computação Paralela Modelo de Programação CUDA Programa Exemplo Arquitetura das GPUs CUDA Hierarquia de Memórias Acesso Coalescido à Memória Global Sincronização em CUDA Sincronização de Threads Streams Escalabilidade do Modelo de Programação CUDA Comparação Sequência-Família Sequência Biológica e Família de Sequências Alinhamento Múltiplo de Sequências e Profile HMM Alinhamento Sequência-Família Ferramenta HMMER Algoritmo de Viterbi Algoritmo MSV Algoritmo SSV Obtenção do Alinhamento Ótimo Trabalhos Relacionados Soluções em FPGA iii
6 4.1.1 Soluções com Perda de Precisão Soluções sem Perda de Precisão Soluções em GPU Soluções para o Algoritmo de Viterbi Soluções com Perda de Precisão Soluções sem Perda de Precisão Soluções para o Algoritmo MSV Outros Trabalhos Soluções que Utilizam Extensões do Conjunto de Instruções Streaming SIMD Extensions Paralelismo de Dados no MSV Considerações Finais Solução em GPU para Comparação Sequência-Família Sequência de Filtros Modelagem dos Kernels Kernel MSV Kernel SSV Estruturas de Dados Probabilidades de Emissão Sequências Otimizações Sequência de Filtros Probabilidades de Emissão Sequências Memória Não-paginada Árvore de Máximos e Redução Loop Unrolling Transferência Host-GPU e Múltiplos Streams Tiling iv
7 6 Resultados Plataformas, Ferramentas e Dados Utilizados Solução Inicial sem Otimizações Solução com Memória Não-paginada Solução com Árvore de Máximos e Redução Otimizada Solução com Scores Representados com Números Inteiros Solução com Scores Representados com Números Naturais Solução com Probabilidades de Emissão na Memória Constante Solução com Probabilidades de Emissão na Memória de Textura Solução com Padding nas Probabilidades de Emissão Solução com Padding nas Sequências Solução com Loop Unrolling Solução com Múltiplos Streams Solução com Acesso Vetorizado às Probabilidades de Emissão Scores Representados com Números Reais e Acesso Vetorizado com Fator Scores Representados com Números Reais e Acesso Vetorizado com Fator Scores Representados com Números Naturais e Acesso Vetorizado com Fator Scores Representados com Números Naturais e Acesso Vetorizado com Fator Scores Representados com Números Naturais e Acesso Vetorizado com Fator Solução com Acesso Vetorizado às Sequências Solução com Tiling Tiling com Fator Tiling com Fator Tiling com Fator Solução Otimizada Final v
8 7 Conclusão Trabalhos Futuros Referências Bibliográficas 93 vi
9 Lista de Figuras 1.1 Crescimento da base de dados de proteínas UniProtKB/TrEMBL entre 1996 e Organização das threads no modelo de programação CUDA: um grid com 2 3 blocos de 2 2 threads Programa em CUDA que executa no host para soma de dois vetores Kernel CUDA para soma de dois vetores Arquitetura de uma GPU CUDA Hierarquia de memórias de uma GPU CUDA e unidades que possuem acesso às memórias (a) Acessos coalescidos à memória global e (b) acessos não coalescidos à memória global (a) Execução de programa com um único stream; (b) execução de programa com dois streams, sobrepondo operações e reduzindo o tempo total de execução Alinhamento múltiplo de cinco sequências de DNA Arquitetura de profile HMM multi-hit local alignment usada pelo HMMER3 no algoritmo de Viterbi Sequência de filtros da ferramenta HMMER Dependências de dados do algoritmo de Viterbi Arquitetura de profile HMM multi-hit ungapped local alignment usada pelo HMMER3 no algoritmo MSV Dependências de dados do algoritmo MSV Dependências de dados do algoritmo SSV Exemplo de alinhamento ótimo MSV obtido pelo algoritmo de Traceback Exploração do paralelismo de tarefas na análise de um conjunto de sequências Arquitetura de profile HMM com duplicação do núcleo do modelo, visando minimizar a perda de precisão causada pela remoção do estado Joining.. 34 vii
10 4.3 Solução em FPGA para algoritmo MSV com caminho crítico na árvore de máximos Técnica de tiling dividindo as linhas da matriz M em duas partições Determinação de âncoras e partições de uma linha da matriz de programação dinâmica para exploração de paralelismo de dados Registradores de 128 bits do SSE2 sendo usados para armazenar: (a) dois dados de 64 bits; (b) quatro dados de 32 bits; (c) oito dados de 16 bits; (d) 16 dados de 8 bits Exemplo de uma operação de soma do tipo SIMD Linhas da matriz M divididas em porções de até 16 células, para exploração de paralelismo de dados com instruções SIMD no algoritmo MSV Sequência de filtros da solução desenvolvida Exploração de paralelismo de tarefas (análise de sequências distintas por diferentes blocos) e exploração de paralelismo de dados (cálculo de células da matriz de programação dinâmica por threads) Código CUDA do kernel MSV Código CUDA do kernel SSV Comparação entre transferências de dados usando memória paginável e memória não-paginada Árvore de máximos para um vetor de tamanho Código CUDA otimizado para cálculo do máximo da linha i (redução) da matriz M do kernel MSV Exemplo de uso de templates no kernel MSV Código para a geração de múltiplos kernels MSV em tempo de compilação Exemplo de kernel SSV com otimização na redução para profile HMMs de tamanho Resultados do CUDA Device Query na GPU usada nos experimentos Tempos gastos com transferências entre host e GPU e execução do kernel MSV com memória paginável Tempos gastos com transferências entre host e GPU e execução do kernel MSV com memória não-paginada viii
11 6.4 Exemplo de caso de uso que usufrui da cache da memória de textura que é otimizada para localidade espacial 2D Determinação do tamanho do batch que proporciona velocidade máxima de transferência Análise visual da sobreposição da execução de kernels e transferências de memória quando usado streams no filtro SSV em GPU ix
12 Lista de Tabelas 2.1 Critério de alinhamento para transações coalescidas na memória global Trabalhos relacionados: plataforma utilizada, algoritmo implementado, formas de paralelismo exploradas e desempenho alcançado Famílias de sequências da base de dados Pfam utilizadas Estatísticas das sequências da base de dados UniProtKB/Swiss-Prot Tempos de execução da solução inicial em GPU e do HMMER Tempos de execução da solução inicial em GPU e do HMMER3.1b Speedups da solução inicial em GPU em comparação com o HMMER Speedups da solução inicial em GPU em comparação com o HMMER3.1b Tempos de execução e speedups da solução com redução otimizada em relação à solução inicial, para os kernels SSV e MSV Tempos de execução e speedups da solução com números inteiros em relação à solução inicial, para os kernels SSV e MSV Tempos de execução e speedups da solução com números naturais em relação à solução inicial, para os kernels SSV e MSV Tempos de execução e speedups da solução com memória constante em relação à solução inicial, para os kernels SSV e MSV Tempos de execução e speedups da solução com memória de textura em relação à solução inicial, para os kernels SSV e MSV Tempos de execução e speedups da solução com padding nas probabilidades de emissão em relação à solução inicial, para os kernels SSV e MSV Tempos de execução e speedups da solução com padding nas sequências em relação à solução inicial, para os kernels SSV e MSV Tempos de execução e speedups da solução com loop unrolling em relação à solução inicial, para os kernels SSV e MSV Tempos de execução e speedups da solução com números reais e acesso vetorizado às probabilidades de emissão com fator 2 em relação à solução inicial, para os kernels SSV e MSV x
13 6.16 Tempos de execução e speedups da solução com números reais e acesso vetorizado às probabilidades de emissão com fator 4 em relação à solução inicial, para os kernels SSV e MSV Tempos de execução e speedups da solução com números naturais e acesso vetorizado às probabilidades de emissão com fator 2 em relação à solução inicial, para os kernels SSV e MSV Tempos de execução e speedups da solução com números naturais e acesso vetorizado às probabilidades de emissão com fator 4 em relação à solução inicial, para os kernels SSV e MSV Tempos de execução e speedups da solução com números naturais e acesso vetorizado às probabilidades de emissão com fator 16 em relação à solução inicial, para os kernels SSV e MSV Tempos de execução e speedups da solução com acesso vetorizado às sequências em relação à solução inicial, para os kernels SSV e MSV Tempos de execução e speedups da solução com tiling com fator 2 em relação à solução inicial, para os kernels SSV e MSV Tempos de execução e speedups da solução com tiling com fator 4 em relação à solução inicial, para os kernels SSV e MSV Tempos de execução e speedups da solução com tiling com fator 8 em relação à solução inicial, para os kernels SSV e MSV Otimizações adotadas na solução final dos filtros SSV e MSV Tempos de execução da solução final em GPU e do HMMER Tempos de execução da solução final em GPU e do HMMER3.1b Speedup da solução final em GPU em comparação com o HMMER Speedup da solução final em GPU em comparação com o HMMER3.1b.. 88 xi
14 Lista de Algoritmos 3.1 Algoritmo de Viterbi Algoritmo MSV Algoritmo SSV Algoritmo de Traceback para obtenção do alinhamento ótimo MSV xii
15 1 Introdução Em Bioinformática, uma sequência biológica é modelada como uma cadeia finita de símbolos que carregam informação biológica. Por sua vez, uma família de sequências é um conjunto de sequências homólogas, isto é, as sequências pertencentes a uma família divergiram durante a evolução de tal forma que compartilham semelhanças estruturais e/ou funcionais [45]. Uma vez que, na natureza, novas sequências evoluem a partir de sequências preexistentes através de mutações, a análise computacional das mesmas é beneficiada disso. Muitas vezes é possível reconhecer semelhanças significativas entre uma nova sequência e uma sequência já conhecida sobre a qual sabe-se algo. Quando isso ocorre é possível transferir informações sobre a estrutura e/ou função da sequência já conhecida para a nova sequência. Assim, as duas sequências relacionadas são então ditas homólogas e a transferência de informações entre elas é feita por homologia [7]. Portanto, o problema da comparação sequência-família, isto é, determinar se uma nova sequência biológica é homóloga a uma família de sequências conhecida é de grande importância na Bioinformática [6]. 1.1 Motivação O volume de dados biológicos vem crescendo de maneira surpreendente. O Projeto Genoma Humano (Human Genome Project) [33], um programa de pesquisa internacional e colaborativo cujo objetivo era o mapeamento completo e a compreensão de todos os genes dos seres humanos, foi declarado completo em 2003, identificando cerca de 3 bilhões de símbolos do código genético humano. Outro exemplo é a base de dados de proteínas UniProtKB/TrEMBL [51], que apresentou um crescimento quase exponencial desde o ano 2000, praticamente dobrando de tamanho a cada dois anos. A Figura 1.1 mostra o crescimento do número de sequências dessa base de dados entre os anos de 1996 e A análise de grandes bases de dados de sequências biológicas e famílias de sequências exige, dentre outras técnicas, uma que permita resolver o problema da comparação sequência-família. O HMMER [9, 19], uma solução em software para computadores de arquitetura convencional, é uma das principais ferramentas utilizadas para este fim e baseia-se em um importante algoritmo denominado algoritmo de Viterbi [14]. Dada a complexidade de tempo quadrática do algoritmo de Viterbi, a análise de grandes bases de sequências e famílias utilizando o HMMER pode demandar longos 1
16 50M Número de sequências 40M 30M 20M 10M Ano Figura 1.1: Crescimento da base de dados de proteínas UniProtKB/TrEMBL entre 1996 e 2013 [51] tempos de execução. Com o rápido crescimento das bases de dados biológicos, esses tempos de execução tornam-se ainda mais críticos [6, 18, 42, 55, 56]. Em vista disso, surge a necessidade de soluções de alto desempenho, capazes de analisar grandes volumes de sequências e/ou famílias em pouco tempo, seja pela adoção de heurísticas ou pela exploração de paralelismo ou até mesmo por ambas [15]. Soluções paralelas baseadas em clusters e aceleradores em hardware implementados em FPGAs (Field-Programmable Gate Arrays) estão entre as propostas desenvolvidas que apresentaram speedups consideráveis se comparados com o HMMER. A ferramenta HMMER, a cada nova versão, também incorpora recursos com o objetivo de melhorar seu desempenho. Tais recursos podem permitir execução em cluster, execução multithread em um processador multi-core e exploração de paralelismo de dados utilizando extensões SIMD do conjunto de instruções de máquina (descritas na Seção 4.3). Mais recentemente, devido aos resultados promissores obtidos em outras áreas utilizando-se GPUs (Graphics Processing Units), pesquisadores voltaram suas atenções para o uso das mesmas como plataforma de computação paralela. Assim, algumas soluções paralelas baseadas em GPU foram propostas para o problema da comparação sequência-família, com o objetivo de obter um bom desempenho. A vantagem da utilização de GPUs está em sua arquitetura many-core que disponibiliza ao desenvolvedor uma grande quantidade de cores que continua a aumentar segundo a Lei de Moore [34]. Assim, problemas que podem ser modelados de forma a tirar vantagem desse crescente aumento na quantidade de cores possuem potencial para ganho de desempenho quando executados explorando paralelismo nesse tipo de plataforma. 2
17 1.2 Objetivos do Trabalho A proposta desse trabalho é a construção de uma solução em GPU para o problema da comparação sequência-família, isto é, determinar se uma nova sequência biológica é homóloga a uma família de sequências conhecida. O objetivo é obter uma solução de alto desempenho que permita a análise de grandes bases de dados biológicos em pouco tempo, utilizando como plataforma de execução uma GPU. Para tal, é necessário estudar os algoritmos utilizados para resolver o problema da comparação sequência-família, suas dependências de dados e como é possível explorar paralelismo na execução desses algoritmos. As abordagens das soluções já desenvolvidas nos variados dispositivos de alto desempenho devem ser estudas e uma possível adaptação das mesmas para GPU deve ser analisada. Também é necessário conhecer a plataforma de computação paralela utilizada, no caso uma GPU, isto é, suas características arquiteturais e seu modelo de programação. Pretende-se compreender como tirar o máximo de proveito dos recursos dessa plataforma e minimizar as desvantagens no uso da mesma, através de boas práticas de programação e otimização de código, com o objetivo de melhorar o desempenho da solução. Por fim, pretende-se realizar uma avaliação experimental de desempenho da solução desenvolvida em GPU, comparando-a com o HMMER e outros trabalhos. O objetivo é avaliar as otimizações aplicadas na solução em GPU e concluir a viabilidade de utilizá-las como plataforma de computação paralela para resolver o problema da comparação sequência-família e problemas similares da área de Bioinformática. 1.3 Organização do Texto O restante desse texto está organizado em seis capítulos. O Capítulo 2 descreve a utilização das GPUs como plataforma de computação paralela, sua arquitetura, hierarquia de memórias, o modelo de programação CUDA e os mecanismos de sincronização deste modelo. No Capítulo 3, conceitos básicos de Bioinformática necessários para a melhor compreensão do problema da comparação sequência-família, são introduzidos. O problema do alinhamento sequência-família é então formalizado. A ferramenta HMMER é descrita e os algoritmos de Viterbi, MSV (Multiple Segment Viterbi) e SSV (Single Segment Viterbi), utilizados para a comparação sequência-família, são apresentados e as dependências de dados presentes nos mesmos são examinadas com o intuito de encontrar possibilidades de exploração de paralelismo. 3
18 Em seguida, no Capítulo 4 são expostos trabalhos relacionados que propõem soluções paralelas para o problema da comparação sequência-família. São apresentadas soluções que utilizam extensões do conjunto de instruções da CPU, aceleradores em hardware implementados em FPGAs e soluções paralelas baseadas em clusters e GPUs, assim como seus resultados. No Capítulo 5, a solução em GPU desenvolvida para o problema da comparação sequência-família é apresentada. As plataformas, ferramentas e dados utilizados para desenvolvimento, execução e experimentos são descritos. A abordagem adotada para a solução é analisada, detalhando-se os kernels desenvolvidos. A organização das estruturas de dados da solução e seu mapeamento na hierarquia de memórias da GPU são apresentados e o consequente impacto no desempenho da solução é analisado. Também são apresentadas as otimizações aplicadas na solução. Os resultados obtidos na avaliação experimental são apresentados no Capítulo 6, sendo analisado o desempenho da solução em GPU e comparado com outras soluções. Por fim, o Capítulo 7 apresenta a conclusão deste trabalho, fazendo uma análise geral sobre o que foi desenvolvido. Outrossim, propostas de melhorias e trabalhos futuros são elicitadas. 4
19 2 Plataforma de Computação Paralela Impulsionadas pela demanda do mercado por gráficos de alta definição em tempo real para aplicações computacionais, as GPUs evoluíram, durante os últimos 30 anos, de dispositivos simples com um único pipeline para dispositivos altamente paralelos formados por vários longos pipelines, capazes de processar imagens interativas complexas de cenas 3D. Paralelamente a isso, muitas das funcionalidades do hardware tornaram-se mais sofisticadas e programáveis pelo programador [25]. Em pipelines gráficos, alguns estágios fazem uma grande quantidade de operações aritméticas com números em ponto flutuante utilizando dados completamente independentes, tais como transformar posições dos vértices de um triângulo ou gerar as cores de pixels. Essa independência de dados, característica dominante em aplicações gráficas, é uma diferença fundamental nas considerações de projetos para GPUs. As funções específicas executadas em alguns poucos estágios dos pipelines gráficos variam dependendo do algoritmo de renderização utilizado. Tal variação motivou o desenvolvimento de GPUs com estágios programáveis. Dois estágios destacaram-se, o vertex shader e o pixel shader. O primeiro mapeia as posições dos vértices de um triângulo na tela, alterando suas posições, cores ou orientação. Usualmente isso consiste em ler a posição de um vértice, que nada mais é que um conjunto de números em ponto flutuante, e calcular uma posição na tela que também é um conjunto de números em ponto flutuante. Já o segundo estágio define a cor de um pixel calculando as cores vermelho, verde, azul e alpha (RGBA), cada uma representada por um número em ponto flutuante, contribuindo assim para a renderização da imagem [25]. Uma vez que as operações ariméticas realizadas nesses dois estágios podiam ser controladas pelo programador, pesquisadores observaram que os dados usados como entrada poderiam na verdade ser quaisquer dados e não necessariamente cores ou posições, bastariam apenas estar no formato de cores ou posições. Em outras palavras, as GPUs poderiam ser enganadas para realizarem tarefas que não eram a de renderização de cenas através da formatação dos dados de entrada de maneira que elas pensassem que estariam renderizando uma cena [48]. Dado ao alto poder aritmético das GPUs, os resultados inciais desses experimentos foram bastante promissores. Entretanto, as limitações e dificuldades de se fazer isso pelo modelo de programação presente até então, as APIs gráficas, dificultavam o desenvolvimento. Esse desejo de usar as GPUs como um dispositivo de computação 5
20 paralela mais genérico motivou a NVIDIA [37], fabricante de GPUs, a desenvolver uma nova arquitetura gráfica unificada e também o modelo de programação CUDA (Compute Unified Device Architecture) [34]. Outro exemplo de modelo de programação paralela é o OpenCL (Open Computing Language) mantido pelo consórcio Khronos Group [23]. O OpenCL é um framework para desenvolvimento e execução de programas em plataformas heterogêneas, consistindo de CPUs, GPUs, FPGAs, entre outros. Assim, as GPUs programáveis evoluíram para um processador many-core, multithread e altamente paralelo. Além disso, esse paralelismo continua a avançar segundo a Lei de Moore [34]. Portanto, o desafio tornou-se desenvolver programas que utilizem esse paralelismo de forma transparente e que ganhem desempenho de forma escalável, para tirar proveito do aumento no número de cores das GPUs que, até o momento, continua a crescer a cada nova geração [11]. 2.1 Modelo de Programação CUDA O modelo CUDA é uma extensão das linguagens de programação C e C++ [25]. O programador escreve programas sequenciais que invocam kernels. Esses, por sua vez, podem ser simples funções ou programas completos. Um kernel é executado em paralelo através de um conjunto de threads paralelas. Essas threads são organizadas em uma hierarquia, na forma de grids de blocos de threads, como pode ser visto na Figura 2.1. Dessa forma, define-se: Thread: unidade básica de processamento; Bloco: um conjunto de threads concorrentes que podem cooperar entre si, sincronizar-se através de uma barreira de sincronização e comunicar-se através da memória compartilhada pertencente ao bloco; Grid: um conjunto de blocos de threads que podem ser executados concorrentemente e assim explorar paralelismo, ou seja, não há dependências entre os mesmos. A Figura 2.1 exemplifica o espaço de execução criado quando um determinado kernel é invocado. Um grid bidimensional é criado, com 3 2 blocos, cada bloco com 2 2 threads. Cada bloco e cada thread possuem uma identificação única. O modelo de programação CUDA permite a criação de grids unidimensionais ou bidimensionais e blocos unidimensionais, bidimensionais ou tridimensionais. O requisito de independência entre os blocos de threads ajuda a evitar possíveis deadlocks e permite que os mesmos possam ser escalonados em qualquer ordem, em 6
21 Grid Bloco (0,0) Bloco (1,0) Bloco (2,0) Bloco (0,1) Bloco (1,1) Bloco (2,1) Bloco (1,1) Thread (0,0) Thread (1,0) Thread (0,1) Thread (1,1) Figura 2.1: Organização das threads no modelo de programação CUDA: um grid com 2 3 blocos de 2 2 threads qualquer quantidade de cores, favorecendo a escalabilidade. O modelo de programação CUDA também não estabelece restrições sobre a localização de cada core, potencialmente permitindo a portabilidade de forma transparente de programas escritos no modelo através de múltiplos dispositivos no futuro [11]. Portanto, um programa CUDA compilado pode ser executado em uma GPU com uma quantidade qualquer de cores e apenas o sistema, em tempo de execução, precisa saber a quantidade de cores físicos [34] Programa Exemplo O problema da soma de dois vetores com elementos cada um, em paralelo, é usado para exemplificar uma solução desenvolvida usando o modelo de programação CUDA. A Figura 2.2 mostra o código CUDA que executa no computador host ao qual a CPU está acoplada (conforme será visto na Seção 2.2), enquanto a Figura 2.3 mostra o código do kernel executado na GPU. Na Figura 2.2, o código que executa no host inicialmente aloca espaço para os vetores na memória RAM do host (linhas 8 a 10) e na memória global da GPU (linhas 13 a 16) 1. Em seguida, os vetores são copiados da memória do host para a memória global 1 A hierarquia de memórias da GPU será descrita na Subseção
22 1 / Função executada no h o s t / 2 int main ( int argc, char argv [ ] ) 3 { 4 / Tamanho dos v e t o r e s / 5 int n = ; 6 7 / Aloca espa ço para v e t o r e s na memó r i a do h o s t / 8 double h_a = ( double ) malloc ( n sizeof ( double ) ) ; 9 double h_b = ( double ) malloc ( n sizeof ( double ) ) ; 10 double h_c = ( double ) malloc ( n sizeof ( double ) ) ; / Aloca espa ço para os v e t o r e s na memó r i a g l o b a l da GPU / 13 double d_a, d_b, d_c ; 14 cudamalloc(&d_a, n sizeof ( double ) ) ; 15 cudamalloc(&d_b, n sizeof ( double ) ) ; 16 cudamalloc(&d_c, n sizeof ( double ) ) ; / I n i c i a l i z a v e t o r e s no h o s t / / Copia v e t o r e s do h o s t para memó r i a g l o b a l da GPU / 22 cudamemcpy (d_a, h_a, n sizeof ( double ), cudamemcpyhosttodevice ) ; 23 cudamemcpy (d_b, h_b, n sizeof ( double ), cudamemcpyhosttodevice ) ; / Invoca k e r n e l com 100 b l o c o s de 1024 t h r e a d s / 26 somavetores <<<100, 1024>>>(d_a, d_b, d_c, n ) ; / Copia r e s u l t a d o s da memó r i a g l o b a l da GPU para o h o s t / 29 cudamemcpy (h_c, d_c, n sizeof ( double ), cudamemcpydevicetohost ) ; / Imprime r e s u l t a d o s / / Libera espa ço alocado para v e t o r e s na memó r i a g l o b a l da GPU / 35 cudafree (d_a ) ; 36 cudafree (d_b ) ; 37 cudafree (d_c ) ; / Libera espa ço alocado para v e t o r e s na memó r i a do h o s t / 40 f r e e (h_a ) ; 41 f r e e (h_b ) ; 42 f r e e (h_c ) ; return ( 0 ) ; 45 } Figura 2.2: Programa em CUDA que executa no host para soma de dois vetores da GPU (linhas 22 e 23). Essa transferência ocorre através do barramento PCI-Express, que conecta o host e a GPU. O host então invoca um kernel, na linha 26, para a execução da soma dos dois vetores em paralelo na GPU. A GPU realiza o processamento e escreve o resultado em sua memória global. O host copia o vetor com o resultado da soma da memória global da GPU para a memória do host (linha 29). Por fim, os espaços alocados nas memórias da 8
23 GPU e do host são liberados. Nesse exemplo o problema é modelado através de um grid unidimensional com 100 blocos de threads cada um. Em outras palavras, o problema é dividido em 100 partições independentes que podem ser executadas concorrentemente. As threads de um mesmo bloco calculam os elementos de cada partição, de forma concorrente às demais partições, visto que são operações independentes. O kernel da Figura 2.3 é executado por cada uma das threads do grid criado. Cada thread fica responsável por realizar a soma de um único par de elementos dos vetores de entrada. Cada thread determina a posição do vetor que ela deve acessar usando a identificação do bloco e da thread para compor o índice da posição no vetor (linha 6 no kernel). 1 / Kernel : cada thread c a l c u l a um elemento do v e t o r c / 2 global void somavetores ( double a, double b, double c, int n ) 3 { 4 / Calcula i d da thread e 5 í ndice do elemento do v e t o r que s e r á c a l c u l a d o p e l a thread / 6 int id = blockidx. x blockdim. x + threadidx. x ; 7 8 / V e r i f i c a se l i m i t e do v e t o r não f o i u l t r a p a s s a d o / 9 i f ( id < n ) 10 c [ id ] = a [ id ] + b [ id ] ; 11 } Figura 2.3: Kernel CUDA para soma de dois vetores 2.2 Arquitetura das GPUs CUDA Atualmente, as GPU são acopladas ao host por um barramento PCI-Express [44] que em sua versão mais recente, a versão 3.0, pode alcançar velocidades de transferência de até 16 GB/s. Apesar de suas altas velocidades, o barramento PCI-Express possui alta latência, o que pode ser um gargalo na utilização de GPUs como dispositivos de computação paralela. Dessa forma, é necessário utilizar técnicas que mascarem essa latência e tornem possível o aproveitamento máximo da execução de tarefas na GPU [25]. A arquitetura desses dispositivos massivamente paralelos baseia-se na replicação de blocos físicos denominados SMs (Streaming Multiprocessors) que por sua vez possuem vários cores denominados SPs (Streaming Processors), isto é, SMs são blocos físicos multithreaded [4, 11, 25]. A Figura 2.4 mostra, de forma simplificada, a arquitetura de uma GPU, baseada no modelo de programação CUDA, com n SMs, cada SM com um determinado número de SPs. Os SPs de um SM podem ser divididos em grupos, onde todos os SPs de um mesmo grupo executam a mesma instrução, contudo, cada SP pode 9
24 operar dados diferentes segundo o modelo SIMD (Single Instruction stream, Multiple Data streams) de execução [11]. Pela taxonomia de Flynn [17], SIMD é a classe de arquitetura de computadores na qual a mesma operação é realizada simultaneamente em múltiplos conjuntos de dados. CPU Host Memória do Sistema GPU CWD SM n Cache de Instruções Escalonador de Warps Escalonador de Warps Unidade de Despacho Unidade de Despacho SM 1 SM n Banco de Registradores... SP SP... SP SP SP... SP Rede de Interconexão Memória Global Rede de Interconexão Memória Compartilhada Figura 2.4: Arquitetura de uma GPU CUDA Warp é a unidade básica para o escalonamento de trabalho dentro de um SM e é, na maioria das GPUs, um grupo de 32 threads de um bloco alocado a um grupo de 32 SPs de um SM, operando de forma SIMD. Como as threads de um warp executam em lock-step, o escalonador não precisa verificar por dependências no fluxo de instruções. Na Figura 2.4 é possível identificar dois escalonadores de warps e duas unidades de despacho, que permitem que dois warps sejam despachados e executados concorrentemente. Visto que cada SM é responsável pelo escalonamento dos seus recursos internos, a cada ciclo de clock os escalonadores de warps do SM decidem quais warps executarão em seguida [11]. Quando threads de um mesmo warp tomam diferentes caminhos de execução dentro de um kernel devido a desvios condicionais, criam-se as chamadas divergências. Esses caminhos são executados sequencialmente até que as threads juntem-se novamente em um mesmo caminho de execução. Portanto, trechos de código com divergências podem causar perda de desempenho, dado a redução do paralelismo na execução dos mesmos [48]. Quando um programa no modelo CUDA executando no host invoca um kernel, a unidade CWD (Compute Work Distribution) realiza o escalonamento dos blocos, isto é, ela enumera os blocos de threads do grid e inicia a distribuição deles aos SMs com 10
25 capacidade de execução disponíveis. As threads de um bloco executam concorrentemente em um SM. Ao fim da execução de um bloco, a unidade CWD atribui um novo bloco ao SM vago. Portanto, a abstração em software de blocos de threads permite o mapeamento natural do kernel em uma GPU com um número arbitrário de SMs. Desse modo, blocos são um conjunto de threads que cooperam entre si, uma vez que apenas threads de um mesmo bloco podem compartilhar dados. Assim, essa abstração se torna uma expressão natural do paralelismo através do particionamento de um problema em múltiplas partições independentes que podem ser executadas concorrentemente. O escalonador pode dar a cada SM um ou mais blocos para serem executados, de acordo com a limitação do hardware da GPU [11]. Essa arquitetura é um fator chave na escalabilidade desses dispositivos, pois adicionando-se mais SMs em uma GPU pode-se executar mais tarefas ao mesmo tempo ou uma mesma tarefa mais rapidamente, se houver paralelismo suficiente para ser explorado nessa tarefa [4] Hierarquia de Memórias Uma GPU possui diferentes memórias, cada uma com características diferentes. Cada SM tem um barramento próprio de acesso à memória global da GPU e, consequentemente, à memória de textura e memória constante que são formas de endereçamento virtual dentro da memória global [4]. Cada bloco possui uma memória compartilhada, visível para todas as threads do bloco, e que tem o mesmo tempo de vida que o bloco. Por fim, as threads podem acessar dados de vários espaços de memória durante o seu tempo de execução. Cada thread possui uma memória local mapeada na memória global. Essa memória local é utilizada, no modelo CUDA, para as variáveis privadas da thread que não cabem nos registradores. Todas as threads de todos os blocos têm acesso à memória global. A Figura 2.5 esboça a hierarquia de memórias de uma GPU segundo o modelo de programação CUDA e quais unidades possuem acesso a essas memórias. Todas as threads de todos os blocos têm acesso às memórias global, de textura e constante da GPU. Todas as threads do bloco 0 têm acesso à memória compartilhada do bloco 0, assim como as threads do bloco 1 têm acesso à memória compartilhada do bloco 1. Cada thread possui também acesso exclusivo à sua memória local. Programas declaram variáveis na memória compartilhada e memória global do dispositivo através dos qualificadores shared e device, respectivamente. Esses espaços de memória correspondem a memórias fisicamente separadas: a memória 11
26 Bloco 0 Bloco 1 Memória Compartilhada Memória Compartilhada Registradores Registradores Registradores Registradores Thread 0 Thread 1 Thread 0 Thread 1 Memória Local Memória Local Memória Local Memória Local Memória Global Memória Constante Memória de Textura Figura 2.5: Hierarquia de memórias de uma GPU CUDA e unidades que possuem acesso às memórias compartilhada por bloco possui baixa latência, enquanto a memória global possui uma latência maior [34]. A função da memória global é ser uma memória com alta largura de banda, embora com uma maior latência que a memória do host [25]. De maneira simplificada, a memória compartilhada da GPU desempenha um papel semelhante à cache L1 de uma CPU (embora tenha que ser programada manualmente pelo desenvolvedor), dada sua baixa latência e por estar próxima aos SPs. Portanto, ela pode prover um acesso rápido e compartilhamento de dados entre as threads de um mesmo bloco. Visto que essa memória tem o mesmo tempo de vida do seu bloco correspondente, o código do kernel normalmente inicializa dados em estruturas alocadas na memória compartilhada, faz cálculos usando essas estruturas e, ao final, copia os resultados da memória compartilhada para a memória global. Por sua vez, blocos de threads pertencentes a grids dependentes comunicam-se através da memória global, usando-a para ler dados de entrada e escrever os resultados [34] Acesso Coalescido à Memória Global Para alcançar um melhor desempenho durante leituras e escritas na memória global, kernels CUDA devem realizar transações coalescidas nessa memória. Acessos não 12
27 coalescidos à memória global geram mais de uma transação de memória, consumindo portanto um maior tempo de execução que os acessos coalescidos. As seguintes restrições devem ser respeitadas para que as leituras ou escritas realizadas pelas threads de um warp sejam coalescidas (considerando GPUs com Compute Capability até 1.1): As palavras acessadas pelas threads do warp devem ter ao menos 32 bits. Ler ou escrever bytes ou palavras de 16 bits são sempre transações não coalescidas; Os endereços acessados pelas threads do warp devem ser contíguos e crescentes, isto é, o offset entre eles deve ser positivo e baseado na identificação da thread; O endereço base do warp, isto é, o endereço acessado pela primeira thread do warp deve ser alinhado de acordo com a Tabela 2.1. Em outras palavras, o endereço base acessado deve ser múltiplo do valor na coluna Alinhamento da Tabela 2.1. Por exemplo, em uma transação em que as threads do warp acessam palavras de 32 bits, o endereço base deve ser múltiplo de 64 bytes. Tabela 2.1: Critério de alinhamento para transações coalescidas na memória global Tamanho da palavra Alinhamento 8 bits Não coalescida 16 bits Não coalescida 32 bits 64 bytes 64 bits 128 bytes 128 bits 256 bytes A Figura 2.6(a) mostra as threads de um warp acessando palavras de 32 bits na memória global. Os acessos geram uma transação coalescida pois todas as restrições são respeitadas: as palavras têm 32 bits os endereços acessados pelas threads são crescentes e contíguos e, por fim, o endereço base do warp está alinhado segundo a Tabela 2.1. Já no caso ilustrado na Figura 2.6(b), os acessos das threads não são coalescidos. Apesar de todas as demais restrições serem atendidas, os endereços acessados pelas threads não são contíguos. Há um intervalo entre o endereço acessado pela thread 0 e a thread Sincronização em CUDA Um programa pode executar vários grids de forma independente ou dependente. Grids independentes podem executar concorrentemente se houver recursos de hardware suficientes. Por sua vez, grids dependentes executam sequencialmente, com uma barreira implícita inserida entre eles, garantindo assim que todos os blocos do primeiro grid 13
28 Warp Thread 0 Thread 1... Thread 31 Memória Global (a) Warp Thread 0 Thread 1... Thread 31 Memória Global (b) Figura 2.6: (a) Acessos coalescidos à memória global e (b) acessos não coalescidos à memória global 14
29 terminarão de executar antes de qualquer bloco do segundo grid, dependente, seja executado. A execução em paralelo e o gerenciamento das threads são feitos de forma automática. Toda criação de thread, escalonamento e finalização são realizados pelo hardware para o programador [4] Sincronização de Threads As threads de um bloco são executadas concorrentemente e podem sincronizar-se em uma barreira através do comando syncthreads(). Isto garante que nenhuma das threads participantes da barreira pode prosseguir até que todas atinjam a barreira, assegurando às mesmas a visualização de todas as escritas realizadas na memória. Dessa forma, threads em um bloco podem se comunicar através da escrita e leitura na memória compartilhada do bloco com uma barreira de sincronização para evitar condições de corrida [4] Streams No modelo CUDA, um stream é uma fila de operações que são executadas em uma ordem definida. Essas operações podem ser a invocação de kernels, transferências de dados entre host e GPU e o início ou parada de eventos. A ordem na qual as operações são adicionadas no stream é a ordem na qual elas serão executadas pela GPU. Assim, cada stream é um conjunto de operações que a GPU deve realizar em uma ordem definida. Uma operação de um stream só é iniciada após todas as operações anteriores daquele stream terem terminado. Entretanto, múltiplos streams podem ser executados em paralelo [48]. O uso de múltiplos streams pode tornar o código mais eficiente. Algumas GPUs mais novas permitem a execução em paralelo de kernels, de até duas transferências (uma do host para a GPU e outra da GPU para o host) e de kernels e transferências entre host e GPU. Sobrepor essas operações reduz o tempo total de execução do programa. A Figura 2.7 contrasta a execução de um conjunto de operações pela GPU sem e com o uso de múltiplos streams. No programa executado na Figura 2.7(a), um único stream é utilizado, logo não há sobreposição na execução das operações. Na Figura 2.7(b), o programa executado utiliza dois streams independentes. Enquanto o kernel de um stream executa na GPU, uma transferência host-gpu do outro stream é realizada. Com a sobreposição na execução de operações dos dois streams, esta solução resulta em um tempo de execução menor que aquele obtido com um único stream. 15
30 Stream 0 Stream 0 Stream 1 Transferência do vetor A do host para GPU Transferência do vetor A do host para GPU Transferência do vetor B do host para GPU Transferência do vetor B do host para GPU Execução do kernel na GPU Execução do kernel na GPU Transferência do vetor A do host para GPU Transferência do vetor C da GPU para host Transferência do vetor B do host para GPU Tempo Transferência do vetor A do host para GPU Transferência do vetor C da GPU para host Execução do kernel na GPU Transferência do vetor B do host para GPU Transferência do vetor C da GPU para host Execução do kernel na GPU Transferência do vetor C da GPU para host (a) (b) Figura 2.7: (a) Execução de programa com um único stream; (b) execução de programa com dois streams, sobrepondo operações e reduzindo o tempo total de execução 2.4 Escalabilidade do Modelo de Programação CUDA O desempenho e a escalabilidade de um programa desenvolvido usando o modelo de programação CUDA depende do particionamento do trabalho realizado pelo programa em blocos de tamanho fixo durante a modelagem do problema. Esses blocos fazem o mapeamento entre o paralelismo do problema em questão e o hardware da GPU [11]. Visto que as threads de um bloco podem compartilhar a memória compartilhada e sincronizar-se através de barreiras, elas residirão no mesmo SM. Contudo, o número de blocos pode exceder a capacidade dos SMs de uma GPU, o que dá ao programador flexibilidade para paralelizar na granularidade que for mais conveniente. Isso permite a decomposição de problemas de forma intuitiva, uma vez que o número de blocos pode ser ditado pelo tamanho dos dados a serem processados e não pela quantidade de SMs no dispositivo. Portanto, um mesmo programa pode ser executado em GPUs com quantidades diferentes de SMs, assim possibilitando o ganho de desempenho na execução desse programa, sem modificações, em gerações de GPUs mais novas e com mais recursos de hardware [34]. 16
31 Para gerenciar o escalonamento de blocos em SMs e fornecer escalabilidade, o modelo CUDA requer que os blocos executem de forma independente. Deve ser possível a execução dos blocos em qualquer ordem, em paralelo ou em série. Diferentes blocos não têm nenhum meio de comunicação direta, embora possam coordenar suas atividades com operações atômicas na memória global, visíveis a todas as threads [34]. Em suma, a abstração de blocos e a replicação de SMs no hardware trabalham em conjunto para fornecer escalabilidade de forma transparente [11]. 17
32 3 Comparação Sequência-Família Novas sequências biológicas evoluem a partir de sequências preexistentes, o que facilita a análise computacional das novas sequências. Muitas vezes é possível reconhecer uma semelhança significativa entre uma nova sequência e uma sequência biológica já conhecida e assim transferir informações sobre a estrutura e/ou função da sequência conhecida para a nova sequência. As duas sequências relacionadas são então ditas homólogas e a transferência de informações entre elas é feita por homologia [7]. À primeira vista, determinar se duas sequências biológicas são semelhantes não é muito diferente de determinar se duas cadeias de caracteres são semelhantes. Um conjunto de métodos para análise de sequências biológicas é, portanto, enraizado na Ciência da Computação, onde há uma extensa literatura sobre métodos de comparação de cadeias de caracteres [7]. 3.1 Sequência Biológica e Família de Sequências Uma sequência biológica pode ser definida como uma cadeia finita de símbolos que carregam informações biológicas. Estes símbolos pertencem a um alfabeto formado por um conjunto finito de símbolos [7]. Em geral, os alfabetos utilizados são [8]: Σ = {A, C, G, T }, para sequências de DNA, portanto composto por quatro símbolos que representam nucleotídeos; Σ = {A, C, G, U}, para sequências de RNA, portanto composto por quatro símbolos que representam nucleotídeos; Σ = {A, C, D, E, F, G, H, I, K, L, M, N, P, Q, R, S, T, V, W, Y }, para sequências de proteínas, portanto composto por vinte símbolos que representam aminoácidos. Assim, uma sequência biológica pode ser formalmente definida como uma palavra S, tal que, S Σ +. Por sua vez, uma família de sequências é um conjunto finito de sequências homólogas. As sequências membros de uma família de sequências divergiram durante a evolução e compartilham semelhanças estruturais e/ou funcionais [45]. 18
33 3.2 Alinhamento Múltiplo de Sequências e Profile HMM O alinhamento múltiplo de sequências é o alinhamento de três ou mais sequências biológicas e consiste na inserção de gaps (representados por hifens ) nas sequências com o objetivo de deixá-las com o mesmo tamanho e criar uma correspondência entre os símbolos das mesmas. Um alinhamento múltiplo pode ser representado por uma matriz, onde cada linha corresponde a uma sequência biológica e cada coluna tem os símbolos alinhados das sequências ou os gaps adicionados para alinhá-las [45]. A Figura 3.1 mostra o alinhamento múltiplo de cinco sequências de DNA. A A G A A A T A A C T G G C C A G C C G G Figura 3.1: Alinhamento múltiplo de cinco sequências de DNA Podem existir vários alinhamentos múltiplos distintos para um mesmo conjunto de sequências. De fato, várias funções e métodos para a obtenção do alinhamento múltiplo ótimo vêm sendo propostos na literatura. Por exemplo, para construir alinhamentos múltiplos de sequências biológicas é usual levar em consideração a relação evolucionária entre as sequências [45]. Dado o problema de determinar se uma nova sequência biológica é homóloga a uma família de sequências conhecida, é necessária uma forma de representar a família de sequências [7]. As variações em uma família de sequências podem ser descritas estatisticamente e essa é a base para a maioria dos métodos de análise de sequências biológicas atualmente [26]. Um profile é um modelo estatístico que mantém informações sobre uma família de sequências e descreve as variações entre os membros da família [31]. Krogh et al. [27] apresentam uma representação que usa Modelos Ocultos de Markov para descrever uma família de sequências e as relações entre seus membros. Um profile HMM (Hidden Markov Model Modelo Oculto de Markov) é um modelo estatístico baseado no alinhamento múltiplo das sequências de uma família e captura informações especificas de cada posição, representando o quão conservada cada coluna do alinhamento múltiplo está e quais símbolos são mais prováveis. Em outras palavras, um profile HMM estende o conceito de profile e modela uma família de sequências representando as similaridades entre seus membros na forma de estados discretos [31]. Uma 19