ads:
Cleomar Pereira da Silva
Computação de Alto Desempenho com Placas
Gráficas para Acelerar o Processamento da Teoria
do Funcional da Densidade
Dissertação de Mestrado
Dissertação apresentada como requisito parcial para
obtenção do grau de Mestre pelo Programa de Pós-
graduação em Engenharia Elétrica do Departamento
de Engenharia Elétrica da PUC-Rio.
Orientador: Prof. Marco Aurélio Cavalcanti Pacheco
Rio de Janeiro
Abril de 2010
PUC-Rio - Certificação Digital Nº 0812691/CA
ads:
Livros Grátis
http://www.livrosgratis.com.br
Milhares de livros grátis para download.
Cleomar Pereira da Silva
Computação de Alto Desempenho com Placas
Gráficas para Acelerar o Processamento da Teoria
do Funcional da Densidade
Dissertação apresentada como requisito parcial para obtenção
do grau de Mestre pelo Programa de Pós-Graduação em
Engenharia Elétrica do Departamento de Engenharia Elétrica do
Centro Técnico Científico da PUC-Rio. Aprovada pela Comissão
Examinadora abaixo assinada.
Prof. Marco Aurélio Cavalcanti Pacheco
Orientador
Departamento de Engenharia Elétrica – PUC-Rio
Prof. Omar Paranaina Vilela Neto
Departamento de Engenharia Elétrica – PUC-Rio
Profa. Cristiana Bentes
UERJ
Prof. Ricardo Cordeiro de Farias
UFRJ
Prof. José Eugenio Leal
Coordenador Setorial do Centro
Técnico Científico – PUC-Rio
Rio de Janeiro, 13 de abril de 2010
PUC-Rio - Certificação Digital Nº 0812691/CA
ads:
Todos os direitos reservados. É proibida a reprodução total ou
parcial do trabalho sem autorização da universidade, do autor e do
orientador.
Cleomar Pereira da Silva
Graduado Engenheiro Eletricista pela Universidade Federal de
Santa Maria (UFSM) em 2008. Mestrado em Engenharia Elétrica,
área de concentração Nanotecnologia, pela Pontifícia Universidade
Católica do Rio de Janeiro (PUC-Rio) em 2010.
Ficha Catalográfica
CDD: 621.3
Silva, Cleomar Pereira da
Computação de alto desempenho com placas gráficas
para acelerar o processamento da teoria do funcional da
densidade / Cleomar Pereira da Silva ; orientador: Marco
Aurélio Cavalcanti Pacheco. – 2010.
87 f. ; 30 cm
Dissertação (mestrado)–Pontifícia Universidade
Católica do Rio de Janeiro, Departamento de Engenharia
Elétrica, 2010.
Inclui bibliografia
1. Engenharia elétrica – Teses. 2. Teoria do funcional
de densidade. 3. SIESTA. 4. Computação de alto
desempenho. 5. GPGPU. 6. CUDA. I. Pacheco, Marco Aurélio
Cavalcanti. II. Pontifícia Universidade Católica do Rio de
Janeiro. Departamento de Engenharia Elétrica. III. Título.
PUC-Rio - Certificação Digital Nº 0812691/CA
Aos Meus Pais Antonio Coraci e Maria Rosa.
PUC-Rio - Certificação Digital Nº 0812691/CA
Agradecimentos
À CAPES e à PUC-Rio pela concessão da bolsa e dos auxílios que tornaram
possível a realização deste trabalho.
Ao meu orientador, Dr. Marco Aurélio Cavalcanti Pacheco, pelo imediato
atendimento das solicitações realizadas durante o período do estudo.
À professora Dra. Marley M. B. R. Vellasco pelos conhecimentos transferidos
através das disciplinas que ministrou.
Aos demais professores pelos ensinamentos.
Ao Dr. Omar Paranaiba V. Neto pela participação na elaboração das idéias
envolvidas com o tema desta dissertação.
Aos amigos Douglas M. Dias, Leandro F. Cupertino, Iury S. O. Bezerra, Daniel S.
Chevitarese, Dr. Juan L. Lazo e Dr. Renato B. Oliveira, pelo apoio e auxilio
prestados no desenvolvimento deste trabalho.
Aos meus pais, pela educação, atenção e carinho.
Aos amigos do ICA por seu contínuo apoio e colaboração.
PUC-Rio - Certificação Digital Nº 0812691/CA
Resumo
Silva, Cleomar Pereira; Pacheco, Marco Aurélio Cavalcanti (Orientador).
Computação de Alto Desempenho com Placas Gráficas para Acelerar o
Processamento da Teoria do Funcional da Densidade. Rio de Janeiro,
2010. 87p. Dissertação de Mestrado - Departamento de Engenharia Elétrica,
Pontifícia Universidade Católica do Rio de Janeiro.
As Unidades de Processamento Gráfico (GPUs), ou Placas Gráficas, são
processadores que foram originalmente projetados para executar tarefas dedicadas
às operações da computação gráfica. Porém, a NVIDIA desenvolveu uma
extensão da linguagem C para programação de GPUs, chamada CUDA (Compute
Unified Device Architecture). Isto permitiu utilizá-las, na Computação de Alto
Desempenho, para processar dados genéricos. Já os sistemas físicos estudados
pela Mecânica Quântica apresentam dimensões próximas da escala atômica, tais
como moléculas, átomos, prótons e elétrons. A Teoria do Funcional da Densidade
(DFT) é um dos métodos iterativos mais usados para encontrar uma solução
aproximada para a equação de Schrödinger. Contudo, os cálculos realizados em
DFT são computacionalmente intensos devido às integrais de troca e correlação
eletrônica, integrais para o cálculo da energia de Hartree e energia cinética dos
elétrons, as quais requerem maior esforço computacional à medida que o número
de elétrons presentes na simulação aumenta. Esta pesquisa teve como objetivo
estudar os cálculos do DFT e identificar partes do algoritmo que, se alteradas,
apresentassem benefícios de desempenho ao serem executadas em GPU. Assim,
funções computacionalmente intensas do método DFT do SIESTA (Spanish
Initiative for Electronic Simulations with Thousands of Atoms) foram
paralelizadas e usadas para calcular propriedades físicas de nanotubos e fulerenos.
Verificou-se que a execução da versão paralela do SIESTA para GPU é capaz de
atingir ganhos em desempenho, em funções individuais, de uma ou até duas
ordens de grandeza, tornando promissor o emprego de GPUs em acelerar o
processamento da Teoria do Funcional da Densidade.
Palavras-chave
Teoria do Funcional da Densidade; SIESTA; Computação de Alto
Desempenho; GPGPU; CUDA.
PUC-Rio - Certificação Digital Nº 0812691/CA
Abstract
Silva, Cleomar Pereira; Pacheco, Marco Aurélio Cavalcanti (Advisor). High
Performance Computing with Graphics Cards to Accelerate Processing
Density Functional Theory. Rio de Janeiro, 2010. 87p. MSc. Dissertation -
Departamento de Engenharia Elétrica, Pontifícia Universidade Católica do
Rio de Janeiro.
The graphics processing units (GPUs), or graphics cards, are processors that
were originally designed to perform dedicated tasks to the computer graphics
operations. However, NVIDIA has developed an extension of the C language for
programming GPUs, called CUDA (Compute Unified Device Architecture). This
allowed the use of GPUs, in High Performance Computing, for processing generic
data. The physical systems studied by quantum mechanics have dimensions close
to atomic scale, such as molecules, atoms, protons and electrons. The Density
Functional Theory (DFT) is one of the most used interactive methods to find an
approximated solution to the Schrödinger equation. However, the calculations in
DFT are computationally intensive because of the exchange and correlation
electronic integrals, integrals to calculate the Hartree energy and electrons kinetic
energy, which requires greater computational effort as the number of electrons
present in the simulation increases. This research aimed to study the DFT
calculations and identify parts of the algorithm that, if changed, experience
performance benefits from execution in GPU. Thus, computationally intensive
DFT functions of the SIESTA method (Spanish Initiative for Electronic
Simulations with Thousands of Atoms) were parallelized and used to calculate the
physical properties of nanotubes and fullerenes. It was found that the
implementation of SIESTA parallel version on the GPU is able to achieve gains in
performance, in individual functions, of one or even two orders of magnitude,
making it promising employment of GPUs to speed up the processing of Density
Functional Theory.
Keywords
Density Functional Theory; SIESTA; High Performance Computing;
GPGPU; CUDA.
PUC-Rio - Certificação Digital Nº 0812691/CA
Sumário
1 Introdução 14
1.1. Motivação 16
1.2. Objetivos 16
1.3. Descrição do Trabalho 16
1.4. Estrutura da Dissertação 18
2 Teoria do Funcional da Densidade 19
2.1. Componentes da Energia Total 20
2.2. Aproximação de Born-Oppenheimer 21
2.3. Teoremas de Hohenberg e Kohn 23
2.4. Equações de Kohn-Sham 24
2.5. Aproximações dos Potenciais de Troca e Correlação 25
2.5.1. Aproximação da Densidade Local 26
2.5.2. Aproximação do Gradiente Generalizado 26
2.6. Pseudopotencial 27
2.7. Hamiltoniano Empregado no SIESTA 29
2.7.1. Funções de Base 30
2.7.2. Orbitais Atômicos Numéricos 31
3 Computação de Propósito Geral em Unidades de Processamento
Gráfico 32
3.1. Programação Paralela 32
3.1.1. Taxonomia de Flynn 33
3.1.2. Lei de Amdahl 35
3.2. Compute Unified Device Architecture (CUDA) 36
3.3. Elementos de Hardware da GPU 39
3.3.1. Processador e Multiprocessador de Fluxo 39
3.3.2. Tipos de Memória da GPU 40
3.4. Elementos de Software da GPU 41
3.4.1. Kernel 42
3.4.2. Grid 43
PUC-Rio - Certificação Digital Nº 0812691/CA
3.4.3. Bloco de Threads 44
3.4.4. Warp 44
3.4.5. Escalabilidade 45
3.5. Técnicas Otimizadas 46
3.5.1. Acesso à Memória Global 46
3.5.2. Acesso à Memória Compartilhada 51
3.5.3. Transferência de Dados Através do Barramento PCI Express 53
3.5.4. Nível de Ocupação do Multiprocessador 54
4 Descrição das Alterações para Aceleração do SIESTA por GPU 56
4.1. Identificação das Partes Adequadas ao Paralelismo de Dados 57
4.1.1. CUFFT 57
4.1.2. CUBLAS 58
4.1.3. MAGMA 58
4.1.4. CULA 59
4.2. Alteração de Partes do SIESTA 59
4.2.1. Potencial de Hartree, Energia de Hartree e Forças de
Deslocamento Atômicas 60
4.2.2. Implementação Paralela para GPUs 61
4.2.3. Cálculo do Dipolo Elétrico 68
4.2.4. Reordenação de Dados 71
5 Estudo de Casos 73
5.1. Caso 1 – Cálculo da Energia e Estrutura de Bandas de Nanotubos 73
5.1.1. Estrutura do Nanotubo 74
5.1.2. Resultados e Testes de Desempenho 75
5.2. Caso 2 – Otimização da Geometria Estrutural de Fulereno 78
5.2.1. Estrutura do Fulereno 78
5.2.2. Resultados e Testes de Desempenho 79
6 Conclusões e Trabalhos Futuros 81
7 Referências Bibliográficas 84
PUC-Rio - Certificação Digital Nº 0812691/CA
Lista de Figuras
Figura 1 – Função de onda de todos os elétrons normalizada (AE)
e a pseudofunção de onda de valência normalizada (PS). 28
Figura 2 – Arquitetura SISD. 34
Figura 3 – Arquitetura SIMD. 34
Figura 4 – Arquitetura MISD. 35
Figura 5 – Arquitetura MIMD. 35
Figura 6 – Efeito da Lei de Amdahl sobre o ganho de velocidade de
algoritmo rodando em múltiplos processadores. 36
Figura 7 – Operações de Ponto Flutuante por Segundo (NVIDIA, 2009b). 37
Figura 8 – Largura de Banda (NVIDIA, 2009b). 38
Figura 9 – A GPU dedica mais transistores para o processamento de
dados (NVIDIA, 2009b). 38
Figura 10 – Tipos de memória das GPUs da NVIDIA (NVIDIA, 2009b). 41
Figura 11 – Fluxo de processamento de um software empregando GPGPU. 42
Figura 12 – Grid e Blocos de Threads (NVIDIA, 2009b). 43
Figura 13 – Modelo de Escalabilidade da GPU (NVIDIA, 2009b). 46
Figura 14 – Segmentos Alinhados de Memória e threads de um half-warp. 47
Figura 15 – As threads do half-warp acessam elementos da memória
global dentro de um mesmo segmento alinhado de 16 palavras (NVIDIA,
2009a). 47
Figura 16 – Elementos não alinhados com o bloco de 16 palavras residem
dentro do mesmo bloco de 128 bytes (NVIDIA, 2009a). 48
Figura 17 – Elementos não alinhados com o bloco de 16 palavras residem
dentro de dois blocos diferentes de 128 bytes (NVIDIA, 2009a). 48
Figura 18 – Largura de banda efetiva para padrões de acesso à memória
com deslocamento (NVIDIA, 2009a). 49
Figura 19 – Acesso a elementos da memória com intervalos entre os
elementos acessados (NVIDIA, 2009a). 49
Figura 20 – Desempenho da largura de banda para padrões de acesso
com intervalos entre os elementos lidos (NVIDIA, 2009a). 50
PUC-Rio - Certificação Digital Nº 0812691/CA
Figura 21 – Leitura, com deslocamentos, na memória global e na memória
de textura (NVIDIA, 2009a). 50
Figura 22 – Padrões de acesso à memória compartilhada sem conflitos de
banco (NVIDIA, 2009b). 52
Figura 23 – Padrões de acesso à memória compartilhada com conflitos de
banco (NVIDIA, 2009b). 53
Figura 24 – Redução do tempo total através da sobreposição dos tempos
de transferência e de execução (NVIDIA, 2009a). 54
Figura 25 – Interface entre Fortran e CUDA. 59
Figura 26 – As equações diferenciais são transformadas em equações
algébricas. 60
Figura 27 – Somatório paralelo realizado na memória compartilhada da GPU
para um único vetor de dados. 61
Figura 28 – Solução no domínio da freqüência. 62
Figura 29 – Uma thread para cada elemento da matriz de densidades. 63
Figura 30 – Nível de Ocupação do Multiprocessador para primeira versão
do kernel da Equação de Poisson. 64
Figura 31 – Uma thread processa um vetor de dados da matriz de densidades 64
Figura 32 – Nível de Ocupação do Multiprocessador para segunda versão
do kernel da Equação de Poisson. 65
Figura 33 – Nível de Ocupação do Multiprocessador para o kernel dos
somatórios da energia de Hartree e das forças de deslocamento atômicas. 66
Figura 34 – Somatório paralelo de sete vetores na memória compartilhada
da GPU. 67
Figura 35 – Padrões usados de acesso à memória compartilhada, sem conflitos
de banco. 68
Figura 36 – Nível de Ocupação do Multiprocessador para o kernel de Cálculo
do Dipolo Elétrico. 69
Figura 37 – Somatório paralelo de três vetores na memória compartilhada
da GPU. 70
Figura 38 – Nível de Ocupação do Multiprocessador para o kernel
da Reordenação de Dados. 72
PUC-Rio - Certificação Digital Nº 0812691/CA
Figura 39 – Estrutura de bandas de nanotubos (a) armchair (5, 5), (b) zigzag
(8, 0) e (c) zigzag (12, 0). O nível de Fermi está deslocado para o zero,
indicado pela linha tracejada (Silva, 2008). 74
Figura 40 – Nanotubo (4,2). a) Vista Frontal, b) Vista Lateral. 75
Figura 41 – Estrutura de bandas para o nanotubo (4,2). O nível de Fermi
está deslocado para o zero, indicado pela linha tracejada. 78
Figura 42 – Fulereno C60. 79
PUC-Rio - Certificação Digital Nº 0812691/CA
Lista de Tabelas
Tabela 1 – Unidades Atômicas 21
Tabela 2 – Funções da Biblioteca CUFFT adequadas para emprego no
SIESTA 58
Tabela 3 – Funções da Biblioteca CUBLAS com possibilidade de emprego
para o SIESTA. 58
Tabela 4 – Funções da Biblioteca MAGMA com possibilidade de emprego
para o SIESTA 58
Tabela 5 – Funções da Biblioteca CULA utilizadas no meio científico 59
Tabela 6 – Operações de adição com warp incompleto por bloco, replicando
sete vezes o padrão ilustrado na Figura 27. 66
Tabela 7 – Operações de adição com warp incompleto por bloco, com o
padrão ilustrado na Figura 34. 67
Tabela 8 – Operações de adição com warp incompleto por bloco, replicando
três vezes o padrão ilustrado na Figura 27. 69
Tabela 9 – Operações de adição com warp incompleto por bloco, com o
padrão ilustrado na Figura 37. 70
Tabela 10 – Resultados de Aceleração com a primeira versão do kernel
da Equação de Poisson, para o nanotubo (4,2). 76
Tabela 11 – Resultados de Aceleração com a segunda versão do kernel
da Equação de Poisson, para o nanotubo (4,2). 76
Tabela 12 – Resultados de Aceleração Reduzindo-se a Transferência de
Dados através do Barramento PCI Express, para o nanotubo (4,2). 77
Tabela 13 – Resultados de Aceleração com a segunda versão do kernel
da Equação de Poisson, para o fulereno C60. 80
Tabela 14 – Resultados de Aceleração Reduzindo-se a Transferência de
Dados através do Barramento PCI Express, para o fulereno C60. 80
PUC-Rio - Certificação Digital Nº 0812691/CA
1
Introdução
As Unidades de Processamento Gráfico (GPUs), ou Placas Gráficas, foram
originalmente desenvolvidas com o propósito de renderização gráfica. Contudo,
nos últimos anos, o desempenho na realização de operações de ponto flutuante das
GPUs superou em muito o desempenho das CPUs (Unidades Centrais de
Processamento). As GPUs possuem um elevado número de elementos de
processamento e uma grande largura de banda de memória. Constituem assim,
uma poderosa ferramenta, tanto para o processamento gráfico, quanto para a
execução de trechos de código de programas de propósito geral.
Recentemente, muitas aplicações científicas foram portadas para a
plataforma de hardware das GPUs. Assim, podemos citar o estudo de Anderson,
Goddard & Schröder (2007) para execução do Método de Monte Carlo Quântico
nesta plataforma. As simulações envolvendo Elementos Finitos também se
apresentam adequadas (Göddeke, Strzodka e Turek, 2007). Stone, Phillips, et al.
(2007) relataram a aceleração de cálculos de mecânica molecular e Meel, Arnold,
et al. (2008) descrevem o emprego de GPUs para simulações de dinâmica
molecular clássica. Outra área que se beneficia da alta capacidade de cálculo das
GPUs é a Simulação de Nanodispositivos, tais como dispositivos nano-fotônicos
(Kelmelis, Durbano, et al., 2006; Price, Humphrey e Kelmelis, 2007).
No campo dos cálculos quânticos da estrutura eletrônica, Yasuda (2008b)
propôs um algoritmo para o cálculo da integral da repulsão de dois elétrons em
GPUs. Mais recentemente, Ufimtsev & Martínez (2008b) revisaram este
algoritmo e propuseram novas estratégias, obtendo resultados promissores. O
campo autoconsistente do algoritmo de Hartree-Fock, utilizado no GAMESS
(General Atomic and Molecular Electronic Structure System), foi acelerado por
GPUs (Ufimtsev e Martínez, 2008a). Métodos de maior acurácia na energia total
calculada, como os que empregam a teoria da perturbação de segunda ordem,
foram parcialmente modificados para execução em GPUs por Vogt, Olivares-
Amaya, et al. (2008).
PUC-Rio - Certificação Digital Nº 0812691/CA
Introdução 15
A Teoria do Funcional da Densidade (DFT) está entre os mais populares e
versáteis métodos disponíveis para estudos de química computacional e física do
estado sólido. É usada para estudar a estrutura eletrônica (principalmente o estado
fundamental) de sistemas de muitos corpos, em particular átomos, moléculas e
sólidos. O seu formalismo foi estabelecido a partir dos dois teoremas de
Hohenberg & Kohn (1964). Eles demonstraram que em princípio a densidade
eletrônica contém toda a informação que pode ser obtida da função de onda de
muitos elétrons. Conforme será detalhado no Capítulo 2, o cálculo da energia de
um sistema físico por DFT geralmente encontra-se subdividido em cinco termos
principais.
Em trabalhos relacionados à DFT, Yasuda (2008a) obteve bons resultados
com o termo de troca e correlação da energia total em GPU. Os termos da energia
relativos ao pseudopotencial do algoritmo BigDFT, o qual usa conjuntos de base
com funções wavelets, foi estudado por Genovese, Ospici, et al. (2009), obtendo-
se acelerações de até 20 vezes em algumas operações. Genovese, Ospici, et al.
(2009) calcularam a solução da Equação de Poisson na CPU e comentaram que a
solução da mesma na GPU teria eliminado a necessidade de duas transferências de
memória para a variável densidade. A transferência de memória entre CPU e
GPU, através do barramento PCI Express, constitui um dos potenciais gargalos
que deve ser evitado para obter bom desempenho.
O SIESTA (Spanish Initiative for Electronic Simulations with Thousands of
Atoms) é um método autoconsistente da Teoria do Funcional da Densidade que
usa o pseudopotencial para representação do caroço iônico, orbitais atômicos
numéricos e funções de base localizada (Soler, Artacho, et al., 2002; Ordejón,
Artacho e Soler, 1996). Este algoritmo representa um esforço para o
desenvolvimento de um método ab initio DFT autoconsistente de ordem N, O(N),
para CPU. A determinação do hamiltoniano autoconsistente em O(N) operações é
difícil de ser obtida usando-se ondas planas como funções de base. Isto levou à
escolha de um conjunto de bases atômicas localizadas para o SIESTA. Contudo,
algumas funções ainda necessitam um número superior de operações, tais como a
Transformada de Fourier empregada na resolução da Equação de Poisson, a qual é
de ordem N log(N).
PUC-Rio - Certificação Digital Nº 0812691/CA
Introdução 16
1.1.
Motivação
Este trabalho envolve o estudo de duas áreas distintas do conhecimento, a
Teoria do Funcional da Densidade (DFT) e a Programação de Propósito Geral em
Unidades de Processamento Gráfico (GPGPU). O poder computacional crescente
das GPUs pode ser empregado para auxiliar na realização dos cálculos envolvidos
na DFT. Assim, o estudo destas duas áreas distintas visa prestar uma contribuição
no sentido de apontar as alterações da DFT para execução acelerada por GPU. A
modificação de partes do código e a realização de simulações de propriedades
físicas empregando o SIESTA são utilizadas para comprovar os ganhos de
desempenho que podem ser alcançados pela colaboração entre estas duas áreas.
1.2.
Objetivos
O objetivo geral desta pesquisa é investigar o apoio das GPUs na aceleração
dos cálculos envolvidos na Teoria do Funcional da Densidade.
Entre os objetivos específicos encontra-se a identificação dos trechos do
código SIESTA que seriam adequados ao paralelismo de dados oferecido pelas
placas gráficas. A alteração de algumas destas partes identificadas como sendo
adequadas para GPUs. Além da realização de cálculos de propriedades físicas em
estruturas, tais como nanotubos de carbono e fulerenos, com a finalidade de
comprovar o ganho de velocidade da execução em GPUs para as partes alteradas.
1.3.
Descrição do Trabalho
Descrição das principais etapas seguidas no desenvolvimento deste trabalho:
1. Estudo da Teria do Funcional da Densidade:
Para isso, uma revisão bibliográfica foi realizada visando o entendimento
das equações fundamentais da Teria do Funcional da Densidade
2. Estudo do Modelo de Programação de Propósito Geral em
Unidades de Processamento Gráfico:
Revisão do conhecimento bibliográfico disponível na literatura evolvendo a
programação CUDA (Compute Unified Device Architecture) para GPUs.
PUC-Rio - Certificação Digital Nº 0812691/CA
Introdução 17
3. Seleção de um Algoritmo DFT:
Foi realizada uma busca por um algoritmo eficiente, de código fonte aberto,
para servir de ponto de partida para estudar a alteração de trechos de código para
GPU. Diversos pacotes de programas, usados em cálculos de química quântica,
foram compilados, optando-se finalmente por estudar mais detalhadamente o
código fonte do SIESTA.
4. Identificação e Alteração de Partes do SIESTA:
O estudo das aproximações DFT e da forma particular usada no SIESTA
para representar os termos da energia total possibilitou a identificação de trechos
de código com potencial de ganhos para execução em GPU. Algumas destas
partes identificadas foram reescritas usando a linguagem CUDA, tais como o
cálculo do potencial de Hartree e da energia de Hartree, usando a resolução da
Equação de Poisson na GPU. A contribuição para as forças de deslocamento
atômicas, devidas ao potencial de Hartree, também são computadas em GPU.
5. Comprovação do Desempenho:
Verificação da correção dos resultados obtidos com a versão das funções
executadas em GPU, quando comparada às funções seqüenciais originais.
Comparação dos tempos das simulações para avaliar a aceleração obtida para
estas funções alteradas para executar em GPUs.
6. Resultados:
Propriedades físicas de nanotubos e de fulerenos foram calculadas na CPU e
na GPU para comparação do desempenho. A solução da Equação de Poisson na
GPU para calcular o potencial de Hartree e a energia de Hartree pode ser
acelerada cerca de 30 vezes. A reordenação de dados em GPU apresenta
importância para reduzir o volume de transferências de memória através do
barramento PCI Express. O cálculo do dipolo elétrico pode ser acelerado em até
duas ordens de grandeza.
7. Sugestões de Trabalhos Futuros:
A continuação do trabalho, com a alteração do código fonte para execução
em GPUs, dos demais termos da energia total do programa SIESTA. Sugestão do
emprego de técnicas inteligentes para redução do esforço necessário para
obtenção de código fonte compilável para execução em GPU.
PUC-Rio - Certificação Digital Nº 0812691/CA
Introdução 18
1.4.
Estrutura da Dissertação
No Capítulo introdutório são listadas aplicações gerais de GPUs no meio
científico, e aplicações voltadas, especificamente, para os métodos empregados na
química quântica.
No Capítulo 2 é feita uma revisão bibliográfica sobre as principais
aproximações empregadas na Teoria do Funcional da Densidade, visando tornar
possível a obtenção de uma solução aproximada para a equação de Schrödinger
em tempo hábil. Ao final deste Capítulo é feita uma ligação da teoria DFT com a
forma particular empregada no SIESTA para fazer a representação dos termos da
energia total no hamiltoniano.
O Capítulo 3 é dedicado a fazer uma revisão do conhecimento disponível na
literatura sobre a Computação de Propósito Geral em Unidades de Processamento
Gráfico. Este conhecimento será empregado para realizar as alterações de partes
do código SIESTA para execução acelerada por GPU.
As partes do código SIESTA selecionadas para alteração são apresentadas
no Capítulo 4. São explicadas as modificações realizadas fazendo-se uma ligação
com os termos da energia do SIESTA explicados no final do Capítulo 3.
O Capítulo 5 descreve os cálculos de propriedades físicas, tais como a
energia e a estrutura de bandas de nanotubos e a otimização da geometria
estrutural de fulerenos. Estas simulações foram realizadas com a finalidade de
comprovar a eficiência das GPUs em acelerar trechos de códigos empregados na
DFT.
Por fim, as conclusões e os trabalhos futuros, identificados durante a
realização desta pesquisa, são apresentados no Capítulo 6.
PUC-Rio - Certificação Digital Nº 0812691/CA
2
Teoria do Funcional da Densidade
A Mecânica Quântica teve seu início com a equação de Schrödinger (1926).
Esta equação determina a função de onda quântica de um sistema, seja ele um
átomo, uma molécula ou um sólido, contendo toda a informação necessária para
determinar o estado do sistema. Contudo, são poucos os sistemas físicos que
possuem solução analítica e mesmo a solução numérica aproximada pode ser
computacionalmente inviável.
Uma aproximação baseada somente na densidade eletrônica foi proposta por
Thomas (1927) e Fermi (1927). Os dois pesquisadores, trabalhando de forma
independente, empregaram um modelo estatístico para aproximar a distribuição
dos elétrons nos átomos, o qual ficou conhecido como modelo de Thomas-Fermi.
Apesar da baixa qualidade das previsões para sistemas reais, este modelo é o
precursor da moderna Teoria do Funcional da Densidade.
O formalismo da Teoria do Funcional da Densidade – DFT (Density
Functional Theory) foi estabelecido a partir dos dois teoremas de Hohenberg &
Kohn (1964). Eles demonstraram que em princípio a densidade eletrônica contém
toda a informação que pode ser obtida da função de onda de muitos elétrons.
O conhecimento da densidade do estado fundamental de um sistema de
muitos elétrons permite deduzir, a menos de uma constante, o potencial externo
no qual os elétrons residem (Hohenberg e Kohn, 1964). Utilizando este teorema e
sabendo-se que todos os outros termos do hamiltoniano que representa um sistema
de muitos elétrons podem ser escritos, em princípio, como um funcional único da
densidade eletrônica, conclui-se que o conhecimento da densidade eletrônica do
estado fundamental determina completamente o problema de um sistema de
muitos corpos. Enquanto a função de onda necessita de 3N
e
variáveis (função de
onda para cada elétron) para a sua descrição, a densidade é uma função real de 3
variáveis (densidade em três dimensões).
A DFT é usada para estudar a estrutura eletrônica (principalmente o estado
fundamental) de sistemas de muitos corpos, em particular átomos, moléculas e
PUC-Rio - Certificação Digital Nº 0812691/CA
Teoria do Funcional da Densidade 20
sólidos. Com esta teoria, as propriedades de um sistema de muitos elétrons podem
ser determinadas usando-se funcionais, isto é, funções de outra função, a qual
neste caso é a densidade eletrônica. A DFT está entre os mais populares e
versáteis métodos disponíveis para estudos de química computacional e física do
estado sólido.
O programa SIESTA (Spanish Initiative for Electronic Simulations with
Thousands of Atoms) faz uso da Teoria do Funcional da Densidade (Soler,
Artacho, et al., 2002; Ordejón, Artacho e Soler, 1996). Ele foi utilizado para a
realização deste trabalho e neste Capítulo serão detalhadas as aproximações
empregadas.
2.1.
Componentes da Energia Total
A resolução da equação de Schrödinger independente do tempo permite
determinar a estrutura do estado fundamental de um sistema de muitos elétrons e
núcleos:
ܪ
ȰൌܧȰ
(2.1)
O Hamiltoniano ܪ
do sistema apresenta o seguinte formato:
ܪ
ൌܶ
ܶ
ܸ
ܸ
ܸ
(2.2)
A energia cinética ܶ
dos M núcleos do sistema é representada pelo termo:
ܶ
ൌെ
ͳ
ʹ
ͳ
ܯ
ఈ
ఈ
ଶ
ெ
ఈୀଵ
(2.3)
Onde
ఈ
ଶ
ൌ
߲
ଶ
߲ݔ
ఈ
ଶ
߲
ଶ
߲ݕ
ఈ
ଶ
߲
ଶ
߲ݖ
ఈ
ଶ
(2.4)
A energia cinética ܶ
dos N elétrons do sistema é dada pelo termo:
ܶ
ൌെ
ͳ
ʹ
ଶ
ே
ୀଵ
(2.5)
O terceiro termo do lado direito da eq. (2.2) é a energia de atração
elétron-núcleo ܸ
dada pela eq. (2.6):
ܸ
ൌെ
ܼ
ఈ
ȁ
ܴ
ఈ
െݎ
ȁ
ே
ୀଵ
ெ
ఈୀଵ
(2.6)
PUC-Rio - Certificação Digital Nº 0812691/CA
Teoria do Funcional da Densidade 21
Onde ܴ
ఈ
e ݎ
representam, respectivamente, os raios nucleares e os raios
eletrônicos.
A energia de repulsão elétron-elétron ܸ
, quarto termo do lado direito da
eq. (2.2), é dada por:
ܸ
ൌ
ͳ
หݎ
െݎ
ห
ே
வ
ே
ୀଵ
(2.7)
O quinto termo da eq. (2.2) representa a energia de repulsão núcleo-núcleo
ܸ
, sendo dado por:
ܸ
ൌ
ܼ
ఈ
ܼ
ఉ
หܴ
ఈ
െܴ
ఉ
ห
ெ
ఉவఈ
ெ
ఈୀଵ
(2.8)
As unidades atômicas foram utilizadas para simplificar as fórmulas. A
Tabela 1, adaptada de Koch & Holthausen (2001), lista as unidades atômicas e
seus valores correspondentes em unidades do sistema internacional (SI).
Tabela 1 – Unidades Atômicas.
Quantidade Unidade Atômica Valor no SI Símbolo
Massa Massa do Elétron 9.1094 x 10
–31
kg
݉
Carga Carga do Elétron 1.6022 x 10
–19
C
݁
Constante de Planck Reduzida
ܥ݊ݏݐܽ݊ݐ݈݁݀݁ܲܽ݊ܿ݇ ʹߨ
Τ
1.0546 x 10
–34
J s
Comprimento
Ͷߨߝ
݉
݁
ଶ
Τ
5.2918 x 10
–11
m
ܽ
ሺ
ܾ݄ݎ
ሻ
Energia
ଶ
݉
ܽ
ଶ
Τ 4.3597 x 10
–18
J
ܧ
ሺ݄ܽݎݐݎ݁݁ሻ
A resolução exata da equação de Schrödinger, eq. (2.1), para o
Hamiltoniano da eq. (2.2) apresenta um nível de complexidade muito elevado.
Assim, são necessárias aproximações para viabilizar o seu emprego em sistemas
reais. Uma das aproximações mais importantes é a de Born-Oppenheimer, descrita
a seguir, que se baseia nos teoremas de Hohenberg e Kohn. Em seguida serão
apresentadas as equações de Kohn-Sham e as duas aproximações que possibilitam
suas soluções: aproximações dos potenciais de troca e correlação, e o
pseudopotencial.
2.2.
Aproximação de Born-Oppenheimer
O núcleo dos átomos tem uma massa muito superior à massa dos elétrons.
Portanto, a velocidade dos elétrons é muito maior do que a velocidade dos
PUC-Rio - Certificação Digital Nº 0812691/CA
Teoria do Funcional da Densidade 22
núcleos, de forma que os elétrons reagem quase que instantaneamente ao
movimento dos núcleos. Assim, a Aproximação de Born-Oppenheimer serve para
desacoplar o movimento dos núcleos e dos elétrons, podendo-se considerar que o
movimento eletrônico ocorre num campo nuclear fixo (Born e Oppenheimer,
1927).
Os elétrons são considerados os responsáveis pela energia cinética do
sistema e estão sujeitos à energia potencial devido às interações elétron-elétron e a
energia potencial externa, devido aos núcleos. Assim, para a solução do problema
eletrônico, utiliza-se o Hamiltoniano eletrônico simplificado com o seguinte
formato:
ܪ
ൌܶ
ܸ
ܸ
(2.9)
Uma vez que ܶ
ൌͲ e ܸ
ൌܿ݊ݏݐܽ݊ݐ݁.
A Aproximação de Born-Oppenheimer permite escrever a função de onda
total com sendo um produto da função de onda dos núcleos, ߰
, e da função de
onda dos elétrons, ߰
:
Ȱൌ߰
߰
(2.10)
A equação de Schrödinger para o problema eletrônico é escrita como:
ܪ
߰
ሺ
ܴ
ଵ
ǡǥǡܴ
Ǣݎ
ଵ
ߪ
ଵ
ǡǥǡݎ
ே
ߪ
ே
ሻ
ൌܧ
ሺ
ܴ
ଵ
ǡǥǡܴ
ሻ
߰
ሺ
ܴ
ଵ
ǡǥǡܴ
Ǣݎ
ଵ
ߪ
ଵ
ǡǥǡݎ
ே
ߪ
ே
ሻ
(2.11)
Onde as coordenadas dos núcleos, ܴ
, entram como parâmetros da função de
onda eletrônica e não como variáveis.
Assim, a aproximação de Born-Oppenheimer reduz o problema de muitos
corpos originalmente representado pelo Hamiltoniano da eq. (2.2). Porém, o
problema representado pelo Hamiltoniano eletrônico da eq. (2.9) ainda é de difícil
solução e entre os métodos mais usados atualmente para torná-lo tratável
computacionalmente estão os métodos baseados na Teoria do Funcional da
Densidade, tais como o empregado no programa SIESTA.
Como estamos usando a Teoria do Funcional da Densidade, a grandeza
fundamental do sistema deixa de ser a função de onda e passa a ser a densidade
eletrônica do sistema. O problema de 3N
e
variáveis, envolvendo a função de onda
para cada elétron, passa a ser um funcional da densidade, a qual é uma função real
das 3 variáveis de dimensão. Assim, podemos escrever a energia como sendo:
ܧ
ሾ
ߩ
ሿ
ൌܶ
ሾ
ߩ
ሿ
ܸ
ሾ
ߩ
ሿ
ܷ
ሾ
ߩ
ሿ
(2.12)
PUC-Rio - Certificação Digital Nº 0812691/CA
Teoria do Funcional da Densidade 23
Onde ܶ
ሾ
ߩ
ሿ
representa a energia cinética, ܸ
ሾ
ߩ
ሿ
o potencial de interação
elétron-núcleo e ܷ
ሾ
ߩ
ሿ
representa a interação elétron-elétron.
O estabelecimento formal da teoria DFT ocorreu em 1964 por meio de dois
teoremas, devidos a Hohenberg & Kohn, descritos a seguir.
2.3.
Teoremas de Hohenberg e Kohn
O trabalho de Hohenberg & Kohn (1964) mostra, através de dois teoremas,
que a partir da densidade eletrônica do sistema no estado fundamental é possível
obter a energia do estado fundamental de maneira exata.
O primeiro teorema estabelece que a densidade de carga ߩ
ሺ
ݎ
ሻ
do estado
fundamental de um sistema de muitos elétrons é determinada de maneira unívoca,
a menos de uma constante aditiva, a partir do potencial externo ݒ
௫௧
ሺ
ݎ
ሻ
.
Para provar este teorema é necessário supor o contrário, que dois potencias
externos ݒ
ሺ
ݎ
ሻ
e ݒ
ᇱ
ሺ
ݎ
ሻ
resultam na mesma densidade ߩ
ሺ
ݎ
ሻ
para o estado
fundamental. Dois hamiltonianos ܪ e ܪ
ᇱ
são definidos pelos dois potenciais ݒ
ሺ
ݎ
ሻ
e ݒ
ᇱ
ሺ
ݎ
ሻ
, associando-se a eles duas funções de onda ߮ e ߮
ᇱ
, respectivamente.
Assim, temos:
ۦ
߮
ȁ
ܪ
ȁ
߮
ۧ
ൌܧ
൏
ۦ
߮
ᇱ
ȁ
ܪ
ȁ
߮
ᇱ
ۧ
ൌ
ۦ
߮
ᇱ
ȁ
ܪ
ᇱ
െܸ
ᇱ
ܸ
ȁ
߮
ᇱ
ۧ
ൌ
ۦ
߮
ᇱ
ȁ
ܪ
ᇱ
ܸെܸ
ᇱ
ȁ
߮
ᇱ
ۧ
ܧ
൏
ۦ
߮
ᇱ
ȁ
ܪ
ᇱ
ܸെܸ
ᇱ
ȁ
߮
ᇱ
ۧ
ൌܧ
ᇱ
නߩ
ሺ
ݎ
ሻሾ
ݒ
ሺ
ݎ
ሻ
െݒ
ᇱ
ሺ
ݎ
ሻሿ
݀ݎ
ܧ
൏ܧ
ᇱ
නߩ
ሺ
ݎ
ሻሾ
ݒ
ሺ
ݎ
ሻ
െݒ
ᇱ
ሺ
ݎ
ሻሿ
݀ݎ
(2.13)
De forma semelhante, podemos fazer:
ۦ
߮
ᇱ
ȁ
ܪ
ᇱ
ȁ
߮
ᇱ
ۧ
ൌܧ
ᇱ
൏
ۦ
߮
ȁ
ܪ
ᇱ
ȁ
߮
ۧ
ൌ
ۦ
߮
ȁ
ܪെܸܸ
ᇱ
ȁ
߮
ۧ
ൌ
ۦ
߮
ȁ
ܪെ
ሺ
ܸെܸ
ᇱ
ሻȁ
߮
ۧ
ܧ
ᇱ
൏
ۦ
߮
ȁ
ܪെ
ሺ
ܸെܸ
ᇱ
ሻȁ
߮
ۧ
ൌܧ
െනߩ
ሺ
ݎ
ሻሾ
ݒ
ሺ
ݎ
ሻ
െݒ
ᇱ
ሺ
ݎ
ሻሿ
݀ݎ
ܧ
ᇱ
൏ܧ
െනߩ
ሺ
ݎ
ሻሾ
ݒ
ሺ
ݎ
ሻ
െݒ
ᇱ
ሺ
ݎ
ሻሿ
݀ݎ
(2.14)
Somando-se a eq. (2.13) com a eq. (2.14):
ܧ
ܧ
ᇱ
൏ܧ
ᇱ
ܧ
(2.15)
Portanto, a suposição inicial de que os potencias externos ݒ
ሺ
ݎ
ሻ
e ݒ
ᇱ
ሺ
ݎ
ሻ
resultavam na mesma densidade ߩ
ሺ
ݎ
ሻ
para o estado fundamental não é verdadeira.
Assim, existe uma correspondência única entre a densidade ߩ
ሺ
ݎ
ሻ
do estado
fundamental e o potencial externo ݒ
௫௧
ሺ
ݎ
ሻ
.
PUC-Rio - Certificação Digital Nº 0812691/CA
Teoria do Funcional da Densidade 24
O segundo Teorema diz que a energia do estado fundamental corresponde
ao mínimo do funcional de energia ܧ
ሾ
ߩ
ሺ
ݎ
ሻሿ
, obtido a partir da densidade exata
do estado fundamental ߩ
ሺ
ݎ
ሻ
. Qualquer densidade diferente ߩ
ᇱ
ሺ
ݎ
ሻ
്ߩ
ሺ
ݎ
ሻ
conduzirá a uma energia ܧ
ሾ
ߩ
ᇱ
ሺ
ݎ
ሻሿ
maior do que a energia do estado fundamental,
ܧ
ሾ
ߩ
ᇱ
ሺ
ݎ
ሻሿ
ܧ
ሾ
ߩ
ሺ
ݎ
ሻሿ
.
Este teorema torna possível o uso do princípio variacional para encontrar a
densidade do estado fundamental. Existem muitas possibilidades para a densidade
eletrônica, porém o problema é resolvido por minimização. Ao encontrar a
densidade para a qual a energia é mínima, encontrou-se a densidade do estado
fundamental.
2.4.
Equações de Kohn-Sham
A equação de Kohn-Sham pode ser vista como sendo a equação de
Schrödinger para um sistema fictício composto de partículas não interagentes, tal
que a densidade eletrônica gerada seja a mesma do sistema real original composto
pelas partículas reais interagentes. A equação de Kohn-Sham é definida por um
potencial externo efetivo ݒ
ሺ
ݎ
ሻ
no qual as partículas não interagentes se
movem.
ܪ
ௌ
߰
ሺ
ݎ
ሻ
ൌെ
ͳ
ʹ
ଶ
ݒ
ሺ
ݎ
ሻ
൨߰
ሺ
ݎ
ሻ
ൌߝ
߰
ሺ
ݎ
ሻ
(2.16)
Onde as funções ߰
ሺ
ݎ
ሻ
são as autofunções da equação de Kohn-Sham e o
potencial externo efetivo ݒ
ሺ
ݎ
ሻ
é dado por:
ݒ
ሺ
ݎ
ሻ
ൌݒ
ሺ
ݎ
ሻ
න
ߩ
ሺ
ݎ
ሻ
ȁ
ݎെݎ
ᇱ
ȁ
݀ݎݒ
௫
ሺ
ݎ
ሻ
(2.17)
A eq. (2.18) foi a solução proposta por Kohn & Sham (1965) para
representar a densidade de carga:
ߩ
ሺ
ݎ
ሻ
ൌ
ȁ
߰
ሺ
ݎ
ሻȁ
ଶ
ே
ୀଵ
(2.18)
A energia total do estado fundamental da eq. (2.12) pode ser reescrita como:
ܧ
ሾ
ߩ
ሿ
ൌܶ
ሾ
ߩ
ሿ
නߩ
ሺ
ݎ
ሻ
ݒ
ሺ
ݎ
ሻ
݀ݎ
ͳ
ʹ
නන
ߩ
ሺ
ݎ
ሻ
ߩ
ሺ
ݎ
ᇱ
ሻ
ȁ
ݎെݎ
ᇱ
ȁ
݀ݎ݀ݎ
ᇱ
(2.19)
Na aproximação de Kohn-Sham, ܶ
ሾ
ߩ
ሿ
é separada em duas componentes.
Uma das componentes corresponde à energia cinética de um gás de elétrons não
PUC-Rio - Certificação Digital Nº 0812691/CA
Teoria do Funcional da Densidade 25
interagentes ܶ
ௌ
ሾ
ߩ
ሿ
, e a outra componente inclui as interações eletrônicas e a
correção da energia cinética no termo de troca e correlação ܧ
௫
ሾ
ߩ
ሿ
. Assim, a
eq. (2.19) se torna:
ܧ
ሾ
ߩ
ሿ
ൌܶ
ௌ
ሾ
ߩ
ሿ
නߩ
ሺ
ݎ
ሻ
ݒ
ሺ
ݎ
ሻ
݀ݎ
ͳ
ʹ
නන
ߩ
ሺ
ݎ
ሻ
ߩ
ሺ
ݎ
ᇱ
ሻ
ȁ
ݎെݎ
ᇱ
ȁ
݀ݎ݀ݎ
ᇱ
ܧ
௫
ሾ
ߩ
ሿ
(2.20)
Isolando ݒ
ሺ
ݎ
ሻ
da eq. (2.17) e substituindo na eq. (2.20), teremos:
ܧ
ሾ
ߩ
ሿ
ൌܶ
ௌ
ሾ
ߩ
ሿ
නߩ
ሺ
ݎ
ሻ
ݒ
ሺ
ݎ
ሻ
݀ݎെනන
ߩ
ሺ
ݎ
ሻ
ߩ
ሺ
ݎ
ᇱ
ሻ
ȁ
ݎെݎ
ᇱ
ȁ
݀ݎ݀ݎ
ᇱ
െනߩ
ሺ
ݎ
ሻ
ݒ
௫
ሺ
ݎ
ሻ
݀ݎ
ͳ
ʹ
නන
ߩ
ሺ
ݎ
ሻ
ߩ
ሺ
ݎ
ᇱ
ሻ
ȁ
ݎെݎ
ᇱ
ȁ
݀ݎ݀ݎ
ᇱ
ܧ
௫
ሾ
ߩ
ሿ
(2.21)
O somatório dos autovalores ߝ
da eq. (2.16) corresponde aos dois primeiros
termos do lado direito da eq. (2.21).
ߝ
ൌܶ
ௌ
ሾ
ߩ
ሿ
නߩ
ሺ
ݎ
ሻ
ݒ
ሺ
ݎ
ሻ
݀ݎ
(2.22)
Portanto, a forma final para a energia total do estado fundamental na
aproximação de Kohn-Sham é dada por:
ܧ
ሾ
ߩ
ሿ
ൌߝ
െ
ͳ
ʹ
නන
ߩ
ሺ
ݎ
ሻ
ߩ
ሺ
ݎ
ᇱ
ሻ
ȁ
ݎെݎ
ᇱ
ȁ
݀ݎ݀ݎ
ᇱ
െනߩ
ሺ
ݎ
ሻ
ݒ
௫
ሺ
ݎ
ሻ
݀ݎܧ
௫
ሾ
ߩ
ሿ
(2.23)
Onde ݒ
௫
ሺ
ݎ
ሻ
é o potencial de troca e correlação e sua definição formal vem
da seguinte derivada funcional:
ݒ
௫
ሺ
ݎ
ሻ
ൌ
߲ܧ
௫
ሾ
ߩ
ሿ
߲ߩ
ሺ
ݎ
ሻ
(2.24)
2.5.
Aproximações dos Potenciais de Troca e Correlação
A expressão exata do funcional da energia de troca e correlação ܧ
௫
ሾ
ߩ
ሿ
não
é conhecida. Assim, para que seja possível utilizar as Equações de Kohn-Sham é
necessário determinar uma boa aproximação para o termo de troca e correlação,
que é o termo de interpretação física mais difícil da DFT. Entre as aproximações
mais utilizadas para este termo desconhecido estão a Aproximação da Densidade
Local (LDA) e um aperfeiçoamento chamado de Aproximação do Gradiente
Generalizado (GGA).
PUC-Rio - Certificação Digital Nº 0812691/CA
Teoria do Funcional da Densidade 26
2.5.1.
Aproximação da Densidade Local
A energia de troca e correlação na Aproximação da Densidade Local (LDA)
é escrita como:
ܧ
௫
ሾ
ߩ
ሺ
ݎ
ሻሿ
ൌනߩ
ሺ
ݎ
ሻ
݁
௫
ሾ
ߩ
ሺ
ݎ
ሻሿ
݀ݎ
(2.25)
Onde ߩ
ሺ
ݎ
ሻ
é a densidade eletrônica no ponto ݎ e ݁
௫
ሾ
ߩ
ሺ
ݎ
ሻሿ
é a energia de
troca e correlação por partícula de um gás homogêneo de elétrons com
densidade ߩ.
Na LDA o termo de troca e correlação é subdivido em:
݁
௫
ሺ
ߩ
ሻ
ൌ݁
௫
ሺ
ߩ
ሻ
݁
ሺ
ߩ
ሻ
(2.26)
Onde a energia de troca ݁
௫
ሺ
ߩ
ሻ
por partícula é dada por:
݁
௫
ሺ
ߩ
ሻ
ൌെ
ͲǤͶͷͺ
ݎ
௦
(2.27)
E o termo referente à correlação eletrônica ݁
ሺ
ߩ
ሻ
foi estimado por Ceperley
e Alder (1980):
݁
ሺ
ߩ
ሻ
ൌെ
ͲǤͶͶ
ݎ
௦
Ǥͺ
(2.28)
Onde ݎ
௦
é o raio da esfera cujo volume é igual ao volume por elétron de
condução:
ݎ
௦
ൌ
൬
Ͷߨ
͵
൰
ߩ൨
ି
ଵ
ଷ
(2.29)
A LDA é adequada para descrever bem sistemas onde a densidade
eletrônica é aproximadamente uniforme. Se a densidade eletrônica for fortemente
não uniforme, será necessário utilizar a Aproximação do Gradiente Generalizado.
2.5.2.
Aproximação do Gradiente Generalizado
A Aproximação do Gradiente Generalizado (GGA) considera, além da
densidade eletrônica ߩ
ሺ
ݎ
ሻ
no ponto ݎ, o gradiente da densidade eletrônica ߩ
ሺ
ݎ
ሻ
neste ponto, onde a densidade de energia de troca e correlação está sendo
calculada. Assim, o termo de troca e correlação é escrito como:
ܧ
௫
ሾ
ߩ
ሺ
ݎ
ሻሿ
ൌන݁
௫
ሾ
ߩ
ሺ
ݎ
ሻ
ǡߩ
ሺ
ݎ
ሻሿ
݀ݎ
(2.30)
PUC-Rio - Certificação Digital Nº 0812691/CA
Teoria do Funcional da Densidade 27
Ao contrário da aproximação LDA, onde existe um único ݁
௫
ሾ
ߩ
ሺ
ݎ
ሻሿ
correto,
na aproximação GGA existem diferentes parametrizações para ݁
௫
ሾ
ߩ
ሺ
ݎ
ሻ
ǡߩ
ሺ
ݎ
ሻሿ
que originam funcionais diferentes. A parametrização GGA mais usual é a
desenvolvida por Perdew, Burke & Ernzerhof (1996).
2.6.
Pseudopotencial
Um átomo pode ser visto como sendo composto por um caroço interno e
pelos elétrons de valência. O caroço interno é composto pelo núcleo atômico e
pelos elétrons mais próximos do núcleo, os quais sentem um forte potencial
atrativo e tem pouca participação nas ligações químicas. Os elétrons de valência
estão fracamente ligados ao núcleo e por isso apresentam grande participação nas
ligações químicas, determinando a maior parte das propriedades físicas de um
sólido ou molécula.
O Pseudopotencial é um potencial efetivo construído pelo programa ATOM
(García, 2006) e usado para substituir o potencial real gerado pelo conjunto dos
prótons e elétrons próximos ao núcleo. Os estados eletrônicos de caroço são
eliminados e os elétrons de valência são descritos por uma pseudofunção de onda
sem nodos. Isto reduz o custo computacional simplificando os cálculos de
estrutura eletrônica.
Pseudopotenciais de norma-conservada são aqueles que satisfazem as
seguintes condições (Hamann, Schlüter e Chiang, 1979):
1. O hamiltoniano real e o pseudo-hamiltoniano devem possuir os
mesmos autovalores.
2. A função de onda de todos os elétrons normalizada (AE) e a
pseudofunção de onda de valência normalizada (PS) devem ser
iguais para ݎݎ
, onde ݎ
é o raio de corte.
3. A integral da densidade de carga das funções de onda AE e PS
devem ser iguais para ݎݎ
.
4. As derivadas logarítmicas das funções de onda AE e PS devem
coincidir para ݎݎ
.
PUC-Rio - Certificação Digital Nº 0812691/CA
Teoria do Funcional da Densidade 28
A função de onda de todos os elétrons normalizada (AE) deve coincidir com
os valores da pseudofunção de onda de valência normalizada (PS) para ݎݎ
,
conforme o ilustrado na Figura 1.
Figura 1 – Função de onda de todos os elétrons normalizada (AE) e a
pseudofunção de onda de valência normalizada (PS).
Para construir um pseudopotencial capaz de representar a ação do caroço
iônico é necessário resolver autoconsistentemente a equação radial de Kohn-Sham
para todos os elétrons do átomo isolado (Kohn e Sham, 1965):
ቈ
െ
ͳ
ʹ
݀
ଶ
݀ݎ
ଶ
ݎ
݈
ሺ
݈ͳ
ሻ
ʹݎ
ଶ
െ
ܼ
ݎ
ܸ
ு
ܸ
௫
ݎܴ
ሺ
ݎ
ሻ
ൌ߳
ܴ
ሺ
ݎ
ሻ
(2.31)
Um pseudopotencial blindado é definido para o átomo isolado e no estado
neutro. Este pseudopotencial pode ser construído usando-se a as soluções
encontradas para ܴ
ሺ
ݎ
ሻ
:
ܸ
ǡ
௦
ሺ
ݎ
ሻ
ൌ߳
െ
݈
ሺ
݈ͳ
ሻ
ʹݎ
ଶ
ͳ
ʹݎܴ
௦
ሺ
ݎ
ሻ
݀
ଶ
݀ݎ
ଶ
ቀݎܴ
௦
ሺ
ݎ
ሻ
ቁ
(2.32)
Um potencial iônico pode ser gerado removendo-se os potenciais ܸ
ு
e ܸ
௫
do potencial blindado ܸ
ǡ
௦
. Este potencial iônico poderá ser usado num
procedimento autoconsistente para determinar a blindagem eletrônica dos elétrons
de valência em ambientes diferentes do considerado para calcular a pseudofunção
de onda ܴ
ሺ
ݎ
ሻ
.
ܸ
ǡ
௦
ሺ
ݎ
ሻ
ൌܸ
ǡ
௦
ሺ
ݎ
ሻ
െܸ
ு
௦
ሺ
ݎ
ሻ
െܸ
௫
௦
ሺ
ݎ
ሻ
(2.33)
O potencial da eq. (2.33) pode ser reescrito como sendo composto por um
termo local e um termo semilocal:
ܸ
௦
ሺ
ݎ
ሻ
ൌܸ
ǡ
௦
ሺ
ݎ
ሻ
ܸ
௦ǡ
ሺ
ݎ
ሻ
ܲ
(2.34)
PUC-Rio - Certificação Digital Nº 0812691/CA
Teoria do Funcional da Densidade 29
Onde ܲ
é um operador de projeção da componente de momento angular
ȁ
݈
ۄۃ
݈
ȁ
da função de onda.
O procedimento descrito por Kleinman & Bylander (1982) pode ser usado
para reescrever o termo de potencial semilocal para uma forma não-local:
ܸ
ିǡ
ሺ
ݎ
ሻ
ൌ
หܸ
௦ǡ
߯
ۄۃ
߯
ܸ
௦ǡ
ห
ൻ
߯
หܸ
௦ǡ
ห
߯
ൿ
(2.35)
A forma final do pseudopotencial é dada pela eq. (2.36):
ܸ
ௌ
ሺ
ݎ
ሻ
ൌܸ
ǡ
௦
ሺ
ݎ
ሻ
ܸ
ିǡ
ሺ
ݎ
ሻ
(2.36)
O pseudopotencial de Troullier & Martins (1991) é um dos modelos mais
usualmente empregados. É considerado como sendo um pseudopotencial suave,
pois a energia total do sistema e suas propriedades físicas apresentam uma
convergência rápida com o aumento das funções de base (Fagan, 2003). A
utilização do programa ATOM é uma forma prática para geração de
pseudopotenciais do tipo Troullier-Martins (García, 2006).
2.7.
Hamiltoniano Empregado no SIESTA
As aproximações descritas para a Teoria do Funcional da Densidade
permitem escrever o hamiltoniano de Kohn-Sham como sendo (Kohn e Sham,
1965):
ܪ
ௌ
ൌ
ቈ
െ
ͳ
ʹ
ଶ
ܸ
ௌ
ሺ
ݎ
ሻ
න
ߩ
ሺ
ݎ
ሻ
ȁ
ݎെݎ
ᇱ
ȁ
݀ݎܸ
௫
ሺ
ݎ
ሻ
(2.37)
Usando o procedimento de Kleinman & Bylander (1982) para separar o
potencial ܸ
ௌ
ሺ
ݎ
ሻ
em componentes local (longo alcance) e não-local (curto
alcance):
ܪ
ௌ
ൌ
ቈ
െ
ͳ
ʹ
ଶ
ܸ
ǡ
௦
ሺ
ݎ
ሻ
ܸ
ିǡ
ሺ
ݎ
ሻ
න
ߩ
ሺ
ݎ
ሻ
ȁ
ݎെݎ
ᇱ
ȁ
݀ݎܸ
௫
ሺ
ݎ
ሻ
(2.38)
Utilizando-se índices atômicos, a eq. (2.38) pode ser reescrita na forma
equivalente:
ܪ
ௌ
ൌ
ܶ
ௌ
ܸ
ூ
ሺ
ݎ
ሻ
୍
ܸ
ூ
୍
ܸ
ு
ሺ
ݎ
ሻ
ܸ
௫
ሺ
ݎ
ሻ
൩
(2.39)
PUC-Rio - Certificação Digital Nº 0812691/CA
Teoria do Funcional da Densidade 30
Onde o índice atômico é representado por ,
σ
ܸ
ூ
ሺ
ݎ
ሻ
୍
e
σ
ܸ
ூ
ሺ
ݎ
ሻ
୍
são,
respectivamente, as contribuições locais e não-locais do pseudopotencial, ܸ
ு
ሺ
ݎ
ሻ
é
o potencial de Hartree e ܸ
௫
ሺ
ݎ
ሻ
é o potencial de troca e correlação.
O termo local do pseudopotencial é um operador de longo alcance que pode
ser custoso em termos computacionais. Assim, o método empregado no SIESTA
faz a substituição do potencial ܸ
ூ
pelo potencial ܸ
ூ
ே
, que é definido por
ܸ
ூ
ே
ൌܸ
ூ
ܸ
ூ
௧Ø
(Soler, Artacho, et al., 2002; Ordejón, Artacho e Soler,
1996). O termo ܸ
ூ
௧Ø
é o potencial eletrônico criado por uma distribuição de
densidade eletrônica ߩ
௧Ø
. Desta forma, a eq. (2.39) pode se escrita de forma
equivalente:
ܪൌ
ܶ
ௌ
ܸ
ூ
୍
ቀܸ
ூ
ே
ሺ
ݎ
ሻ
െܸ
ூ
௧Ø
ሺ
ݎ
ሻ
ቁ
୍
ܸ
ு
ሺ
ݎ
ሻ
ܸ
௫
ሺ
ݎ
ሻ
൩
(2.40)
O potencial de Hartree ܸ
ு
ሺ
ݎ
ሻ
é reescrito como sendo:
ܸ
ு
ሺ
ݎ
ሻ
ൌߜܸ
ு
ሺ
ݎ
ሻ
ܸ
ூ
௧Ø
ሺ
ݎ
ሻ
(2.41)
Assim, o hamiltoniano usado no SIESTA assume a sua forma final:
ܪൌ
ܶ
ௌ
ܸ
ூ
୍
ܸ
ூ
ே
ሺ
ݎ
ሻ
୍
ߜܸ
ு
ሺ
ݎ
ሻ
ܸ
௫
ሺ
ݎ
ሻ
൩
(2.42)
O termo ܸ
ூ
ே
é nulo para ݎݎ
. Desta forma, conseguiu-se uma formulação
onde os termos de longo alcance foram eliminados. Os dois primeiros termos da
eq. (2.42) envolvem integrais de dois centros dos orbitais atômicos calculados e
armazenados em função da distância interatômica. Os três últimos termos são
calculados através de somas discretas numa matriz de densidades tridimensional
do espaço real.
2.7.1.
Funções de Base
Os orbitais ߰
ሺ
ݎ
ሻ
da equação de Kohn-Sham, eq. (2.16), precisam ser
expressos através de uma base para a realização dos cálculos. No programa
SIESTA são empregadas funções de base atômicas localizadas, as quais possuem
dois parâmetros importantes, número de orbitais por átomo e o alcance destes
orbitais (Soler, Artacho, et al., 2002; Anglada, Soler, et al., 2002). Com este tipo
de base é possível obter boa precisão de cálculos com baixo custo computacional.
PUC-Rio - Certificação Digital Nº 0812691/CA
Teoria do Funcional da Densidade 31
A base mínima SZ (single-ߞ), possui somente uma função radial para cada
canal de momento angular e somente para os estados de valência do átomo
isolado. Ela permite a realização de cálculos rápidos com um grande número de
átomos para obtenção de tendências qualitativas das ligações químicas e uma boa
descrição da banda de valência.
A base DZ (double-ߞ), possui duas funções radiais por canal de momento
angular permitindo uma melhor descrição da parte radial. É possível adicionar
flexibilização angular através das chamadas funções de polarização. Assim
podemos ter as bases SZP e DZP, representadas a partir das duas anteriormente
descritas (Artacho, Gale, et al., 2008).
Os raios das funções de base confinadas são definidos de maneira
sistemática através de um único parâmetro, a correção de energia (energy shift). O
seu significado é a quantidade de energia incrementada a um orbital devido ao
fato de estar confinado. Todos os raios de corte são escolhidos de maneira que o
incremento de energia seja sempre o mesmo, gerando uma base equilibrada.
2.7.2.
Orbitais Atômicos Numéricos
O programa SIESTA utiliza Orbitais Atômicos Numéricos (NAO), os quais
são obtidos pela resolução da equação de Kohn-Sham para os pseudo-átomos
isolados (Junquera, Paz, et al., 2001; Artacho, Sánchez-Portal, et al., 1999). O
confinamento das funções de base é obtido estabelecendo-se um raio de corte no
ponto onde os valores são suficientemente pequenos para serem desprezados.
Os orbitais atômicos confinados apresentam como vantagem o fato de que
as interações se estendem numa região finita de átomos vizinhos. Desta forma,
obtêm-se melhor eficiência computacional, pois não é necessário calcular a
interação de um átomo com os seus vizinhos que estejam além do raio de corte
dos seus orbitais. Isso também dá origem a matrizes esparsas de hamiltoniano e de
sobreposição, as quais podem ser mais eficientemente armazenadas na memória.
PUC-Rio - Certificação Digital Nº 0812691/CA
3
Computação de Propósito Geral em Unidades de
Processamento Gráfico
As Unidades de Processamento Gráfico (GPUs) foram originalmente
desenvolvidas para o processamento de gráficos e eram difíceis de programar.
Contudo, o desenvolvimento de interfaces de programação e o suporte a
linguagens como o C permitiram que as GPUs fossem utilizadas para Computação
de Propósito Geral (GPGPU) (Harris, 2005). Atualmente as GPUs são
processadores massivamente paralelos, compostos por um grande número de
elementos de processamento num único chip, com alto desempenho em operações
de ponto flutuante.
Neste Capítulo descreveremos inicialmente alguns conceitos fundamentais
da programação paralela, passando em seguida para o modelo de programação
paralela desenvolvido especificamente para placas gráficas da NVIDIA,
denominada de CUDA (Compute Unified Device Architecture).
3.1.
Programação Paralela
A programação paralela consiste em subdividir um problema em partes
independentes de maneira que seja possível utilizar múltiplos elementos de
processamento para executá-las simultaneamente (Barney, 2009). O interesse pela
computação paralela tem aumentado devido às restrições físicas que impedem o
aumento crescente da freqüência.
A Lei de Moore é uma observação empírica de que a densidade de
transistores duplica a cada 18 ou 24 meses (Moore, 1998). A dificuldade de
aumentar a freqüência faz com que os transistores adicionais sejam usados para
adicionar hardware extra para computação paralela.
O paralelismo pode ocorrer em diversos níveis. O Paralelismo em Nível de
Bit significa aumentar o tamanho da palavra, reduzindo a quantidade de instruções
que um processador deve executar para realizar uma operação em variáveis cujo
PUC-Rio - Certificação Digital Nº 0812691/CB
Computação de Propósito Geral em Unidades de Processamento Gráfico 33
tamanho é maior do que o da palavra. Um processador de 32 bits teria que
executar duas operações para realizar a soma de dois números inteiros com
precisão de 64 bits, enquanto um processador de 64 bits poderia fazê-lo em uma
única operação.
O Paralelismo em Nível de Instrução consiste em reordenar o fluxo de
instruções e reagrupá-las, de modo que estes grupos possam ser executados em
paralelo sem alterar o resultado do programa. O pipeline dos processadores
modernos possui múltiplos estágios. Cada estágio corresponde a uma ação
diferente que o processador executa em determinada instrução; um processador
com um pipeline de N estágios pode ter até N diferentes instruções em diferentes
estágios de execução.
O Paralelismo de Dados requer a distribuição dos dados por diferentes nós
computacionais de maneira a viabilizar a execução de uma mesma instrução sobre
diferentes conjuntos de dados simultaneamente (Hillis e Steele, 1986). O
paralelismo de dados é inerente aos laços de repetição e só pode ser empregado se
não houver a dependência de uma iteração do laço com a saída de uma ou mais
iterações anteriores.
Diferentemente do Paralelismo de Dados, no Paralelismo de Tarefas
operações diferentes podem ser executadas sobre um mesmo conjunto de dados ou
sobre conjuntos de dado diferentes.
3.1.1.
Taxonomia de Flynn
A Taxonomia de Flynn é um modelo de classificação de arquitetura de
computadores baseado no fluxo de instruções e dados (Flynn, 1972; Duncan,
1990). Esta classificação é dividida em quatro categorias: SISD, SIMD, MISD e
MIMD.
A classificação SISD (Single Instruction, Single Data) equivale a um
programa puramente seqüencial, sendo também conhecida como arquitetura Von
Neumann. Um fluxo único de instruções é aplicado sobre um conjunto único de
dados conforme o ilustrado na Figura 2. Na ilustração UC é utilizado para
representar a unidade de controle e UP é usado para representar a unidade de
processamento.
PUC-Rio - Certificação Digital Nº 0812691/CB
Computação de Propósito Geral em Unidades de Proces
Figura 2 –
Arquitetura SISD
O modelo das
arquiteturas vetoriais
múltiplos operando
s simultaneamente, é denominada
Multiple Data). A Figura 3
Figura 3 –
Arquitetura SIMD
Na classificação
MISD (
unidades de processamento realiz
conjunto de dados. Conforme o ilustrado na
dados
, uma unidade de processamento realiza uma operação
passa à unidade de processamento seguinte que execu
sobre o mesmo dado.
Computação de Propósito Geral em Unidades de Processamento G
ráfico
Arquitetura SISD
.
arquiteturas vetoriais
, onde uma
operação é executada sobre
s simultaneamente, é denominada
SIMD (
Single Instruction,
ilustra esta arquitetura.
Arquitetura SIMD
.
MISD (
Multiple Instruction, Single Data)
, múltiplas
unidades de processamento realiz
am operações diferentes sobre
um mesmo
conjunto de dados. Conforme o ilustrado na
Figura 4, numa m
áquina de fluxo de
, uma unidade de processamento realiza uma operação sobre o dado e o
passa à unidade de processamento seguinte que executará uma operação diferente
34
operação é executada sobre
Single Instruction,
, múltiplas
um mesmo
áquina de fluxo de
sobre o dado e o
tará uma operação diferente
PUC-Rio - Certificação Digital Nº 0812691/CB
Computação de Propósito Geral em Unidades de Proces
Figura 4 –
Arquitetura MISD
A arquitetura
MIMD (
multiprocessadores é apresentada na
instruções são execu
tadas sobre diferentes conjuntos de dados simultane
usando
unidades de processamento diferentes controladas po
controle independentes.
Figura 5 –
Arquitetura MIMD
3.1.2.
Lei de Amdahl
O ganho de velocidade de um a
processadores é limitado pelo tempo consumido pela
Computação de Propósito Geral em Unidades de Processamento G
ráfico
Arquitetura MISD
.
MIMD (
Multiple Instruction, Multiple Data)
empregada em
multiprocessadores é apresentada na
Figura 5
. Nesta arquitetura, diferentes
tadas sobre diferentes conjuntos de dados simultane
unidades de processamento diferentes controladas por unidades de
Arquitetura MIMD
.
O ganho de velocidade de um a
lgoritmo rodando em múltiplos
processadores é limitado pelo tempo consumido pela fração seqüencial do
35
empregada em
. Nesta arquitetura, diferentes
tadas sobre diferentes conjuntos de dados simultaneamente,
r unidades de
lgoritmo rodando em múltiplos
fração seqüencial do
PUC-Rio - Certificação Digital Nº 0812691/CB
Computação de Propósito Geral em Unidades de Processamento Gráfico 36
programa. Assim, o máximo ganho de velocidade esperado quando um sistema
não pode ser totalmente paralelizado é dado pela Lei de Amdahl (1967):
ܵ
ൌ
ͳ
ͳ
െ
ܲ
ܲ
ܰ
(3.1)
Onde ܵ é o ganho de velocidade do programa, ܲ é a fração do programa que
pode ser paralelizada e ܰ é o número de processadores.
A Figura 6 ilustra o efeito da Lei de Amdahl para os casos em que uma
fração de 50%, 75%, 90% ou 95% do programa podem ser efetivamente
paralelizados.
Figura 6 – Efeito da Lei de Amdahl sobre o ganho de velocidade de algoritmo
rodando em múltiplos processadores.
3.2.
Compute Unified Device Architecture (CUDA)
Em Novembro de 2006, a NVIDIA apresentou uma extensão da linguagem
C para programação em GPUs, denominada CUDA (Compute Unified Device
Architecture). Esta linguagem simplificada possibilitou a utilização da alta
capacidade de cálculos das GPUs para solucionar problemas de propósito geral no
meio científico (NVIDIA, 2010).
10
0
10
1
10
2
10
3
10
4
10
5
0
2
4
6
8
10
12
14
16
18
20
Lei de Amdahl
Número de Processadores
Ganho de Velocidade
95 %
90 %
75 %
50 %
PUC-Rio - Certificação Digital Nº 0812691/CB
Computação de Propósito Geral em Unidades de Processamento Gráfico 37
A demanda, no mercado de jogos, por gráficos 3D de alta definição em
tempo real estimulou o desenvolvimento das GPUs. As GPUs se tornaram
processadores de múltiplos elementos paralelos (NVIDIA, 2009b). A Figura 7
ilustra a evolução do desempenho em operações de ponto flutuante das GPUs
comparada com CPUs nos últimos anos. A taxa de transferência de dados através
do barramento da memória, ou largura de banda da memória, também teve um
crescimento mais acelerado nas GPUs, conforme é apresentado na Figura 8.
Figura 7 – Operações de Ponto Flutuante por Segundo (NVIDIA, 2009b).
PUC-Rio - Certificação Digital Nº 0812691/CB
Computação de Propósito Geral em Unidades de Processamento Gráfico 38
Figura 8 – Largura de Banda (NVIDIA, 2009b).
As GPUs, que eram originalmente projetadas para resolver o problema da
renderização gráfica, seguiram um caminho de evolução que as tornaram
especializadas para o processamento paralelo de dados. Uma maior quantidade de
transistores é utilizada para processamento, reduzindo o espaço para controle de
fluxo e memória cache. A CPU, por outro lado, apresenta um mecanismo de
controle de fluxo mais sofisticado e maior área do chip dedicada a cache. Estas
diferenças são ilustradas na Figura 9.
Figura 9 – A GPU dedica mais transistores para o processamento de dados
(NVIDIA, 2009b).
Algumas das limitações do CUDA são: (i) impossibilidade de utilizar
recursão e ponteiro para função; (ii) a latência de acesso à memória, ocasionada
PUC-Rio - Certificação Digital Nº 0812691/CB
Computação de Propósito Geral em Unidades de Processamento Gráfico 39
pela falta de cache, pode ser minimizada pela divisão do problema de forma a
aproveitar a capacidade de administrar um elevado número de threads. Enquanto
um grupo de 32 threads, denominado warp, espera pelo dado a ser lido da
memória, os recursos do hardware são usados para realizar operações em outros
warps.
3.3.
Elementos de Hardware da GPU
Os principais elementos do hardware da GPU são o Processador de Fluxo
(SP), o Multiprocessador de Fluxo (SM) e os diferentes tipos de memórias
disponíveis, os quais serão descritos em maiores detalhes nos próximos itens.
3.3.1.
Processador e Multiprocessador de Fluxo
Um multiprocessador é constituído por oito processadores escalares, uma
unidade de instrução, memória compartilhada interna ao chip e duas unidades
especiais para execução de operações complexas, como seno, cosseno,
exponencial e logaritmo.
Cada um dos oito processadores possui uma unidade de ponto flutuante de
32 bits, capaz de realizar operações de ponto flutuante e de inteiros, tais como
multiplicação, adição, multiplicação-adição combinadas numa operação,
subtração, comparação, etc. Geralmente uma placa de elevado nível possui vários
multiprocessadores. Podemos citar como exemplo a Tesla C1060, a qual possui
30 multiprocessadores, totalizando 240 processadores.
O multiprocessador foi projetado para criar, administrar e executar threads
concorrentes com um custo mínimo devido à overhead e sincronização. As
barreiras de sincronização consomem uma única instrução. Uma instrução emitida
pela unidade de instruções é executada por até 32 threads. Isso possibilita um
paralelismo fino, sendo possível atribuir uma thread para cada elemento de dado,
como por exemplo, atribuir uma thread para cada pixel de uma imagem
(NVIDIA, 2009b).
PUC-Rio - Certificação Digital Nº 0812691/CB
Computação de Propósito Geral em Unidades de Processamento Gráfico 40
3.3.2.
Tipos de Memória da GPU
A memória de maior capacidade e também a de maior latência é a memória
global, ou memória da placa. Esta memória aceita operações de leitura e escrita
provenientes tanto da CPU, quanto da GPU, sendo a principal memória de troca
de dados entre CPU e GPU. A latência de acesso a esta memória varia entre 400 e
600 ciclos de clock. Uma GPU Tesla C1060 possui 4 GB de memória RAM
GDDR3.
A memória compartilhada é uma memória interna ao chip, o que a torna
rápida, porém possui um tamanho bastante limitado de apenas 16 KB por
multiprocessador. É útil para a troca de dados entre threads de um mesmo bloco.
A latência pode ser de apenas um ou dois ciclos se forem empregadas técnicas
adequadas de programação, evitando conflitos de bancos de memória.
A memória de textura é uma memória somente leitura para a GPU, devendo
ser gravada pela CPU, que possui cerca de 6 a 8 KB de memória cache por
multiprocessador. Ela é implementada no mesmo espaço físico da memória global
e, portanto, em casos de falta de dados na memória cache, o tempo de latência é o
mesmo da memória global. Possui ainda modos de endereçamento diferentes ou
filtros para alguns formatos específicos de dados.
A memória de constantes é uma memória somente leitura com capacidade
total de 64 KB e possui cache. A maneira mais rápida de acessar este tipo de
memória ocorre quando todas as threads de um half-warp (grupo de 16 threads)
precisam ler um mesmo elemento. Neste caso a latência será de um ou dois ciclos.
Se for necessário fazer a leitura de múltiplos elementos, a latência será
proporcionalmente maior, em relação ao número de elementos diferentes lidos
pelas threads de um mesmo half-warp (NVIDIA, 2009b).
Uma GPU com capacidade de cálculo 1.3 possui 16384 registradores de 32
bits por multiprocessador. Considerando que neste tipo de GPU o número máximo
permitido de threads ativas é de 1024 por multiprocessador, então teremos um
mínimo de 16 registradores por thread. Os diferentes tipos de memória
empregados numa GPU são apresentados na Figura 10.
PUC-Rio - Certificação Digital Nº 0812691/CB
Computação de Propósito Geral em Unidades de Processamento Gráfico 41
Figura 10 – Tipos de memória das GPUs da NVIDIA (NVIDIA, 2009b).
3.4.
Elementos de Software da GPU
Um algoritmo empregando GPGPU consiste de um fluxo de controle
rodando seqüencialmente numa CPU, operações de transferência de dados entre a
memória da CPU e da GPU e chamadas de kernel (função escrita em CUDA para
GPU) que executarão a tarefa usando o paralelismo de dados na GPU.
O primeiro passo do algoritmo de controle, rodando na CPU, é a chamada
de funções para realizar alocações e transferências dos dados para a memória da
GPU. As instruções CUDA são então enviadas para a GPU e executadas sobre os
dados. Quando o kernel termina sua execução, é necessário realizar a transferência
dos resultados para a memória da CPU, completando o processo ilustrado na
Figura 11.
PUC-Rio - Certificação Digital Nº 0812691/CB
Computação de Propósito Geral em Unidades de Proces
Figura 11 –
Fluxo de processamento de um software empregando GP
A execução do
kernel
hierarquia de três níveis
. O
de todas as threads
executadas numa chamada de
blocos de threads
, que por sua vez contém as
3.4.1.
Kernel
Uma função escrita para rodar na GPU é chamada de
compilação __global__
é indispensável na declaração de um
será realizada a partir do programa principal rodan
__device__
torna a função invisível
kernel
chamado a partir da GPU.
O programa principal faz a chamada de um
qualquer funç
ão, exceto pela necessidade de incluir os parâmetro
configuração de paralelismo antes de passar as vari
݂ݑ̴݊ܿ݇݁ݎ݈݊݁ ൏൏൏ܽǡܾǡ
ܰݏ
das três
dimensões é especificado através de
Computação de Propósito Geral em Unidades de Processamento G
ráfico
Fluxo de processamento de um software empregando GPGPU
kernel
na GPU usando múltiplas threads
é baseada numa
. O
grid é o nível mais alto e é
constituído pelo conjunto
executadas numa chamada de
kernel. O grid
é subdividido em
, que por sua vez contém as
threads
no nível mais baixo.
Uma função escrita para rodar na GPU é chamada de
kernel
. A diretiva de
é indispensável na declaração de um
kernel
cuj
será realizada a partir do programa principal rodando na CPU.
A diretiva
torna a função invisível
para
a CPU e somente poderá ser usada
chamado a partir da GPU.
O programa principal faz a chamada de um
kernel
de maneira semelhante a
ão, exceto pela necessidade de incluir os parâmetro
configuração de paralelismo antes de passar as variáveis usando a notação
ܰݏ
ሺ
ܸͳǡܸʹǥ
ሻ
. O tamanho do grid
em cada uma
dimensões é especificado através de
ܽ, que
é uma estrutura com
42
GPU
.
é baseada numa
constituído pelo conjunto
é subdividido em
no nível mais baixo.
. A diretiva de
cuj
a chamada
A diretiva
a CPU e somente poderá ser usada
num
de maneira semelhante a
ão, exceto pela necessidade de incluir os parâmetros para
áveis usando a notação
em cada uma
é uma estrutura com
os
PUC-Rio - Certificação Digital Nº 0812691/CB
Computação de Propósito Geral em Unidades de Processamento Gráfico 43
campos ܽǤݔ, ܽǤݕ e ܽǤݖ. De forma semelhante, a configuração do tamanho dos
blocos é realizada através de ܾ. ܰݏ é reservado para os casos em que se deseja
especificar dinamicamente a quantidade de memória compartilhada a ser usada
em cada bloco.
3.4.2.
Grid
O nível mais alto da hierarquia de paralelismo de uma GPU é o grid, o qual
é composto por todas as threads executando numa chamada de kernel. O grid é
subdividido em blocos de threads, sendo limitado ao tamanho máximo de 65535
X 65535 X 1 blocos. O programa principal define explicitamente a quantidade de
blocos e a sua distribuição no grid através da passagem de parâmetros ao kernel.
A organização das threads em blocos e dos blocos no grid é ilustrada na Figura
12.
Figura 12 – Grid e Blocos de Threads (NVIDIA, 2009b).
PUC-Rio - Certificação Digital Nº 0812691/CB
Computação de Propósito Geral em Unidades de Processamento Gráfico 44
3.4.3.
Bloco de Threads
O segundo nível da hierarquia de threads é composto pela organização
destas em blocos. Cada bloco é limitado a um total de 512 threads que podem ser
distribuídas em três dimensões. De forma semelhante ao grid, o programa
principal define dinamicamente a configuração dos blocos de threads através da
passagem de parâmetros na chamada do kernel.
Os registradores de um multiprocessador são alocados para todas as threads
de um bloco simultaneamente e somente serão liberados quando o bloco terminar
de executar. Por isso é necessário um número relativamente grande de
registradores por multiprocessador para que cada thread tenha um número
razoável de registradores para sua execução. Um multiprocessador pode ativar e
executar até oito blocos de threads simultaneamente, desde que o limite de 1024
threads ativas por multiprocessador não seja ultrapassado.
Quando as threads de um bloco são sincronizadas, o poder computacional
do hardware da GPU pode ser empregado para executar operações em outro bloco
independente. Assim, o custo de sincronização das threads dentro de um bloco é
pequeno. A memória compartilhada pode ser acessada para leitura e escrita por
todas as threads de um mesmo bloco, permitindo a colaboração e a troca de dados
entre elas de forma rápida e eficiente. As threads de blocos diferentes só podem
ser sincronizadas com segurança no término da execução do kernel e qualquer
troca de dados só pode ser realizada através da memória global, a qual possui uma
latência elevada.
3.4.4.
Warp
A unidade de instruções do multiprocessador pode emitir uma instrução a
cada dois ciclos. Os oito processadores que compõem o multiprocessador rodam
no dobro da freqüência do resto da placa e podem iniciar uma instrução a cada
ciclo. Assim, cada instrução emitida pela unidade de instruções pode ser
executada por 32 threads, antes que a próxima instrução seja emitida, definindo o
tamanho do warp. Para obter-se o máximo desempenho possível é importante
PUC-Rio - Certificação Digital Nº 0812691/CB
Computação de Propósito Geral em Unidades de Processamento Gráfico 45
manter o número de threads de um bloco como sendo múltiplo do tamanho do
warp, ou seja, múltiplo de 32.
O multiprocessador emprega um modelo de arquitetura chamado SIMT
(Single Instruction, Multiple Thread). Esse modelo é muito semelhante ao modelo
SIMD. Quando um bloco de threads é atribuído para ser executado num
multiprocessador, o bloco é inicialmente dividido em warps. Todas as threads de
um warp são inicializadas juntas, contudo elas são livres para divergir e executar
independentemente. Assim, uma das principais diferenças do SIMT para o SIMD
é que o controle da execução e dos desvios é especificado para cada thread.
Contudo, para obter melhor desempenho é importante manter todas as threads de
um mesmo warp executando a mesma instrução, evitando desvios de fluxo dentro
do warp. Comandos condicionais que desviam o fluxo de todas as threads de um
warp em relação a outro warp podem ser empregados sem redução de
desempenho.
3.4.5.
Escalabilidade
O modelo de escalabilidade empregado facilita a adaptação do código
CUDA para a execução em GPUs diferentes. Quando o programa principal
rodando na CPU chama um kernel que executará na GPU, os blocos de um grid
são enumerados e distribuídos para os multiprocessadores com capacidade de
execução disponíveis. Todas as threads de um mesmo bloco executam
concorrentemente num único multiprocessador. Quando os blocos terminam,
novos são atribuídos aos multiprocessadores que ficaram vazios. Isto permite que
o mesmo código possa ser executado em duas GPUs, com números de
multiprocessadores diferentes, sem ter que ser compilado duas vezes. O código
automaticamente se adapta ao hardware disponível para extrair o máximo de
desempenho, conforme o ilustrado na Figura 13.
PUC-Rio - Certificação Digital Nº 0812691/CB
Computação de Propósito Geral em Unidades de Processamento Gráfico 46
Figura 13 – Modelo de Escalabilidade da GPU (NVIDIA, 2009b).
3.5.
Técnicas Otimizadas
Alguns cuidados devem ser levados em consideração ao realizarem-se
leituras e escritas nos diferentes tipos de memória para que seja obtido o melhor
desempenho possível com GPUs. As transferências de memória através do
barramento PCI Express constituem um dos possíveis gargalos. A seguir serão
descritos as formas mais eficientes de acesso às memórias global e compartilhada,
bem como técnicas para redução do tempo total despendido com transferências de
memória.
3.5.1.
Acesso à Memória Global
A memória global de uma GPU pode ser vista em termos de segmentos
alinhados de 16 e de 32 palavras de 32 bits. A Figura 14, adaptada de NVIDIA
(2009a), ilustra a divisão da memória da GPU. Cada linha representa um
segmento alinhado de 64 bytes, ou 16 números com precisão simples de ponto
flutuante. Duas linhas da mesma cor representam um segmento alinhado de 128
PUC-Rio - Certificação Digital Nº 0812691/CB
Computação de Propósito Geral em Unidades de Proces
bytes, ou 32 números com precisão simples
são representados as 16
threads
Figura 14 –
Segmentos Alinhados de Memória e
Uma GPU com capacidade de cálculo
pode acessar a memória global realizando a leitura
32 palavras de 32 bits numa única operação. A opera
todas as threads de um
half
será o menor, d
entro das opções 32, 64 ou 128 bytes
as threads de um
half-
warp
memória.
O melhor desempenho é obtido quando as
elementos de dados dentro de um mesmo segmento alin
todas as threads
precisam participar
nenhuma thread do
half-
warp
palavras. A Figura 15
ilustra este padrão de acesso à memória global.
Figura 15 – As threads
d
dentro de um mesmo segmento alinhado de 16 palavras
Para o padrão de acesso apresentado na
transferência de um bloco de dados de 64 bytes. Se
elementos acessados pelas
cálculo 1.1 terão que realizar 16 operações de tran
Computação de Propósito Geral em Unidades de Processamento G
ráfico
bytes, ou 32 números com precisão simples
de ponto flutuante. Na parte inferior
threads
de um half-warp.
Segmentos Alinhados de Memória e
threads de um half-
warp
Uma GPU com capacidade de cálculo
superior à 1.2, como a
Tesla C1060
pode acessar a memória global realizando a leitura de blocos de dados de 8, 16 ou
32 palavras de 32 bits numa única operação. A operação de leitura pode atender
half
-warp
em paralelo. O tamanho do bloco
de dados
entro das opções 32, 64 ou 128 bytes
, com capacidade
para atender
warp
. Assim, evita-se
desperdício de largura de banda
O melhor desempenho é obtido quando as
threads do
half-
warp
elementos de dados dentro de um mesmo segmento alinhado de 64 bytes. Nem
precisam participar
n
a operação de leitura, mas é importante que
warp
acesse elementos fora do segmento alinhado de 16
ilustra este padrão de acesso à memória global.
d
o half-warp
acessam elementos da memória global
dentro de um mesmo segmento alinhado de 16 palavras
(NVIDIA, 2009a)
Para o padrão de acesso apresentado na
Figura 15
, será realizada uma única
transferência de um bloco de dados de 64 bytes. Se houver permutação dos
elementos acessados pelas
threads do
half-warp
, dispositivos com capacidade de
cálculo 1.1 terão que realizar 16 operações de transferênci
a de dados. Porém, os
47
de ponto flutuante. Na parte inferior
warp
.
Tesla C1060
,
de blocos de dados de 8, 16 ou
ção de leitura pode atender
de dados
lido
para atender
desperdício de largura de banda
de
warp
acessam
hado de 64 bytes. Nem
a operação de leitura, mas é importante que
acesse elementos fora do segmento alinhado de 16
acessam elementos da memória global
(NVIDIA, 2009a)
.
, será realizada uma única
houver permutação dos
, dispositivos com capacidade de
a de dados. Porém, os
PUC-Rio - Certificação Digital Nº 0812691/CB
Computação de Propósito Geral em Unidades de Processamento Gráfico 48
dispositivos com capacidade de cálculo 1.2 ou superior não terão seu desempenho
afetado, continuando a realizar uma única operação de transferência de dados.
As threads do
half-warp
podem acessar 16 palavras de 32 bits que não
estejam alinhadas com um bloco de 64 bytes da memória global, conforme o
ilustrado na Figura 16. Isto gerará 16 transferências de dados para dispositivos
com capacidade de cálculo 1.1. Para dispositivos com capacidade de cálculo 1.2
ou superior, poderá ser realizada uma única operação de transferência de dados de
128 bytes. Ou seja, 32 palavras são lidas da memória e somente metade delas será
efetivamente usada. Contudo, o desempenho é melhor do que realizar 16
operações de leitura de blocos de 8 palavras, onde somente uma palavra por bloco
será aproveitada.
Figura 16 – Elementos não alinhados com o bloco de 16 palavras residem dentro
do mesmo bloco de 128 bytes (NVIDIA, 2009a).
Os elementos não alinhados com o bloco de 16 palavras podem estar dentro
de dois blocos diferentes de 128 bytes, conforme o ilustrado na Figura 17. Neste
caso, os dispositivos com capacidade de cálculo 1.2 ou superior terão que realizar
2 operações de transferência de dados. Uma operação fará a leitura de 16 palavras
das quais serão utilizadas 15 e a outra operação fará a leitura de 8 palavras, das
quais serão utilizadas apenas uma.
Figura 17 – Elementos não alinhados com o bloco de 16 palavras residem dentro
de dois blocos diferentes de 128 bytes (NVIDIA, 2009a).
A Figura 18 apresenta a influencia das diversas formas de acesso à memória
sobre o desempenho da largura de banda efetiva. São apresentados os casos para
dois dispositivos, a NVIDIA GeForce GTX 280 com capacidade de cálculo 1.3 e
PUC-Rio - Certificação Digital Nº 0812691/CB
Computação de Propósito Geral em Unidades de Processamento Gráfico 49
o NVIDIA Quadro® FX 5600 com capacidade de cálculo 1.0. Deslocamento
unitário é o tipo de acesso ilustrado na Figura 16.
Figura 18 – Largura de banda efetiva para padrões de acesso à memória com
deslocamento (NVIDIA, 2009a).
Pelo ilustrado na Figura 18, podemos verificar que, dispositivos com
capacidade de cálculo 1.1 ou inferiores, tem seu desempenho de largura de banda
drasticamente afetado por padrões de acesso com deslocamento diferente de zero
ou 16. O desempenho da largura de banda para este caso ficou reduzido em cerca
de oito vezes. Por outro lado, dispositivos com capacidade de cálculo 1.2 ou
superior sofrem uma menor redução de desempenho. A redução na largura de
banda para deslocamentos diferentes de zero ou 16 foi de apenas duas vezes.
As threads do
half-warp
podem acessar 16 elementos da memória deixando
intervalos entre os elementos acessados conforme o ilustrado na Figura 19. Neste
caso específico, todos os elementos lidos residem dentro do mesmo bloco de 128
bytes alinhados e um dispositivo com capacidade de cálculo 1.3 ou superior
poderia realizar uma única operação de transferência de dados. Contudo, 32
palavras são lidas, sendo que somente 16 são efetivamente utilizadas.
Figura 19 – Acesso a elementos da memória com intervalos entre os elementos
acessados (NVIDIA, 2009a).
PUC-Rio - Certificação Digital Nº 0812691/CB
Computação de Propósito Geral em Unidades de Processamento Gráfico 50
O padrão de acesso ilustrado na Figura 19 deve ser evitado. À medida que o
número de intervalos entre os elementos lidos aumenta, a largura de banda
efetivamente utilizada é drasticamente reduzida, conforme o ilustrado na Figura
20.
Figura 20 – Desempenho da largura de banda para padrões de acesso com
intervalos entre os elementos lidos (NVIDIA, 2009a).
Um dos fatores responsáveis pela redução da largura de banda em acessos à
memória com deslocamentos é o fato de que a memória global da GPU não possui
cache. A memória de textura, apesar de ser somente leitura, possui cache. Assim,
para o caso de variáveis que serão somente lidas, a memória de textura poderá ser
uma alternativa para evitar a redução de desempenho em acessos com
deslocamento diferente de zero ou 16. A Figura 21 compara o desempenho da
largura de banda em leituras, com deslocamentos, para a memória global e para a
memória de textura.
Figura 21 – Leitura, com deslocamentos, na memória global e na memória de
textura (NVIDIA, 2009a).
PUC-Rio - Certificação Digital Nº 0812691/CB
Computação de Propósito Geral em Unidades de Processamento Gráfico 51
3.5.2.
Acesso à Memória Compartilhada
A memória compartilhada é interna ao chip e por isso é muito mais rápida
do que a memória global. O acesso a memória compartilhada pode ser tão rápido
quanto o acesso aos registradores, desde que não existam conflitos de banco entre
as threads pertencentes a um
half-warp
.
A memória compartilhada é dividida em 16 módulos exatamente iguais, os
quais são chamados de bancos. Uma operação de transferência de dados pode
atender as 16 threads de um
half-warp
em paralelo. Para isso, é necessário que
cada thread acesse um banco diferente, não importando à ordem dos acessos. Se
duas threads acessarem um mesmo banco, então teremos um conflito de banco e
serão necessárias duas operações de transferência para atender todas as threads de
um
half-warp
. Assim, o tempo de acesso aumenta devido aos conflitos de banco.
Na Figura 22 são apresentados padrões de acesso à memória compartilhada
sem conflitos de banco. No lado esquerdo temos um padrão de acesso linear e no
lado direito temos um padrão de acesso com permutações randômicas entre as
threads e os bancos. Em nenhum dos casos ocorre de duas threads usarem o
mesmo banco e, portanto, não há conflitos de bancos, podendo todas as threads do
half-warp
serem servidas numa única operação de transferência.
PUC-Rio - Certificação Digital Nº 0812691/CB
Computação de Propósito Geral em Unidades de Processamento Gráfico 52
Figura 22 – Padrões de acesso à memória compartilhada sem conflitos de banco
(NVIDIA, 2009b).
Padrões de acesso a memória compartilhada com conflitos de banco são
ilustrados na Figura 23. No padrão da esquerda as 16 threads de um
half-warp
acessam 16 palavras de 32 bits com um intervalo entre cada elemento de dado.
Dessa forma, todos os acessos recaem sobre somente oito bancos, gerando
conflitos que necessitarão de duas operações de leitura. No padrão de acesso da
direita, existe um intervalo de sete elementos vazios entre cada elemento
acessado. Desta forma as leituras são concentradas em apenas dois bancos e serão
necessárias oito operações de transferência de dados para atender todas as threads.
PUC-Rio - Certificação Digital Nº 0812691/CB
Computação de Propósito Geral em Unidades de Processamento Gráfico 53
Figura 23 – Padrões de acesso à memória compartilhada com conflitos de banco
(NVIDIA, 2009b).
3.5.3.
Transferência de Dados Através do Barramento PCI Express
A maneira tradicional de realização de cálculos usando uma GPU consiste
em transferir um vetor de dados da memória da CPU para a memória da GPU
através do barramento PCI Express. Quando a transferência termina, um kernel é
invocado para realizar operações sobre estes dados. Ao final da execução do
kernel, os dados são novamente copiados da CPU para a GPU.
Contudo, em alguns casos, se as dependências entre os dados permitirem,
poderá ser empregado um modelo onde os dados são copiados assincronamente.
Assim, não é necessário esperar que todo o vetor de dados tenha sido transferido
para a memória da GPU para que os cálculos comessem a ser realizados. Neste
PUC-Rio - Certificação Digital Nº 0812691/CB
Computação de Propósito Geral em Unidades de Proces
caso são empregados múltiplos fluxos de processamen
chegam à
memória da GPU eles são processados e podem ser cop
para CPU de forma
assíncrona com a tra
dados (NVIDIA, 2009a)
. A
possibilidade de redução do tempo total pela sobrep
transferência e do tempo de execução.
Figura 24 – Redução
do tempo total através da sobreposição dos tempos d
transferência e de execução
3.5.4.
Nível de Ocupação do Multiprocessador
Um multiprocessador de um dispositivo com capacidad
superior pode manter até 32
instrução executada por um Multiprocessador sobre u
apresenta um grande período de latência, os re
aproveitados para executar outras instruções de out
possuem maior latência são os acessos de leitura e
O Nível de Ocupação do Multiprocessador é a razão e
warps
ativos e o número máximo possível de
Ocupação significa que o Multiprocessador não terá
instruções em outros
warps
Desta forma, o desempenho
cálculo 1.2 ou superior recomenda
18.75
%. Contudo, melhores resultados podem ser obtidos c
elevados de Ocupação
(NVIDIA, 200
Por outro lado, um elevado Nível de Ocupação não ne
traduz num elevado nível de desempenho. Por exemplo
Ocupação de 66 % para 100 % não significa que o des
Computação de Propósito Geral em Unidades de Processamento G
ráfico
caso são empregados múltiplos fluxos de processamento e à medida que os dados
memória da GPU eles são processados e podem ser copiados de volta
assíncrona com a tra
nsferência de outras partes do vetor de
. A
Figura 24
ilustra os dois casos, mostrando a
possibilidade de redução do tempo total pela sobreposição dos tempos de
transferência e do tempo de execução.
do tempo total através da sobreposição dos tempos d
transferência e de execução
(NVIDIA, 2009a).
Nível de Ocupação do Multiprocessador
Um multiprocessador de um dispositivo com capacidade de cálculo 1.2 ou
superior pode manter até 32
warps
ativos simultaneamente. Quando uma
instrução executada por um Multiprocessador sobre um determinado
apresenta um grande período de latência, os re
cursos de hardware podem ser
aproveitados para executar outras instruções de outros
warps.
As instruções que
possuem maior latência são os acessos de leitura e escrita na memória global.
O Nível de Ocupação do Multiprocessador é a razão entre o número de
ativos e o número máximo possível de
warps
ativos. Um baixo Nível de
Ocupação significa que o Multiprocessador não terá a possibilidade de executar
warps
, quando ocorrer uma instrução com grande latência.
Desta forma, o desempenho
é reduzido. Para dispositivos com capacidade de
cálculo 1.2 ou superior recomenda
-se manter um Ní
vel mínimo de Ocupação de
%. Contudo, melhores resultados podem ser obtidos com níveis mais
(NVIDIA, 200
9a).
Por outro lado, um elevado Nível de Ocupação não necessariamente se
traduz num elevado nível de desempenho. Por exemplo, aumentar o Nível de
Ocupação de 66 % para 100 % não significa que o desempenho vá aumentar na
54
to e à medida que os dados
iados de volta
nsferência de outras partes do vetor de
ilustra os dois casos, mostrando a
osição dos tempos de
do tempo total através da sobreposição dos tempos de
e de cálculo 1.2 ou
ativos simultaneamente. Quando uma
m determinado
warp
cursos de hardware podem ser
As instruções que
escrita na memória global.
ntre o número de
ativos. Um baixo Nível de
a possibilidade de executar
, quando ocorrer uma instrução com grande latência.
é reduzido. Para dispositivos com capacidade de
vel mínimo de Ocupação de
om níveis mais
cessariamente se
, aumentar o Nível de
empenho vá aumentar na
PUC-Rio - Certificação Digital Nº 0812691/CB
Computação de Propósito Geral em Unidades de Processamento Gráfico 55
mesma proporção. Na prática, uma vez que o nível de 50 % tenha sido atingido,
continuar aumentado o Nível de Ocupação não trará benefícios adicionais
(NVIDIA, 2009a).
PUC-Rio - Certificação Digital Nº 0812691/CB
4
Descrição das Alterações para Aceleração do SIESTA por
GPU
Inicialmente avaliou-se a possibilidade de trabalhar com os principais
algoritmos de código fonte aberto empregados na química quântica. Entre eles
esta o MPQC (Massively Parallel Quantum Chemistry Program), que emprega
métodos tais como o Hartree-Fock e a Teoria do Funcional da Densidade
(Janssen, Nielsen, et al., 2010). Outro pacote considerado foi o Quantum
ESPRESSO (opEn-Source Package for Research in Electronic Structure,
Simulation, and Optimization), que possui cálculos de estrutura eletrônica
realizados com a Teoria do Funcional da Densidade, empregando ondas planas e
pseudopotenciais (DEMOCRITOS, 2009). O YAeHMOP (Yet Another extended
Huckel Molecular Orbital Package), que é um método tight-binding usado para
obter uma rápida descrição qualitativa (Landrum e Glassey, 2001). O escolhido
para utilização no trabalho foi o SIESTA (Spanish Initiative for Electronic
Simulations with Thousands of Atoms), o qual usa DFT, pseudopotenciais, orbitais
atômicos numéricos e funções de base localizada para aumentar a sua eficiência
(Soler, Artacho, et al., 2002). Este programa empregado na química quântica
possui 281 arquivos, totalizando 142023 linhas de código fonte em Fortran.
Devido a sua complexidade, apenas partes com maior custo computacional serão
portadas para execução em GPU.
Conforme o descrito no item 2.7, o SIESTA possui uma maneira particular
de representar os termos da energia total. Uma densidade atômica é definida para
reescrever os termos do pseudopotencial local e do potencial de Hartree. Isto é
feito com o objetivo de eliminar os termos de longo alcance, aumentando assim a
eficiência computacional. As funções de bases atômicas localizadas, usadas no
SIESTA foram escolhidas para viabilizar o desenvolvimento de um método ab
initio DFT autoconsistente de ordem N. Pois, o emprego de ondas planas como
funções de base impossibilitaria a determinação do hamiltoniano autoconsistente
PUC-Rio - Certificação Digital Nº 0812691/CB
Descrição das Alterações para Aceleração do SIESTA por GPU 57
em O(N) operações. Por outro lado, funções tais como a Transformada de Fourier
ainda necessitam de ordem N log(N) (Soler, Artacho, et al., 2002).
Neste Capítulo, as partes do código SIESTA com potencial para execução
em GPU são listadas e são apresentadas algumas alterações realizadas para testes
e comprovação de desempenho.
4.1.
Identificação das Partes Adequadas ao Paralelismo de Dados
O termo da eq. (2.42), visto no Capitulo 2, relacionado às contribuições para
a energia total da parte não-local do pseudopotencial e o termo de cálculo da
energia cinética dos elétrons, são operações repetitivas sobre dados diferentes e
poderiam ser estudados para execução acelerada por GPU. Uma matriz de
densidade é usada para as demais integrais da eq. (2.42), como a energia de troca
e correlação, o potencial e a energia de Hartree e a parte local do pseudopotencial
do hamiltoniano do SIESTA.
Inicialmente foi estuda a melhor forma de se reescrever em CUDA, os
termos referentes ao cálculo do potencial de Hartree e da energia de Hartree, com
a resolução da Equação de Poisson na GPU, e a contribuição para as forças de
deslocamento atômicas, devidas ao potencial de Hartree.
A próxima etapa será o estudo do termo da energia de troca e correlação e a
parte local do pseudopotencial do hamiltoniano do SIESTA. Isto permitiria
trabalhar a matriz de densidade na memória da GPU. Somente seriam necessárias
transferências de memória para os parâmetros de inicialização e os resultados das
integrais. É importante reduzir as transferências de memória entre CPU e GPU
devido à latência da interface PCI Express.
4.1.1.
CUFFT
A biblioteca CUFFT (CUDA Fast Fourier Transform) possui funções para a
realização de Transformadas de Fourier, no sentido direto ou no sentido inverso,
sobre matrizes de dados reais ou complexas com uma, duas ou três dimensões
(NVIDIA, 2007b). Para emprego no SIESTA, são adequadas as funções
apresentadas na Tabela 2.
PUC-Rio - Certificação Digital Nº 0812691/CB
Descrição das Alterações para Aceleração do SIESTA por GPU 58
Tabela 2 – Funções da Biblioteca CUFFT adequadas para emprego no SIESTA.
1 cufftExecR2C Transformada de Fourier no sentido direto para matrizes tridimensionais.
2 cufftExecC2R Transformada Inversa de Fourier para matrizes tridimensionais.
4.1.2.
CUBLAS
A biblioteca CUBLAS é uma implementação das funções básicas de álgebra
linear BLAS (Basic Linear Algebra Subprograms) para execução paralela, em
GPUs, usando CUDA (NVIDIA, 2007a). A Tabela 3 lista e descreve brevemente
algumas das funções CUBLAS que poderiam substituir funções de álgebra linear
no código SIESTA.
Tabela 3 – Funções da Biblioteca CUBLAS com possibilidade de emprego para o SIESTA.
1 cublasDgemm
calcula o produto de duas matrizes A e B, multiplicando o resultado por um
escalar alpha. A seguir soma com o produto de uma matriz C por um escalar
beta.
2 cublasIsamax
encontra o máximo de um vetor. Se o máximo não for um ponto único, o
ponto retornado é o que possui o menor índice.
3 cublasSasum calcula a soma dos valores absolutos de um vetor.
4 cublasSaxpy multiplica um vetor por um escalar.
5 cublasScopy copia os elementos de um vetor x para um vetor y.
6 cublasSdot calcula o produto escalar de dois vetores.
7 cublasSscal multiplica um vetor x por um escalar alpha.
8 cublasCaxpy
multiplica um vetor de números complexos x por um escalar alpha e em
seguida adiciona um segundo vetor de números complexos y.
9 cublasCscal multiplica um vetor de números complexos x por um escalar alpha.
10
cublasCswap
troca os elementos de um vetor de números complexos x, com os elementos
de outro vetor de números complexos y.
4.1.3.
MAGMA
A biblioteca de álgebra linear MAGMA (Matrix Algebra on GPU and
Multicore Architectures) foi projetada para ser similar ao LAPACK (Linear
Algebra PACKage) em nível de funções suportadas (Tomov, Nath, et al., 2009).
Algumas das suas funções, listadas na Tabela 4, apresentam compatibilidade com
funções empregadas no SIESTA.
Tabela 4 – Funções da Biblioteca MAGMA com possibilidade de emprego para o SIESTA.
1
magma_dpotrf_gpu calcula a fatoração de Cholesky de uma matriz real simétrica positiva.
2
magma_zpotrf_gpu
calcula a fatoração de Cholesky de uma matriz Hermitiana complexa
positiva.
3
magmablas_dtrsm
resolve equações matriciais em GPU, com matrizes X e B, m por n,
triangulares superior ou inferior e sendo possível aplicar o operador de
transposição à matriz A.
ሺ
ܣ
ሻ
כ
ܺ
ൌ
݈ܽ
݄
ܽ
כ
ܤ
ܺ
כ
ሺ
ܣ
ሻ
ൌ
݈ܽ
݄
ܽ
כ
ܤ
PUC-Rio - Certificação Digital Nº 0812691/CB
Descrição das Alterações
para
A biblioteca MAGMA
meio científico, tais como
lineares baseados nas decomposições LU, QR e
2009).
4.1.4.
CULA
Outra biblioteca
de álgebra linear, acelerada por GPUs,
meio científico é a CULA
(EM Photonics, 2010)
disponibilizadas
, para números reais e para números complexos,
Tabela 5.
Tabela 5 –
Funções da Biblioteca CULA
1
SGESV e CGESV
resolve um sistema geral de equações
2
SGETRF e CGETRF
calcula a fatoração LU de uma matriz geral.
3
SGEQRF e CGEQRF
calcula a fatoração QR de uma matriz retangular ger
4
SGELS e CGELS
resolve um sistema de equações lineares por mínimos
5
SGGLSE e CGGLSE
resolve um sistema linear com restrições por mínimo
6
SGESVD e CGESVD
calcula a decomposição em valor singular de uma mat
geral.
As funções de resolução de sistemas lineares das bi
CULA poderiam ser empregadas para
complexas, usadas no código
4.2.
Alteração de Partes do SIESTA
A última versão estável do código fonte do SIESTA,
em Fortran, foi utilizado como
código para GPU.
O paralelismo de dados da
linguagem CUDA
, e a linguagem C é usada como uma interface entre
CUDA
, conforme o ilustrado na
Linux CentOS 5.5.
Figura 25 –
Interface entre
para
Aceleração do SIESTA por GPU
A biblioteca MAGMA
também possui outras funções de grande utili
um conjunto de funções CUDA para resolver sistemas
lineares baseados nas decomposições LU, QR e
Cholesky
(Tomov, Nath,
de álgebra linear, acelerada por GPUs,
com aplicação no
(EM Photonics, 2010)
. Algumas das principais funções
, para números reais e para números complexos,
são
listadas na
Funções da Biblioteca CULA
utilizadas no meio científico.
resolve um sistema geral de equações
lineares AX=B.
calcula a fatoração LU de uma matriz geral.
calcula a fatoração QR de uma matriz retangular geral.
resolve um sistema de equações lineares por mínimos quadrados.
resolve um sistema linear com restrições por mínimos quadrados.
calcula a decomposição em valor singular de uma matriz retangular
geral.
As funções de resolução de sistemas lineares das bibliotecas
MAGMA
CULA poderiam ser empregadas para
calcular
a inversa de matrizes reais ou
complexas, usadas no código
SIESTA.
Alteração de Partes do SIESTA
A última versão estável do código fonte do SIESTA, originalmente
em Fortran, foi utilizado como
ponto de partida para estudar
e alterar
O paralelismo de dados da
s
GPUs é acessado através da
, e a linguagem C é usada como uma interface entre
, conforme o ilustrado na
Figura 25.
O sistema operacional utilizado foi o
Interface entre
Fortran e CUDA.
59
também possui outras funções de grande utili
dade no
um conjunto de funções CUDA para resolver sistemas
(Tomov, Nath,
et al.,
com aplicação no
. Algumas das principais funções
listadas na
quadrados.
s quadrados.
riz retangular
MAGMA
e
a inversa de matrizes reais ou
originalmente
escrito
e alterar
trechos de
GPUs é acessado através da
, e a linguagem C é usada como uma interface entre Fortran e
O sistema operacional utilizado foi o
PUC-Rio - Certificação Digital Nº 0812691/CB
Descrição das Alterações
para
4.2.1.
Potencial de Hartree,
Energia
Atômicas
Inicialmente, o termo referente
reescrito
usando a linguagem CUDA, empregando a resolução da
Poisson na GPU.
O quarto termo do lado direito da eq.
Hartree, ߜܸ
ு
ሺ
ݎ
ሻ
,
é produzido
consistente ߜߩ
ሺ
ݎ
ሻ
. A solução
implementada no SIESTA permite
com a variação da densidade
A
Equação de Poisson
eletrostática, na engenharia e na física teórica. N
é freqüentemente escrita
como:
Considerando o sistema de coordenadas cartesianas p
equação pode ser reescrita como:
ቆ
߲
ଶ
߲ݔ
ଶ
߲ݕ
Para problemas que envolvem
Transformada de Fourier
equações algébricas
, as quais são resolvidas com
Em seguida a transformada inversa é utilizada para
para o espaço real. A
Figura
de equações diferencias.
Figura 26 –
As equações diferenciais são transformadas em equaç
A
energia de Hartree
atômicas, devido ao potencial de Hartree, estão
código da
Equação de Poisson
elemento da matriz de densidades
para
Aceleração do SIESTA por GPU
Energia
de Hartree e Forças de De
slocamento
Inicialmente, o termo referente
ao cálculo do potencial de Hartre
usando a linguagem CUDA, empregando a resolução da Equação de
O quarto termo do lado direito da eq.
(2.42),
potencial de
é produzido
pela diferença d
as densidades atômica
. A solução
da equação de Poisson
na forma original
implementada no SIESTA permite
encontrar o potencial de Hartree
com a variação da densidade
.
Equação de Poisson
é uma equação diferencial parcial muito utilizada n
eletrostática, na engenharia e na física teórica. No espaço euclidiano esta equação
como:
ଶ
ɔ
ൌ
Considerando o sistema de coordenadas cartesianas por simplicidade, esta
equação pode ser reescrita como:
߲
ଶ
߲ݕ
ଶ
߲
ଶ
߲ݖ
ଶ
ቇ
ɔ
ሺ
ǡ
ǡ
ሻ
ൌ
ሺ
ǡ
ǡ
ሻ
Para problemas que envolvem
condições periódicas de contorno
Transformada de Fourier
permite transformar as
equações diferenciais em
, as quais são resolvidas com
operações algébricas normais
Em seguida a transformada inversa é utilizada para trazer a solução novamente
Figura
26
ilustra o emprego deste processo para a resolução
As equações diferenciais são transformadas em equações algébricas
energia de Hartree
e a contribuição para as forças de deslocamento
atômicas, devido ao potencial de Hartree, estão
combinadas
na mesma seção do
Equação de Poisson
e também são computadas
em GPU.
elemento da matriz de densidades
ߜߩ
ሺ
ݎ
ሻ
contribui individualmente para a
60
slocamento
ao cálculo do potencial de Hartre
e foi
Equação de
potencial de
as densidades atômica
s e auto-
na forma original
associado
é uma equação diferencial parcial muito utilizada na
o espaço euclidiano esta equação
(4.1)
or simplicidade, esta
(4.2)
condições periódicas de contorno
, o uso da
equações diferenciais em
operações algébricas normais
.
trazer a solução novamente
ilustra o emprego deste processo para a resolução
ões algébricas
.
e a contribuição para as forças de deslocamento
na mesma seção do
em GPU.
Cada
contribui individualmente para a
energia
PUC-Rio - Certificação Digital Nº 0812691/CB
Descrição das Alterações
para
de Hartree
e para as forças de deslocamento atômicas. Assim, o
de Hartree
e das forças de deslocamento atômicas em GPU envolv
variável calculada, somatório
O somatório
paralelo realizado na memória compartilhada de uma
consiste na redução do vetor de dados pela metade e
elementos do hardware paralelo. A
paralelo em GPU
para um único veto de dados
Figura 27 –
Somatório paralelo realizado na memória compartilha
um único vetor de dados.
4.2.2.
Implementação Paralela para GPUs
A
Equação de Poisson
periódicas de contorno e
é resolvida no domínio da freqüência. Inicialmente
encontrados os componentes da Transformada de Fouri
ߜߩ
ሺ
ݎ
ሻ
. A solução é realizada no
contribuição de cada elemento da matriz de densidad
Hartree
e as forças de deslocamento atômicas.
em seguida transformado para
A
Transformada de Fourier
o algoritmo GPFA (
Generalized Prime Factor FFT
(1992)
. Este algoritmo era originalmente empregado
densidades ߜߩ
ሺ
ݎ
ሻ
do espaço real para o domínio da freqüência e para
para
Aceleração do SIESTA por GPU
e para as forças de deslocamento atômicas. Assim, o cálculo da
e das forças de deslocamento atômicas em GPU envolve, para cada
variável calculada, somatório
s sobre todos os elementos d
a matriz de densidades.
paralelo realizado na memória compartilhada de uma
consiste na redução do vetor de dados pela metade em cada operação dos
elementos do hardware paralelo. A
Figura 27
ilustra a realização de um somatório
para um único veto de dados
.
Somatório paralelo realizado na memória compartilhada da GPU
Implementação Paralela para GPUs
Equação de Poisson
usada no SIESTA
está baseada em condições
é resolvida no domínio da freqüência. Inicialmente
encontrados os componentes da Transformada de Fourier da densidade eletrônica
. A solução é realizada no
domínio da freqüência
e um somatório
contribuição de cada elemento da matriz de densidades produz a
energia de
e as forças de deslocamento atômicas.
O potencial de Hartree
,
em seguida transformado para
o espaço real conforme o ilustrado na
Figura
Transformada de Fourier
originalmente empregada
no código SIESTA
Generalized Prime Factor FFT
)
desenvolvido por Temperton
. Este algoritmo era originalmente empregado
para transformar a matriz de
do espaço real para o domínio da freqüência e para transformar
61
cálculo da
energia
e, para cada
a matriz de densidades.
paralelo realizado na memória compartilhada de uma GPU
m cada operação dos
ilustra a realização de um somatório
da da GPU
para
está baseada em condições
é resolvida no domínio da freqüência. Inicialmente são
er da densidade eletrônica
e um somatório
da
energia de
,
ߜܸ
ு
ሺ
ݎ
ሻ
, é
Figura
28.
no código SIESTA
é
desenvolvido por Temperton
para transformar a matriz de
do espaço real para o domínio da freqüência e para transformar
PUC-Rio - Certificação Digital Nº 0812691/CB
Descrição das Alterações
para
o potencial de Hartree,
ߜܸ
espaço real.
Na
Equação de Poisson
GPFA pelo algoritmo
CUFFT
NVIDIA (2007b).
Assim, a Transformada de Fourier da matriz de densi
ߜߩ
ሺ
ݎ
ሻ
e a Transformada Inversa de Fourier do potencial de
realizadas na placa gráfica.
de Fourier poderiam contribuir p
código SIESTA
, tal como o algoritmo proposto por
Para o cálculo do potencial de Hartre
escrito um kernel baseado
etapa consistiu n
a tradução do algoritmo da linguagem Fortran para a
C. Em seguida estudou-
se a melhor forma de tirar proveito do paralelismo
dados oferecido pela GPU. As funções de transfer
GPU foram introduzidas e os laços originais foram s
de hardware da GPU. Foram realizados testes para av
otimizações de acesso à memória global, de acesso à
transferência através do barramento
3, nos itens 3.5.1, 3.5.2 e
3.5.3
Figura 28 – So
lução no domínio da freqüência
A CPU tem uma unidade de instruç
a utilização de 240
thread
poderia ser adequada.
Nas GPUs o paradigma é diferente, pois a GPU possui
unidade de instruções para cada 8 elementos de proc
multiprocessador),
o tempo necessário
aplicá-la sobre 4 threads
instrução carregada na unidade de instruções
para
Aceleração do SIESTA por GPU
ߜܸ
ு
ሺ
ݎ
ሻ
, c
alculado no domínio da freqüência, para o
Equação de Poisson
para execução em GPU foi substituído
o algoritmo
CUFFT
(CUDA Fast Fourier Transform)
desenvolvido pela
Assim, a Transformada de Fourier da matriz de densi
e a Transformada Inversa de Fourier do potencial de
Hartree,
ߜܸ
realizadas na placa gráfica.
Algoritmos GPU mais eficientes para a Transformada
de Fourier poderiam contribuir p
ara aumentar o desempenho desta etapa do
, tal como o algoritmo proposto por
Nukada, Ogata,
et al.
Para o cálculo do potencial de Hartre
e, ߜܸ
ு
ሺ
ݎ
ሻ
, propriamente dito, foi
na função seqüencial original do SIESTA. A primeira
a tradução do algoritmo da linguagem Fortran para a linguagem
se a melhor forma de tirar proveito do paralelismo
dados oferecido pela GPU. As funções de transfer
ência de dados entre CPU e
GPU foram introduzidas e os laços originais foram substituídos pelo par
de hardware da GPU. Foram realizados testes para avaliar os efeitos das
otimizações de acesso à memória global, de acesso à memória compartilhada e d
transferência através do barramento
PCI Express
, conforme o descrito no Capítulo
3.5.3
, respectivamente.
lução no domínio da freqüência
.
A CPU tem uma unidade de instruç
ões para cada
núcleo de processamento e
thread
s num cluster
com 240 núcleos de processamento
Nas GPUs o paradigma é diferente, pois a GPU possui
somente
unidade de instruções para cada 8 elementos de processamento
o tempo necessário
para mudar de in
strução é suficiente para
de cada elemento de processamento.
Ou seja,
instrução carregada na unidade de instruções
pode ser aplicada em 32
thread
62
alculado no domínio da freqüência, para o
o algoritmo
desenvolvido pela
Assim, a Transformada de Fourier da matriz de densidades
ߜܸ
ு
ሺ
ݎ
ሻ
, são
Algoritmos GPU mais eficientes para a Transformada
ara aumentar o desempenho desta etapa do
et al.
(2008).
, propriamente dito, foi
na função seqüencial original do SIESTA. A primeira
a tradução do algoritmo da linguagem Fortran para a linguagem
se a melhor forma de tirar proveito do paralelismo de
ência de dados entre CPU e
ubstituídos pelo par
alelismo
aliar os efeitos das
memória compartilhada e d
e
, conforme o descrito no Capítulo
núcleo de processamento e
com 240 núcleos de processamento
somente
uma
essamento
(um
strução é suficiente para
Ou seja,
cada
thread
s, um
PUC-Rio - Certificação Digital Nº 0812691/CB
Descrição das Alterações
para
warp, executand
o em 8 elementos de processamento.
foram
ser projetadas para serem
Além disso, existem os tempos de latência no acesso
memória compartilhada
e aos registradores
a possibilidade de executar instruções em outros
instrução com grande latência. Assim
dado
ser lido ou escrito nas memórias
mínimo de 192 threads
por bloco para que o multiprocessador possa ter sem
instruções disponíveis
, reduzindo
2009b)
. Desta forma, um número mínimo de
192 threads / bloco x 30
desempenho será obtido se o problema for pensado de
número maior de threads
na GPU.
Na primeira versão do
SIESTA na GPU
foi atribuída uma
densidades ߜߩ
ሺ
ݎ
ሻ
(Figura
29
registradores e 3616 bytes
Multiprocessador
é de 75 %, conforme pode ser observado na
Figura 29 – Uma thread
para cada elemento da matriz de densidades
para
Aceleração do SIESTA por GPU
o em 8 elementos de processamento.
Assim, as
thread
s
ser projetadas para serem
muito leves.
Além disso, existem os tempos de latência no acesso à memória global
e aos registradores
. Com 32 threads por bloco
a possibilidade de executar instruções em outros
warps
, quando ocorre uma
instrução com grande latência. Assim
,
o multiprocessador fica ocioso esperando o
ser lido ou escrito nas memórias
. Por isso, a NVIDIA r
ecomenda manter um
por bloco para que o multiprocessador possa ter sem
, reduzindo
o tempo ocioso do multiprocessador
. Desta forma, um número mínimo de
threads
deveria ficar em torno de
multiprocessadores = 5760 threads
. Contudo, melhor
desempenho será obtido se o problema for pensado de forma
a dividi
-
na GPU.
Na primeira versão do
kernel para a resolução da
Equação de Poisson
foi atribuída uma
thread
para cada elemento da matriz de
29
). Neste caso, cada bloco com 256
threads
de memória compartilhada. O
Nível de Ocupação do
é de 75 %, conforme pode ser observado na
Figura
30
para cada elemento da matriz de densidades
.
63
thread
s
de GPUs
à memória global
, à
não existe
, quando ocorre uma
o multiprocessador fica ocioso esperando o
ecomenda manter um
por bloco para que o multiprocessador possa ter sempre
(NVIDIA,
deveria ficar em torno de
. Contudo, melhor
-
lo por um
Equação de Poisson
do
para cada elemento da matriz de
threads
utiliza 18
Nível de Ocupação do
30
.
PUC-Rio - Certificação Digital Nº 0812691/CB
Descrição das Alterações
para
Figura 30 –
Nível de Ocupação do Multiprocessador
kernel
da Equação de Poisson
Foi desenvolvida u
ma segunda versão d
de Poisson
na GPU, onde cada
densidades ߜߩ
ሺ
ݎ
ሻ
(
Figura
registradores e 7200 bytes
de memória compartilhada. A intensidade aritmética
cada thread
é maior, por outro lado o
reduzido para 50 %
, conforme o apresentado na
Figura 31 – Uma thread
processa um vetor de dados da matriz de densidades
para
Aceleração do SIESTA por GPU
Nível de Ocupação do Multiprocessador
para
primeira versão do
da Equação de Poisson
.
ma segunda versão d
o kernel
para resolução da Equação
na GPU, onde cada
thread processa um vetor de dados
da matriz de
Figura
31). Cada bloco com 256 threads
de memória compartilhada. A intensidade aritmética
é maior, por outro lado o
Nível de Ocupação do Multiprocessador
, conforme o apresentado na
Figura 32.
processa um vetor de dados da matriz de densidades
64
primeira versão do
para resolução da Equação
da matriz de
utiliza 28
de memória compartilhada. A intensidade aritmética de
Nível de Ocupação do Multiprocessador
é
processa um vetor de dados da matriz de densidades
.
PUC-Rio - Certificação Digital Nº 0812691/CB
Descrição das Alterações
para
Figura 32 –
Nível de Ocupação do Multiprocessador
kernel
da Equação de Poisson
O kernel usado
para resolução da Equação de Poisson
das contribuições dos elementos da matriz de densid
da energia de Hartree
(um somatório) e das forças de desloc
somatórios). Os pré-
somatórios são realizados para cada bloco. Assim, p
obterem-
se os valores finais, foi utilizado um segundo
registradores e 3620
bytes
threads. O
Nível de Ocupação do Multiprocessador
ilustrado na Figura 33.
para
Aceleração do SIESTA por GPU
Nível de Ocupação do Multiprocessador
para segunda
da Equação de Poisson
.
para resolução da Equação de Poisson
realiza pré-
somatórios
das contribuições dos elementos da matriz de densidades
ߜߩ
ሺ
ݎ
ሻ
para os cálculos
(um somatório) e das forças de desloc
amento atômicas (seis
somatórios são realizados para cada bloco. Assim, p
se os valores finais, foi utilizado um segundo
kernel
, o qual utiliza
bytes
de memória compartilhada para cada bloco com 256
Nível de Ocupação do Multiprocessador
é de 100 %, conforme o
65
versão do
somatórios
para os cálculos
amento atômicas (seis
somatórios são realizados para cada bloco. Assim, para
, o qual utiliza
13
de memória compartilhada para cada bloco com 256
é de 100 %, conforme o
PUC-Rio - Certificação Digital Nº 0812691/CB
Descrição das Alterações
para
Figura 33 –
Nível de Ocupação do
da
energia de Hartree e das forças de deslocamento atô
A parte paralela dos sete somatórios
kernel
da Equação de Poisson
poderia ser constituída de sete replicas do padrão
se isto fosse feito, teríamos um elevado número de
warps incompletos. A
Tabela
35 operações de adição realizadas com o
Tabela 6 – Operações
de adição
padrão ilustrado na Figura 27.
N Total Adições / Bloco
Threads
896
448
224
112
56
28
14
7
para
Aceleração do SIESTA por GPU
Nível de Ocupação do
Multiprocessador para o kernel dos
somatório
energia de Hartree e das forças de deslocamento atômicas
.
A parte paralela dos sete somatórios
,
realizados de forma combinada no
da Equação de Poisson
, ou no kernel
para obtenção dos valores finais,
poderia ser constituída de sete replicas do padrão ilustrado na
Figura
27
se isto fosse feito, teríamos um elevado número de operações realizadas
Tabela
6 mostra que, para cada bloco de
threads
35 operações de adição realizadas com o
warp incompleto.
de adição
com warp incompleto por bloco,
replicando se
Operação com
Threads
Ativas Adições / Thread
Completo
Incompleto
128 7 28
64 7 14
32 7 7
16 7
8 7
4 7
2 7
1 7
49
66
somatório
s
realizados de forma combinada no
para obtenção dos valores finais,
27
. Contudo,
operações realizadas
com
threads
, teríamos
replicando se
te vezes o
Operação com
Warp
Incompleto
7
7
7
7
7
35
PUC-Rio - Certificação Digital Nº 0812691/CB
Descrição das Alterações
para
Assim, com o objetivo de reduzir as operações de ad
warps
incompletos, os sete somatórios são realizados segu
na Figura 34.
Neste caso, ao invés de declarar sete espaços na me
compartilhada, um único espaço de memória é declara
pertencentes a um mesmo somatório são armazenados m
não usadas na memória.
operações de adição com
warp
Figura 34 – So
matório paralelo de sete vetor
GPU.
Tabela 7 – Operações
de adição com
Figura 34.
N Total Adições / Bloco
Threads
896
448
224
112
56
28
14
7
O s
omatório paralelo d
34, realiza um
padrão de acesso linear
portanto, não existem conflitos de banco. Na etapa
as threads
que calculam os valores a
compartilhada segundo um padrão de seis intervalos
acessado. Conforme o ilustrado na esquerda da
também está livre de conflitos de bancos.
para
Aceleração do SIESTA por GPU
Assim, com o objetivo de reduzir as operações de adição realizadas com
incompletos, os sete somatórios são realizados seguindo o padrão ilustrado
Neste caso, ao invés de declarar sete espaços na me
compartilhada, um único espaço de memória é declarado. Os elementos
pertencentes a um mesmo somatório são armazenados mantendo
-
se seis posições
Conforme o apresentado na Tabela 7
, o número de
warp
s incompletos foi reduzido para 5.
matório paralelo de sete vetor
es na memória compartilhada da
de adição com
warp incompleto por bloco,
com o padrão ilustrado na
Operação com
Threads
Ativas Adições / Thread
Completo
Incompleto
256 3 24
128 1 4
256 1 8
192 1 6
224 1 7
112 1 3
56 1 1
28 1
14 1
7 1
53
omatório paralelo d
os sete vetores
, da forma como o ilustrado na
padrão de acesso linear
na memória compartilhada da GPU
portanto, não existem conflitos de banco. Na etapa anterior ao somatório paralelo,
que calculam os valores a
serem somados, acessam a memória
compartilhada segundo um padrão de seis intervalos entre cada elemento
acessado. Conforme o ilustrado na esquerda da
Figura 35, este t
ipo de acesso
também está livre de conflitos de bancos.
67
ição realizadas com
indo o padrão ilustrado
Neste caso, ao invés de declarar sete espaços na memória
do. Os elementos
se seis posições
, o número de
es na memória compartilhada da
com o padrão ilustrado na
Operação com
Warp
Incompleto
1
1
1
1
1
5
, da forma como o ilustrado na
Figura
na memória compartilhada da GPU
e,
anterior ao somatório paralelo,
serem somados, acessam a memória
entre cada elemento
ipo de acesso
PUC-Rio - Certificação Digital Nº 0812691/CB
Descrição das Alterações
para
Figura 35 –
Padrões usados de acesso à memória compartilhada, s
banco.
4.2.3.
Cálculo do Dipolo Elétrico
Uma versão CUDA
d
do SIESTA
, onde os laços foram substituídos pelo paralelismo
As contribuições dos elementos da matriz de densida
Elétrico são calculadas no
kernel
realizado.
Cada bloco deste
3104 bytes
de memória compartilhada. O
Multiprocessador
, apresentado na
para
Aceleração do SIESTA por GPU
Padrões usados de acesso à memória compartilhada, sem conflitos de
Cálculo do Dipolo Elétrico
d
esta função foi escrita para substituir a
função original
, onde os laços foram substituídos pelo paralelismo de dados da GPU.
As contribuições dos elementos da matriz de densidades
ߜߩ
ሺ
ݎ
ሻ
para o
kernel
e em seguida, um somatório
paralelo em GPU
Cada bloco deste
kernel, com 256 threads
, utiliza 18 registradores e
de memória compartilhada. O
Nível de Ocupação do
, apresentado na
Figura 36, é de 75 % para este caso.
68
em conflitos de
função original
de dados da GPU.
para o
Dipolo
paralelo em GPU
é
, utiliza 18 registradores e
Nível de Ocupação do
PUC-Rio - Certificação Digital Nº 0812691/CB
Descrição das Alterações
para
Figura 36 –
Nível de Ocupação do Multiprocessador
Dipolo Elétrico.
Para o cálculo do
Dipolo Elétrico
somatórios. O padrão ilustrado na
contudo, teríamos 15 operações de adição realizadas
(Tabela 8).
Tabela 8 – Operações
de adição com
padrão ilustrado na Figura 27.
N Total Adições / Bloco
Threads
384
192
96
48
24
12
6
3
A realização dos três somatórios de forma combinada
ilustrado na Figura 37
, reduz o número operações de adição realizadas com
incompletos para 5, conforme os cálculos apresentad
para
Aceleração do SIESTA por GPU
Nível de Ocupação do Multiprocessador
para o kernel de
Cálculo do
Dipolo Elétrico
são necessários
a realização de três
somatórios. O padrão ilustrado na
Figura 27 poderia
ser repetido três vezes,
contudo, teríamos 15 operações de adição realizadas com
warps
incomp
de adição com
warp incompleto por bloco,
replicando
Operação com
Threads
Ativas Adições / Thread
Completo
Incompleto
128 3 12
64 3 6
32 3 3
16 3
8 3
4 3
2 3
1 3
21
A realização dos três somatórios de forma combinada, usando o padrão
, reduz o número operações de adição realizadas com
incompletos para 5, conforme os cálculos apresentados na
Tabela 9.
Neste caso,
69
Cálculo do
a realização de três
ser repetido três vezes,
incomp
letos
replicando
três vezes o
Operação com
Warp
Incompleto
3
3
3
3
3
15
, usando o padrão
, reduz o número operações de adição realizadas com
warps
Neste caso,
PUC-Rio - Certificação Digital Nº 0812691/CB
Descrição das Alterações
para
os
elementos pertencentes a um mesmo somatório
duas
posições não usadas na memória
Figura 37 – So
matório paralelo de três vetores na memória compart
GPU.
Tabela 9 – Operações
de adição com
Figura 37.
N Total Adições / Bloco
Threads
384
192
96
48
24
12
6
3
Durante a etapa do somatório paralelo dos
memória compartilhada é realizado de forma linear,
cálculo dos elementos a serem somados, na etapa ant
acesso com dois interva
compartilhada. Este padrão está ilustrado na direit
visto, também está livre
de conflitos de
para
Aceleração do SIESTA por GPU
elementos pertencentes a um mesmo somatório
são
armazenados mantendo
posições não usadas na memória
.
matório paralelo de três vetores na memória compart
de adição com
warp incompleto por bloco,
com o padrão ilustrado na
Operação com
Threads
Ativas Adições / Thread
Completo
Incompleto
256 1 8
128 1 4
192 1 6
96 1 3
48 1 1
24 1
12 1
6 1
3 1
22
Durante a etapa do somatório paralelo dos
três vetores, Figura 37
, o acesso à
memória compartilhada é realizado de forma linear, sem conflitos de banco. O
cálculo dos elementos a serem somados, na etapa anterior, emprega um padrão de
acesso com dois interva
los entre cada posição acessada na memória
compartilhada. Este padrão está ilustrado na direita da
Figura 35 e,
como pode ser
de conflitos de
bancos.
70
armazenados mantendo
-se
matório paralelo de três vetores na memória compartilhada da
com o padrão ilustrado na
Operação com
Warp
Incompleto
1
1
1
1
1
5
, o acesso à
sem conflitos de banco. O
erior, emprega um padrão de
los entre cada posição acessada na memória
como pode ser
PUC-Rio - Certificação Digital Nº 0812691/CB
Descrição das Alterações para Aceleração do SIESTA por GPU 71
4.2.4.
Reordenação de Dados
O SIESTA emprega duas formas de organizar os elementos da matriz de
densidades,
ߜߩ
ሺ
ݎ
ሻ
, e do potencial de Hartree, ߜܸ
ு
ሺ
ݎ
ሻ
. Assim, a função de
reordenação de dados faz o mapeamento entre estas duas organizações diferentes.
Antes da execução da Equação de Poisson é feita a reordenação da matriz de
densidades e após é feita a reordenação do potencial de Hartree. A reordenação é
realizada também antes e depois do Cálculo do Dipolo Elétrico. O principal
motivo de executar a reordenação de dados em GPU é a economia de
transferências de memória através do barramento PCI Express. Pois ela é
executada sobre os mesmos conjuntos de dados empregados para as demais
funções GPU anteriormente descritas.
Nesta função foi utilizada a otimização de transferência de memória,
conforme o descrito no item 3.5.3. Assim, são empregados múltiplos fluxos de
processamento e a transferência assíncrona de dados através do barramento PCI
Express. Isto é feito com a finalidade de redução do tempo total de execução pela
sobreposição dos tempos de transferência de memória com o tempo de
processamento da reordenação.
O kernel da Reordenação de Dados necessita de 18 registradores para cada
bloco com 256 threads. O Nível de Ocupação do Multiprocessador para este caso
é de 75 %, conforme o apresentado na Figura 38.
PUC-Rio - Certificação Digital Nº 0812691/CB
Descrição das Alterações
para
Figura 38 –
Nível de Ocupação do Multiprocessador
Reordenação de Dados.
Os Níveis
de Ocupação do Multiprocessador
descritos e utilizados neste trabalho estão dentro
NVIDIA
. Pois é sugerido um valor mínimo
capacidade de cálculo 1.2 ou superior.
Ocupação do Multiprocess
práticos para o aumento de desempenho
para
Aceleração do SIESTA por GPU
Nível de Ocupação do Multiprocessador
para
o
de Ocupação do Multiprocessador
empregados nos
descritos e utilizados neste trabalho estão dentro dos limites recomendados pela
. Pois é sugerido um valor mínimo
de 18.75 %
para dispositivos com
capacidade de cálculo 1.2 ou superior.
Por outro lado, o aumento do
Ocupação do Multiprocess
ador
além do limite de 50 % não traz benefícios
práticos para o aumento de desempenho
(NVIDIA, 2009a).
72
o
kernel da
empregados nos
kernels
dos limites recomendados pela
para dispositivos com
Por outro lado, o aumento do
Nível de
além do limite de 50 % não traz benefícios
PUC-Rio - Certificação Digital Nº 0812691/CB
5
Estudo de Casos
Neste capítulo, serão descritos os sistemas físicos e as suas propriedades
calculadas com o objetivo comparar o desempenho computacional de uma GPU
com o desempenho obtido com uma CPU.
Na realização deste estudo foi utilizada uma GPU Tesla C1060, a qual
possui 240 elementos de processamento (1.3 GHz) e 4 GB de memória RAM com
largura de banda de acesso à memória de 102 GB/segundo através de um
barramento de 512 bits (NVIDIA, 2010). O computador utilizado possuía quatro
destas GPUs e um processador AMD Phenom II (4 núcleos), trabalhando na
freqüência de 3 GHz, com 16 GB de memória RAM. Todas as comparações foram
realizadas entre uma GPU e um núcleo da CPU. A largura de banda do
barramento PCI Express (Peripheral Component Interconnect) pode ser
considerada baixa. Pois enquanto a largura de banda das transferências internas,
entre memória GDDR3 e elementos de processamento, de uma Tesla C1060 é de
102 GB/segundo, a transferência pela PCI Express fica em torno de 2
GB/segundo.
5.1.
Caso 1 – Cálculo da Energia e Estrutura de Bandas de Nanotubos
Os nanotubos de carbono são estruturas cristalinas (cilíndricas) formadas
por átomos de carbono (estruturas alotrópicas de carbono). Possuem alta
resistência a tensão mecânica, podendo ser usados como aditivos em compostos
para melhorar suas características. Dependendo do diâmetro ou da quiralidade do
tubo, os nanotubos de carbono são condutores ou semi-condutores, podendo ter
aplicações em circuitos micro e nano-eletrônicos.
PUC-Rio - Certificação Digital Nº 0812691/CC
Estudos de Casos 74
5.1.1.
Estrutura do Nanotubo
O vetor quiral define a estrutura atômica dos nanotubos e são usadas
denominações especiais para cada caso: os nanotubos (n,n) são denominados
armchair, enquanto os nanotubos (n,0) são denominados zigzag. Os nanotubos
(n,m), com n
m e m 0, são denominados genericamente quirais, em oposição
aos nanotubos armchair e zigzag, que são aquirais.
A ilustração da Figura 39 (a) mostra o cruzamento das bandas de valência e
condução no nível de Fermi, sendo este comportamento característico dos
nanotubos metálicos ou nanotubos armchair (n,n). A Figura 39 (b) ilustra a
estrutura de banda de um nanotubo (8,0) semicondutor. Os nanotubos (n,m) com
n - m 3i, onde i é um número inteiro, apresentam um comportamento
semicondutor. Já os nanotubos (n,m) com n - m = 3i são classificados como
semicondutores com gap quase nulo. Para este caso, o nanotubo (12,0) é ilustrado
na Figura 39 (c).
Figura 39 – Estrutura de bandas de nanotubos (a) armchair (5, 5), (b) zigzag (8, 0)
e (c) zigzag (12, 0). O nível de Fermi está deslocado para o zero, indicado pela
linha tracejada (Silva, 2008).
Uma ilustração da estrutura geométrica da célula unitária do nanotubo (4,2)
com 56 átomos de carbono, utilizado para os cálculos de comparação de
desempenho entre CPU e GPU realizados neste trabalho, é apresentada na Figura
PUC-Rio - Certificação Digital Nº 0812691/CC
Estudos de Casos
40
para a vista frontal e a lateral
comportamento semicondutor é esperado.
Figura 40 – Nanotubo
(4,2)
5.1.2.
Resultados
e Testes de Desempenho
As propriedades físicas, energia e estrutura de ban
foram calculadas por primeiros princípios fazendo
emprego do programa SIESTA.
utilizada foi a LDA, com a parametrização desenvolv
(1980).
Os cálculos de energia para o nanotubo
para cada teste. Uma vez com a versão seqüencial or
em seguida com a versão que emprega funções alterad
Isto foi feit
o para verificar a acurácia dos resultados produzid
funções CUDA e para comparação de desempenho. Confo
Tabela 10, Tabela 11 e
Tabela 12
da Energia Total encontra-
se na sexta casa decimal.
A
resolução da Equação de Poisson na GPU é utilizada
potencial de Hartree. A energia de Hartree e a cont
atômicas de deslocamento,
e encontram-
se combinadas na mesma porção do código para resolu
Equação de Poisson
, sendo também calculadas na GPU. A
resultados obtidos com a utilização da
Poisson, a qual
apresentava
foi descrito no Capítulo 4.
Estudos de Casos
para a vista frontal e a lateral
. Uma vez que 4 - 2 3i
, então um
comportamento semicondutor é esperado.
(4,2)
. a) Vista Frontal, b) Vista Lateral.
e Testes de Desempenho
As propriedades físicas, energia e estrutura de bandas, do nanotubo
foram calculadas por primeiros princípios fazendo
-
se o uso da DFT, com o
emprego do programa SIESTA.
A aproximação para o termo de troca e correlação
utilizada foi a LDA, com a parametrização desenvolvida por
Ceperley
Os cálculos de energia para o nanotubo
(4,2)
foram realizados duas vezes
para cada teste. Uma vez com a versão seqüencial original do SIESTA na CPU e
em seguida com a versão que emprega funções alteradas para execução em GPU.
o para verificar a acurácia dos resultados produzidos com emprego de
funções CUDA e para comparação de desempenho. Conforme o apresentado nas
Tabela 12
, é possível verificar
, que para este caso, o erro
se na sexta casa decimal.
resolução da Equação de Poisson na GPU é utilizada para calcular o
potencial de Hartree. A energia de Hartree e a contribuição para as forças
devidas ao potencial de Hartree,
envolvem somatórios
se combinadas na mesma porção do código para resolu
, sendo também calculadas na GPU. A
Tabela 10
resultados obtidos com a utilização da
primeira versão do kernel
da Equação de
apresentava
um nível de ocupação do multiprocessador
75
, então um
das, do nanotubo
(4,2)
se o uso da DFT, com o
A aproximação para o termo de troca e correlação
Ceperley
& Alder
foram realizados duas vezes
iginal do SIESTA na CPU e
as para execução em GPU.
os com emprego de
rme o apresentado nas
, que para este caso, o erro
para calcular o
ribuição para as forças
envolvem somatórios
se combinadas na mesma porção do código para resolução
da
mostra os
da Equação de
de 75 % e
PUC-Rio - Certificação Digital Nº 0812691/CC
Estudos de Casos 76
Tabela 10 – Resultados de Aceleração com a primeira versão do kernel da Equação de
Poisson, para o nanotubo (4,2).
CPU GPU Aceleração
Energia Total
-8592.057279 eV -8592.057281 eV
Tempo da Equação de Poisson
988.0372 s 34.7373 s
28.4431 X
Um segundo kernel foi desenvolvido para resolução da Equação de Poisson
na GPU. Conforme o descrito no Capítulo 4, este kernel possui uma maior
intensidade aritmética e um menor nível de ocupação do multiprocessador, 50 %.
Os resultados obtidos para este caso, Tabela 11, são ligeiramente superiores aos
do primeiro kernel.
Tabela 11 – Resultados de Aceleração com a segunda versão do kernel da Equação de
Poisson, para o nanotubo (4,2).
CPU GPU Aceleração
Energia Total
-8592.057279 eV -8592.057282 eV
Tempo da Equação de Poisson
988.0372 s 33.796 s
29.2353 X
Ao portar novas funções de cálculo do dipolo elétrico e de reordenação de
dados para execução em GPU, foi reduzida a necessidade de transferência de
memória através do barramento PCI Express. A matriz de densidade, necessária
para a resolução da Equação de Poisson, já se encontrava na memória da GPU
devido à execução das outras funções. Isso contribui para aumentar a aceleração
da Equação de Poisson, conforme o observado na Tabela 12.
O Cálculo do Dipolo Elétrico é uma pequena porção do código, quando
comparada a Equação de Poisson, envolvendo apenas três somatórios. Não há a
necessidade de operações mais dispendiosas computacionalmente, tais como a
Transformada de Fourier, e não são realizadas transferências de memória. Esta
função é empregada sobre a matriz de densidade, a qual é transferida para a GPU
durante a Reordenação de Dados. Isso resultou num elevado valor de aceleração,
de duas ordens de grandeza, Tabela 12.
A Reordenação de Dados é responsável pelas transferências de memória
relacionadas com a matriz de densidade. Apesar de a transferência ser realizada
empregando transferência assíncrona e execução concorrentes, a baixa intensidade
aritmética desta função resulta num valor menor de aceleração de 2.9 quando
comparado às demais funções, Tabela 12. Esta função foi portada para GPU por
atuar sobre os mesmos conjuntos de dados que a Equação de Poisson. Portar todas
as funções relacionadas com a matriz de densidades para execução em GPU
PUC-Rio - Certificação Digital Nº 0812691/CC
Estudos de Casos 77
poderá eliminar a necessidade de sua transferência através do barramento PCI
Express. Seria necessário transferir somente os parâmetros de inicialização e o
resultado das integrais.
Tabela 12 – Resultados de Aceleração Reduzindo-se a Transferência de Dados através do
Barramento PCI Express, para o nanotubo (4,2).
CPU GPU Aceleração
Energia Total
-8592.057279 eV -8592.057282 eV
Tempo de Cálculo do Dipolo Elétrico
51.6835 s 0.2093 s
246.9350 X
Tempo de Reordenação de Dados
105.5779 s 39.0563 s
2.7032 X
Tempo da Equação de Poisson
988.0372 s 26.6871 s
37.0230 X
O tempo total de execução do SIESTA, para o nanotubo (4,2), com a
Equação de Poisson, o Cálculo do Dipolo Elétrico e a Reordenação de Dados em
GPU, foi de 353.2323 segundos. A versão seqüencial original do SIESTA
consumiu 1390.1599 segundos. Assim, a aceleração global do código levando em
consideração as chamadas de função C a partir do Fortran, para este caso, foi de
3.9355.
A estrutura de bandas do nanotubo (4,2), obtida nos cálculos realizados com
o SIESTA, é apresentada na Figura 41. Pode se verificar que está coerente com o
fato do nanotubo (4,2) ser um semicondutor (4 - 2 3i,), pois não ocorre o
cruzamento das bandas de valência e condução no nível de Fermi, de forma
semelhante ao discutido anteriormente para a Figura 39 (b).
PUC-Rio - Certificação Digital Nº 0812691/CC
Estudos de Casos 78
Figura 41 – Estrutura de bandas para o nanotubo (4,2). O nível de Fermi está
deslocado para o zero, indicado pela linha tracejada.
5.2.
Caso 2 – Otimização da Geometria Estrutural de Fulereno
O fulereno é uma estrutura formada por átomos de carbono organizados nos
vértices de um icosaedro truncado que tem a forma de uma bola de futebol (com
pentágonos e hexágonos).
5.2.1.
Estrutura do Fulereno
A estrutura do fulereno é esférica, formada por hexágonos interligados por
pentágonos, sendo estes últimos responsáveis pela sua curvatura e,
conseqüentemente, por sua forma tridimensional (Kroto, Allaf e Balm, 1991). O
representante mais conhecido da família dos fulerenos é o C60 (com 60 carbonos),
com um diâmetro de aproximadamente 1 nm. Este fulereno é apresentado na
Figura 42 e é nosso segundo estudo de caso.
PUC-Rio - Certificação Digital Nº 0812691/CC
Estudos de Casos 79
Figura 42 – Fulereno C60.
5.2.2.
Resultados e Testes de Desempenho
Uma estrutura muito próxima da estrutura final foi utilizada como ponto de
partida para otimização da geometria estrutural do fulereno C60. A energia total
da estrutura é calculada e pequenas variações das coordenadas atômicas são
realizadas, de acordo com os resultados obtidos para as forças de deslocamento
atômicas. Foi utilizada a DFT do programa SIESTA, com a aproximação GGA
para o termo de troca e correlação da energia. Os cálculos são repetidos duas
vezes, em CPU e em CPU-GPU, para comparação de desempenho. Com relação à
acurácia, para este caso, o Erro da Energia Total encontra-se na quarta casa
decimal, conforme pode ser observado nas Tabela 13 e Tabela 14. A otimização
da geometria estrutural requer a realização de mais de um cálculo de energia. Isto
justifica o aumento do erro, neste caso, pois o erro do cálculo de energia é
acumulado entre os passos do processo iterativo de otimização da geometria.
Juntamente com o código para o cálculo do potencial de Hartree pela
resolução da Equação de Poisson na GPU, encontra-se, de forma combinada, o
código para o cálculo da energia de Hartree e a contribuição para as forças de
deslocamento atômicas, devidas ao potencial de Hartree. A aceleração obtida pela
resolução da Equação de Poisson na GPU é apresentada na Tabela 13. Neste caso
foi utilizada somente a segunda versão do kernel da Equação de Poisson,
conforme o descrito no Capítulo 4.
PUC-Rio - Certificação Digital Nº 0812691/CC
Estudos de Casos 80
Tabela 13 – Resultados de Aceleração com a segunda versão do kernel da Equação de
Poisson, para o fulereno C60.
CPU GPU Aceleração
Energia Total
-7853.604264 eV -7853.604698 eV
Tempo da Equação de Poisson
187.1953 s 10.7989 s
17.3346 X
Portar um maior volume de operações, realizadas sobre a mesma matriz de
densidades, para a GPU contribui para melhorar o desempenho. Comparando a
Tabela 13 com a Tabela 14, podemos observar o aumento de desempenho da
Equação de Poisson na GPU, pois a matriz de densidade havia sido usada por
outras funções GPU, não necessitando a transferência através do barramento PCI
Express.
Tabela 14 – Resultados de Aceleração Reduzindo-se a Transferência de Dados através do
Barramento PCI Express, para o fulereno C60.
CPU GPU Aceleração
Energia Total
-7853.604264 eV -7853.604698 eV
Tempo de Cálculo do Dipolo Elétrico
14.4850 s 0.0774 s
187.1447 X
Tempo da Equação de Poisson
187.1953 s 8.9205 s
20.9848 X
PUC-Rio - Certificação Digital Nº 0812691/CC
6
Conclusões e Trabalhos Futuros
Os métodos iterativos empregados na química quântica, como a Teoria do
Funcional da Densidade, são computacionalmente dispendiosos devido ao grande
número de integrais realizadas para obter-se o valor da energia total do sistema
físico em estudo. A dificuldade aumenta com o número de átomos e elétrons
presentes na simulação. A energia de um sistema físico é, em geral, dividida em
cinco termos principais. Com o objetivo de investigar o apoio das GPUs na
aceleração de partes dos cálculos envolvidos na DFT, o termo do cálculo do
potencial e da energia de Hartree foi alterado para execução em GPU. Métodos de
mecânica quântica geralmente possuem um código fonte muito extenso e a
codificação usando CUDA é trabalhosa e necessita de cuidados especiais. Entre
estes cuidados, temos a memória compartilhada que está sujeita a conflitos de
bancos que precisam ser resolvidos manualmente pelo programador para evitar a
redução de desempenho. O acesso à memória global, que não possui cache,
precisa ser projetado para tirar proveito da capacidade de leitura de blocos de
dados de 8, 16 ou 32 palavras, em paralelo, em uma única operação. As threads de
GPUs empregam paralelismo finamente granulado e, preferencialmente, devem
ser criadas em múltiplos de 32. Outro ponto importante é minimizar as
transferências de memória através do barramento PCI Express, portando partes do
código que atuam sobre um mesmo conjunto de dados.
Neste sentido, foram estudadas as equações fundamentais da Teoria do
Funcional da Densidade, permitindo identificar trechos de código com
características adequadas para o emprego do paralelismo de dados oferecido pelas
Unidades de Processamento Gráfico. Entre os algoritmos de código fonte aberto
disponíveis, optou-se por fazer uso do SIESTA, o qual foi desenvolvido visando à
eficiência computacional.
As integrais para cálculo da energia cinética dos elétrons e o cálculo da
contribuição para a energia total relativo ao termo não-local do pseudopotencial
são adequados para aceleração por GPUs. Os termos da energia de troca e
PUC-Rio - Certificação Digital Nº 0812691/CA
Conclusões e Trabalhos Futuros 82
correlação, o potencial e a energia de Hartree e a parte local do pseudopotencial
do hamiltoniano do SIESTA são calculados sobre uma mesma matriz de
densidades.
O termo referente ao potencial e a energia de Hartree foi acelerado através
da execução em GPUs, conseguindo-se um ganho em torno de 30 vezes. A
reordenação de dados em GPU apresentou um ganho inferior às demais funções
para GPU, contudo sua importância é devida à redução do volume de
transferências de memória através do barramento PCI Express. O cálculo do
dipolo elétrico pode ser acelerado até duas ordens de grandeza na GPU. A
continuação do trabalho deveria se concentrar prioritariamente nas integrais para o
cálculo da energia relativo ao termo local do pseudopotencial do SIESTA e no
cálculo da energia de troca e correlação. Assim, a matriz de densidade poderia ser
inicializada e trabalhada na memória do dispositivo, reduzindo drasticamente as
transferências de memória entre CPU e GPU, pois seriam transferidos somente os
parâmetros de inicialização e os resultados das integrais.
Cálculos da energia e da estrutura de bandas de nanotubos e da otimização
da geometria estrutural de fulerenos foram realizados usando a versão seqüencial
original e a versão alterada para a execução de trechos de código do SIESTA em
GPUs. Estes cálculos foram realizados para comprovar a eficiência das GPUs em
calcular o termo do potencial de Hartree, a energia de Hartree e o dipolo elétrico.
A aplicação de técnicas inteligentes (redes neurais, algoritmos genéticos e
programação genética), em trabalhos futuros, para redução do esforço necessário
para obtenção de código fonte eficiente para GPUs contribuiria para facilitar o
desenvolvimento de aplicações. Neste caso, uma ferramenta deveria ser projetada
para administrar o uso das memórias, evitando conflitos de banco na memória
compartilhada e minimizando o desperdício de largura de banda da memória
global. Um desafio a ser enfrentado seria a coordenação da execução concorrente
de um grande número de threads leves de forma automática, sem perder em
eficiência. Sua utilização para a evolução de uma Transformada de Fourier mais
eficiente que a atualmente disponibilizada pela NVIDIA, poderia ser empregada
na resolução da equação de Poisson do SIESTA, para melhorar a aceleração
obtida no termo do potencial e da energia de Hartree.
A continuação do trabalho, alterando o código fonte para execução em
GPUs dos demais termos sugeridos, poderia trazer reduções significativas do
PUC-Rio - Certificação Digital Nº 0812691/CA
Conclusões e Trabalhos Futuros 83
tempo total de cálculo do programa SIESTA. Isto contribuiria para as pesquisas
realizadas na química computacional e na física dos materiais, reduzindo o tempo
de espera por previsões precisas baseadas em cálculos teóricos ab initio.
PUC-Rio - Certificação Digital Nº 0812691/CA
7
Referências Bibliográficas
AMDAHL, G. M. Validity of the Single Processor Approach to Achieving Large-
Scale Computing Capabilities. Spring Joint Computer Conference, Sunnyvale,
p. 483-485, Abril 1967.
ANDERSON, A. G.; GODDARD, W. A.; SCHRÖDER, P. Quantum Monte Carlo
on graphical processing units. Computer Physics Communications, v. 177, n. 3,
p. 298-306, 30 Março 2007.
ANGLADA, E. et al. Systematic generation of finite-range atomic basis sets for
linear-scaling calculations. Physical Review B, v. 66, p. 205101, 7 Novembro
2002.
ARTACHO, E. et al. Linear-Scaling ab-initio Calculations for Large and Complex
Systems. physica status solidi (b), v. 215, p. 809-817, 3 Setembro 1999.
ARTACHO, E. et al. SIESTA User's Guide 2.0.2. Fundación General
Universidad Autónoma de Madrid. Madrid, p. 89. 2008.
BARNEY, B. Lawrence Livermore National Laboratory. Introduction to
Parallel Computing, 2009. Disponivel em:
<https://computing.llnl.gov/tutorials/parallel_comp/>. Acesso em: 5 Janeiro 2010.
BORN, M.; OPPENHEIMER, R. On the quantum theory of molecules. Annalen
der Physik, v. 84, n. 20, p. 457-484, 25 Agosto 1927.
CEPERLEY, D. M.; ALDER, B. J. Ground State of the Electron Gas by a
Stochastic Method. Physical Review Letters, v. 45, n. 7, p. 566-569, 18 Agosto
1980.
DEMOCRITOS. Quantum ESPRESSO, 2009. Disponivel em:
<http://www.quantum-espresso.org/index.php>. Acesso em: 20 Janeiro 2010.
DUNCAN, R. A Survey of Parallel Computer Architectures. Computer, v. 23, n.
2, p. 5-16, Fevereiro 1990.
EM PHOTONICS. CULA tools. Function List, 2010. Disponivel em:
<http://www.culatools.com/functions>. Acesso em: 15 Fevereiro 2010.
FAGAN, S. B. Funcionalização de Nanotubos de Carbono: uma Abordagem
de Primeiros Princípios. Universidade Federal de Santa Maria. Santa Maria, p.
207. 2003.
PUC-Rio - Certificação Digital Nº 0812691/CA
Referências Bibliográficas 85
FERMI, E. Un Metodo Statistico per la Determinazione di alcune Priorieta
dell'Atome. Rend. Accad. Naz. Lincei, v. 6, p. 602-607, 1927.
FLYNN, M. J. Some Computer Organizations and Their Effectiveness. IEEE
Transactions on Computers, v. 21, n. 9, p. 948-960, Setembro 1972. ISSN 0018-
9340.
GARCÍA, A. ATOM User Manual. Universidad del País Vasco. Bilbao,
Espanha, p. 17. 2006.
GENOVESE, L. et al. Density Functional Theory calculation on many-cores
hybrid CPU-GPU architectures. ArXiv e-prints, n. 0904.1543, 9 Abril 2009.
GÖDDEKE, D.; STRZODKA, R.; TUREK, S. Performance and accuracy of
hardware-oriented native-, emulated-and mixed-precision solvers in FEM
simulations. International Journal of Parallel, Emergent and Distributed
Systems, v. 22, n. 4, p. 221-256, 4 Janeiro 2007.
HAMANN, D. R.; SCHLÜTER, M.; CHIANG, C. Norm-Conserving
Pseudopotentials. Physical Review Letters, v. 43, n. 20, p. 1494-1497, 12
Novembro 1979.
HARRIS, M. Mapping computational concepts to GPUs. ACM SIGGRAPH
2005 Courses, New York, v. 31, p. 493-508, Agosto 2005.
HILLIS, W. D.; STEELE, G. L. Data Parallel Algorithms. Communications of
the ACM, v. 29, n. 12, p. 1170-1183, Dezembro 1986.
HOHENBERG, P.; KOHN, W. Inhomogeneous Electron Gas. Physical Review,
v. 136, n. 3B, p. B864-B871, 9 Novembro 1964.
JANSSEN, C. et al. The Massively Parallel Quantum Chemistry Program,
2010. Disponivel em: <http://www.mpqc.org/index.php>. Acesso em: 5 Janeiro
2010.
JUNQUERA, J. et al. Numerical atomic orbitals for linear-scaling calculations.
Physical Review B, v. 64, p. 235111, 28 Novembro 2001.
KELMELIS, E. J. et al. Modeling and simulation of nanoscale devices with a
desktop supercomputer. Nanomodeling II, v. 6328, 8 Setembro 2006.
KLEINMAN, L.; BYLANDER, D. M. Efficacious Form for Model
Pseudopotentials. Physical Review Letters, v. 48, n. 20, p. 1425-1428, 17 Maio
1982.
KOCH, W.; HOLTHAUSEN, M. C. A Chemist's Guide to Density Functional
Theory. 2ª Edição. ed. Weinheim: Wiley-VCH, 2001. ISBN 3-527-30372-3.
KOHN, W.; SHAM, L. J. Self-Consistent Equations Including Exchange and
Correlation Effects. Physical Review, v. 140, n. 4A, p. A1133-A1138, 15
Novembro 1965.
PUC-Rio - Certificação Digital Nº 0812691/CA
Referências Bibliográficas 86
KROTO, H. W.; ALLAF, A. W.; BALM, S. P. C60: Buckminsterfullerene.
Chemical Reviews, v. 91, n. 6, p. 1213-1235, Julho 1991.
LANDRUM, G.; GLASSEY, W. The YAeHMOP Home Page, 2001. Disponivel
em: <http://yaehmop.sourceforge.net/>. Acesso em: 15 Dezembro 2009.
MEEL, J. A. V. et al. Harvesting graphics power for MD simulations. Molecular
Simulation, v. 34, n. 3, p. 259-266, 2 Fevereiro 2008.
MOORE, G. E. Cramming More Components Onto Integrated Circuits.
Proceedings of the IEEE, v. 86, n. 1, p. 82-85, Janeiro 1998.
NUKADA, A. et al. Bandwidth Intensive 3-D FFT kernel for GPUs using CUDA.
Proceedings of the 2008 ACM/IEEE conference on Supercomputing, Austin,
n. 5, p. 11, Novembro 2008. ISSN 978-1-4244-2835-9.
NVIDIA. CUBLAS Library. NVIDIA. Santa Clara, CA, p. 80. 2007a.
NVIDIA. CUDA CUFFT Library. NVIDIA. Santa Clara, CA, p. 17. 2007b.
NVIDIA. NVIDIA CUDA C Programming Best Practices Guide. NVIDIA.
Santa Clara, CA, p. 66. 2009a.
NVIDIA. NVIDIA CUDA Programming Guide. NVIDIA. Santa Clara, CA, p.
145. 2009b.
NVIDIA. CUDA Zone, 2010. Disponivel em:
<http://www.nvidia.com/object/cuda_home.html#>. Acesso em: 25 Janeiro 2010.
ORDEJÓN, P.; ARTACHO, E.; SOLER, J. M. Self-consistent order-N density-
functional calculations for very large systems. Physical Review B, v. 53, p.
R10441-R10444, 15 Abril 1996.
PERDEW, J. P.; BURKE, K.; ERNZERHOF, M. Generalized Gradient
Approximation Made Simple. Physical Review Letters, v. 77, n. 18, p. 3865-
3868, 28 Outubro 1996.
PRICE, D. K.; HUMPHREY, J. R.; KELMELIS, E. J. Accelerated Simulators for
Nano-Photonic Devices. Numerical Simulation of Optoelectronic Devices, p.
103-104, Setembro 2007.
SCHRÖDINGER, E. An Undulatory Theory of the Mechanics of Atoms and
Molecules. Physical Review, v. 28, n. 6, p. 1049-1070, Dezembro 1926.
SILVA, L. B. D. Campos Elétricos Transversais Sobre Nanotubos De
Carbono: Um Estudo De Primeiros Princípios. Universidade Federal de Santa
Maria. Santa Maria, p. 103. 2008.
SOLER, J. M. et al. The SIESTA method for ab initio order-N materials
simulation. Journal of Physics: Condensed Matter, v. 14, p. 2745-2779, 8
Março 2002.
PUC-Rio - Certificação Digital Nº 0812691/CA
Referências Bibliográficas 87
STONE, J. E. et al. Accelerating Molecular Modeling Applications with Graphics
Processors. Journal of Computational Chemistry, v. 28, n. 16, p. 2618-2640, 25
Setembro 2007.
TEMPERTON, C. A Generalized Prime Factor FFT Algorithm for any N=2p3q5r.
SIAM Journal on Scientific and Statistical Computing, v. 13, n. 3, p. 676-686,
Maio 1992.
THOMAS, L. H. The calculation of atomic fields. Proceedings of the
Cambridge Philosophical Society, v. 23, n. 05, p. 542-548, Janeiro 1927.
TOMOV, S. et al. MAGMA Users’ Guide. Univ. of Tennessee; Univ. of
California; Univ. of Colorado. Knoxville; Berkeley; Denver, p. 83. 2009.
TROULLIER, N.; MARTINS, J. L. Efficient pseudopotentials for plane-wave
calculations. Physical Review B, v. 43, n. 3, p. 1993-2006, 15 Janeiro 1991.
UFIMTSEV, I. S.; MARTÍNEZ, T. J. Graphical Processing Units for Quantum
Chemistry. Computing in Science & Engineering, v. 10, n. 6, p. 26-34,
Dezembro 2008a.
UFIMTSEV, I. S.; MARTÍNEZ, T. J. Quantum Chemistry on Graphical
Processing Units. 1. Strategies for Two-Electron Integral Evaluation. Journal of
Chemical Theory and Computation, v. 4, n. 2, p. 222-231, 2008b.
VOGT, L. et al. Accelerating Resolution-of-the-Identity Second-Order Møller-
Plesset Quantum Chemistry Calculations with Graphical Processing Units.
Journal of Physical Chemistry, v. 112, n. 10, p. 2049-2057, 2008.
YASUDA, K. Accelerating Density Functional Calculations with Graphics
Processing Unit. Journal of Chemical Theory and Computation, v. 4, n. 8, p.
1230-1236, 2008a.
YASUDA, K. Two-Electron Integral Evaluation on the Graphics Processor Unit.
Journal of Computational Chemistry, v. 29, n. 3, p. 334-342, 2008b.
PUC-Rio - Certificação Digital Nº 0812691/CA
Livros Grátis
( http://www.livrosgratis.com.br )
Milhares de Livros para Download:
Baixar livros de Administração
Baixar livros de Agronomia
Baixar livros de Arquitetura
Baixar livros de Artes
Baixar livros de Astronomia
Baixar livros de Biologia Geral
Baixar livros de Ciência da Computação
Baixar livros de Ciência da Informação
Baixar livros de Ciência Política
Baixar livros de Ciências da Saúde
Baixar livros de Comunicação
Baixar livros do Conselho Nacional de Educação - CNE
Baixar livros de Defesa civil
Baixar livros de Direito
Baixar livros de Direitos humanos
Baixar livros de Economia
Baixar livros de Economia Doméstica
Baixar livros de Educação
Baixar livros de Educação - Trânsito
Baixar livros de Educação Física
Baixar livros de Engenharia Aeroespacial
Baixar livros de Farmácia
Baixar livros de Filosofia
Baixar livros de Física
Baixar livros de Geociências
Baixar livros de Geografia
Baixar livros de História
Baixar livros de Línguas
Baixar livros de Literatura
Baixar livros de Literatura de Cordel
Baixar livros de Literatura Infantil
Baixar livros de Matemática
Baixar livros de Medicina
Baixar livros de Medicina Veterinária
Baixar livros de Meio Ambiente
Baixar livros de Meteorologia
Baixar Monografias e TCC
Baixar livros Multidisciplinar
Baixar livros de Música
Baixar livros de Psicologia
Baixar livros de Química
Baixar livros de Saúde Coletiva
Baixar livros de Serviço Social
Baixar livros de Sociologia
Baixar livros de Teologia
Baixar livros de Trabalho
Baixar livros de Turismo
