Upload
geslane-schepers
View
24
Download
1
Embed Size (px)
DESCRIPTION
Programaçao
Citation preview
Paralelismo em Mecânica Computacional
Renato Elias, Marcos Martins, Alvaro Coutinho{renato, marcos, alvaro}@nacad.ufrj.br
Núcleo de Atendimento em Computação de Alto Desempenho E Departamento de Engenharia Civil
Universidade Federal do Rio de Janeiro, Brasilwww.nacad.ufrj.br
2
Conteúdo
Conceitos básicos envolvidos na paralelização de códigos de mecânica
computacional
– Revisão de métodos discretos (elementos, volumes, diferenças, finitas)
– Representação dos dados e esquemas de armazenamento
Paralelização de programas com o OpenMP
– Revisão de OpenMP;
– Blocagem de dados Mesh coloring;
– Operações básicas de códigos de mecânica computacional em OpenMP.
Paralelização de programas com o MPI
– Particionamento e distribuição de dados Metis;
– Operações básicas de códigos de mecânica computacional em MPI.
3
Revisão de métodos discretos em 6 slides (1/6)
O que é uma discretização?
infinitoanalíticolimitado
geralmente não
incógnitas tratamento
aplicabilidade aproximação
finitonuméricoamplageralmente sim
Métodos discretos: Diferenças finitas, volumes finitos e elementos finitos
4
Diferenças Finitas (2/6)
As equações diferenciais são convertidas em “stencils” nos nós de uma grade
Simples, rápido
Pouco flexível
5
Volumes Finitos (3/6)
As equações diferenciais são reformuladas em termos de fluxo através dos
contornos dos elementos
Flexível
Geração de malha em 3D difícil
6
Elementos Finitos (4/6)
Equação diferencial reformulada em termos de formulação variacional e funções de
teste
Flexível
Geração de malha em 3D difícil
7
O que todos tem em comum? (5/6)
No final sempre teremos um sistema de equações
Frequentemente muito grande para ser resolvido com um método direto
Métodos iterativos são utilizados, introduzindo um outro nível de aproximação
O paralelismo viabiliza que problemas cada vez maiores sejam resolvidos
8
Representação computacional de uma malha não estruturada (6/6)
ndim: número de dimensões espaciais
nnos: número de nós
nnoel: número de nós por elemento
nel: número de elementos
9
Formas de representação de dados:Malhas, Grades, Grafos e Matrizes
1,1 1,2 1,3 1,4 1,5 1,6 1,7
2,1 22 2,3 2,6 2,9 2,10
3,1 3,2 3,3 3,4 3,8 3,9 3,11
4,1 4,3 4,4 4,5 4,8
5,1 5,4 5,5 5,7
6,1 6,2 6,6 6,7 6,10
7,1 7,5 7,6 7,7
8,3 8,4 8,8 8,11
9,2 9,3 9,9 9,10 9,11
10,2 1
a a a a a a a
a a a a a a
a a a a a a a
a a a a a
a a a a
a a a a a
a a a a
a a a a
a a a a a
a a 0,6 10,9 10,10
11,3 11,8 11,9 11,11
a a
a a a a
Matriz esparsa
representação geométrica representação relacional representação algébrica
10
Reordenação Nodal
O perfil da matriz resultante está diretamente relacionado à ordem com que os nós da
malha estão numerados;
Algoritmos de ordenação tal como o Reverse Cuthill Mckee reordenam a malha com
intuito de reduzir a largura de banda
exemplo básico
problema real
11
Efeito da reordenação de dados
0
20
40
60
80
0 100 200 300 400 500 600
Edges x 1000
No
des
x 1
000
Node i
Node j
0
20
40
60
80
0 100 200 300 400 500 600
Edges x 1000
No
de
s x
100
0
Node i
Node j
12
Outras formas de reordenação
0
20
40
60
80
100
0 100 200 300 400 500 600 700
x 10
00
x 1000
Edges
Nod
es
Node i
Node j
Original edge order
0
20
40
60
80
100
0 100 200 300 400 500 600 700
x 10
00
x 1000
Edges
Nod
es
Node i
Node j
Improving data locality: RCM vertex reordering and edge
reordering in ascending vertex order
Improving data locality with NO memory dependencies
0
20
40
60
80
100
0 100 200 300 400 500 600 700
x 10
00
x 1000
Edges
Nod
es
Node i
Node j
0
20
40
60
80
100
0 100 200 300 400 500 600 700
x 10
00
x 1000
Edges
Nod
es
Node i
Node j
Improving data locality for each superedge type.
Superedge improves the use of CPU registers
0
20
40
60
80
100
0 100 200 300 400 500 600 700
x 10
00
x 1000
Edges
Nod
es
Node i
Node j
Reordering for minimizing i/a
Data reordering schemes provided by
EdgePack®
13
Armazenamento de matrizes esparsas
Os métodos descritos anteriormente dão origem à matrizes com um
alto grau de esparsidade e consequentemente uma elevada
quantidade de termos nulos;
Existem diversos esquemas destinados ao armazenamento de
matrizes esparsas;
Estes esquemas tiram proveito da esparsidade da matriz e
armazenam somente os termos não-nulos de uma forma compacta;
Mostraremos alguns desses esquemas e algumas implicações
relacionadas à eles;
14
Estrutura de dados: considerações gerais
O uso de uma estrutura de dados apropriada é crucial para se
atingir uma boa performance x espaço de armazenamento;
Núcleos de álgebra linear básicos (p. ex., produtos matriz-vetor)
dependem da estrutura de dados;
Estruturas compactas para o armazenamento de matrizes esparsas
são convenientes somente para métodos iterativos de solução;
Pacotes de solução de problemas de álgebra linear (p. ex. Petsc
Trilinos, etc...) geralmente utilizam estruturas padronizadas e
reconhecidas (p. ex., CSR, BSR).
15
Exemplos de estruturas de dados
DNS: Dense
BND: Linpack Banded;
COO: Coordinate;
CSR: Compressed Sparse Row;
CSC: Compressed Sparse
Column;
MSR: Modified CSR;
ELL: Ellpack-Itpack;
DIA: Diagonal;
BSR: Block Sparse Row
SSK: Symmetric Skyline;
NSK: Nonsymmetric Skyline;
JAD: Jagged Diagonal;
EBE: Element-by-Element;
EDS: Edge-by-Edge
16
1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12.
1 4 1 2 4 1 3 4 5 3 4 5
1 3 6 10 12 13
1 2 3 4 5 6 7 8 9 10 11 12
AA(nnz)=
JA(nnz)=
IA(n+1)= 1 2 3 4 5 6
CSR: Compressed Sparse Row
Armazena a matriz na forma global (“assemblada”)
Muito utilizada em pacotes de álgebra linear;
Difícil de programar;
Favorece a construção de pré-condicionadores;
AA: coeficientes não-nulos de A
JA: Índices das colunas dos não-nulos
IA: número de não-nulos por linha
A
0. 0.
0. 0. 0.
6. 0. 7.
0. 12.
3. 4. 0.
10. 11. 0
1. 0. 0. 2. 0.
.
9.
. 0.
8.
5
17
Produto matriz-vetor CSR (SERIAL)
do i = 1, n
k1 = ia(i)
k2 = ia(i+1)-1
y(i) = dotproduct(a(k1:k2),x(ja(k1:k2)))
End do
Onde n é o número de linhas da matriz e dotproduct é o produto escalar entre a linha da matriz (a(k1:k2)) e o vetor de entrada.
Para realizar a operação dotproduct normalmente utilizamos a rotina ddot da BLAS
1. 0. 0. 2. 0.
3. 4. 0. 5. 0.
6. 0. 7. 8. 9.
0. 0. 10. 11. 0.
0. 0. 0. 0. 12.
A
1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12.
1 4 1 2 4 1 3 4 5 3 4 5
1 3 6 10 12 13
1 2 3 4 5 6 7 8 9 10 11 12
AA(nnz)=
JA(nnz)=
IA(n+1)= 1 2 3 4 5 6
18
BSR: Block Sparse Row
Adequada para armazenar coeficientes construídos aos blocos (normalmente quando se
trabalha com vários graus de liberdade acoplados);
Eventualmente armazenará termos nulos;
Mais complexo de se trabalhar do que o CSR
1. 3. 9. 11. 17. 20.5. 7. 13. 15. 22. 24.2. 4. 10. 12. 18. 21.6. 8. 14. 16. 23. 25.
AA =
1 5 3 5 1 5
1 3 5 7
JA =
IA =
1 2 3 4 5 6
1 2 3 4
A
0. 0.
0. 0.
0. 0. 9. 10. 11. 12.
0. 0. 13. 14. 15. 16.
17. 18. 0. 0. 20. 21.
22. 23. 0. 0
6.
. 24.
5
2. 4.1. 3.
.7
5.
. 8.
2
19
EBE: Element-by-Element
Esquema onde as contribuições de cada elemento são armazenadas
individualmente;
Como as contribuições são armazenadas individualmente a matriz global
não necessita ser montada, ou seja, esse esquema é guarda coeficientes
não-assemblados;
Dificulta o pré-condicionamento;
20
EBE: Element-by-Element
Sendo:nel : número de elementosnnoel: número de nós por elementongl : número de graus de liberdade
t(ngl*nnoel,ngl*nnoel,nel)
Considere a seguinte malha de elementos finitos formada por triângulos (3
nós por elemento) contendo 2 graus de liberdade por nó (ngl):
21
Produto matriz-vetor EBE (SERIAL)
para cada elemento fazer:
Localize os coeficientes no vetor global
Realize a multiplicação
Espalhe o resultado na posição global correspondente
fim
do ie = 1, nel
neq1 = lm(1,ie)
neq2 = lm(2,ie)
neq3 = lm(3,ie)
neq4 = lm(4,ie)
...
retrieve and multiply 16 coefs
...
p(neq1) = p(neq1) + ap1
p(neq2) = p(neq2) + ap2
p(neq3) = p(neq3) + ap3
p(neq4) = p(neq4) + ap4
end do
22
EDS: Edges-by-Edges
Similar à estrutura EBE porém ao invés de armazenarmos os
coeficientes dos elementos, armazenamos somente as
contribuições de suas arestas;
A globalização (“assembling”) das contribuições das arestas só é
realizada no momento da multiplicação matriz-vetor (assim como
ocorre com o esquema EBE);
Dificuldade na construção de pré-condicionadores;
A estrutura de dados EDS oferece uma série de benefícios em
relação às estruturas EBE e CSR
23
Comparativo de estruturas de dados (CSR, EBE, EDS)
10
100
1000
10000
1 2 3 4 5 6ndof
by
tes
EBEEDSCSR 10
100
1000
10000
1 2 3 4 5 6ndof
flo
ps
EBE
EDS
CSR
10
100
1000
1 2 3 4 5 6
ndof
i/a
EBEEDS
CSR
24
EDS: Edges-by-Edge (cont.)
11 12 13 14
21 22 23 24
31 32 33 34
41 42 43 44
e
J J J J
J J J JJ
J J J J
J J J J
1 111 121 121 22
1
0 0
0 0
0 0 0 0
0 0 0 0
e
T T
T TT
2 222 23
2 2 232 33
0 0 0 0
0 0
0 0
0 0 0 0
e
T TT
T T
3 311 13
3 3 331 33
0 0
0 0 0 0
0 0
0 0 0 0
e
T T
TT T
4 411 14
4
4 441 44
0 0
0 0 0 0
0 0 0 0
0 0
e
T T
T
T T
5 522 24
5
5 542 44
0 0 0 0
0 0
0 0 0 0
0 0
e
T TT
T T
6 6 633 346 643 44
0 0 0 0
0 0 0 0
0 0
0 0
e
TT T
T T2 6541 3e eeee ee J T TT TTT
“…desmontando algebricamente um tetraedro…”
Na prática, os bloco diagonais nodais são utilizados como pré-condicionadores do tipo block-Jacobi.
25
EDS: Edge-by-Edge
Compare os esquemas de armazenamento EBE e EDS:
t(ngl*3,ngl*3,nel) t(ngl*2,ngl*2,nedges)
Elemento tetraédrico linear
u-p fully coupled incompressible flow (4 dofs)
nel 5.5nnodes, nedges 7nnodes.
Custo computacional
26
Produto matriz-vetor EDS (SERIAL)
do ie = 1, nel
neq1 = lm(1,ie)
neq2 = lm(2,ie)
neq3 = lm(3,ie)
neq4 = lm(4,ie)
...
retrieve and multiply 16 coefs
...
p(neq1) = p(neq1) + ap1
p(neq2) = p(neq2) + ap2
p(neq3) = p(neq3) + ap3
p(neq4) = p(neq4) + ap4
end do
Elements (1 degree of freedom)
do ie = 1, nedges
neq1 = lm(1,ie)
neq2 = lm(2,ie)
...
retrieve and multiply 4 coefs.
...
p(neq1) = p(neq1) + ap1
p(neq2) = p(neq2) + ap2
end do
Edges
STORAGE REQUIREMENTS:
For 4 degrees of freedom:
Elements: (nel x 192 coefs.) + BDiag
Edges: (nedges x 32 coefs) + BDiag
(1 degree of freedom)
27
Paralelizando códigos de mecânica computacional
Considerações básicas:
– O que deverá ser paralelizado?
– Que modelo de paralelismo devo adotar?
– Quais as implicações do modelo de paralelismo escolhido?
• Arquitetura que o programa rodará
• Estrutura de dados a ser adotada
• Eficiência computacional
• Facilidade de implementação e manutenção
– Adaptarei um programa ou iniciarei o trabalho do zero?
“Botando a mão na massa”
– Por onde eu começo?
– Como me certifico que o programa está rodando corretamente?
28
O que deverá ser paralelizado?
REGRA BÁSICA: O mesmo código deverá rodar em paralelo e serial! “...Um programa serial nada mais é do que um programa paralelo rodando em 1
processador, processo ou tarefa...”
Devemos considerar:
– O programa será integralmente paralelo ou somente algumas partes?
– Quais partes do programa merecem ser paralelizados? “...considere os requisitos
de esforço computacional e volume de dados manipulados...”
– Os algoritmos e métodos escolhidos podem ser eficientemente paralelizados
com o modelo escolhido?
– Leitura e escrita de dados será em paralelo? “...esse questionamento se aplica mais
ao MPI do que ao OpenMP...”
29
Que modelo de paralelismo adotar?
OpenMP: Paralelismo exclusivamente de memória OpenMP: Paralelismo exclusivamente de memória compartilhadacompartilhada
– “fácil” implementação; baixa escalabilidade; restrições de hardware;
dependência de memória (mesh coloring);
MPI: Paralelismo de memória MPI: Paralelismo de memória distribuídadistribuída E/OU E/OU compartilhadacompartilhada
– Implementação “difícil”; alta escalabilidade; pequena restrição de hardware;
particionamento de dados, balanceamento e distribuição de tarefas.
Hibrido: MPI+OpenMPHibrido: MPI+OpenMP
– Implementação dificil, alta escalabilidade e flexibilidade, dados devem tratar
dependência de memória e particionamento ao mesmo tempo.
30
Adaptar ou começar do zero?
Depende:
– Qual o tamanho do programa?
– Qual a complexidade do programa? “...programas que concentram esforço computacional em pequenos
trechos (rotinas) podem tirar proveito do paralelismo através de pequenas intervenções estratégicas sem
modificações drásticas...”
– Está documentado?
– A documentação é realmente satisfatória?
– Conversão será feita pelo mesmo programador que criou a versão original serial?
Considere:
– Vale a pena iniciar do zero?! “...as vezes a adaptação leva mais tempo do que a recriação...”
– Iniciar um programa do zero acrescentando os componentes do antigo;
– Iniciar um programa do zero aproveitando para implementar novas tecnologias. “
O paralelismo deve ser um fator fortemente considerado no projeto de softwares
novos. “...mesmo que a primeira versão não seja paralela o programador deverá vislumbrar que isso será
necessário futuramente e prover mecanismos que facilitem as modificações que serão necessárias...”
31
Por onde eu começo?
Considere a seguinte anatomia básica de programas de mecânica computacional
Pré-processamento
Geração de malha
Condições iniciais
Laço de integração no tempo
Laço de iterações não-lineares
Formação do sistema de equações
Solução do sistema de equações
fim do laço
Fim do laço
Pós-processamento
Saída de dados
Visualização
32
Por onde eu começo?
O que os programas de mecânica computacional normalmente fazem?
– Para códigos de elementos, volumes e diferenças finitas, basicamente utilizamos uma
representação discreta (uma malha) de algo real para resolvermos equações
diferenciais de uma forma aproximada;
– Na parte mais interna desse tipo de programa geralmente teremos a solução de
sistemas de equações lineares;
– Esses sistemas de equações lineares são resolvidos diversas vezes
• CONCLUSÃO: Uma boa forma de obter ganhos reais de processamento seria acelerar a
solução do sistema de equações ou reduzir a quantidade de vezes que ele é resolvido.
– O esforço em resolver esse sistema de equações é diretamente proporcional a:
• Qualidade da discretização número de elementos, arestas, nós
• Física do problema número de graus de liberdade
• Acoplamentos condicionamento do sistema
• Outros fatores localidade de dados
ou seja, encontramos alguns pontos que podem ser explorados.
33
Algumas formas de pensarmos o paralelismo
2
1
, , ,
,
, ,
, , ..
. .
.
, .
...
n
f
f x
f
x
y
x
y
y z t
t
z
z
t
CPU 0
CPU 1
CPU n
Funções
, , , ,x y zvf v pv t
CPU 0
CPU 1
CPU 2
CPU 3
CPU 4
Graus de liberdade
, ...f x y
, ...f x y
Domínio
Note que nos 2 primeiros exemplos o paralelismo fica limitado pelo número de funções ou graus
de liberdade que estão sendo resolvidos;
Sabemos que o tamanho da malha está diretamente relacionada ao tamanho do sistema de
equações;
Uma boa forma de abordar o paralelismo pode ser obtida tratando-se diretamente de nossas
discretizações;
NOTA: Podemos paralelizar somente o esforço computacional (tarefas), somente os dados do
problema (memória) ou ambos simultaneamente. “... buscaremos sempre tratar de ambos...”
Introdução ao OpenMP
35
O que é OpenMP?
OpenMP é uma API para o desenvolvimento de aplicações multi-
tarefa (multithread) em ambientes de memória compartilhada;
Consiste de um conjunto de diretivas do compilador e rotinas da
biblioteca;
Relativamente fácil de criar aplicações multi-tarefa em Fortran, C e
C++
Maiores informações em http://www.openmp.org
36
Modelo de programação OpenMP
Uma tarefa “mestre” cria outras equipes de tarefas quando requisitado através de
diretivas;
O paralelismo é acrescentado incrementalmente ao programa;
O programa serial possui diversos trechos paralelos:
Tarefa mestre(serial)
Tarefa mestre(serial)
Regiões paralelas
37
Programando em OpenMP
O programador insere diretivas OpenMP (comentários em Fortran e
#pragma em C) em pontos específicos do código-fonte;
O compilador interpreta essas diretivas e gera chamadas à biblioteca para
paralelizar essas regiões
Serial:
void main() { double x[1000]; for (int i=0; i<1000; i++) { big_calc(x[i]); }}
Paralelo:
void main() { double x[1000];#pragma omp parallel for for (int i=0; i<1000; i++) { big_calc(x[i]); }}
Divide as iterações do laço entre as tarefas
38
Programando em OpenMP
O número de tarefas pode ser controlado de dentro do programa ou
definindo-se a variável de ambiente OMP_NUM_THREADS
O programador é responsável por sincronizar os dados e resolver
as dependências!!!!
39
Como as tarefas interagem entre si?
Memória compartilhada
– As tarefas se comunicam compartilhando os dados da memória (PERIGO!!!)
O compartilhamento não intencional de dados pode criar condições
de competição entre as tarefas;
– Para evitar tais situações utilize sincronização ou remova TODAS as
dependências de dados
A sincronização de dados é computacionalmente cara!
– Ou seja, REMOVA AS DEPENDÊNCIAS DE DADOS!!!
40
Importância da (IN)dependência de memória
Resolver as dependências de memória são importantes para:
– OpenMP: paralelismo real;
– Auto-paralelismo: paralelismo gerado pelo compilador;
– Vetorização;
– Pipelining;
– OOO: Execução de instruções fora de ordem;
41
Tipos de dependência
Dependência de fluxo:
– Variáveis lidas e escritas em
iterações diferentes
for (j=1; j<MAX; j++) {
A[j]=A[j-1];
}
A[1]=A[0];
A[2]=A[1];
Anti-dependência:
– Variáveis escritas e lidas em
iterações diferentes
for (j=1; j<MAX; j++) {
A[j]=A[j+1];
}
A[1]=A[2];
A[2]=A[3];
42
Tipos de dependência
Dependência externa:
– Variável sobrescrita em iterações
diferentes
for (j=1; j<MAX; j++) {
A[j]=B[j];
A[j+1]=C[j];
}
A[1]=B[1];
A[2]=C[1];
A[2]=B[1];
A[3]=C[1];
Dependência interna:
– Dependência no corpo do laço
k = 1;
for (j=1; j<MAX; j++) {
A[j]=A[j]+1;
B[k]=A[k]+1;
k = k + 2;
}
A[1]=A[1]+1;
B[1]=A[1]+1;
43
Suponha o seguinte laço sendo
executado em 2 threads:
!$OMP PARALLEL DO
do i=1,1000
a(i) = a(i+1)**2
enddo
!$OMP END PARALLEL DO
OpenMP: Como funciona na prática?
Em OpenMP o laço será dividido de tal
maneira que as iterações pares sejam
executadas pela thread 1 e as impares
pela thread 2 (veremos mais tarde
como controlar a divisão dessas
tarefas);
Qual o problema desse laço?
– Note que a thread 1 dependerá de valores
calculados pela thread 2 e vice-versa, ou
seja, as threads serão interdependentes e
competirão por informação (dados).
44
OpenMP: Regras gerais da sintaxe
As construções em OpenMP são em sua maioria diretivas de processamento:
– C e C++:
#pragma omp construct [clausula [clausula]...]
– Fortran
c#omp construct [clausula [clausula]...]
!#omp construct [clausula [clausula]...]
*#omp construct [clausula [clausula]...]
Como são utilizadas diretivas, os compiladores que não suportam o OpenMP
conseguem compilar programas em OpenMP sem maiores problemas (obviamente
que serão produzidos programas seriais...);
O código-fonte sofre somente pequenas modificações pela adição das diretivas
OpenMP.
45
OpenMP: Regras gerais
As diretivas de OpenMP só se aplicam à blocos estruturados, ou
seja, com início e fim bem definidos:
C$OMP PARALLEL
10 wrk(id) = junk(id)
res(id) = wrk(id)**2
if (conv(res)) goto 10
C$OMP END PARALLEL
print *, id
C$OMP PARALLEL
10 wrk(id) = junk(id)
30 res(id) = wrk(id)**2
if (conv(res)) goto 20
goto 10
C$OMP END PARALLEL
if (not_done) goto 30
20 print *, idBloco estruturadoBloco não-estruturado
46
Regiões paralelas
Construções fundamentais que iniciam a execução paralela:
Fortran:
c$omp parallel
c$omp& shared(var1, var2, ...)
c$omp& private(var1, var2, ...)
c$omp& firstprivate(var1, var2, ...)
c$omp& reduction(operador|intrinseco: var1, var2, ...)
c$omp& if(expression)
c$omp& default (private|shared|none)
! algum bloco estruturado
c$omp end parallel
47
Diretivas de OpenMP
shared(var1, var2, ...)
– Variáveis compartilháveis entre as tarefas (utilizam a mesma área de memória)
private(var1, var2, ...)
– Cada tarefa mantém sua própria cópia da variável dentro da região paralela
firstprivate(var1, var2, ...)
– Variáveis privadas que são iniciadas fora da região paralela
lastprivate(var1, var2, ...)
– Variáveis privadas que são salvas ao sair da região paralela
if(expressão)
– Somente paraleliza se a expressão for verdadeira
default(shared|private|none)
– Define o escopo padrão das variáveis no interior de região paralela
schedule(type [,chunk])
– Controla como as iterações dos laços serão distribuídas entre as tarefas
reduction(operator|intrinsic:var1, var2, ...)
– Garante que a operação de redução (p. ex., soma global) seja executada de modo seguro
48
Exemplo de private, default
C$OMP PARALLEL DO SHARED(A)
C$OMP& PRIVATE(MYID,X)
myid=omp_get_thread_num()
x = work(myid)
if (x<1.0) then
a(myid) = x
end if
C$OMP END PARALLEL
Equivale a:
C$OMP PARALLEL DO DEFAULT(PRIVATE)
C$OMP& SHARED(A)
Cada tarefa tem sua própria cópia de x
e myid
A menos que x seja feito privado, seu
valor é indeterminado dentro da região
paralela
Valores de variáveis privadas são
indefinidos no inicio e no fim da região
paralela
A clausula default automaticamente
torna x e myid privados
49
Exemplo de firstprivate
O firstprivate faz com que as
variáveis privadas (local para cada
thread), sejam iniciadas com o último
valor da região serial
program first
integer myid, c
c=98
C$OMP PARALLEL PRIVATE(MYID) &
C$OMP FIRSTPRIVATE(C)
myid = omp_get_thread_num()
print *, ‘T:’,myid,’c=‘,c
C$OMP END PARALLEL
end program
Executando-se 4 threads o código
produzirá:
T:1 c=98
T:3 c=98
T:2 c=98
T:0 c=98
50
A cláusula schedule(type, [chunk])
Controla como o esforço é distribuído entre as threads
chunk é usado para especificar o tamanho de cada parcela de esforço
(número de iterações)
type pode ser uma das seguintes opções:
– static
– dynamic
– guided
– runtime
O argumento chunk é opcional.
51
schedule(static)
As iterações são divididas igualmente entre as tarefas.
C$OMP PARALLEL DO SHARED(X) PRIVATE(I)C$OMP& SCHEDULE(STATIC)do i=1,1000 x(i) = aenddo
thread 3 (i=751,1000)
thread 2 (i=501,750)
thread 1 (i=251,500)
thread 0 (i=1,250)
thread 0 thread 0
52
schedule(static,chunk)
Divide o esforço em parcelas com o tamanho de chunk.
C$OMP PARALLEL DO SHARED(X) PRIVATE(I)C$OMP& SCHEDULE(STATIC,1000)do i=1,12000 x(i) = aenddo
thread 3 (3001,4000),(7001,8000),(11001,12000)
thread 2 (2001,3000),(6001,7000),(10001,11000)
thread 1 (1001,2000),(5001,6000),(9001,10000)
thread 0 (1,1000),(4001,5000),(8001,9000)
thread 0 thread 0
53
schedule(dynamic,chunk)
Divide o esforço em parcelas de tamanho
chunk
Quando a thread finaliza seu trabalho ela
assume a próxima
Valor padrão de chunk é 1
Maior overhead, mas promove um melhor
balanço de carga
C$OMP DO SHARED(X) PRIVATE(I)
C$OMP& SCHEDULE(DYNAMIC,1000)
do i=1,10000
... trabalho ...
end do
thread 0 thread 1
thread 3 thread 2
iteração 1
iteração 2
iteração 3
iteração 4
iteração 5
iteração 6
iteração 7
54
schedule(guided,chunk)
Similar ao agendamento dinâmico porém o tamanho do chunk varia dinamicamente;
O tamanho do chunk depende do número de iterações órfãs (íterações restantes)
Alcança um bom balanço de carga com um overhead relativamente baixo
Garante que nenhuma thread ficará presa fazendo um alto número de trabalho
restante enquanto outras permaneçam ociosas
C$OMP DO SHARED(X) PRIVATE(I)
C$OMP& SCHEDULE(GUIDED,55)
do i=1,12000
...trabalho...
end do
55
reduction(operator|intrinsic:var1[,var2])
Realiza um cálculo global ou uma
comparação segura;
Uma cópia privada de cada variável listada é
criada e iniciada dependendo do operator
ou intrinsic (p.ex., 0 para +)
Somas parciais e demais operações são
determinadas pelas threads em paralelo
As somas parciais são adicionadas para se
obter o valor global
Os mínimos parciais são comparados para
obter o valor global
C$OMP DO SHARED(X) PRIVATE(I)
C$OMP& REDUCTION(+:SUM)
do i=1,n
sum = sum + x(i)
enddo
C$OMP DO SHARED(X) PRIVATE(I)
C$OMP& REDUCTION(MIN:GMIN)
do i=1,n
gmin = min(gmin,x(i))
enddo
56
Operações de reduction
As variáveis listadas devem ser shared e fechadas dentro da região
paralela:
– Em Fortran:
• operator pode ser +, *, -, .and., .or., .eqv., neqv.
• intrinsic pode ser max, min, iand, ior, ieor
– Em C:
• operator pode ser +, *, -, &, ^, |, &&, ||
• ponteiros não são permitidos em operações de redução!
57
schedule(runtime)
A distribuição das tarefas é definida em tempo de execução
Depende do valor da variável de ambiente OMP_SCHEDULE
Esta variável de ambiente é verificada em tempo de execução, e a divisão
de tarefas é consequentemente ajustada;
O método de agendamento padrão é o static
O tamanho do chunk é opcional
Útil para experimentar os diferente métodos de distribuição sem
necessidade de recompilação
setenv OMP_SCHEDULE static,1000
setenv OMP_SCHEDULE dynamic
58
OpenMP em métodos discretos
O OpenMP pode ser facilmente adotado por qualquer código já
existente desde que sejam tomados os devidos cuidados para
solucionar problemas de dependência de memória;
Qualquer laço onde não haja dependência de memória pode ser
paralelizado com o OpenMP;
Aconselha-se a inserção das diretivas de paralelismo OpenMP de
uma maneira incremental ao programa, do laço mais intenso para o
menos intenso, avaliando-se sempre os efeitos no desempenho do
código;
59
Dependência de memória em OpenMP (e vetorização)
Considere o seguinte laço nos elementos da malha sendo executado por 2 threads:
C$OMP PARALLEL DOdo i=1,nel ! recupere os nós do elemento x(no) = x(no) + aenddoC$OMP END PARALLEL DO
Note no exemplo acima que as threads estarão COMPETINDO para modificar o valor da variável x
para o nó 3 pois cada thread estará trabalhando com o seu elemento e os seus dados
correspondentes (nós, vetores, etc...).
CONCLUSÃO: No exemplo acima as threads não poderiam compartilhar a informação referente
ao nó 3
60
Resolvendo a dependência de memória em malhas
Para resolvemos o problema de dependência de memória em malhas, basta separarmos os
elementos da malha em blocos de elementos onde não haja o compartilhamento dos nós;
Esse procedimento chama-se coloração ou blocagem de malha
ielm = 0do icor = 1, ncores nvec = ielblk(icor)C$OMP PARALLEL DOC$OMP& FIRSTPRIVATE(NVEC) do i = ielm+1, ielm+nvec ! recupere os nós do elemento x(no) = x(no) + a enddoC$OMP END PARALLEL DO ielm = ielm+nvecenddo
Utilizamos um algoritmo guloso (greedy) para colorir a malha;
A coloração de malha perturba qualquer ordenação feita previamente realizada (RCM)
61
Produto Matriz vetor EBE ou EDS (OpenMP)
Exemplo para esquema EDS com 1 grau de liberdade
iside = 0 DO iblk = 1, nedblknvec = iedblk(iblk)!DIR$ IVDEP!$OMP PARALLEL DO DO ka = iside+1, iside+nvec, 1 neq1 = lm(1,ka) neq2 = lm(2,ka) ... ! retrieve and multiply 4 coefs. ... p(neq1) = p(neq1) + ap p(neq2) = p(neq2) + ap ENDDO!$OMP END PARALLEL DOiside = iside + nvecENDDO
Diretiva de vetorização. Afirma ao compilador que no bloco a seguir não há nenhuma dependência de dados
62
Produto matriz-vetor CSR (OpenMP)
!$OMP PARALLEL DO
DO i = 1, n
k1 = ia(i)
k2 = ia(i+1)-1
y(i) = dotproduct(a(k1:k2),x(ja(k1:k2)))
ENDDO
!$OMP END PARALLEL DO
Onde n é o número de linhas da matriz e dotproduct é o produto escalar entre a linha da matriz (a(k1:k2)) e o vetor de entrada.
Para realizar a operação dotproduct normalmente utilizamos a rotina ddot da BLAS
1. 0. 0. 2. 0.
3. 4. 0. 5. 0.
6. 0. 7. 8. 9.
0. 0. 10. 11. 0.
0. 0. 0. 0. 12.
A
1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12.
1 4 1 2 4 1 3 4 5 3 4 5
1 3 6 10 12 13
1 2 3 4 5 6 7 8 9 10 11 12
AA(nnz)=
JA(nnz)=
IA(n+1)= 1 2 3 4 5 6
PARALELISMO POR TROCA DE MENSAGENS
Message Passing Interface
64
Aspectos Gerais
Fácil de entender difícil de implementar:
Bastante versátil/portátil pois funciona em sistemas de memória distribuída
e/ou compartilhada;
Idéia básica:
– convertemos um problema grande em vários problemas menores e atribuímos
cada sub-problema a um processador (ou processo);
Várias cópias idênticas do programa são executadas simultaneamente;
Cada cópia atua em sua porção do problema independentemente;
Partes em comum são “sincronizadas” através de operações de trocas
de mensagens utilizando as rotinas do MPI;
Paralelismo baseado na chamada de rotinas da biblioteca MPI para a troca
de informações entre os processos (MPI_BROADCAST, MPI_SEND,
MPI_ALLREDUCE, ...)
65
O que e como paralelizar?
A forma mais intuitiva, e amplamente difundida, de se aplicar o MPI à problemas de
mecânica computacional é através da divisão do domínio em diversos subdomínios
menores e mais fáceis de serem calculados. Desta forma, cada processador pode atuar de
forma colaborativa na solução mais rápida do problema global.
conclusão: temos de aprender como particionar um problema adequadamente
66
Peculiaridades do MPI
Enquanto alguns algoritmos podem ser diretamente implementados
no modelo de paralelismo proposto pelo MPI, outros são
particularmente desafiantes, p. ex.:
– métodos baseados em avanço de fronteira ou algoritmos no qual o método
segue a fronteira de solução são bons exemplos no qual o paralelismo com o
MPI requer cuidados especiais.
– Avalie o seguinte exemplo:
67
Formas de representação de dados:Malhas, Grades, Grafos e Matrizes
1,1 1,2 1,3 1,4 1,5 1,6 1,7
2,1 22 2,3 2,6 2,9 2,10
3,1 3,2 3,3 3,4 3,8 3,9 3,11
4,1 4,3 4,4 4,5 4,8
5,1 5,4 5,5 5,7
6,1 6,2 6,6 6,7 6,10
7,1 7,5 7,6 7,7
8,3 8,4 8,8 8,11
9,2 9,3 9,9 9,10 9,11
10,2 1
a a a a a a a
a a a a a a
a a a a a a a
a a a a a
a a a a
a a a a a
a a a a
a a a a
a a a a a
a a 0,6 10,9 10,10
11,3 11,8 11,9 11,11
a a
a a a a
Matriz esparsa
68
Particionamento de malhas e grafos
Dada a seguinte malha e suas possíveis representações, avaliar o seu
particionamento:
Malha
Grafo nodal
Grafo dual
matriz
69
Partição Dual (por elemento)
MalhaGrafo dual
Note que os nós 3 e 4 pertencem as duas sub-malhas simultaneamente, ou seja, esses nós são nós de interface entre as partições;
Na prática, teremos de manter as informações dos nós de interface devidamente compatibilizadas (iguais ou sincronizadas)
Conclusão: O volume de comunicação entre as partições é diretamente proporcional ao número de nós de interface
Note que os nós 3 e 4 pertencem as duas sub-malhas simultaneamente, ou seja, esses nós são nós de interface entre as partições;
Na prática, teremos de manter as informações dos nós de interface devidamente compatibilizadas (iguais ou sincronizadas)
Conclusão: O volume de comunicação entre as partições é diretamente proporcional ao número de nós de interface
Malha particionada
70
Partição nodal
Alternativamente, poderíamos utilizar o grafo nodal para realizarmos
o particionamento. Na prática, o particionamento dual ou nodal é
escolhido de acordo com o número de partições que se deseja
produzir.
Malha Grafo nodal
71
Particionamento de dados(Aspectos Importantes)
O bom particionamento deverá:
– Minimizar o volume de comunicação mantendo uma carga (esforço)
computacional balanceado;
• Minimizar o volume de comunicação significa minimizar o número de cortes de
arestas do grafo (nodal ou dual);
• Balancear carga significa distribuir uniformemente o número de vértices do grafo
por partição;
– NOTA: sistemas heterogêneos, formados por máquinas com capacidades de
processamento distintas, requerem particionamentos desbalanceados
estrategicamente.
72
Qualidade do Particionamento
POUCA comunicação
Carga BALANCEADA
MUITA comunicação
Carga BALANCEADA
POUCA comunicação
Carga (DES)BALANCEADA
Ruim ou boa?! Depende do sistema...Core I3 Core I7
Core I7 Core I7
73
Particionando uma malha na prática
10.264.863 elementos 1.858.246 nós
Como eu vou partir isso?
Por onde eu começo?
74
Softwares para particionamento de dados
Chaco
– http://www.cs.sandia.gov/~bahendr/chaco.html
Jostle
– http://staffweb.cms.gre.ac.uk/~c.walshaw/jostle/
Metis/ParMetis
– http://glaros.dtc.umn.edu/gkhome/views/metis/
index.html
PARTY
– http://wwwcs.uni-paderborn.de/fachbereich/AG/
monien/RESEARCH/PART/party.html
Scotch
– http://www.labri.fr/perso/pelegrin/scotch/
S-Harp
75
Comparação entre os algoritmos de particionamento
76
Comparativo de particionadoresMetis Jostle Party/DB
Metis Jostle Party/DB
Metis Jostle Party/DB
77
Metis
O que é a Metis?
– Particionador de grafosgrafos e malhasmalhas;
– Reduz largura de banda de matrizesmatrizes esparsas;
– Extremamente rápido;
– Produz partições balanceadas ou não (particionamento ponderado);
– Baseado em métodos de minimização de corte de arestas;
– Desenvolvido em CC com interface para FortranFortran;
– Possui versão em paralelo (ParMetisParMetis);
– Gratuita e de código aberto.Gratuita e de código aberto.
– Maiores informações em: http://www-users.cs.umn.edu/~karypis/metis/
– Pode ser usada como um programa stand-alone ou biblioteca
78
Malha particionada!
METIS
79
Biblioteca Metis
Particionamento de grafos:
– METIS_PartGraphRecursive;
– METIS_PartGraphKway;
– METIS_PartGraphVKway;
– METIS_mCPartGraphRecursive;
– METIS_mCPartGraphKway;
– METIS_WPartGraphRecursive;
– METIS_WPartGraphKway;
– METIS_WPartGraphVKway;
Particionamento de malhas:
– METIS_PartMeshNodal;
– METIS_PartMeshDual;
Reordenação de matrizes
esparsas
– METIS_EdgeND;
– METIS_NodeND;
– METIS_NodeWND
Rotinas auxiliares:
– METIS_MeshToNodal;
– METIS_MeshToDual;
– METIS_EstimateMemory.
80
Usando a Metis(por onde eu começo?)
Baixe os arquivos em:
– http://glaros.dtc.umn.edu/gkhome/metis/metis/download
Metis-4.0.tar.gz (codigo fonte)
Metis-4.0.zip (precompilações para Win32)
81
Usando a Metis stand-alone(Exemplo prático)
10 21 2 6 3 1 8 4 3 1 5 8 6 3 1 8 6 8 6 7 3 5 6 9 8 8 6 11 7 6 10 9 11 8 9 11 6 8 11 9 12
arq.msh
C:\partdmesh.exe arq.msh 2
0010010111
000010011111
arq.msh.epart.2 arq.msh.npart.2
82
Metis para “Fortranianos”(usando a biblioteca)
“linkar” seu programa com a biblioteca
ifort –o <meuprog> *.o libmetis.a
Arquivo libmetis.a
83
FAQ (Frequently Asked Questions)
A Metis particiona malhas híbridas (com mais de um tipo de
elemento)?
– Sim e não: A Metis particiona preferencialmente GRAFOS, ou seja, a Metis
conseguirá particionar a malha desde que seja feita a conversão da mesma
para um grafo nodal ou dual.
– OBS.: A versão paralela da Metis (ParMetis) possui suporte nativo para
particionamento de malhas hibridas.
Como eu descubro quais são os nós de interface?
– A Metis não retorna essa informação explicitamente pois ela é bastante
simples de ser obtida. Toda vez que um nó estiver contido em uma partição
de numeração diferente daquela referente ao elemento que o nó pertence
significa que esse nó é de interface.
84
Metis: Descobrindo os nós de interface
1 1 2 6 3 0 2 1 8 4 3 0 3 1 5 8 6 1 4 3 1 8 6 0 5 8 6 7 3 0 6 5 6 9 8 1 7 8 6 11 7 0 8 6 10 9 11 1 9 8 9 11 6 1 10 8 11 9 12 1
1 02 03 04 05 16 07 08 19 110 111 112 1
nó1 nó2 nó3 nó4 partie nó part
Note que o nó 1 pertence a partição 0 e faz parte do elemento 3 que está na partição 1, ou seja, o nó 1 é um nó de interface!
85
Indo além do particionamento
Após realizarmos o particionamento da malha devemos tentar transformar o
problema original em vários sub-problemas menores e independentes (a
não ser pelos nós de interface);
Uma boa alternativa é renumerar independentemente cada entidade de
malha nas partições criadas, assim teríamos uma forma natural de construir
problemas independentes;
Existem diversas técnicas de se realizar essa reordenação. Aqui trataremos
da forma mais simples;
Note que os dados referentes à interface serão replicados para cada
partição;
Em alguns tratamentos mais criteriosos são criadas camadas de elementos
fantasmas (ghosts) com intuito se de minimizar a quantidade de
comunicação entre as partições (não faremos esse tratamento)
86
Reordenação das partições
Após a reordenação nós teremos como tratar cada partição como se fosse um problema completamente independente e combinar as informações da interface somente quando necessário
87
Tratamento de interface
O que seria melhor:
– trocar mensagens grandes ou muitas pequenas?
88
Uma outra alternativa: Utilização de ghosts
89
Raciocinando em Paralelo
Programas em MPI normalmente possuem um “maestro”: O rank 0
– O rank 0 normalmente controla o fluxo principal do programa
– O rank 0 também deve ser utilizado em processamento – “não desperdice
nenhum processador!!!”
Programas em MPI devem funcionar também em modo serial
(somente 1 processador)
O resultado de uma rodada paralela deve ser igual ao de uma
rodada serial (mas nem sempre é idêntico!)
DICA DE PROGRAMADOR: Utilize macros de pré-processamento
para produzir versões 100% seriais de seu programa (veremos
exemplos)
90
Construindo programas em MPI
Quem faz o que?
Quem lê os dados?
Cada processo lê o seu dado ou o rank 0 lê e distribui?
O que deve ser comunicado/sincronizado?
Quando eu devo comunicar/sincronizar os processos?
Quem imprimirá mensagens na tela?
Quem gravará os arquivos de saída?
Os arquivos de saída devem ser unificados?
“...Depurar em MPI é a arte de saber escrever mensagens na tela em
trechos estratégicos do programa...”
91
program foo
#ifdef MPICODE
include 'mpif.h'
character*60 message
dimension :: istatus(MPI_STATUS_SIZE)
#endif
integer :: nprocs, myrank
irank=0; nprocs=1
#ifdef MPICODE
CALL MPI_INIT(ierr)
CALL MPI_COMM_SIZE(MPI_COMM_WORLD, nprocs, ierr)
CALL MPI_COMM_RANK(MPI_COMM_WORLD, myrank, ierr)
#endif
if (myrank.eq.0) then
! Inicie algo
#ifdef MPICODE
else
! Termine algo iniciado no rank 0
#endif
endif
#ifdef MPICODE
call MPI_BARRIER(MPI_COMM_WORLD,ierr)
call MPI_FINALIZE(ierr)
#endif
stop
end program
program foo
integer :: nprocs, myrank
irank=0; nprocs=1
if (myrank.eq.0) then
! faça algo
endif
stop
end program
program foo
include 'mpif.h'
character*60 message
dimension :: istatus(MPI_STATUS_SIZE)
integer :: nprocs, myrank
irank=0; nprocs=1
CALL MPI_INIT(ierr)
CALL MPI_COMM_SIZE(MPI_COMM_WORLD, nprocs, ierr)
CALL MPI_COMM_RANK(MPI_COMM_WORLD, myrank, ierr)
if (myrank.eq.0) then
! Inicie algo
else
! Termine algo iniciado no rank 0
endif
call MPI_BARRIER(MPI_COMM_WORLD,ierr)
call MPI_FINALIZE(ierr)
stop
end program
Código fonte original
ifort –c –fpp foo.f90
ifort –c –fpp MPICODE foo.f90
Código serial
Código paralelo
Garantindo portabilidade
92
Produto escalar (MPI)
.8
.1 .3
.1 .3
.2
P0
P1
.01
.09
.68
++
.78
2
2
2
.8 .8 .64
.1 .
x0 x0 x0
1 .01
.3 .3 .09
2
2
2
.1 .1 .01
.3 .
.2 .
x1 = x1 x1
3
2 0
9
. 4
.0
MPI_ALLREDUCE
.68
.01
.09
++
.78
2
2
2
2.8 .8
.1 .1
.3
.2 .2
.
x x x
3
= .78
Cálculo serial
Cálculo paralelo
93
Produto escalar (MPI) - Algoritmo
function pddot(n,nit,dx,dy)use mpidefs! n=numero de valores internos, nit=numero de valores de interfaceinteger :: n, nit, ddot1, ddot2, drec, pddotreal*8 :: dx(1), dy(1)
#ifdef MPICODEif (nprocs.ne.1) then ddot1 = ddot( n,dx( 1),1,dy( 1),1) ! BLAS ddot local ddot2 = ddot(nit,dx(n+1),1,dy(n+1),1) ! BLAS ddot da interface CALL MPI_ALLREDUCE( ddot1, drec, 1, MPI_DOUBLE_PRECISION, & MPI_SUM, MPI_COMM_WORLD, ierr ) pddot = drec + ddot2 ! global ddotelse#endif pddot = ddot(n,dx(1),1,dy(1),1)#ifdef MPICODEendif#endifend function
94
Produto Matriz-vetor (MPI)
Matri-Vetor global
Matri-Vetor (MPI)
95
Produto Matriz-Vetor EDSAlgoritmo (MPI)
do ie = 1, nedges
neq1 = lm(1,ie)
neq2 = lm(2,ie)
...
retrieve and multiply 4 coefs.
...
p(neq1) = p(neq1) + ap
p(neq2) = p(neq2) + ap
enddo
#ifdef MPICODE
! Adicionando influencia da interface
call MPI_AllReduce
#endif
96
Paralelismo híbrido: Tratamento dos dados
How do we work with data partitioning and memory dependency?
Particionamento de dados (Metis)
Distributed memory (MPI)
Data partitioning + Mesh coloring and vectorization
Mesh coloring
Vectorization (ivdep) and/or shared memory parallelism (OpenMP)
Pipelines in Itanium-2 (ivdep)
97
Produto matriz-vetor híbrido (OpenMP+MPI)
iside = 0
DO iblk = 1, nedblk
nvec = ia_edblk(iblk)
!dir$ ivdep
!$OMP PARALLEL DO
DO ka = iside+1, iside+nvec, 1
...MATVEC computations...
ENDDO
!$OMP END PARALLEL DO
ENDDO
...over interface nodes...
#ifdef MPICODE
call MPI_AllReduce
#endif
Element-by-Element Edge-by-Edge
98
Outras operações importantes
Devemos tratar qualquer operação que envolva a aglutinação de
valores nodais (entidade mais primitiva e indivisível de uma malha e
que é compartilhada entre as partições)
Exemplo de outras operações importantes:
– Montagem de resíduo;
– Avaliação de pré-condicionadores bloco-diagonal nodal.
99
Pré-condicionamento Bloco Diagonal
Os blocos diagonais nodais são naturalmente obtidos em esquemas de
armazenamento EBE e EDS e podem ser utilizados como pré-condicionador (serial
ou paralelo).
A estrutura da matriz bloco diagonal nodal pode ser representada por:
Na prática os blocos diagonais (geralmente matrizes de dimensão ngl x ngl) são
invertidos, bloco-a-bloco e pré-multiplicados pelo vetor RHS antes da chamada do
método iterativo;
No interior do método iterativo a multiplicação do lado esquerdo é realizada logo
após o produto matriz-vetor.
W =Ax = b
-1 -1W Ax = W b
100
Pré-condicionamento Bloco Diagonal (cont.)
subroutine LinearSolution (A, W, x, b)
! Resolve W-1 A x = W-1 b
#ifdef MPICODE
! Globalize blocos de interface
#endif
do i=1,nblocos
! inverta Wi
enddo
! Multiplique W-1 b
call GMRES(A, W, x, b)
! Produto matriz-vetor z A x
! Multiplique W-1 z
end solution