Embed Size (px)
Citation preview
Pontifícia Universidade Católica do Rio Grande do Sul
Faculdade de Informática
Pós-Graduação em Ciência da Computação
Estudo sobre a aplicação da Computação Paralela na
resolução de sistemas lineares
Aluno: Cleber Roberto Milani
Orientador: Prof. Dr. Luiz Gustavo Leão Fernandes
Introdução à Pesquisa I
Porto Alegre, junho de 2008
Conteúdo
LISTA DE FIGURAS iv
LISTA DE SÍMBOLOS E ABREVIATURAS v
Capítulo 1: Introdução 1
1.1 Objetivos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2 Organização do trabalho . . . . . . . . . . . . . . . . . . . . . . . . . . .. . 3
Capítulo 2: Álgebra Linear e Equações Diferenciais 4
2.1 Equações Diferenciais . . . . . . . . . . . . . . . . . . . . . . . . . . . .. . 5
2.2 Sistemas de Equações Lineares Algébricas . . . . . . . . . . . .. . . . . . . . 6
2.3 Solução de Sistemas de Equações Lineares . . . . . . . . . . . . .. . . . . . . 8
2.3.1 Métodos Diretos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.3.2 Métodos Iterativos . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.4 Aplicações práticas de modelagem baseada em EDs e SELAs .. . . . . . . . . 15
2.4.1 Método dos Elementos Finitos aplicado à Odontologia .. . . . . . . . 16
2.4.2 EDs e Método dos Elementos Finitos aplicados à Engenharia . . . . . . 16
2.4.3 EDs aplicadas à Física e Matemática . . . . . . . . . . . . . . . .. . . 17
2.5 Conclusão . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
i
Capítulo 3: Computação paralela 19
3.1 Contextualização Geral . . . . . . . . . . . . . . . . . . . . . . . . . . .. . . 20
3.2 Modelagem de programas paralelos . . . . . . . . . . . . . . . . . . .. . . . 20
3.3 Medidas em Computação Paralela . . . . . . . . . . . . . . . . . . . . .. . . 24
3.3.1 Complexidade de tempo em algoritmos especificados porum DAG . . 25
3.3.2 Encontrando um DAG ótimo . . . . . . . . . . . . . . . . . . . . . . . 26
3.3.3 SpeedUPe Eficiência . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.3.4 Granulosidade e Escalabilidade . . . . . . . . . . . . . . . . . .. . . 29
3.4 Classificação das arquiteturas paralelas . . . . . . . . . . . .. . . . . . . . . . 30
3.4.1 Classificações Acadêmicas . . . . . . . . . . . . . . . . . . . . . . .. 30
3.4.2 Classificações Comerciais . . . . . . . . . . . . . . . . . . . . . . .. 34
3.5 Paradigmas de programação paralela . . . . . . . . . . . . . . . . .. . . . . . 37
3.5.1 Task-Farming (ou Master/Slave) . . . . . . . . . . . . . . . . . .. . . 37
3.5.2 Single Program Multiple Data (SPMD) . . . . . . . . . . . . . . .. . 38
3.5.3 Data Pipelining . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.5.4 Divide and Conquer . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.5.5 Speculative Parallelism . . . . . . . . . . . . . . . . . . . . . . . .. . 39
3.6 Conclusão . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
Capítulo 4: Computação Paralela aplicada os Métodos Numéricos 41
4.1 Paralelização de algoritmos numéricos simples . . . . . . .. . . . . . . . . . 42
4.1.1 Adição de escalares . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
4.1.2 Produto Interno . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
4.1.3 Adição e Multiplicação de Matrizes . . . . . . . . . . . . . . . .. . . 44
4.1.4 Potência de uma Matriz . . . . . . . . . . . . . . . . . . . . . . . . . 45
4.2 Paralelização dos Métodos Iterativos para resolução deSELAs . . . . . . . . . 46
4.2.1 Paralelização do Método de Jacobi . . . . . . . . . . . . . . . . .. . . 46
4.2.2 Paralelização de Gauss-Seidel . . . . . . . . . . . . . . . . . . .. . . 49
4.3 Implementação da paralelização dos métodos iterativosclássicos . . . . . . . . 54
ii
4.4 Conclusão . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
Capítulo 5: Bibliotecas e pacotes para resolução numérica de SELAs com alto
desempenho 57
5.1 BLAS – Basic Linear Algebra Subprograms . . . . . . . . . . . . . .. . . . . 57
5.2 LINPACK . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
5.3 LAPACK – Linear Algebra PACKage . . . . . . . . . . . . . . . . . . . . .. 61
5.4 scaLAPACK - Scalable LAPACK . . . . . . . . . . . . . . . . . . . . . . . .62
5.5 HPL - High Performance Linpack . . . . . . . . . . . . . . . . . . . . . .. . 63
Capítulo 6: Conclusão 64
BIBLIOGRAFIA 65
iii
Lista de Figuras
3.1 DAG de um algoritmo para avaliação da expressão(x1 + x2)(x2 + x3). Fonte:
[8]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.2 DAG de outro algoritmo para avaliação da expressão(x1+x2)(x2+x3). Fonte:
[8]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.3 Classificação segundo o compartilhamento de memória.Fonte: [Hwa98] apud
[19]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
4.1 DAG para soma de 16 escalares utilizando 8 processadores. Fonte: [8]. . . . . 43
4.2 DAG para soma de 16 escalares com 4 processadores.Fonte: [8]. . . . . . . . 44
4.3 Computação paralela da potência de uma matriz n x n.Fonte: [8]. . . . . . . . 46
4.4 Em (a) Grafo de dependências e (b) DAG para função com o grafo de dependên-
cias definido por (a).Fonte: Adaptação de [8]. . . . . . . . . . . . . . . . . . 48
4.5 DAG para possíveis iterações de Gauss-SeidelFonte: Adaptação de [8]. . . . . 50
4.6 Grafo para iterações Gauss-Seidel com coloração.Fonte: Adaptação de [8]. . . 54
iv
Lista de Símbolos e Abreviaturas
ED Equação Diferencial 4
SELA Sistema de Equações Lineares Algébricas 4
EDO Equação Diferencial Ordinária 4
EDP Equação Diferencial Parcial 4
MDF Método das Diferenças Finitas 6
MEF Método dos Elementos Finitos 6
SOR Sucessive Overrelaxation 14
ED Equação Diferencial 16
SISD Single Instruction Stream Single Data Stream 30
SIMD Single Instruction Stream Multiple Data Streams 31
MISD Multiple Instruction Streams Single Data Stream 31
MIMD Multiple Instruction Streams Multiple Data Streams 31
PVM Parallel Virtual Machine 32
MPI Message Passing Interface 32
UMA Uniform Memory Access 32
NUMA Non-Uniform Memory Access 33
NORMA Non-Remote Memory Access 33
RPC Remote Procedure Call 37
SPMD Single Program Multiple Data 37
v
SPMD Single Program Multiple Data 38
BLAS Basic Linear Algebra Subprograms 57
vi
Resumo
Diversos problemas das mais variadas áreas podem ser modelados através de Equações
Diferenciais (EDs) e Sistemas de Equações Lineares. Uma equação diferencial é definida em
um conjunto infinito de pontos (espaço contínuo), entretanto, para que ela possa ser resolvida
computacionalmente, é necessário discretizá-la (criar umespaço discreto e finito), processo esse
que implica na solução de sistemas de equações lineares. Um sistema linear é tipicamente com-
posto por "n"equações com "n"variáveis, sendo que quando empregados para discretização de
EDs ou modelagem de problemas, os sistemas resultantes costumam apresentar um número el-
evado de variáveis, o que torna seu processo de solução bastante complexo. Nesse contexto, a
computação de alto desempenho surge como ferramenta essencial para viabilizar a execução em
tempo viável dos algoritmos que resolvem tais sistemas. Neste trabalho é feita, inicialmente,
uma revisão dos conceitos da Álgebra Linear envolvidos na solução de sistemas lineares. Na
seqüência é apresentada uma revisão dos diversos aspectos da computação paralela, desde as
arquiteturas até os paradigmas de programação existentes.Feita essa contextualização discute-
se a aplicação da computação paralela nos métodos numéricose, ao final, são apresentadas as
bibliotecas e pacotes desoftwaremundialmente mais utilizados e que representam o estado da
arte na solução de sistemas lineares por métodos numéricos.
1 Introdução
Nos últimos anos tem-se observado que os processadores disponíveis no mercado atingiram
freqüências declock bastante próximas do limite suportado pelas atuais tecnologias, o que tem
dificultado crescimentos maiores de sua velocidade. Os fatores que levam a esse cenário de
dificuldade são diversos, destacando-se, dentre eles, o aumento do consumo de energia e seus
respectivos efeitos colaterais, como dissipação térmica efenômenos de capacitância e indutância
parasitas, por exemplo. Tal situação faz com que o aumento dafreqüência dosclocksseja um es-
forço com custos muito elevados e ganhos, proporcionalmente, muito baixos [20]. Além disso,
é importante observar que o crescimento da freqüência doclock, isoladamente, não proporciona
aumento considerável no desempenho do computador. Outros aspectos, como a eficiência no
acesso à memória, devem acompanhar o crescimento doclock. O comportamento conhecido
comoGargalo de von Neumanndefine que o poder de processamento disponibilizado é limi-
tado em função da taxa de transferência de dados e instruçõesentre a memória e o processador.
Assim, o uso de paralelismo nas arquiteturas é uma estratégia para superar essas limitações
[34]. Nesse contexto, a solução encontrada pelas fabricantes tem sido empacotar diversos pro-
cessadores em um único processador, criando os chamados processadoresmulti-core[20]. Esse
processo deu início à popularização das arquiteturas paralelas entre os usuários domésticos.
Em complemento a isso, sabe-se que os fenômenos naturais sãoinerentemente paralelos,
logo, nada mais natural que expressar as computações pertinentes aos mesmos de forma par-
alela. Em algumas situações a ordem de execução é importantepara o melhor entendimento do
problema real, mesmo que em outras ela seja irrelevante [32]. Todavia, na computação seqüen-
cial não é possível trabalhar com essa flexibilidade, pois faz-se sempre necessário definir uma
1
ordem para a execução das ações. De fato, o nicho original da computação paralela é a com-
putação científica, dado que, além das características de paralelismo inerente, o estudo desses
fenômenos exige grande poder de processamento para resolução dos mais diversos algoritmos
numéricos. A necessidade original da computação de alto desempenho na computação numérica
veio de uma série de contextos envolvendo equações diferenciais parciais, como dinâmicas de
fluídos, previsão do tempo, processamento de imagens entre outros. Nessas aplicações existe
uma quantidade massiva cálculos a serem realizados, os quais podem ser facilmente decompos-
tos, tornando-os, portanto, candidatos naturais à paralelização [8].
Nos últimos anos tem crescido também o interesse em outros tipos de computação de larga
escala, como, por exemplo, análise, simulação e otimizaçãode sistemas de larga escala inter-
conectados, simuladores de filas, dentre outros. Outra categoria de problemas abrange a solução
de sistemas de equações gerais, programação matemática e problemas de otimização. Em alguns
casos desses problemas, eles podem ser decompostos, porém,a quantidade de subtarefas obti-
das na decomposição tende a ser menor e as tarefas mais complexas do que aquelas obtidas no
contexto das equações diferenciais parciais. Conseqüentemente, tende-se a utilizar um número
menor de processadores mais poderosos, coordenados por um mecanismo de controle mais com-
plexo. Em ambas as classes de aplicação, a preocupação principal é a mesma, a relação entre
custo e velocidade. Ohardwarenão deve ser caro a ponto de se tornar proibitivo, entretantoa
computação deve terminar em um tempo que seja aceitável parauma determinada aplicação [8].
Cabe observar que, mesmo com o crescimento da velocidade dosprocessadores, a computação
paralela possibilita muitas vezes uma relação custo/benefício melhor que a obtida ao utilizar
equipamentos com processadores de última geração. Isso porque, agrupar em um equipamento
paralelo processadores mais antigos prove uma alternativacomputacional com custo menor e su-
porte natural a técnicas de tolerância a falhas (característica inerente à redundância dehardware)
[34].
Atualmente, decorridos diversos anos de ampla utilização dos sistemas de computação par-
alela, pode-se ter uma melhor compreensão do potencial e limitações da área. Em muitos casos,
a computação paralela representa as melhores esperanças deresolver problemas maiores do que
aqueles que se consegue resolver nos dias de hoje [8]. Nesse contexto, este trabalho propõe um
2
estudo sobre o paralelismo aplicado aos algoritmos numéricos para resolução de sistemas de
equações lineares, problema esse que está presente na grande maioria dos cenários envolvendo
computação em larga escala.
1.1 Objetivos
Os objetivos desse trabalho consistem em: revisar os conceitos de álgebra linear e com-
putação numérica para resolução de sistemas lineares; revisar tópicos em geral da computação
paralela; relacionar a aplicação da computação paralela com sua aplicação para solução de prob-
lemas da computação numérica; investigar pacotes e bibliotecas largamente utilizados para res-
olução de sistemas lineares com computação paralela.
1.2 Organização do trabalho
Este trabalho está dividido em seis capítulos, sendo o primeiro deles a presente introdução.
O Capítulo 2 apresenta uma revisão sobre equações diferenciais e álgebra linear, focando, prin-
cipalmente, na solução de sistemas de equações lineares. NoCapítulo 3 é feita uma contextu-
alização geral sobre a computação paralela, apresentando modelos dehardwaree softwarepara
sua exploração e algumas medidas que auxiliam na avaliação do paralelismo.
O Capítulo 4 faz a relação entre os dois anteriores, apresentando a aplicação dos conceitos
da computação paralela para otimização dos algoritmos numéricos. No Capítulo 5 é feito um
estudo de caso de bibliotecas e aplicações que utilizam a computação de alto desempenho para
resolução de sistemas lineares. Por fim, o Capítulo 6 apresenta as conclusões resultantes do
trabalho.
3
2 Álgebra Linear e Equações
Diferenciais
A Computação Científica sempre teve como uma de suas grandes aplicações possibilitar,
através de técnicas de modelagem e simulação, um melhor entendimento dos problemas e fenô-
menos físicos. Tais fenômenos dependem de variáveis contínuas, ou seja, são funções dessas
variáveis, e podem ser representados por Equações Diferenciais (ED). A variável contínua é
independente e geralmente corresponde a parâmetros como tempo, velocidade, distância, entre
outros. Já a função é uma variável dependente e deve descrever da maneira mais fiel possível o
fenômeno que se propunha reproduzir. Essa representação matemática recebe o nome de mode-
lagem [11].
As equações diferenciais são definidas em um conjunto infinito de pontos (espaço contínuo),
entretanto, para que possam ser resolvidas computacionalmente, é necessário discretizá-las (criar
um espaço discreto e finito), processo esse que implica na solução de sistemas de equações
lineares algébricas (SELAs) [17, 12]. Esses sistemas, alémde serem utilizados para discretizar
as EDs, também podem ser utilizados para modelagem de problemas variados. Um sistema
linear é tipicamente composto porn equações comn variáveis, sendo que na modelagem de
problemas reais o número de variáveis costuma ser bastante elevado, resultando em sistemas
bastante complexos.
Este capítulo apresenta inicialmente conceitos básicos sobre as Equações Diferenciais Or-
dinárias (EDOs) e Equações Diferenciais Parciais (EDPs), que são, em geral, aquelas que mod-
elam os fenômenos ou problemas. Na seqüência, é feita uma revisão sobre os sistemas de
4
equações lineares e alguns métodos para resolução dos mesmos. Ao final, são apresentados
alguns casos que exemplificam o uso da modelagem e simulação computacional em problemas
de áreas variadas do conhecimento.
2.1 Equações Diferenciais
Uma Equação Diferencial é aquela cujas incógnitas são funções que aparecem na equação
sob a forma das respectivas derivadas. As equações diferenciais podem ser de dois tipos [11, 17]:
• Equações Diferenciais Ordinárias: envolvem uma função desconhecida e algumas de suas
derivadas. As incógnitas dependem de apenas uma variável independente e a equação
contém apenas derivadas totais. Podem ser escritas na forma:
F (x, y(x), y′(x), y′′(x), y3(x), ..., yn−1(x), yn(x)) = 0, n ≥ 1 (2.1)
• Equações Diferenciais Parciais: envolvem funções com maisde uma variável indepen-
dente, ou seja, estabelecem relações entre uma variável dependente de duas ou mais var-
iáveis independentes. Nesse caso as derivadas são parciais. Uma equação que descreve
uma variável dependenteu em termos de duas variáveisx e t é escrita na forma:
(∂u
∂x,∂u
∂x, u,D) = A
∂2u
∂x2+ B
∂2u
∂x∂t+ C
∂2u
∂t2(2.2)
ondeA,B,C e D podem ser constantes ou funções das variáveis independentes. Uma
derivada parcial de uma função com dois ou mais argumentos é asua derivada em relação
a uma daquelas variáveis, mantendo as demais constantes. Considerando como exemplo
a fórmula para cálculo do volume de um cone [11], tem-se que ela depende da alturah e
do raior:
V =r2hπ
3(2.3)
a derivada parcial deV com relação ar define a taxa com que o volume do cone cresce
conforme seu raio aumenta mantendo-se a altura constante. Pode ser escrita como:
5
∂V
∂r=
2rhπ
3(2.4)
A solução de uma ED é uma função que, além de satisfazer a equação dada, não contém
derivadas nem diferenciais e pode existir ou não e, ainda, caso exista, pode ser única ou não.
Existem três tipos de soluções para EDs: a geral (solução quecontém tantas constantes arbi-
trárias quantas forem as unidades da ordem de integração), aparticular (solução deduzida da
solução geral atribuindo-se valores particulares às constantes) e a solução singular (solução não
deduzida da solução geral e que só existe em alguns casos). A equação 2.1 possui uma família de
soluções, para especificar uma das curvas que formam a família de soluções precisamos impor
condições adicionais na função Y, as quais são chamadas condições iniciais. O problema com
as condições iniciais é chamado de problema de valor inicialou problema de condições iniciais.
Além dos problemas de valores iniciais, podemos ter problemas com condições de contorno,
isto é, além da condição no início do fenômeno, temos também uma condição a atingir no fim
do fenômeno [17].
A solução analítica de equações diferenciais é, em geral, impraticável, o que leva à grande
utilização de métodos numéricos que, por sua vez, substituem a solução exata analítica por uma
solução aproximada. Porém, para que as EDs possam ser resolvidas computacionalmente, faz-se
necessário discretizá-las, uma vez que elas são definidas emum conjunto infinito de pontos. Ex-
istem vários métodos que podem ser utilizados para esse propósito, dentre os quais se destacam
o Método das Diferenças Finitas (MDF) e Método dos ElementosFinitos (MEF) , cada um ap-
resenta vantagens e desvantagens, cabendo ao usuário definir qual será utilizado. Entretanto, os
métodos têm em comum o fato de que implicam na resolução de sistemas de equações algébricas
lineares [12, 18, 30], que é o foco do presente trabalho.
2.2 Sistemas de Equações Lineares Algébricas
Diversos problemas da Matemática Numérica são modelados emtermos de um Sistema de
Equações Lineares Algébricas. Isso vale em geral para o tratamento numérico de equações
funcionais lineares que ocorrem, entre outras, como equações diferenciais parciais ou ordinárias
6
e equações integrais. Um sistema linear é tipicamente composto porn equações comn variáveis
e pode ser escrito da seguinte forma [17]:
A~x = ~y (2.5)
onde A é uma matriz formada pelos coeficientes,~x é o vetor das incógnitas do sistema, as
quais deverão ser encontradas, e~y é o vetor dos termos independentes. A teoria sobre SELAs é
oriunda da álgebra linear e não será aprofundada neste trabalho, uma vez que o foco do mesmo é
o tratamento numérico que a resolução desses sistemas requer. Dado que a solução numérica dos
SELAs baseia-se em conceitos de matrizes, algumas definições básicas, as quais serão usadas
no decorrer do trabalho, são rapidamente revisadas a seguir[17, 33]:
1. Uma matriz que possui o mesmo número de linhasm e colunasn é dita matriz quadrada;
2. Uma matriz em que todos os elementos que não compõem a diagonal principal valem
zero (∀(i, j), i 6= j → aij = 0) e os elementos da diagonal principal têm valor diferente
dezero (∀(i, j), i = j → aij 6= 0) é dita matriz diagonal;
3. Uma matriz onde os elementos da diagonal principal e os acima dela são diferentes de
zero (∀(i, j), i < j → aij 6= 0) e aqueles abaixo da diagonal principal são iguais azero
(∀(i, j), i > j → aij = 0) é dita matriz triangular superior;
4. Uma matriz onde os elementos da diagonal principal e os abaixo dela são diferentes de
zero (∀(i, j), i > j → aij 6= 0) e aqueles acima da diagonal principal são iguais azero
(∀(i, j), i < j → aij = 0) é dita matriz triangular inferior;
5. Se a diagonal principal tiver seus elementos nulos(∀(i, j), i = j → aij = 0), as matrizes
definidas nos itens3 e4 são ditas estritamente superior e estritamente inferior, respectiva-
mente;
6. As matrizes quadradas têm um determinante definido por umafórmula recursiva:det(A)
ou |A| = a11M11 − a12M12 + a12M13 . . . ondeM1k é o determinante da matriz obtida
deA1 eliminando sua linha1 e colunak;
7
7. Uma matriz cujo determinante é diferente de0 é dita matriz regular ou não singular;
8. Chama-se tranposta deA à matrizAT obtida convertendo-se as linhas deA em colunas e
as colunas em linhas;
9. SeA = AT , entãoA é dita matriz simétrica;
10. Uma matriz que possui aproximadamente 80% ou mais de seuselementos nulos é dita
matriz esparsa. Esse tipo de matriz ocorre freqüentemente nos problemas práticos, princi-
palmente naqueles modelados por EDPs;
11. Uma matriz na qual todos os elementos não nulos concentram-se na vizinha da diagonal
principal é dita matriz banda;
2.3 Solução de Sistemas de Equações Lineares
Os métodos para resolver sistemas de equações lineares são divididos em duas classes: méto-
dos diretos e indiretos. Um método é dito direto quando, através da aplicação de um número
finito de operações aritméticas, chega-se a uma solução exata, ou seja, com precisão infinita. Já
os métodos indiretos, mais conhecidos como iterativos, baseiam-se na repetição de uma seqüên-
cia finita de repetição de operações que objetivam aproximara solução.
A solução de sistemas lineares em uma máquina apresenta errotanto ao utilizar os métodos
diretos como utilizando os iterativos. Nos algoritmos diretos a solução apresenta erros causados
pelo arredondamento, o erro total é composto pela soma do erro de entrada, ou seja, aquele
resultante da representação do SELA no sistema de ponto flutuante da máquina, com o erro de
aritmética, também devido à representação no sistema de ponto flutuante. No caso dos sistemas
iterativos, além desses dois erros, soma-se ainda o erro de discretização que advém do fato de
truncarmos (cortarmos) a seqüência infinita de repetições após um número finito de iterações
[17].
8
2.3.1 Métodos Diretos
Os métodos diretos são aqueles que após a aplicação de um número finito de operações arit-
méticas elementares obtêm uma solução exata. Existem diversos métodos desse tipo, e em geral
eles aplicam-se mais eficientemente a matrizes (ou sistemaslineares) densas. Serão abordados
apenas apenas os métodos de Gauss, com algumas alterações domesmo, e a Fatoração LU. Os
algoritmos que implementam computacionalmente esses métodos não serão apresentados neste
trabalho e podem ser encontrados em [17] e [7].
Eliminação de Gauss
É o método direto mais utilizado para resolver SELAs densos de pequeno (matrizes de ordem
até 30) e médio (ordem até 50) porte. Para os sistemas densos de grande porte, bem como para
os sistemas esparsos de qualquer porte, existem outras técnicas mais eficientes. O algoritmo
básico de Gauss, em essência, nada mais é do que a aplicação esquemática das propriedades
fundamentais da álgebra linear. Pode ser dividido em duas etapas: triangularização da matriz e
retrosubstituição das variáveis [17].
A triangularização consiste em transformar a matrizA em uma matriz triangular superior
U aplicando para isso as seguintes operações: adição de uma linha com um múltiplo de outra
linha e substituição de uma dessas linhas pelo resultado da adição; multiplicação de uma linha
por uma constante; troca de linhas. A cada passo o elemento dediagonal da primeira linha
da matriz-resto é chamado pivô. O segundo passo, a retrosubstituição, consiste no cálculo do
vetor~x, solução deA~x = ~y através da substituição do último componente de~x nas equações
anteriores [17, 33].
Eliminação de Gauss com Pivotamento
Dependendo do sistema, o algoritmo básico de Gauss pode apresentar problemas no que diz
respeito aos resultados. Esses problemas consistem em erros causados pelos erros de arredonda-
mento, sendo que, nos casos em que o elemento pivô for muito pequeno o erro aumenta. Visando
a minimizá-lo, surge a idéia do algoritmo de Gauss com pivotamento, o qual nada mais é do que
9
o mesmo algoritmo anterior, porém com uma troca sistemáticadas linhas [17].
A idéia do algoritmo com pivotamento é buscar sempre colocarna posição de pivô o el-
emento com maior valor absoluto. Para isso, deve-se freqüentemente trocar linhas e colunas
não somente quando o pivô tem valor zero, mas também quando seu valor for muito pequeno.
Quando a escolha do pivô considera apenas determinada coluna, diz-se que o algoritmo é de
Pivotamento Parcial, e, nos casos em que escolha do pivô considera toda a matriz resto, diz-se
que o algoritmo é de Pivotamento Total. Cabe observar, porém, que a melhora de precisão pro-
porcionada por este método costuma não ser vantajosa quandoconsiderada em conjunto com
o aumento que o método introduz na complexidade do algoritmo[17], ficando a cargo de que
desenvolve a aplicação decidir o que é mais vantajoso em cadacaso.
Condicionamento de uma matriz
Um problema é dito mal condicionado quando pequenas alterações nos valores de entrada
acarretam em grandes erros no resultado final. No caso de SELAs de ordem 2 esse condiciona-
mento pode ser visualizado graficamente de maneira simples,entretanto, em sistemas de ordem
maior o procedimento de traçar e analisar o gráfico é muito difícil. Assim, é necessário um meio
de medir o condicionamento [17]. O condicionamento é dado por:
cond(A) = ‖A‖.‖A−1‖ (2.6)
onde‖A‖ é a norma da matrizA. Quanto maior o valor decond(A), mais sensível é o sis-
tema. O condicionamento é uma medida interessante porque métodos como o de Gauss, mesmo
com pivotamento, não produzem nenhuma estimativa sobre a exatidão da solução quando ex-
ecutados em uma máquina. Com o condicionamento pode-se determinar o quão confiável é a
solução do SELA. O condicionamento é um número do qual interessa a ordem de grandeza e,
além disso, está relacionado com a precisão de máquina usada. Second(A) = 105 mas a pre-
cisão decimal da máquina for de 17 casas, então não há problema em se perder os últimos cinco
dígitos. Porém, se a precisão decimal for de apenas 6 casas decimais, então o resultado será
confiável apenas até a primeira casa decimal [17].
10
Existem alguns métodos que são particularmente úteis no caso de matrizes mal condi-
cionadas, porém não serão abordados neste trabalho. Um exemplo é o Método de Kaczmarz, o
qual pode ser encontrado em [9].
Método de Gauss-Jordan
É um método direto complementar ao método de Gauss. Consisteem transformar a matriz
A em uma matriz diagonal, o que é feito em 3 passos: transformação deA~y em uma matriz tri-
angular superior; transformar o resultado em uma matriz triangular inferior; solução do sistema
diagonal resultante. Esse método requer um número maior de operações que a eliminação de
Gauss sendo, portanto, menos eficiente para resolver sistemas do tipoA~x = ~y [17].
Decomposição LU
Também conhecido como Fatoração LU, esse método consiste emtransformar a matrizA
em um produto de duas outras matrizesL eU , ondeL é a matriz triangular inferior eU a matriz
triangular superior. É um dos métodos mais empregados para resolução de SELAs, principal-
mente para sistemas densos de grande porte. Algumas vantagens dos processos de fatoração
é que eles permitem resolver qualquer sistema que tenhaA como matriz de coeficientes e no
caso de o vetor~y ser alterado, a resolução do novo sistema é quase imediata. Em princípio
seria necessário realizar cálculos com matrizes inversas para encontrar as matrizesL eU , entre-
tanto, na prática, pode-se calcularL e U com a simples aplicação da definiçâo de produto e de
igualdade de matrizes, ou seja, impondo queLU = A.
Para resolução de SELAsA~x = ~y de ordemn, seA satisfaz as condições para fatoração
LU , então o sistema pode ser escrito na formaLU~x = ~y. Faz-se entãoU~x = ~b e a equação fica
reduzida aL~b = ~y. Resolve-se o sistema triangular inferiorL~b = ~y para obter o vetor~b. Na
seqüência~b é substituído emU~x = ~b e obtém-se o sistema triangular superior cuja solução é o
vetor~x. Logo, a solução de SELAs através de fatoração LU consiste nasolução de dois sistemas
triangulares [7].
11
2.3.2 Métodos Iterativos
Os métodos iterativos consistem em arbitrar valores iniciais para~x0 e, a partir deles, repetir
uma seqüência de passos em que, a cada passo, a aproximação inicial é modificada para se
aproximar da solução. Essas modificações são chamadas passos de relaxação. O processo é
interrompido quando o vetor de soluções chegar a uma vizinhança suficientemente próxima da
solução exata, o que é definido através de um parâmetro chamado critério de parada. Assim, o
número de iterações que serão executadas é desconhecido e a solução do sistema será dada pelo
vetor encontrado na última iteração [30, 17]. Os métodos iterativos são mais indicados do que
os diretos para resolver sistemas esparsos, pois os algoritmos diretos não levam em conta uma
propriedade importante desses sistemas que é o fato de que muitas operações não precisam ser
feitas ao resolvê-los, fator esse que é considerado nos métodos iterativos [17].
A idéia mais simples de um algoritmo iterativo para resolverA~x = ~y é obtida, resolvendo-se
a primeira equação parax1, a segunda parax2 e assim sucessivamente. Generalizando temos
[17]:
xi =1
aii(yi − ai1x1 − ai2x2 − · · · − ai,i−1xi−1 − ai,i+1xi+1 − · · · − ainxn) (2.7)
onde são conhecidos os valoresx1, x2, . . . , xi−1, xi+1, . . . , xn usados no lado direito.
Assim, os métodos iterativos podem ser facilmente deduzidos dividindo-se a matrizA da
equaçãoA~x = ~y em uma soma de duas matrizes mais simples. SeA = S − T entãoA~x = ~y
será o mesmo queS~x = T~x + ~y. Dessa forma podemos experimentar a seguinte iteração [33]:
S~xk+1 = T~xk + ~y (2.8)
Naturalmente apenas isso não garante que o método seja adequado. Para o resultado ser
correto, a divisão da matrizA deve satisfazer dois requisitos: o novo vetor~xk+1 deve ser fácil de
computar, dessa forma,S deve ser uma matriz simples e invertível e deve também ser diagonal
ou triangular; a seqüência~xk deve convergir para a solução verdadeira~x. A matrizS é também
chamada de pré-condicionador e sua escolha é crucial na análise numérica, pois dessa escolha
dependem os critérios de convergência e a eficiência do método [33]. Simplificando, pode-se
12
dizer que um pré-condicionador é qualquer forma de modificação de um sistema linear original,
seja ela implícita ou explícita, que faz com que o sistema possa ser mais facilmente resolvido
por um método iterativo [30].
A escolha da matrizS determina o método iterativo que teremos. Existem diversasalterna-
tivas para decomposição, as mais populares são os tiposL, D eU , ondeL é a matriz triangular
estritamente inferior;U a matriz triangular estritamente superior;D a matriz diagonal. Cabe
observar que a matrizD, por definição, não pode conter zeros na diagonal principal,caso isso
aconteça, deve-se trocar linhas ou colunas na matriz. Além disso, o erro contido na solução
será diretamente proporcional ao tamanho dos coeficientes da equação de recorrência, ou seja,
organizando-se a matriz de modo a deixar os maiores elementos na diagonal, tende-se a ter um
erro menor [33]. Na seqüência são apresentados os métodos iterativos mais populares, os quais
surgem de acordo com a escolha dos pré-condicionadores citados. Assim como nos métodos
diretos, não serão apresentados os algoritmos que implementam esses métodos computacional-
mente, os quais podem ser encontrados em [17, 30]. Cabe aindaobservar que esses métodos são
chamados esquemas de relaxação pontual, pois atualizam um componente por vez, e podem ser
generalizados para esquemas de relaxação de blocos, onde atualiza-se um conjunto inteiro de
valores a cada passo [30].
Método de Jacobi
É obtido considerando-se a decomposiçãoA = D−L−U . ComoA = S−T pode-se obter
S = D eT = L + U que substituindo em (2.8) implica emD~x = (L + U)~x + ~y. Isolando-se o
vetor das incógnitas e o reescrevendo de acordo com (2.7) tem-se a equação iterativa do Método
de Jacobi, a qual, na forma matricial é dada por [33, 30]:
~xk+1 = D−1(L + U)~xk + D−1~y (2.9)
Reescrevendo na forma de componentes tem-se o passo:
xk+1i =
1
aii(yi −
n∑
j=1;j 6=i
aijxkj ) i = 1, . . . , n. (2.10)
13
Um inconveniente do método de Jacobi é que ele requer que se mantenham armazenados
todos os componentes de~xk até que o cálculo dexk+1 esteja completo, o que dependendo
do tamanho das matrizes pode ser um problema. Uma abordagem mais natural e que requer
apenas a metade do armazenamento é começar a utilizar cada componente do novo vetor~xk+1
conforme for computado, deixando então~xk livre para ser destruído. É o que ocorre no método
de Gauss-Seidel [33, 17].
Método de Gauss-Seidel
Neste método considera-se S = D + L e T = -U. Portanto, a equaçãona forma matricial do
método é definida como [33, 30]:
~xk+1 = (D − L)−1U~xk + (D − L)−1~y (2.11)
Cujo passo iterativo na forma de componentes é:
xk+1i =
1
aii(−
i−1∑
j=1
aijxk+1i −
n∑
j=i+1
aijxkj + yi) i = 1, . . . , n. (2.12)
Em geral, o método de Gauss-Seidel converge mais rapidamente do que o método de Jacobi,
entretanto, há casos em que isso não acontece, bem como também há casos em que um método
converge e outro não [33].
Método de Relaxamento Sucessivo
O método de Gauss-Seidel obtém o valor de uma incógnita na interaçãok somando o valor
da mesma na interaçãok−1 com uma correção calculada com base em todas as variáveis. Após
anos de utilização desse método percebeu-se que a convergência pode ser acelerada com intro-
dução de um fatorω de relaxamento, ou seja, atribuindo-se um peso a essa correção. Atribuindo-
se o valor1 paraω continuamos com Gauss-Seidel, já com valores em queω > 1 o método passa
a utilizar essa aceleração e é conhecido como Método de Relaxamento Successivo (Sucessive
Overrelaxation - SOR). A escolha do melhor valor paraω depende de cada problema, embora
nunca seja maior do que2 e, em geral, fica próximo de1.9 [33].
14
Considerando-se novamente a decomposiçãoA = L−D−U e ponderando a correção com
ω, a equação matricial do método SOR é dada por [33, 30]:
(D − ωL)~xk+1 = [ωU + (1 − ω)D]~xk + ω~y (2.13)
Cujo passo iterativo na forma de componentes é:
xk+1i = ω(
1
aii(−
i−1∑
j=1
aijxk+1i −
n∑
j=i+1
aijxkj + yi)) + 1(1 − ω)xk
i i = 1, . . . , n. (2.14)
A busca do valor ótimo paraω foi o foco das pesquisas por diversos anos, em [33] são
apresentados os passos que levaram à conclusão de que a escolha ótima é aquela que faz com que
o maior autovalor deL seja o menor o possível. Em uma matriz de ordem21, por exemplo, pode-
se obter com um passo SOR o resultado que o método de Jacobi levaria30 passos para alcançar
e, ainda, com uma redução do erro em 25%. Esse é considerado umresultado muito bom gerado
a partir de uma idéia bastante simples. Sua aplicação real, entretanto, não é em problemas de
uma dimensão, pois os sistemasA~x = ~y nesses casos são simples de serem resolvidos, mas
naqueles com mais dimensões, como por exemplo na solução de equações diferenciais parciais
[33].
2.4 Aplicações práticas de modelagem baseada em EDs e SELAs
Conforme colocado anteriormente, diversas são as áreas queutilizam a modelagem e simu-
lação para resolver e obter uma melhor compreensão dos mais variados problemas. As equações
diferenciais têm inúmeras aplicações práticas em medicina, engenharia, química, biologia etc,
onde as soluções dessas equações são empregadas no processamento de imagens, projeto de
pontes, automóveis, aviões, dentre outros. Esta seção exemplifica brevemente algumas dessas
aplicações.
15
2.4.1 Método dos Elementos Finitos aplicado à Odontologia
O desenvolvimento do Método dos Elementos Finitos iniciou-se no final do século XVIII,
quando Gauss propôs a utilização de funções de aproximação para a solução de problemas
matemáticos. Conforme visto anteriormente, consiste na discretização de um meio contínuo
(EDs), mantendo as mesmas propriedades do meio original e seu uso só se tornou viável com o
advento dos computadores para solucionar os enormes SELAs gerados na solução do MEF [29].
O MEF vem sendo utilizado há alguns anos em experimentos relacionados às mais diver-
sas especialidades da Odontologia, em [29] os autores apontam uma série desses trabalhos. A
utilização do MEF nesses casos traz diversas vantagens se comparada aos métodos tradicionais,
dentre as quais cabe destacar: não há necessidade de laboratórios super equipados com instru-
mentos específicos, o que reduz os custos dos experimentos; as formas geométricas podem ser
representadas mais próximas das reais, enquanto métodos tradicionais costumam trabalhar com
formas ideais; os resultados numéricos apresentam maior exatidão, ao contrário das abordagens
analíticas que costumam trazer resultados simplificados; os experimentos são mais fáceis de
conduzir e reproduzir do que aqueles que utilizam seres vivos e tendem a ter resultados mais
exatos devido ao maior controle da situação [29].
O MEF é aplicado na Ortodontia para: estudo da distribuição das tensões durante o movi-
mento dentário; estudo do efeito de forças ortopédicas no complexo craniofacial; avaliação da
resistência de base de braquetes e do desempenho de molas para fechamento de espaços, dentre
outros. Esses estudos obtêm informações, por exemplo, sobre o comportamento dos tecidos e
movimentos do complexo maxilar, o que permite otimizações no uso de aparelhos e em im-
plantes dentários [29].
2.4.2 EDs e Método dos Elementos Finitos aplicados à Engenharia
Em todas as áreas da engenharia uma série de situações podem ser simuladas através de
equações diferenciais. Um problema freqüente na engenharia civil e engenharia mecânica é
análise das tensões e deformações de materiais, o qual é facilmente modelado por EDPs. Na
engenharia geotécnica um fenômeno muito conhecido é o fluxo de água em meios porosos, que é
16
analisado, por exemplo, em casos de barragens de terra com fluxo e enchimento estacionários, ou
seja, em que a carga hidráulica não varia com o tempo. Outra aplicação é a análise das barragens
de terra com fluxo transiente, ou seja, em que o nível de água pode variar em condições rápidas,
levando a uma variação na umidade do material e, conseqüentemente, da condutividade elétrica
do mesmo. Nesse caso, a modelagem é dada por uma equação diferencial em que são definidas
condições de contorno e condições iniciais [12].
A engenharia elétrica também aplica equações diferenciaisem diversos estudos, dentre eles,
pode-se citar a análise de capacitância, indutância e diferenças de potencial dos circuitos elétri-
cos. Já na Engenharia Química, um procedimento bastante comum é efetuar balanceamentos de
massa, volume ou energia em um determinado sistema, como um reator químico, por exemplo,
modelagem que pode ser feita através de uma EDO [11].
2.4.3 EDs aplicadas à Física e Matemática
Diversas fórmulas bastante conhecidas da Física e Matemática são oriundas de ou podem ser
vistas como equações diferenciais, a segunda Lei de Newton (~F = m~a) é um exemplo. Outra lei
de Newton que consiste em uma ED é a Lei de Resfriamento, a qualdiz que a taxa de variação da
temperatura de um corpo em resfriamento é proporcional à diferença entre sua temperatura atual
e a temperatura constante do meio ambiente. Cálculos de trajetórias, escoamento de fluídos, a
Lei de Hooke, a Equação da Onda, Equações de Maxwell, Equaçãodo Calor, Lei de Torricelli,
enfim, grande parte dos problemas cotidianos da física podemser modelados por EDs [11].
No contexto da matemática essa situação se repete. O estudo das dinâmicas de populações
baseia-se em equações diferenciais que modelam o crescimento ou decaimento da população ao
longo do tempo (modelos de Malthus e Verhulst). Um exemplo damatemática financeira mod-
elada por EDPs é o cálculo de juros e investimentos. Considerando uma caderneta de poupança
que recebe um investimento inicial, o crédito dos juros mensais e também depósitos adicionais
contínuos, podem-se relacionar esses parâmetros por uma EDpara obter informações do tipo:
qual deve ser o valor dos depósitos adicionais para atingir oobjetivo de acumular um determi-
nado valor em um período de tempo [11].
17
2.5 Conclusão
As equações diferenciais são uma das principais ferramentas dos cientistas nas mais variadas
áreas do conhecimento. Entretanto, a solução analítica dessas equações é bastante complicada.
A alternativa, que é a solução numérica, era impraticável nopassado mas tornou-se viável graças
aos avanços da computação, principalmente, nas últimas duas décadas, o que tem possibilitado
diversos avanços na simulação e otimização de problemas. A solução numérica das EDs exige
uma discretização a qual é feita, em geral, pelos Métodos dasDiferenças Finitas e Métodos dos
Elementos Finitos, em ambos os casos o processo resulta na necessidade de resolver sistemas de
equações lineares algébricas. Além disso, muitos problemas podem ser modelados diretamente
na forma de sistemas lineares, assim como a resolução de outros tipos de equações e modelos
teóricos passa por passos em que é necessário resolver SELAs. Essas realidades demonstram a
importância de pesquisar constantemente maneiras para otimizar os métodos de solução tanto
em termos de eficiência como de precisão.
18
3 Computação paralela
A perspectiva da exploração do paralelismo nos sistemas computacionais é tão antiga quanto
os primeiros computadores. Trabalhos desenvolvidos por von Neumann, na década de 40, já
ponderavam a possibilidade de utilizar algoritmos paralelos. Registros apontam que no período
de 1944 a 1947 Stibitz e Williams teriam desenvolvido, nos laboratórios da Bell Telephone,
aquele que é considerado um dos primeiros trabalhos envolvendo paralelismo com implemen-
tação emhardware, o sistema MODEL V. Esse sistema era formado por dois processadores
e três canais de entrada e saída, possibilitando a execução simultânea de dois programas dis-
tintos [34]. Durante essa trajetória, o emprego do paralelismo já passou por fases de extremo
otimismo, quando era visto como única e promissora solução para as mais complexas aplicações
computacionais, bem como por fases pessimistas, nas quais foi encarado como um nicho de
mercado decadente [32]. Nos últimos anos a computação paralela avançou muito e, atualmente,
está presente em todos os níveis da computação indo desde o paralelismo de instruções em um
pipeline, passando pelos processadores multi-core, colocados comosolução para os problemas
de limitação no crescimento doclock, até sistemas distribuídos de larga escala como as grades
computacionais.
Neste capítulo serão inicialmente abordados aspectos gerais da computação paralela e definidos
alguns conceitos fundamentais da área. Na seqüência, discute-se a necessidade de modelos
para desenvolvimento de software e algumas medidas utilizadas na computação paralela são
descritas. Feito isso, são apresentadas as classificações dos sistemas paralelos em relação à ar-
quitetura dehardware. Por fim, são abordados os paradigmas de programação utilizados para
exploração desses sistemas.
19
3.1 Contextualização Geral
Muitas vezes os termoscomputação paralelae computação distribuídasão utilizados em
conjunto, é importante, no entanto, saber que podem tratar-se de sistemas diferentes. Uma dis-
tinção inicial é dizer que os sistemas de computação paralela são aqueles compostos por um
conjunto de processadores fisicamente próximos, os quais têm como objetivo trabalhar conjun-
tamente para resolver um problema e se comunicam de forma confiável e previsível. Os sistemas
distribuídos, por outro lado, seriam aqueles cujos processadores podem estar geograficamente
distribuídos e entre os quais a comunicação é uma questão mais problemática, dado que men-
sagens podem ser perdidas ou serem entregas com atraso desconhecido. A topologia do sistema
distribuído pode ser alterada em tempo de execução, seja devido a problemas de comunicação
como links quebrados ou nós que deixam de funcionar ou pela adição/remoção de processadores.
Além disso, os sistemas distribuídos apresentam pouco ou nenhum controle centralizado, cada
processador pode estar executando suas próprias atividades ao mesmo tempo que coopera com
o SD [8].
Assim, pode-se dizer que os sistemas paralelos são aqueles cujo hardwareé dedicado e
tende a ser homogêneo, enquanto na computação distribuída ohardwareé mais heterogêneo
e a computação distribuída em si pode ser apenas uma atividade secundária do sistema. Essa
diferenciação é importante, pois, embora não exista uma linha clara dividindo os dois, ajuda a
perceber que cada caso possui questões próprias, mesmo que muitas delas sejam compartilhadas
e possuam, inclusive, soluções semelhantes [8].
3.2 Modelagem de programas paralelos
A compreensão do inter-relacionamento entre as características de um algoritmo paralelo e
as da arquitetura em que ele irá executar é fundamental, poispara que o processamento possa
ter desempenho otimizado em uma arquitetura paralela é necessário haver equilíbrio entre as
características dohardwaree software. Entretanto, as arquiteturas para processamento paralelo
estão em constante evolução e, portanto, nenhuma arquitetura pode ter a expectativa de se manter
como solução ideal por mais do que alguns anos. Assim, buscando contornar o problema da falta
20
de portabilidade e curto tempo de vida útil dosoftwareparalelo, surgiu a idéia da criação de
modelos suficientemente abstratos que possibilitem ocultar os aspectos da arquitetura e manter
o desempenho da aplicação quando houver mudanças dohardware[34, 32].
Os autores de [32] argumentam que um modelo desse tipo pode ser construído da seguinte
forma: adota-se uma máquina abstrata para a qual osoftwareé desenvolvido e a máquina, por
sua vez, é emulada pelas diferentes arquiteturas paralelas. Dessa forma, os aspectos do projeto de
softwarepodem ser desassociados daqueles de sua implementação. Argumentam ainda que um
modelo ideal deve atender algumas necessidades mínimas: oferecer facilidade e metodologias
para o desenvolvimento desoftware, ser independente de arquitetura, ser de fácil compreensão,
apresentar garantias de desempenho e medidas de custo.
Uma classificação dos modelos para desenvolvimento desoftwareparalelo, utilizando como
critério o nível de abstração com que o paralelismo pode ser explorado, é sugerida em [32].
As categorias apresentadas são as seguintes: Modelos nos quais o paralelismo é explorado de
forma totalmente implícita; Modelos com assinalamento do paralelismo explícito; Modelos com
assinalamento e decomposição do paralelismo explícitos; Modelos com assinalamento, decom-
posição e mapeamento do paralelismo explícitos; Modelos com assinalamento, decomposição,
mapeamento e comunicação explícitos; Modelos em que o paralelismo é explorado de forma
totalmente explícita [32];
Os esforços da comunidade científica têm sido direcionados para os modelos de mais alto
nível, os quais buscam facilitar a programação [34]. O OpenMP, API para desenvolvimento de
aplicações paralelas em arquiteturas de memória compartilhada, embora não seja um modelo,
pode também ser visto como um esforço no sentido de facilitara programação. O programador
OpenMP precisa apenas delimitar, através de diretivas de compilação, os trechos paralelizáveis
do código, sem se preocupar com aspectos de sincronização e comunicação. Feito isso, os
compiladores geram automaticamente o código paralelo [20]. Entretanto, até os dias de hoje
nenhum modelo que atendesse as necessidades propostas por [32] foi desenvolvido. Neste tra-
balho os algoritmos serão apresentados graficamente através de DAGs (Directed Acyclic Graphs
ou Grafos Dirigidos Acíclicos), uma modelagem independente dehardwaree das questões de
comunicação que permite a visualização do problema em um alto grau de abstração.
21
Um grafo acíclico dirigido é um grafo dirigido sem ciclos, isto é, para qualquer vérticev,
não há nenhum caminho dirigido começando e acabando emv. SejaG = (N,A) um DAG,
ondeN = 1, . . . , N é o conjunto de nós eA é o conjunto dos arcos dirigidos, cada nó
representa uma operação que pode ser feita pelo algoritmo e os arcos as dependências de dados.
Em particular, um arco(i, j) ∈ A indica que a operação correspondente ao nój usa o resultado
da operação correspondente ao nói. As operações podem ser elementares, como a adição de
dois escalares, ou operações de alto nível, como a execução de uma subrotina. A Figura 3.1
exemplifica essa modelagem, o algoritmo representado corresponde à expressão(x1+x2)(x2+
x3). O símbolo dentro de cada nó ilustra a operação por ele realizada, sendo queS corresponde
a operação elevação ao quadrado (square) [8].
Figura 3.1: DAG de um algoritmo para avaliação da expressão(x1 + x2)(x2 + x3). Fonte: [8].
Um nó i ∈ N é dito predecessor de outro nój ∈ N se (i, j) ∈ A. Define-se como
grau de entrada de um nói ∈ N o número de predecessores desse nó e como grau de saída o
número de nós dos quaisi é predecessor. Os nós com grau zero de entrada são chamados nós
de entrada e os nós com grau zero de saída são ditos nós de saída. O conjunto de todos os nós
de entrada é representado porN0. Um caminho positivo é uma seqüênciai0, . . . , ik de nós tais
que(ik, ik+1) ∈ A parak = 0, . . . ,K − 1, ondeK é dito o tamanho do caminho. Denota-se
por D a profundidade de um DAG, a qual é dada pelo tamanho do maior caminho positivo.
Considerandoxi o resultado da operação correspondente aoi-ésimo nó do DAG, então o DAG
22
pode ser visto como uma representação de dependências funcionais na forma [8]:
xi = fi(xj/j precedei) (3.1)
Nessa representação,fi é uma função descrevendo a operação correspondente aoi-ésimo nó.
Sei é um nó de entrada, entãoxi não depende de outras variáveis e é visto como uma variável
de entrada externa. Assim, a operação correspondente a um nóde entradai é essencialmente
o tempo gasto para ler o valor da variável de entradaxi, o qual assume-se como insignificante.
Para um nói que não é de entrada (i ∈ N0), assume-se que a operação correspondente, ou seja,
a avaliação da funçãofi, demora uma unidade de tempo. Entretanto, essa premissa somente
é razoável quando cada nó representa uma operação aritmética, pois em algoritmos numéricos
mais complicados o tempo de execução em cada nó pode ser bastante diferente. Por isso, utiliza-
se a definição auxiliar de escalonamento descrita a seguir [8].
Um DAG é uma representação apenas parcial do algoritmo, a qual especifica que operações
acontecem sob quais operandos e impõe restrições de precedência entre elas. Para representar
completamente o algoritmo paralelo é preciso ainda especificar qual processador será respon-
sável por cada operação e em que tempo isso acontecerá. Para isso, assume-sep como sendo o
conjunto dos processadores disponíveis e também que cada processador é capaz de fazer qual-
quer uma das operações desejadas. Para cada nói que não é um nó de entrada (i ∈ N0),
representa-se porPi o processador para o qual foi designada a responsabilidade de realizar a op-
eraçãoi. Também parai ∈ N0, assume-seti como uma variável inteira positiva que especifica
o tempo que leva para a operação correspondente ao nói ser completada. Nenhum processador
é designado para os nós de entrada e utiliza-se a convençãoti = 0 para cada nó de entradai.
Existem duas restrições que precisam ser impostas [8]:
• Um processador pode realizar no máximo uma operação por vez,logo (i ∈ N0, j ∈
N0, i 6= j, ti = tj) → Pi 6= Pj .
• (i, j) ∈ A → tj ≥ ti +1. Ou seja, se existe uma aresta entre um nói e outroj, a operação
correspondente aoj somente pode iniciar após a operação do nói estar completa.
23
DefinidosPi e ti sujeitos às restrições acima, diz-se que o DAG estáescalonadopara exe-
cução paralela e define-se como escalonamento o conjunto(i, Pi, ti)/i ∈ N0. Essa configu-
ração pode ser aplicada a uma grande variedade de implementações. Por exemplo, o processador
Pi pode armazenar o resultadoxi de sua operação em uma memória compartilhada e esse ser
recuperado por outros processadores. Alternativamente, em uma implementação por troca de
mensagens, o processadorPi envia uma mensagem com o valor dexi para qualquer processador
Pj que precise desse valor, ou seja,(i, j) ∈ A. Na prática, os acessos à memória ou trans-
missões de mensagens demoraram um certo tempo, o qual foi ignorado até então. Assim, se a
transmissão de uma mensagem, por exemplo, requerτ unidades de tempo, e se(i, j) ∈ A, então
a restriçãotj ≥ ti + 1 poderia ser modificada parati ≥ ti + 1 sePi = Pj e ti ≥ ti + τ + 1 se
Pi 6= Pj [8].
3.3 Medidas em Computação Paralela
Existem diversas medidas no contexto da computação paralela, sendo as principais as me-
didas de complexidade, as quais objetivam quantificar o total de recursos computacionais uti-
lizados por um algoritmo paralelo, e as medidas de desempenho e eficiência. Algumas medidas
consideram, por exemplo, o número de processadores, o tempogasto até o algoritmo terminar
(complexidade de tempo) ou o número de mensagens transmitidas na execução do algoritmo
(complexidade da comunicação). As medidas de complexidadesão freqüentemente expres-
sas como funções do tamanho do problema a ser resolvido e informalmente definidas como
o número de entradas da computação. Assim, um algoritmo parasomarn inteiros, por exemplo,
é dito de tamanhon. No caso dos algoritmos distribuídos existe uma preocupação adicional
para definição da complexidade de tempo: é possível que um algoritmo termine e nem todos
os processadores tenham conhecimento disso. Nesses casos énatural somar o tempo adicional
necessário para os processadores tomarem conhecimento da terminação [8]. A seguir são discu-
tidas algumas dessas medidas.
24
3.3.1 Complexidade de tempo em algoritmos especificados porum DAG
SejaG = (N,A) um DAG representando um algoritmo paralelo e seja(i, Pi, ti)/i ∈ N0
um escalonamento deste DAG que usap processadores, o tempo gasto por esse escalonamento
é igual aMAXi∈N ti. Define-seTp como sendo o mínimo deMAXi∈N ti, o qual é pego dentre
todos os possíveis escalonamentos que usamp processadores. Tem-se assimTp como a complex-
idade de tempo do algoritmo descrito porG, a qual é uma função do número de processadores
disponíveis [8]. Assume-se:
T∞ = minTp p ≥ 1 (3.2)
Dado queTp é um valor inteiro, existe um valor mínimop∗ tal queTp = T∞∀(p ≥ p∗).
Considera-seT∞ como a complexidade de tempo do algoritmo especificado porG quando um
número suficientemente grande de processadores (pelo menosp∗) está disponível, sendo que
T∞ é igual a profundidadeD do DAG. Define-seT1 como sendo o tempo necessário para uma
execução seqüencial do algoritmo modelado, o qual, naturalmente, será igual ao número de nós
do DAG que não são nós de entrada [8].
Dado um valor arbitrário dep, tem-seT1 ≥ Tp ≥ T∞. Em geral, O valor exato deTp não é
fácil de ser determinado, porém, esse problema pode ser contornado com a definição de alguns
limites paraTp [8]:
1. Supondo que para um nó de saídai qualquer existe um caminho positivo de cada nó de
entrada e que o grau de entrada de cada nó é no máximo 2, então:T∞ ≥ longn, onden é
o número de nós de entada;
2. Sejac um inteiro positivo,q = cp → Tp ≤ cTq;
3. Para cadap tem-seTp < T∞ + T1
p ;
4. (a)p ≥ T1
T∞
→ Tp < 2T∞. Generalizando:p = Ω( T1
T∞
) → Tp = O(T∞).
(b) p ≤ T1
T∞
→ T1
p ≤ Tp < 2T1
p . Generalizando:p = O( T1
T∞
) → Tp = θ(T1
p )
A propriedade (1) define uma limitação fundamental da velocidade dos algoritmos paralelos
enquanto a propriedade (2) expressa o fato de que, se o númerode processadores é reduzido
25
em certo fator, então o tempo de execução é aumentado até, no máximo, esse mesmo fator. Os
dois últimos resultados são de fundamental importância, pois estabelecem que, apesar deT∞ ser
definido com o pressuposto de um número ilimitado de processadores,Ω(T1/T∞) processadores
são suficientes para ter um fator constante deT∞ (4a). Dessa forma, pode-se obter um escalon-
amento simplesmente modificando o escalonamento ótimo parao caso de um número ilimitado
de processadores, ao invés de resolver o problema, geralmente difícil, de escalonamento. Isso
sugere uma metodologia na qual primeiro desenvolve-se um algoritmo paralelo como se um
número de processadores ilimitados estivesse disponível e, então, faz-se a adaptação do mesmo
para o número de processadores realmente disponível. Já a propriedade (4b) significa que, con-
tanto quep = O( T1T∞
), a disponibilidade dep processadores permite acelerar a computação em
um fator diretamente proporcional ap, que é o melhor possível. Assim, para um número de
processadores próximo deT1
T∞
obtém-se ao mesmo tempo a execução e aceleração ótimas do
algoritmo [8].
3.3.2 Encontrando um DAG ótimo
Em geral, para um mesmo problema computacional existem diferentes algoritmos e, ainda,
para cada algoritmo diferentes DAGs. Assim, o objetivo é encontrar o DAG parap proces-
sadores com menorTp, o qual é chamado DAG ótimo. SejaT ∗p o valor deTp correspondente ao
DAG ótimo, entãoT ∗p corresponde ao menor tempo de processamento paralelo, usando p pro-
cessadores, para o problema em questão. Assim,T ∗p é a medida de complexidade do problema,
diferentemente deTp que é a medida de complexidade para um algoritmo em particular. Entre-
tanto, a avaliação explícita deT ∗p é, normalmente, muito difícil [8]. A Figura 3.2 apresenta um
DAG otimizado para o mesmo problema proposto na Figura 3.1. Para o DAG otimizado tem-se
T1 = 3 eT∞ = D = 2, enquanto o DAG da Figura 3.1 resolvia o mesmo problema comT1 = 7
eT∞ = D = 3.
3.3.3 SpeedUP e Eficiência
Uma vez que um modelo particular de computação paralela tenha sido escolhido, seja o mod-
elo DAG anteriormente descrito ou outro qualquer, dado um problema computacional parametrizado
26
Figura 3.2: DAG de outro algoritmo para avaliação da expressão (x1 + x2)(x2 + x3). Fonte:
[8].
por uma variáveln que representa o tamanho do problema, a complexidade de tempo é geral-
mente dependente den [8], e, por isso,n é incorporada na notação a seguir.
Seja um algoritmo paralelo que usap processadores (p pode depender den) e termina no
tempoTp(n) e sejaT ∗(n) o tempo requerido pelo algoritmo seqüencial ótimo para resolver
o mesmo problema, define-se comospeedupa medida relativa entre esses dois tempos. Em
outras palavras,speedupé a medida, dada pela equação 3.3, que indica o número de vezesque
o programa paralelo é mais rápido que o seqüencial [8, 34].
Sp(n) =T ∗(n)
Tp(n)(3.3)
Eficiência é uma medida utilizada em conjunto com ospeedup. Demonstra a fração de
tempo que os processador estão sendo utilizados e é expressaem valores percentuais. O cálculo
da eficiência é dado pela equação 3.4 [8, 34].
Ep(n) =Sp(n)
p=
T ∗(n)
pTp(n)(3.4)
Valores ideais despeedupe eficiência sãoSp(n) = p e Ep(n) = 1, o que é equivalente
a dizer que disponibilidade dep processadores aumenta a velocidade da computaçãop vezes e
que todos os processadores são totalmente aproveitados durante todo o tempo. Entretanto, para
27
isso ocorrer, o algoritmo paralelo teria que ser tal que nenhum processador em momento algum
fizesse trabalho desnecessário ou ficasse parado à espera de trabalho. Por motivos diversos, essa
situação não é possível na prática e um objetivo mais realista é tentar manter o valor deEp tão
distante quanto possível de zero conformep aumenta.
Porém, as definições despeedupe eficiência são relativamente controversas, dado que, no
geral, o tempo seqüêncial ótimoT ∗(n) é desconhecido, mesmo para problemas aparentemente
simples. Por essa razão, encontra-se diferentes definiçõesparaT ∗(n) , sendo as alternativas
mais comuns as seguintes [8]:
1. T ∗(n) é o tempo gasto pelo melhor algoritmo seqüencial existente;
2. T ∗(n) é o tempo gasto por um algoritmo seqüencial de teste. Por exemplo, θ(n3) é um
valor razoável para o problema de multiplicação de duas matrizes densasn x n, mesmo
que existam algoritmos para resolver o mesmo problema em menor tempo;
3. T ∗(n) é o tempo gasto por um único processador para executar o algoritmo paralelo que
está sendo analisado em particular, ou seja, um único processador simula a operação dep
processadores paralelos.
Utilizando a definição (3), a eficiência corresponde ao quão bem o algoritmo em particular
foi paralelizado, porém não provê informações sobre os méritos absolutos do algoritmo, como
acontece nas opções (a) e (b). Além disso, seT ∗(n) é definido de acordo com (c) e o algoritmo
especificado com o modelo DAG, entãoT ∗(n) equivale aT1(n) [8].
Uma observação importante ao tratar despeedupé que mesmo programas com trechos de
fácil paralelização possuem também seções inerentemente seqüenciais. Dessa forma, em um
ambiente com um grande número de processadores a parte paralela do código será rapidamente
executada e as seções seqüenciais serão o gargalo da aplicação. Essa situação é caracterizada
pela Lei de Amdahl, a qual define que todo programa paralelo possui uma fração "f"seqüencial
que limita seuspeedupde acordo com a equação 3.5. Por outro lado, existem numerosos proble-
mas em que, proporcionalmente, o trecho "f"do código tende azero conforme o problema cresce
e, portanto, a Lei de Amdahl deixa de ser uma preocupação [8].
28
Sp(n) ≤1
f + (1−f)p
∀p (3.5)
3.3.4 Granulosidade e Escalabilidade
Granularidade, ou granulosidade, é a medida que diz respeito à quantidade de computação
de uma atividade. Diz-se de granulosidade fina os problemas em que a quantidade de cálculo
é pequena, ou seja, a comunicação e sincronização entre os fluxos de execução é freqüente; de
granulosidade grossa aqueles em que a quantidade de cálculoé grande e, portanto, a freqüência
de sincronização é menor; de granularidade média aqueles emque a relação é mais equilibrada.
Quanto maior a granularidade menor serão os custos com criação de processos e com a comuni-
cação entre eles. É comum existir abundância de paralelismocom baixa granulosidade, o qual,
via de regra, não pode ser explorado com eficiência. O desempenho da arquitetura paralela é
obtido através de um equilíbrio entre granulosidade e comunicação [34].
Na implementação de uma solução paralela, a opção por explorar granularidade fina ou
grossa pode depender de vários aspectos. Por exemplo, na programação por troca de mensagens
é interessante se reduzir a comunicação entre processos devido ao tempo gasto na comunicação
entre os computadores, o que pode sugerir o aumento do númerode instruções por processo. Já
na programação comthreadsos custos com a comunicação podem não ser uma preocupação[34].
Escalabilidade, por sua vez, é um termo que pode ser aplicadotanto em relação aohardware
como aosoftware. No caso dohardware, é usado para indicar que a arquitetura permite expansão
no seu tamanho físico, sendo que, isso faz com que seu desempenho também aumente. Ou seja,
diz-se que uma arquitetura apresenta boa escalabilidade quando é possível adicionar-lhe mais
nós de processamento e esse incremento faz aumentar, também, a performance. No contexto do
softwareo termo é usado para referenciar a capacidade de o mesmo se adaptar às mudanças da
arquitetura, ou seja, umsoftwareé dito escalável quando consegue explorar eficientemente novos
recursos que venham a ser adicionados à arquitetura. Uma maneira de fazer isso é decompor o
algoritmo em mais tarefas do que o número de processadores disponíveis, assim, à medida que
o número de processadores cresce, as tarefas passam a ser distribuídas também para os novos
processadores, reduzido o número de tarefas por processador sem a necessidade de modificar o
29
algoritmo.
3.4 Classificação das arquiteturas paralelas
Em [31] o autor apresenta três motivos pelos quais é importante classificar as arquiteturas:
o primeiro é entender o que já foi desenvolvido, o segundo é que a classificação pode ajudar os
projetistas a visualizarem novas configurações ainda não pensadas e o terceiro é que a classi-
ficação pode ajudar na construção de modelos de performance.Modelos esses que podem ser
utilizados para avaliação de desempenho e também ajudar a revelar porque uma dada arquitetura
possui maior flexibilidade para aumento do desempenho em relação a outras.
Além disso, uma classificação deve agrupar as máquinas existentes de modo que, sabendo
a classificação de uma máquina, seja possível inferir informações diversas sobre ela. Assim, a
escolha do critério é uma tarefa difícil e que acaba por determinar o sucesso ou não da clas-
sificação. Um critério demasiadamente genérico pode não sersuficiente para que se obtenha
qualquer informação útil da classificação, já um critério muito específico pode torná-la depen-
dente de tecnologias e, portanto, de sua própria época, o quecertamente a tornará obsoleta
precocemente. Classificações atuais consideram critérioscomo compartilhamento de memória
e o tratamento da coerência nas memóriascache.
3.4.1 Classificações Acadêmicas
Uma das primeiras classificações, largamente utilizada atéhoje, para arquiteturas paralelas
é a Taxonomia de Flynn, proposta em [21, 22]. Essa proposta baseia-se na relação entre o fluxo
de instruções e o fluxo de dados, que podem ser únicos ou paralelos. Paralelismo de instruções
(ou funcional) ocorre quando cada tarefa executa cálculos diferentes para resolver um problema,
sejam eles sobre um mesmo conjunto de dados ou sobre dados diferentes. Paralelismo de dados
acontece quando um processo executa uma mesma série de operações sobre diferentes dados.
Como esses fluxos são independentes, Flynn combinou-os e propôs quatro classes possíveis
[22, 34]:
• SISD -Single Instruction Stream Single Data Stream(Fluxo Único de Instruções Fluxo
30
Único de Dados): são os computadores seqüenciais convencionais, onde uma única unidade
de controle decodifica seqüencialmente as instruções que operam sobre um único conjunto
de dados;
• SIMD - Single Instruction Stream Multiple Data Streams(Fluxo Único de Instruções
Múltiplos Fluxos de Dados): são os processadores matriciais, nos quais uma única unidade
de controle ativa diversos elementos processadores. A unidade de controle está submetida
a um único programa e repassa suas instruções para os processadores que, por sua vez, as
executam concorrentemente sobre os dados em suas memórias locais;
• MISD - Multiple Instruction Streams Single Data Stream(Múltiplos Fluxos de Instruções
Fluxo Único de Dados): costuma-se dizer que essa classe é considerada para fins didáticos
mas que, porém, não existem computadores reais enquadradosnesta categoria. O autor,
entretanto, argumenta em [22] que é fácil prever e desenvolver processadores desse tipo
e que o desinteresse histórico nos mesmos deve-se a fatores relacionados à programação.
Coloca ainda que um MISD é, na verdade, umpipelinecom múltiplas unidades funcionais
executando independentemente e com adiantamento em um mesmo fluxo de dados. Ar-
gumenta também que no nível microarquitetural é dessa mesmaforma que funciona um
processador vetorial. No entanto, pode-se também considerar que, após passar por um
estágio dopipeline, o dado não é mais o mesmo e, portanto, o argumento do autor ficaria
comprometido. Alguns autores consideram que os arranjos sistólicos de processadores se
enquadram nesta classe;
• MIMD - Multiple Instruction Streams Multiple Data Streams(Múltiplos Fluxos de In-
struções Múltiplos Fluxos de Dados): é a categoria na qual seenquadra a grande maioria
dos sistemas paralelos. Os elementos processadores executam de maneira independente
instruções independentes sobre dados também independentes.
Exemplos de classificações propostas que não obtiveram êxito podem ser encontrados em
[31] e [27]. No primeiro caso a classificação considera oito critérios: número de processadores
de instruções; número de memórias de instruções; tipo de chaveamento conectando os proces-
sadores de instruções às memórias de instruções; número de processadores de dados; número
31
de memórias de dados; tipo de chaveamento conectando processadores de dados e memórias de
dados; tipo de chaveamento conectando processadores de instruções e processadores de dados;
e o tipo de chaveamento de dados conectando processadores dedados entre si. Já em [27], a
proposta é subdividir a classe MIMD da Taxonomia de Flynn em quatro novas subclasses de
acordo com com a estrutura de memória (global ou distribuída) e com os mecanismos de co-
municação/sincronização (variáveis compartilhadas ou troca de mensagem) entre os processos.
Entretanto, embora essas propostas tenham caído em desuso ou não tenham sido adotadas, elas
serviram para atentar ao fato de que a Taxonomia de Flynn não era suficiente para diferenciar
adequadamente a diversidade de máquinas paralelas que surgiam já naquela época.
A abordagem atualmente mais usada para classificação das arquiteturas paralelas subdivide
a categoria das arquiteturas MIMD considerando critérios relacionados à memória. O primeiro
nível da hierarquia é baseado no compartilhamento de memória, resultando em duas classes
[19]:
• Multiprocessadores: são as arquiteturas cujos processadores compartilham memória (shared
memory). Nesse caso, independentemente da memória estar implementada de forma
distribuída (distributed memory) ou centralizada (centralized memory), todos os proces-
sadores compartilham um mesmo espaço de endereçamento e se comunicam através de
operações de escrita e leitura em variáveis compartilhadas(programaçãomultithread);
• Multicomputadores: arquiteturas que possuem múltiplos espaços de endereçamento pri-
vado (multiple private address spaces), ou seja, cada processador possui um espaço próprio.
A implementação da memória é distribuída e os processadoresse comunicam por trocas
de mensagens. As duas principais bibliotecas para implementação de aplicações com tro-
cas de mensagens são a PVM (Parallel Virtual Machine) [5] e MPI (Message Passing
Interface) [4] .
O segundo nível da hierarquia classifica as máquinas de acordo com a latência de acesso à
memória, subdividindo os multiprocessadores em UMA e NUMA eos multicomputadores em
NORMA [19]:
• UMA - Uniform Memory Access(Acesso Uniforme à Memória): computadores nos quais
32
a latência de acesso à memória é igual (uniforme) para todos os processadores que com-
põem o sistema, independentemente da rede de interconexão utilizada ou de as transações
serem únicas ou paralelas. A implementação da memória é centralizada, ou seja, encontra-
se à mesma distância de todos os processadores.
• NUMA - Non-Uniform Memory Access(Acesso Não Uniforme à Memória): arquiteturas
nas quais a latência de acesso à memória não é uniforme. Cada processador é diretamente
ligado a um módulo de memória e acessa os demais módulos do sistema por meio dos
outros procesadores. Logo, a implementação da memória é distribuída e a latência de
acesso dependerá de o endereço referenciado estar ou não na memória diretamente ligada
ao processador que a referenciou.
• NORMA - Non-Remote Memory Access(Sem acesso a Variáveis Remotas): arquiteturas
nas quais cada processador é, na verdade, uma arquitetura inteira replicada e possui, por-
tanto, registradorres de endereçamento próprios que referenciam apenas sua memória lo-
cal.
Existe ainda um terceiro nível na hierarquia que classifica as máquinas de acordo com o trata-
mento da coerência decache. São quatro as variações descritas [19]:
• NCC-NUMA - Non-Cache-Coherent Non-Uniform Memory Access (Acesso Não Uni-
forme à Memória Sem Coerência decache): máquinas em que a coerência decachenão
é garantida porhardware;
• CC-NUMA - Cache-Coherent Non-Uniform Memory Access (Acesso Não Uniforme à
Memória Sem Coerência deCache): arquiteturas onde a coerência decacheé garantida
por hardware;
• SC-NUMA - Software-Coherent Non-Uniform Memory Access (Acesso Não Uniforme à
Memória Com Coerência deCacheem Software): É a alternativa utilizada para garantir
a coerência naquelas máquinas em que ohardwarenão garante, ou seja, nas arquiteturas
NCC-NUMA e NORMA. Essa camada de software, mais conhecida como DSM (Dis-
tributed Shared Memory), é transparente ao usuário;
33
• COMA - Cache Only Memory Architeture (Arquiteturas de memória somente comcache):
as memórias locais de cada processador são todas do tipocache, porém com mais capaci-
dade do que as tradicionais. Existe umhardwarede suporte que integra a gerência das
cachese de memória.
A Figura 3.3 apresenta uma visão geral dessa classificação baseada em compartilhamento de
memória, onde a linha tracejada indica que as classes NCC-NUMA e NORMA podem ser trans-
formadas em máquinas SC-NUMA através da inclusão da camada de software que implementa
coerência decache[19].
Figura 3.3: Classificação segundo o compartilhamento de memória. Fonte: [Hwa98] apud [19].
3.4.2 Classificações Comerciais
Historicamente, as fabricantes de computadores têm criadonomenclaturas próprias para
diferenciar suas tecnologias das concorrentes ou descrever diferentes abordagens adotadas para
tratar problemas antigos. Uma vez que essas nomenclaturas tornam-se populares, pode-se dizer
que formam uma classificação alternativa, embora freqüentemente os critérios utilizados não se-
jam totalmente claros. Em [15] os autores reapresentam critérios, previamente publicados dois
anos antes, para diferenciar os tipos de clusters existentes de acordo com a classificação ado-
tada na lista dos 500 computadores mais rápidos [6]. Além disso, em [15] os autores propõem
34
também uma classificação acadêmica alternativa, a qual, entrentanto, não foi adotada até o mo-
mento. Dentre as classificações comerciais, as mais utilizadas são:
• PVP - Parallel Vector Processing (Processadores VetoriaisParalelos): são arquiteturas
com poucos processadores, os quais são do tipo vetorial e sãodesenvolvidos especifica-
mente para estas máquinas. Em geral não possuem memóriacache, os blocos de memória
são entrelaçados e a interconexão é feita com redes de alta vazão não bloqueantes, per-
mitindo acessos paralelos à memória. O espaço de endereçamento é único, ou seja, esses
processadores podem ser enquadrados na categoria UMA da classificação anterior [19];
• SMP - Symmetric Multiprocessors (Multiprocessadores Simétricos): máquinas constituí-
das por processadores de prateleira (of the shelf). O termo simétrico refere-se ao fato
de que todos possuem privilégios equivalentes de acesso à memória e ao barramento,
diferentemente dos sistemas com processadores mestre-escravo. Os processadores de
prateleira costumam utilizar amplamente memóriascachese possuem memória compar-
tilhada, à qual se conectam através de um barramento de alta velocidade. Assim, o espaço
de endereçamento é único e essas máquinas podem ser classificadas como UMA [19].
Atualmente, essas arquiteturas tem se popularizado tambémentre os usuários domésticos
devido aos processadores multi-core;
• MPP - Massively Parallel Processors (Processadores Maciçamente Paralelos): podem ser
vistas como arquiteturas alternativas àquelas do tipo PVP.A diferença é nos MPPs são
utilizados milhares de processadores (daí o termo maciçamente) de prateleira conectados
por uma rede proprietária, ao invés de processadores específicos. Os espaços de endereça-
mento são próprios para cada nó, logo, se enquandram na classe das máquinas NORMA
[19];
• NOW - Network of Workstations (Rede de Estações de Trabalho): são sistemas forma-
dos por várias estações de trabalhos interligadas através de uma rede tradicional. Nessas
arquiteturas a execução da aplicação paralela não são o objetivo principal e diferem em
relação aos MPPs em três fatores: apresentam hierarquia de barramento nas estações, pre-
35
sença de um disco local nos nós e utilizam redes de interconexão de menor custo e latência
mais alta. São máquinas NORMA [19];
• COW - Cluster of Workstations (Máquinas Agregadas): arquiteturas projetadas especi-
ficamente para aplicações paralelas, diferentemente das NOWs, porém, também consti-
tuídas por várias estações de trabalho. Nesse caso as estações utilizadas como nós de
processamento não possuem alguns periféricos (monitor, mouse, teclado etc.) e o sistema
operacional é otimizado para esse fim. Pode-se dizer que COWssão NOWs dedicadas ao
processamento paralelo, uma evolução dessas. Duas estratégias diferentes são adotadas
no desenvolvimento de COWs, a primeira é interligar milhares de estações de trabalho
através de redes de baixo custo, enquanto a segunda utiliza redes de baixa latência para
ligação de um número menor de processadores mais poderosos.A primeira resulta em
máquinas NORMA, enquanto a outra abordagem pode resultar emmáquinas NUMA se
utilizadas placas de rede que implementam espaço único de endereçamento, embora, no
geral, também sejam do tipo NORMA [19].
Em [15] os autores classificam as COWs construídas com a primeira abordagem como Com-
modity Cluster, subdividindo-as em Constellations, quando o número de processadores em um
nó é maior do que o número de nós da arquitetura, e clusters-NOW, quando o número de nós
é maior do que o número de processadores em cada nó. O argumento utilizado é de que essa
diferenciação influencia na forma como o sistema será programado, pois uma Constellation pode
ser melhor aproveitada com o uso de múltiplos threads, enquanto clusters-NOW necessitam de
programação por troca de mensagens. Um caso especial de commodity cluster são os chamados
Clusters Beowulf [3] em que, além de ser construído somente com equipamentos de prateleira,
a arquitetura opera com sistema operacional Linux [15, 3]. Oconceito de Cluster Beowulf teve
grande importância no sentido de abrir novos caminhos, porém, atualmente é uma nomenclatura
em desuso. As classificações Constellation e cluster-NOW embora utilizadas na lista dos 500
computadores mais rápidos [6] também não são muito comuns.
36
3.5 Paradigmas de programação paralela
As aplicações paralelas podem ser classificadas em alguns paradigmas bem definidos de
programação, sendo que existe um número pequeno desses paradigmas que é usado na maioria
dos programas. Cada paradigma é uma classe de algoritmos quepossui a mesma estrutura de
controle. A escolha do paradigma a ser utilizado deve basear-se no tipo de recursos disponíveis
e de paralelismo inerente ao problema que será tratado. No que trata dos recursos, deve-se,
basicamente, definir o nível de granularidade que pode ser explorada eficientemente no sistema.
Já em relação ao tipo de paralelismo, deve-se observar a estrutura dos dados e da aplicação. Se
existir paralelismo em relação aos dados, pode-se executarem paralelo as mesmas operações
sobre diferentes partes dos dados, enquanto no caso de o paralelismo ocorrer nas operações, mas
não nos dados, pode-se implementar paralelamente a ação de trechos diferentes da aplicação sob
os mesmos dados [10].
As aplicações de Computação Distribuída em geral se baseiamno paradigma cliente/servidor,
no qual os processos comunicam-se entre si, por RPCs (Remote Procedure Calls), por exemplo,
mas não há paralelismo inerente à aplicação. Dado que o foco deste trabalho é em aplicações
paralelas, esse tipo de paradigma não será abordado. Não existe um consenso absoluto na liter-
atura para classificação dos paradigmas paralelos. Em [10] são referenciadas algumas das prin-
cipais classificações e os paradigmas mais utilizados são brevemente descritos. Tais paradigmas
são:Task-Farming, ouMaster/Slave; Single Program Multiple Data(SPMD);Data Pipelining;
Divide and Conquer; Speculative Parallelism[10].
3.5.1 Task-Farming (ou Master/Slave)
O paradigamaTask-farming, também conhecido comoMaster/Slave, ou, ainda, Mestre Es-
cravo, baseia-se na existência de duas entidades, um mestree diversos escravos. O processo
mestre é responsável por decompor o problema em tarefas menores e distribuí-las entre uma
"fazenda"de processos escravos e no final recolher os resultados. Os processos escravos exe-
cutam um ciclo bastante simples: recebem uma mensagem com a tarefa, processam a tarefa e
enviam o resultado para o mestre. Normalmente, a comunicação se dá apenas entre o mestre e
37
os escravos [10].
Existem duas abordagens para esse paradigma, uma delas é fazer a carga das tarefas de
forma estática e a outra dinâmica. Na estática, a distribuição das tarefas é feita no início da
execução, deixando o mestre livre para participar da computação. Já a abordagem dinâmica é
mais adequada para casos em que o número de tarefas é maior do que o número de processadores,
quando o número de tarefas é desconhecido no início da aplicação, quando o tempo de execução
não é previsível ou quando tratando de problemas não balanceados. Essa abordagem implementa
mais facilmente a tolerância a falhas, pois possibilita quena queda de algum processo sua carga
de trabalho seja distribuído para outro [10].
O paradigmaMaster/Slavepode alcançar altos índices despeedupe um bom grau de escal-
abilidade. Entretanto, para um número grande de processos ofato de utilizar o controle central-
izado em um único processo mestre pode torná-lo um gargalo nodesempenho da aplicação. Por
outro lado, é possível contornar essa limitação com a utilização de múltiplos processos mestres
[10].
3.5.2 Single Program Multiple Data (SPMD)
O SPMD é o paradigma mais utilizado, nele cada processo executa o mesmo trecho de
código porém em partes diferentes dos dados, o que envolve a decomposição da aplicação para os
processadores disponíveis. Esse paralelismo é também encontrado na literatura com os nomes:
paralelismo geométrico, decomposição do domínio ou paralelismo de dados.
Problemas como os apresentados na Seção 2.4 possuem uma estrutura geométrica regular
com interações limitadas espacialmente. Essa homogeneidade permite que esses problemas
tenham seus decompostos e distribuídos uniformemente entre os processadores, ficando cada
um responsável por uma área espacial definida. Os processos comunicam com os processos
vizinhos e a carga de comunicação é proporcional ao limite dos elementos, enquanto a carga da
computação é proporcional ao volume dos elementos. Em alguns casos é necessário fazer uma
sincronização global periódica dos processos, sendo que o padrão de comunicações em geral é
bastante regular e previsível.
Aplicações construídas em SPMD são bastante eficientes quando os dados são bem distribuí-
38
dos e o sistema é homogêneo. Uma desvantagem do paradigma é que ele é bastante suscetível a
deadlocks, dependendo da situação, a queda de um processo pode ser suficiente para ocorrência
dedeadlock[10]
3.5.3 Data Pipelining
É um paradigma de granularidade mais fina que os dois anteriores. Baseia-se em uma abor-
dagem de decomposição funcional, sendo um dos paradigmas dedecomposição funcional mais
simples e populares. É bastante utilizado em aplicações de processamento de imagens. Seu
funcionamento é análogo aopipelinede um processador, cada processador da máquina paralela
executa pequenas tarefas paralelizáveis do algoritmo [10].
O padrão de comunicações doData Pipelinigem geral é bastante simples, e pode ser total-
mente assíncrona, pois o fluxo da execução se dá através de estágios adjacentes dopipeline. A
eficiência do paradigma é diretamente proporcional à possibilidade de balanceamento da carga
de execução dentre os estágios dopipeline[10].
3.5.4 Divide and Conquer
O paradigmaDivide and Conquer, ou Dividir e Conquistar, é bastante conhecido na pro-
gramação seqüencial. Consiste em dividir um problema em dois ou mais problemas menores,
computá-los e ao final reunir o resultado dessas computaçõespara obter o resultado do problema
inicial. Na abordagem paralela doDivide and Conquer, os subproblemas podem ser resolvidos
simultaneamente atribuindo-se cada um a um processador e não há necessidade de comunicação
entre eles. As tarefas de dividir o problema e reunir os resultados também podem fazer uso
de algum nível de paralelismo, embora nesse caso seja necessária comunicação. O paradigma
Master/Slave pode ser visto como uma modificação desse do Divide and Conquer [10].
3.5.5 Speculative Parallelism
É um paradigma empregado quando o uso dos anteriores não é possível. Alguns proble-
mas possuem dependências de dados complexas, o que reduz as possibilidades de execução
39
paralela. Nesses casos, uma solução é executar o problema empequenas partes, porém, uti-
lizando algum nível de especulação e supondo-se uma execução otimista, ou seja, executar as
operações em paralelo esperando que não haja violação da consistência dos dados. Outro uso
desse paradigma é na execução de diferentes algoritmos pararesolver um mesmo problema, de
modo que o primeiro a terminar é escolhido como solução [10].
3.6 Conclusão
A computação paralela tem alcançado grandes avanços nas últimas décadas. Porém, existe
ainda uma série de questões abertas em relação a padronizações de modo geral. Dessa forma,
é importante observar que as escolhas da arquitetura, modelo e paradigma de programação etc.,
no processo de implementação de uma solução baseada em paralelismo devem ser feitas com
base em uma análise dos requisitos particulares da aplicação em questão. Nessa análise devem
ser considerados fatores como granularidade que o algoritmo permite, custos da arquitetura,
necessidades de escalabilidade, se o paralelismo tem intuito apenas de oferecer alto desempenho
ou também tolerância a falhas, dentre outros.
40
4 Computação Paralela aplicada os
Métodos Numéricos
Ao embarcar no estudo dos métodos numéricos paralelos e distribuídos é importante refle-
tir sobre suas diferenças em relação aos métodos seqüenciais. Existem muitas questões rela-
cionadas à paralelização que não aparecem quando tratando os mesmos algoritmos de forma se-
qüencial. A primeira delas é a alocação de tarefas que consiste na quebra do total do trabalho em
tarefas menores, as quais serão atribuídas a diferentes processadores, e o correto seqüenciamento
das mesmas quando houver interdependência entre elas. A segunda questão é a comunicação
dos resultados da computação entre os processadores, deve-se buscar manter a comunicação tão
eficiente quanto possível e estimar o impacto da mesma na performance do sistema. A terceira
é a sincronização entre os processadores. Em alguns métodos, chamados de métodos síncronos,
os processadores devem esperar pontos pré-determinados para completar certas operações ou
utilizar certos dados, sendo que o mecanismo utilizado paratal pode ter efeitos importantes na
performance. Na outra abordagem, os chamados métodos assíncronos, não há espera em pontos
pré-determinados e nesse caso as implicações correspondentes devem ser cuidadosamente avali-
adas para assegurar a validade do método. Existem ainda outras questões, as quais relacionam-se
aos assuntos tratados no capítulo anterior, como o desenvolvimento de medidas apropriadas de
performance para os métodos paralelos e os efeitos da arquitetura do sistema sobre elas [8].
Na Seção 2.3 foram apresentadas as diferenças entre os métodos diretos e iterativos para
solução dos sistemas lineares. Entretanto, uma análise mais detalhada revela que, na reali-
dade, não existe uma linha totalmente clara dividindo essasduas classes. Isso porque, por ex-
41
emplo, uma técnica aplicada para melhorar os resultados obtidos com a eliminação de Gauss
pura é aplicar um refinamento iterativo, o que pode ser visto como um método iterativo que
utiliza Eliminação de Gauss como pré-condicionador. Por outro lado, costuma-se utilizar os
pré-condicionadores para acelerar os métodos iterativos,o que, em geral, envolve a solução
direta de sistemas próximos do original. Normalmente, esses sistemas próximos são tais que
permitam uma solução mais econômica ou que o processo iterativo possa explorar melhor
o paralelismo. Assim, a diferença maior entre a solução do sistema original e de um pré-
condicionador com métodos diretos é que, no original é necessário reorganizar a computação
ou o algoritmo para obter um melhor paralelismo, enquanto com os pré-condicionadores pode-
se simplesmente descartar as entradas do sistema aproximado se isso melhorar o paralelismo
[16].
Este capítulo apresenta a aplicação dos conceitos de paralelismo, anteriormente discutidos,
nos métodos iterativos clássicos para solução de sistemas lineares. Visando introduzir as no-
tações, inicialmente é apresentada a paralelização de alguns algoritmos simples numéricos sim-
ples. Na seqüência é abordada a paralelização dos métodos iterativos e, por fim, a implemen-
tação voltada às arquiteturas que se comunicam por troca de mensagens dessas paralelizações.
4.1 Paralelização de algoritmos numéricos simples
As paralelizações a seguir tratam de tarefas bastante elementares, porém, muito comuns em
todos os métodos numéricos. A representação dos algoritmosé feita através de DAGs e na
discussão dos mesmos é assumido que as operações de adição e multiplicação demoram uma
unidade de tempo e que os processadores são capazes de trocarinstantaneamente resultados in-
termediários, ou seja, os atrasos relativos a comunicação são ignorados. Contudo, os algoritmos
apresentados são bastante simples e, portanto, podem ser implementados em diversas arquite-
turas assumindo-se ooverheadde comunicação como sendo insignificante [8].
42
4.1.1 Adição de escalares
A mais simples das tarefas computacionais é a adição den escalares. Evidentemente o
melhor algoritmo seqüencial requern − 1 operações. Assim,T ∗(n) = n − 1. Assumindo
n como uma potência de2, para uma implementação paralela, pode-se dividir osn escalares
emn/2 pares disjuntos e utilizarn/2 diferentes processadores para somar cada um dos pares.
Assim, após uma unidade de tempo, a tarefa de somarn/2 escalares estaria completa e após log
n estágios, a computação estaria completa. A Figura 4.1 ilusta o DAG desse algoritmo para soma
de16 escalares com8 processadores. O algoritmo é facilmente generalizado paraos casos onde
n não é potência de2, sendo que o tempo de execução passa a ser logn paran/2 processadores
[8].
Figura 4.1: DAG para soma de 16 escalares utilizando 8 processadores.Fonte: [8].
A eficiência desse algoritmo é dada porn−1(n/2)(logn) , a qual tende azero conformen cresce.
Um algoritmo paralelo alternativo pode ser obtido da seguinte forma: assumindo, por simplici-
dade, que logn é um inteiro e quen/logn é um inteiro potência de2, divide-se osn números em
n/logn grupos de logn números cada. Utiliza-sen/logn processadores e oi-ésimo processador
efetua a soma dos números noi-ésimo grupo, operação executada em tempo logn − 1. Dessa
forma, efetua-se a soma den/logn números, o que poderia ser alcançado no algoritmo anterior
com tempolog(n/logn) ≤ logn, usandon/(2logn) processadores. Esse algoritmo de duas
fases requer um tempo aproximadamente igual a2logn, ou seja, a velocidade é reduzida em um
fator de 2, porém utiliza apenasn/logn processadores e, dessa forma, sua eficiência aproxima-
43
se de1/2. O número de processadores escolhido é aproximadamente igual aT1(n)/T∞(n) pelo
fato de que essa escolha tende a levar a algoritmos eficientes. Assim, percebe-se que uma pe-
quena redução na velocidade pode melhorar consideravelmente a eficiência, pois a diminuição
no número dos processadores pode decrementar substancialmente a quantidade de comunicação
do algoritmo. O DAG do algoritmo alternativo é apresentado na Figura 4.2 [8].
Figura 4.2: DAG para soma de 16 escalares com 4 processadores. Fonte: [8].
4.1.2 Produto Interno
O produto interno∑n
i=1 xiyi de dois vetoresn-dimensionais pode ser calculado em tempo
log n+1 empregandon processadores da seguinte forma: em um primeiro passo, cadaproces-
sadori calcula o produtoxiyi e, então, aplica-se o algoritmo de adição de escalares de com
tempo logn anteriormente descrito [8].
4.1.3 Adição e Multiplicação de Matrizes
A soma den matrizes de dimensõesm x mpode ser calculada em tempo logn usando
m2[n/2] processadores empregando-se um grupo diferente den/2 processadores calculando
uma diferente entrada da soma. De maneira análoga, a multiplicação de duas matrizesm x nen
x l consiste na avaliação deml produtos internos de vetoresn-dimensionais e pode ser, portanto,
alcançada no tempo logn+1 empregando-senml processadores. No caso em quen = m = l,
44
tem-se o requisito den3 processadores. O número correspondenteT1(n) de operações éθ(n3)
[8].
4.1.4 Potência de uma Matriz
Supondo queA é uma matrizn x na qual deseja-se calcularAk para algum inteirok. Sek
é uma potência de2, o cálculo pode ser realizado através da repetição de operações de elevação
ao quadrado, computando-se a primeiraA2, em seguidaA2A2 = A4 etc. Após logk estágios,
tem-se o valor deAk. Esse procedimento envolve logk multiplicações consecutivas de matrizes
e pode, dessa forma, ser obtido no tempo logk([log n]+1) empregando-sen3 processadores.
Uma modificação simples desse procedimento permite calcular Ak em tempoθ(log k.log n)
mesmo quandok não é potência de2 [8].
Portanto, todas as potênciasA2, . . . , An de uma matrizn x npodem ser calculadas em tempo
θ(log2n) usandon4 processadores bastanto para isso empregar um grupo diferente den3 pro-
cessadores para a computação de cada potênciaAk. Um método alternativo para o cálculo das
potênciasA2, . . . , An, o qual elimina duplicações desnecessárias de esforços computacionais,
é mostrado na Figura 4.3. Os nós rotulados com "S"representam uma operação de elevação ao
quadrado da matriz. No primeiro estágioA2 é calculada, no segundoA3 e A4 são obtidas pela
multiplicação das matrizes anteriores. Generalizando-se, nok-ésimo estágio, tem-se as matrizes
A2k−1+1, . . . , A2k
. Assim,θ(log n) estágios são suficientes para calcularA2, . . . , An, sendo que
cada estágio envolve no máximoθ(n) multiplicações simultâneas de matrizes e, empregando-se
n4 processadores, pode ser alcançado em tempoθ(log n), levando a um tempo médioθ(log2n)
[8].
Os algoritmos previamente discutidos foram desenvolvidoscom intuito de obter os tempos
de execução mais rápido possíveis, porém, essa abordagem freqüentemente conduz a requisitos
excessivos de processadores e implica em baixa eficiência. Por outro lado, conforme discutido
anteriormente, pode-se tornar os mesmos algoritmos eficientes escolhendo-se um número de
processadores tal quep = O(T1(n)/T∞(n)). O produto duas matrizesn x n, por exemplo, pode
ser computado em tempoθ(log n) empregandoθ(n3/log n) processadores, sendo a eficiência
correspondenteθ(1). Reduzindo-se o número de processadores ainda mais, o tempode exe-
45
Figura 4.3: Computação paralela da potência de uma matriz n xn. Fonte: [8].
cução passa a serθ(T1(n)/p). Assim, duas matrizesn x npodem ser multiplicadas em tempo
θ(n) quandon2 processadores são empregados e em tempoθ(n2) quandon processadores são
utilizados [8].
4.2 Paralelização dos Métodos Iterativos para resolução deSELAs
A seguir é discutida a paralelização de alguns dos métodos clássicos para resolução de SE-
LAs. A notação empregada neste capítulo é mais próxima das linguagens de programação, e,
portanto, ligeiramente diferente daquela utilizada em 2.3.2. A equação 2.8 é então reescrita da
seguinte forma [8]:
x(t + 1) = f(x(t)), t = 0, 1, . . . , (4.1)
onde cada~x(t) é um vetor den dimensões ef é uma função do tipoR → R que equivale ao
lado direito da Equação 2.8.
4.2.1 Paralelização do Método de Jacobi
Sejaxi(t) o i-ésimo componente dex(t) e sejafi o i-ésimo componente da funçãof , então
x(t + 1) = f(x(t)) pode ser escrito na forma [8]:
46
xi(t + 1) = fi(x1(t), . . . , xn(t)), i = 1, . . . , n (4.2)
O algoritmox := f(x) pode ser paralelizado alocando para cada um dosn processadores
a tarefa de atualizar um componente diferente dex, conforme na Equação 4.2. A cada estágio,
o i-ésimo processador conhece o valor de todos os componentes de x(t) dos quaisfi depende,
calcula o novo valorxi(t + 1), e comunica-o para os outros processadores de modo a começar
a próxima iteração. A comunicação necessária para a execução de uma iteração pode ser com-
pactamente descrita em termos de um grafo dirigidoG = (N,A) chamadografo de dependên-
cia. O conjunto de nósN é 1, . . . , n, correspondente aos componentes dex. Para quaisquer
dois nós distintosi e j, existe um arco(i, j) do grafo de dependência se e somente se a função
fi depende dexi, ou seja, se e somente se o processadori precisa comunicar os valores dexi(t)
para o processadorj [8].
A Figura 4.4 (a) exemplifica um grafo de dependência para as iterações:x1(t + 1) =
f1(x1(t), x3(t)); x2(t+1) = f2(x1(t), x2(t)); x3(t+1) = f3(x2(t), x3(t), x4(t)) ex4(t+1) =
f4(x2(t), x4(t)). Assumindo-se que a iteração deve ser realizada apenas parat = 0, 1, . . . , T ,
ondeT é um número positivo inteiro, a estrutura do algoritmo pode ser representada em ter-
mos de um DAG, o qual é ilustrado em 4.4 (b) e corresponde a duasiterações quando a função
f tem o grafo de dependência definido em 4.4 (a). Os nós do DAG sãona forma(i, t), onde
i ∈ 1, . . . , n e t é o contador de iterações. Os arcos são da forma((i, t), (j, t+1)), onde(i, j)
é um arco do grafo de dependência oui = j [8].
Uma paralelização de granularidade grossa da iteraçãox := f(x), pode ser obtida decompondo-
se o espaço vetorialRn em um produto cartesiano de subespaços de dimensões menoresRnj ,
j = 1, . . . , p onde∑p
j=1 nj = n. Da mesma forma, qualquer vetorx ∈ Rn é decomposto
comox = (x1, . . . , xj , . . . , xp), onde cadaxj é um vetornj-dimensional, o qual é chamado
de bloco componentede x, ou simplesmentecomponente. De maneira análoga, a iteração
x(t + 1) = f(x(t)) pode ser escrita como [8]:
xj(t + 1) = fj(x(t)), j = 1, . . . , p (4.3)
onde cadafj é uma função vetorialRn → Rnj . Atribuindo-se cada um dos diferentes
47
Figura 4.4: Em (a) Grafo de dependências e (b) DAG para funçãocom o grafo de dependências
definido por (a).Fonte: Adaptação de [8].
blocos-componentes para um processadorp, de acordo com a equação 4.3, diz-se que o algo-
ritmo está paralelizado em bloco. Um grafo de dependênciaG = (N,A) pode ser novamente
definido comN = 1, . . . , p e A = (i, j)/fjdepende dexi. Existem diversas razões pelas
quais a paralelização em blocos é interessante. A primeira éque podem haver poucos proces-
sadores disponíveis, nesse caso, basta que se atribua mais de um componente para cada pro-
cessador. A segunda é que algumas funções escalaresfi podem envolver cálculos em comum,
nesse caso é natural agrupa-las. Por fim, a paralelização em blocos reduz a necessidade de
comunicação do algoritmo [8].
No geral, uma paralelização que atribui a atualização de diferentes componentes para difer-
entes processadores é significativa quando os cálculos envolvidos na atualização de cadaxi são
diferentes, porém, pode ser ruim quando forem o semelhantes. Por exemplo, supondo que cada
fi tem a forma da equação a seguir [8]:
fi(x1, . . . , xn) = xi + (
n∑
j=1
x2j)
1/2 (4.4)
Nesse exemplo, é desnecessário que todos processadores calculem simultaneamente o valor
de(∑n
j=1 x2j )
1/2. Isso pode ser feito por um único processador que então, na seqüência, comu-
nica o resultado para os demais. Uma abordagem ainda mais eficiente é que todos processadores
48
cooperem no cálculo do somatório, fazendo, por exemplo, cálculos na forma apresentada na
Seção 4.1.1. Entretanto, a avaliação de cadafi em muitos casos envolve muito pouca ou nen-
huma duplicação de esforços, sendo essas situações as que apresentam maior interesse [8].
4.2.2 Paralelização de Gauss-Seidel
Conforme visto na Seção 2.3.2, o algoritmo de Gauss-Seidel ao invés de atualizar todos os
componentes de~x simultaneamente, atualiza cada um em separado e utiliza os novos valores dos
outros componentes conforme os mesmos vão sendo calculados. A equação dos componentes
(Equação (2.12)) pode ser reescrita como [8]:
xi(t + 1) = fi(x1(t + 1), . . . , xi−1(t + 1), xi(t), . . . , xn(t)), i = 1, . . . , n (4.5)
Podem haver casos em que a iteração de Gauss-Seidel é completamente não paralelizável,
pois se cada funçãofi depender de todos os componentesxj , então somente um componente
pode ser atualizado por vez. Por outro lado, quando o grafo dedependência é esparso, é pos-
sível que a atualização de certos componentes seja feita em paralelo. Cabe observar que existem
diversos algoritmos correspondentes à mesma funçãof de Gauss-Seidel, dado que é possível
escolher a ordem com que os componentes são atualizados. Porexemplo, pode-se querer atu-
alizar os componentes de~x iniciando comxn e retroceder até quex1 seja atualizado no final.
Diferentes ordens de atualização correspondem a diferentes algoritmos e, em geral, a resultados
diferentes. Por outro lado, em diversas aplicações, um algoritmo de Gauss-Seidel converge em
um número de iterações para um mesmo valor, independente da ordem de atualização. Como a
velocidade de convergência do método correspondente às diferentes ordens de atualização não
apresenta grandes diferenças, torna-se natural escolher uma ordenação tal que o paralelismo em
cada iteração seja maximizado [8].
A Figura 4.5 apresenta duas paralelizações possíveis para Gauss-Seidel para o grafo de de-
pendência def da Figura 4.4 (a). No exemplo 4.5 (a), o DAG ilustra as dependências de dados
em uma iteração de Gauss-Seidel e a a forma do algoritmo éx1(t + 1) = f1(x1(t), x3(t));
x2(t + 1) = f2(x1(t + 1), x2(t)); x3(t + 1) = f3(x2(t + 1), x3(t), x4(t)); x4(t + 1) =
49
f4(x2(t + 1), x4(t)). Há quatro atualizações a serem feitas, mas a profundidade do grafo é
3, percebe-se assim quex3 e x4 podem ser atualizados em paralelo. Em 4.5 (b) é apresen-
tada uma opção em que o paralelismo é aumentado modificando-se a ordem das atualizações.
A ordem nesse caso é:x1(t + 1) = f1(x1(t), x3(t)); x3(t + 1) = f3(x2(t), x3(t), x4(t));
x4(t + 1) = f4(x2(t), x4(t)) ex2(t + 1) = f2(x1(t + 1), x2(t)). A profundidade passa a ser2
e a iteração pode ser feita em dois passos usando dois processadores [8].
Figura 4.5: DAG para possíveis iterações de Gauss-SeidelFonte: Adaptação de [8].
A seguir é apresentada uma formulação de teoria dos grafos para o problema de encontrar
uma ordem de atualização que minimize o tempo paralelo necessário para o algoritmo: Dado
um grafo de dependênciaG = (N,A), umacoloraçãode G usandoK cores, define-seh : N →
1, . . . ,K como um mapeamento que atribui uma "cor"k = h(i) para cada nói ∈ N . A idéia
é que variáveis de mesma cor sejam atualizados em paralelo. Oresultado a seguir demonstra
que maximizar o paralelismo é equivalente ao problema de "coloração ótima"[8].
Proposição 1: Os dois enunciados são equivalentes:
1. Existe uma ordenação das variáveis tal que uma iteração doalgoritmo de Gauss-Seidel
correspondente pode ser feita emK passos paralelos.
2. Existe uma coloração do grafo de dependências que usaK cores e com a propriedade de
que não existem ciclos positivos em que todos os nós do ciclo tenham a mesma cor.
50
Prova: Primeiro demonstra-se que (1) implica (2). Considerando-se uma ordenação das var-
iáveis com a qual uma iteração de Gauss-Seidel demoraK passos paralelos. Define-seh(i), a cor
do nói, sendo igual ak se a variávelxi é atualizada nok-ésimo estágio paralelo. Considerando
um ciclo positivoi1, i2, . . . , im com im = i1 e escolhendo o nóil no ciclo (1 ≤ l ≤ m) que
é o primeiro na ordem assumida, desde queil+1 é ordenado apósil, e que(il, il+1) ∈ A, a
variávelxil+1(t + 1) depende dexil(t + 1). Segue quexil e xil+1
não podem ser atualizados
simultaneamente e tem-seh(il) 6= h(il+1). Isso mostra que em cada ciclo positivo, existem dois
nós com cores diferentes, o que prova (2) [8].
Será provado agora um resultado auxiliar. Mostrar-se-á queseG é um DAG, então seus nós
podem ser ordenados de maneira que if(i, j) ∈ A, entãoj precedei. A prova é como segue.
Para cada nói, supõe-sedi como o número de arcos maior possível em um caminho positivo
que inicia emi (di = 0 se o nói não tem arestas de saída). Sabe-se quedi é finito como uma
conseqüência da aciclicidade. Ordena-se então os nós em ordem crescente de valoresdi (laços
entre nós com o mesmo valor dedi são quebrados arbitrariamente). Vê-se que se(i, j) ∈ A
entãodi > dj . Isso porque é possível selecionar o maior caminho iniciando emj e colocar o
arco(i, j) para obter um caminho ainda mais longo iniciando emi. Conclui-se quej precedei
sempre que(i, j) ∈ A, conforme desejado [8].
Assume-se agora a validade de (2). Sejah uma coloração comk cores e sem ciclos positivos
no qual todos os nós tenham a mesma cor. Para cada cork, sejaGk um subgrafo deG obtido
deixando somente os nós com a cork e os arcos juntando eles. CadaGk é acíclico e, de acordo
com o resultado do parágrafo anterior, os nós emGk podem ser ordenados de forma quej
vem antes dei sempre que(i, j) ∈ A. Ordena-se os nós emG em ordem crescente de cor;
laços com entre nós com a mesma cork são quebrados utilizando a ordenação do grafoGk.
Consderando-se a iteração de Gauss-Seidel correspondentea esta ordenação, sejami e j dois
nós distintos com a mesma cork, se(i, j) /∈ A e(j, i) /∈ A, entãoxi exj podem ser atualizados
em paralelo. O caso em que(i, j) ∈ A e (j, i) ∈ A é impossível porqueGk é acíclico. Se
(i, j) ∈ A e (j, i) /∈ A, entãoj precedei na ordem construída e, dessa forma, a computação de
xj(t + 1) somente precisa do valor dexi(t) e não do valor dexi(t + 1). Finalmente, o caso em
que(j, i) ∈ A e (i, j) /∈ A é similar. Conclui-se que todos osxi com a mesma cor podem ser
51
atualizados em paralelo, provando assim (1) [8].
Para o grafo de dependência da Figura 4.4 (a) duas cores são suficientes. Em particular,
é possível terh(1) = h(3) = h(4) = 1 e h(2) = 2. Desde que todos os ciclos possíveis
passem pelo nó2, não existem ciclos positivos com todos os nós com a mesma cor, conforme
requerido. Em particular, o subgrafoG1 no qual somente os nós com a cor1 são mantidos é
acíclico. Comdi definido como na prova da propriedade, anteriormente, tem-sed1 = 0, d3 = 1
e d4 = 2. A ordenação das variáveis construída na prova é(1, 3, 4, 2), e a iteração de Gauss-
Seidel correspondente é mostrada na Figura 4.2.2(b). De acordo com a propriedade acima, esta
ordenação requer o menor número possível de estágios paralelos por iteração. Além disso, a
proposição 1 pode ser simplificada nos casos em que o grafo de dependência tem a propriedade
de simetria [8].
Proposição 2: Supondo que(i, j) ∈ A se e somente se(j, i) ∈ A. Então, os enunciados a
seguir são equivalentes [8]:
1. Existe uma ordem das variáveis tal que uma iteração do algoritmo correspondente Gauss-
Seidel pode ser feita emK passos paralelos.
2. Existe uma coloração do grafo de dependências que usa no máximo K cores e tal que os
nós adjacentes tem cores diferentes, ou seja, se(i, j) ∈ A entãoh(i) 6= h(j).
Prova: é suficiente dizer que a condição (2) dessa proposiçãoé equivalente à condição (2) da
proposição 1. Supondo-se a existência de um ciclo positivo com todos os nós daquele ciclo com
a mesma cor, então existem dois nós adjacentes com a mesma cor. Inversamente, se(i, j) ∈ A,
então(j, i) ∈ A e os dois arcos(i, j) e (j, i) formam um ciclo positivo. Assim, se dois nós
vizinhos tem a mesma cor, então existe um ciclo positivo em que todos os nós do ciclo têm a
mesma cor, o que prova a equivalência das duas condições e conclui a prova [8].
Infelizmente, o problema da coloração ótima acima descritoé computacionalmente intratável
(NP-completo). Entretanto, na prática, os problemas freqüentemente têm uma estrutura especial
e coloração com relativamente poucas cores, pode ser as vezes encontrada por observação. Por
exemplo, se todos os nós emG tem no máximoD vizinhos, entãoD + 1 cores são suficientes.
Assumindo-se que o primeiro nói é colorido, seD + 1 cores são utilizadas, então existe uma
52
cor passível de ser usada para o nói + 1 e que assegura que sua coloração é diferente de seus
vizinhos já coloridos [8].
Quando um esquema de coloração é empregado para implementação paralela de um algo-
ritmo de Gauss-Seidel, não há necessidade de atribuir um processador diferente para cada com-
ponente de~x, pois cada processador se mantém parado enquanto as variáveis de cores diferentes
são atualizadas. Uma medida para contornar essa limitação éutilizar menos processadores, com
cada processador sendo responsável por diversas variáveisque possuem cores diferentes associ-
adas [8].
Exemplo: Coloração Rubro-negra de uma matriz
Nos algoritmos iterativos para solução de equações diferenciais parciais, por exemplo, cada
componente de~x é associado com um ponto particular de uma certa região do espaço bi-
dimensional. SejaN o conjunto de todos os pontos(i, j) ∈ R2, tal quei e j são inteiros
que satisfazem0 ≤ i ≤ M e 0 ≤ j ≤ M . Sejaxij o componente do vetor~x correspondente
ao ponto(i, j). Conectando-se os vizinhos mais próximos forma-se um grafoG = (N,A),
como ilustrado na Figura 4.2.2(a). Pode-se verG como um grafo dirigido fazendo os arcos bi-
direcionais, assumindo que esse é o grafo de dependência associado a iteraçãox := f(x). A
execução paralela dessa iteração no método de Jacobi é direta, basta atribuir-se um processador
diferente para cada ponto(i, j). Esse processador é responsável por atualizarxij e, para fazê-lo,
somente precisa saber os valores dos componentes de~x associados com pontos vizinhos. Assim,
é mais natural assumir que os processadores associados são ligados por links de comunicação
diretos [8].
Para implementação do algoritmo relacionado ao Método de Gauss-Seidel, pode-se colorir
o grafo da Figura 4.2.2(a) com duas cores, conforme indicadona Figura 4.2.2(b). Atribuindo-
se um processador para cada componentexij , cada processador permanece parado metade do
tempo. Assim, é razoável atribuir dois componentes de coresdiferentes para cada processador.
Como pode ser visto na Figura 4.2.2(b), isso pode ser feito sem perder a propriedade de que
somente os vizinhos mais próximos se comunicam uns com os outros. Na prática, o número
de pontos envolvido é freqüentemente grande o suficiente de modo que para cada processador
53
sejam atribuídos mais do que dois componentes de~x. A coloração indicada na Figura 4.2.2(b)
é comumente conhecida como coloração rubro-negra e o algoritmo Gauss-Seidel associado é
conhecido como Gauss-Seidel com ordenação rubro-negra [8].
Figura 4.6: Grafo para iterações Gauss-Seidel com coloração. Fonte: Adaptação de [8].
4.3 Implementação da paralelização dos métodos iterativosclássi-
cos
A seguir, discute-se a implementação das paralelizações apresentadas na Seção?? para os
algoritmos iterativos clássicos (Seção 2.3.2). Embora atéeste momento os algoritmos tenham
sido abordados de maneira independente da arquitetura utilizada, nesta seção o foco é voltado
às arquiteturas que se comunicam por troca de mensagens uma vez, naquelas em que a comuni-
cação ocorre por memória compartilhada, a implementação tende a ser mais simples [8].
O método de Jacobi para solução deA~x = ~y é diretamente paralelizável, pois cada in-
teração envolve uma multiplicação de matriz e vetor. Supondo-se que hajamn processadores
disponíveis, e que oi-ésimo processador é responsável por calcularxi(t) a cada iteraçãot (o
caso em que o número de processadores é menor do que o número devariáveis é semelhante).
Supondo que oi-ésimo processador conhece as entradas dai-ésima linha deA. Para calcular
xi(t + 1), o processadori precisa saber os valores dexj(t) obtidos na iteração prévia pelos
54
processadoresj para os quaisaij 6= 0. Se a maioria das entradas deA são diferentes dezero, é
fácil transmitirxj para todos os processadoresi, mesmo queaij seja igual azero, pois transmis-
sões seletivas podem introduzir umoverheaddespropositado. Trata-se, portanto de um broad-
cast multi-nós. Na prática, essa abordagem pode resultar emum tempo de execução dominado
pelo tempo de comunicação, nesse caso, deveriam haver menosprocessadores com mais com-
ponentes atribuídos para cada um. Uma implementação alternativa, é considerar que existem
novamenten processadores, com oi-ésimo processador responsável peloi-ésimo componente
xi. Entretanto, assume-se que oi-ésimo processador tem acesso ai-ésima coluna deA, ao con-
trário dai-ésima coluna deA. Assumindo queA é uma matriz densa, a computação acontece
da seguinte maneira: cada processadori avaliaaijxi paraj = 1, . . . , n. Então para cadaj, as
quantidadesajixi parai = 1, . . . , n são propagadas para o processadorj, com somas parciais
formadas ao longo do caminho, o qual é uma acumulação multi-nós. Os requisitos de comuni-
cação desse método e do anterior são os mesmos, para o caso de sistemas densos. A preferência
por uma das duas abordagens pode depender da arquitetura do computador que será utilizado.
Considerando a matrizA como sendo uma matriz esparsa, a estrutura esparsa deA determina
que um grafo dirigidoG = (N,A), ondeN = 1, . . . , n e o conjunto de arcosA é o conjunto
(i, j)/i 6= jandaji 6= 0 de todos os processos pares(i, j) tais quei precisa comunicar comj.
Utilizando o grafo de dependência introduzido em 4.2. Dada uma arquitetura especial, obtém-se
implementações paralelas eficientes se o grafo de dependência puder ser encaixado no grafo que
descreve a topologia da interconexão de rede. Nesse caso, todas as comunicações acontecem
entre processadores vizinhos e as penalidades de comunicação são mínimas. Quando tal encaixe
não pode ser encontrado, os requisitos de comunicação do método de armazenamento da matriz
são bastante diferentes, o que é uma diferença importante para o caso de matrizes densas [8].
Os algoritmos de Gauss-Seidel e SOR não são naturalmente ideais para paralelização, no
geral. Isso pode ser constatado, por exemplo, considerandoo algoritmo de Gauss-Seidel e
supondo uma matrizA totalmente densa, ou seja,aij 6= 0 para todoi e j. Logo, para que o
processadori calculexi(t + 1), o valor dexj(t + 1) será necessário para todoj < i. Isso torna
o algoritmo inerentemente seqüencial dado que dois componentes dex não podem ser atualiza-
dos simultaneamente. Tal situação é particularmente inconveniente, pois os algoritmos SOR,
55
quando utilizando um valor adequado paraω, convergem, freqüentemente, de maneira muito
mais rápida que os demais. Por outro lado, como apresentado em 4.2, essa dificuldade pode ser
contornada aplicando-se o esquema de coloração nos casos emque a matrizA é esparsa. Feliz-
mente, quando a matrizA é obtida a partir da discretização de uma equação diferencial parcial,
ela é sempre esparsa. Além disso, na prática, o númerop de processadores é freqüentemente
muito menor do que o númeron de variáveis. Logo, diversas variáveis podem ser atribuídas
para cada processador [8].
4.4 Conclusão
Neste capítulo foram consideradas possíveis paralelizações para os métodos iterativos clás-
sicos para resolução de SELAs. Fica evidente que, mesmo havendo restrições e limitações
teóricas, esses métodos podem ser eficientemente implementados para execução em máquinas
de alto desempenho diversas dado que, na prática, os sistemas lineares resultantes das modela-
gens e discretizações de equações não se enquadram nos requisitos que tornam a paralelização
inviável. A granularidade que arquitetura permite explorar com esses algoritmos influencia di-
retamente na eficiência dos métodos paralelos, pois é ela o principal fator para determinar se a
movimentação dos dados, ou seja, operações de leitura e escrita, requerida por um dado método
terá um custo compatível com o método e SELA em questão.
56
5 Bibliotecas e pacotes para
resolução numérica de SELAs com alto
desempenho
Neste Capítulo serão 0apresentadas as bibliotecas, pacotes e softwares mundialmente mais
utilizados para problemas de álgebra linear, dentre eles, aresolução de sistemas lineares.
5.1 BLAS – Basic Linear Algebra Subprograms
BLAS (Basic Linear Algebra Subprograms) é uma coleção de rotinas que proveêm as bases
para implementar programas que façam operações com matrizes e vetores. A especificação orig-
inal da BLAS foi desenvolvida na forma de subprogramas em Fortran 66 e na seqüência Fortran
77 [23, 1]. As funcionalidades da biblioteca são divididas em três níveis: no primeiro estão as
funcionalidades que executam operações de escalares com escalares, escalares com vetores e
vetores com vetores; no nível dois encontram-se as operações entre matrizes e vetores; e o nível
3 implementa operações de matrizes com matrizes. O nível 1 foi o primeiro a ser desenvolvido,
foi o resultado de um projeto colaborativo entre 1973 e 1977.Com o surgimento das máquinas
vetoriais, máquinas com hierarquias de memória e máquinas paralelas com memória compartil-
hada, as especificações para os níveis 2 e 3 da BLAS foram desenvolvidos de 1984 à 1986 e de
1987 à 1988, respectivamente. Originalmente, a BLAS oferecia rotinas apenas para trabalhar
com matrizes densas e do tipo banda, atualmente as matrizes esparsas são também suportadas
57
[1].
Muitos dos algoritmos freqüentemente utilizados na álgebra linear numérica podem ser cod-
ificados de forma que o volume da computação seja feito por chamadas a rotinas do nível 2. Nas
máquinas com processadores vetoriais, por exemplo, um dos objetivos na implementação é man-
ter os tamanhos dos vetores utilizados o maior possível e na maioria dos algoritmos o resultado
é computado como um vetor por vez, o qual pode ser uma linha ou uma coluna. Nas máquinas
de registradores vetoriais a performance pode ser ainda maior fazendo-se reuso dos resultados
de um registrador vetorial ao invés de guardar o vetor de volta na memória. Entretanto, essa
abordagem não é, em geral, adequada para máquinas que fazem uso de hierarquias de memória
e máquinas com processadores paralelos [13].
Nessas arquiteturas é preferível particionar as matrizes em blocos e fazer as operações de
matriz com matriz nos blocos, pois, organizando a computação dessa forma, é possível ter um
reuso completo dos dados enquanto o bloco continua na memória cacheou memória local.
Essa abordagem evita a movimentação excessiva de dados da memória. Adicionalmente, nas
arquiteturas que proveêm processamento paralelo, o paralelismo pode ser explorado de duas
formas. A primeira é através de operações em paralelo nos diferentes blocos, a segunda é fazer
em paralelo as operações de escalares e vetores de um mesmo bloco [13]. Em [13] os autores
apresentam de maneira mais detalhada o funcionamento do nível 3 da BLAS e o porquê de cada
escolha durante sua criação.
A BLAS oferce, por padrão, duas formas de manipulação de erros: através de um manip-
ulador de erros, o BLAS_ERROR, e por códigos de retorno de erro. A BLAS para matrizes
densas utiliza o manipulador de erros, enquanto a BLAS para matrizes esparsas trabalha com os
códigos de retorno de erro. Cada função padrão da BLAS determina se e quando um mecanismo
manipulador deve erro ser chamado e a especificação da funçãodeve conter as condições, se
existirem, que disparam o manipulador de erros [1].
Por ser uma biblioteca eficiente, portável e amplamente utilizada, a BLAS costuma ser
empregada no desenvolvimento dos softwares de álgebra linear de alta qualidade. A BLAS é
disponibilizada gratuitamente e pode ser livremente incorporada em softwares comerciais. Além
disso, existem versões otimizadas para diversas máquinas,ou seja, versões otimizadas depen-
58
dentes de máquina, desenvolvidas pelas fabricantes e fornecedores de software independentes
(ISV - independent software vendors). Alguns exemplos dessas versões são [23]:
• ACML - Desenvolvido pela AMD, é uma versão multithread que implementa os níveis 1,
2 e 3 da BLAS otimizados para processadores opteron;
• Velocity Engine - Implementação da Apple para processadores G4 e G5;
• MLIB - Versão implementada pela HP que inclui todas as rotinas dos 3 níveis além de
outras facilidades para a resolução de sistemas densos e esparsos;
• MKL - Desenvolvida pela Intel, é uma implementação multithread que pode ser adquirida
como um produto em separado, como adicional ao Intel ClusterToolkit ou com as versões
profissionais dos compiladores da Intel.
É possível também utilizar a biblioteca ATLAS (Automatically Tuned Linear Algebra Soft-
ware) [2], a qual permite gerar automaticamente uma versão otimizada da BLAS para a arquite-
tura escolhida [24]. Exemplos de funcionamento e código-fonte com chamadas para a BLAS
podem ser obtidos em [13].
5.2 LINPACK
O LINPACK é um pacote que contém uma coleção de programas paramétodos diretos para
resolução de sistemas lineares na formaA~x = ~y onde a matrizA pode ser de vários tipos:
genérica (GE), banda (GB), positiva definida (PO), banda positiva definida (HI), triangular (TR),
tridiagonal positiva definida (PT) que ainda, com exceção das duas primeiras, podem estar com-
pactadas ou não [17]. As rotinas do LINPACK são escritas em Fortran e a abordagem utilizada
é resolver as matrizes através da abordagem decomposicional da álgebra linear, ou seja, dada
uma matrizA, decompõe-se a mesma em um produto de outras matrizes mais simples e bem
estruturadas, as quais podem ser facilmente manipuladas para resolver o problema original. Seu
desenvolvimento foi baseado no nível 1 da BLAS e a maioria dasoperações de ponto flutuante
dos algoritmos do LINPACK são realizadas através de chamadas as rotinas da BLAS [26].
59
Atualmente, o LINPACK é mais conhecido enquanto benchmark do que biblioteca. A versão
original do benchmark LINPACK é vista como um acidente, tendo sido originalmente concebido
para ajudar os usuários da biblioteca LINPACK no fornecimento de informações sobre os tem-
pos de execução necessários para resolver um sistema de equações lineares. A primeira aparição
do LINPACK como benchmark foi em 1979, com os anos ele foi recebendo incrementos e, nos
dias de hoje, é, na verdade, composto por três benchmarks, sendo o Highly Parallel Computing
Benchmark (HPL) o mais importante deles, o qual serve como medida para elaboração da lista
dos 500 computadores mais rápidos do mundo [6, 14]. O método usado no benchmark é decom-
posição LU com pivotamento parcial, a matriz é do tipo densa,composta por elementos inteiros
distribuídos aleatoriamente entre -1 e 1. A resolução do sistema de equações requerO(n3) op-
erações de ponto flutuante, mais especificamente,2/3n3 +2n2 +O(n) adições e multiplicações
de ponto flutuante [26].
Dentre as técnicas para melhora da performance na resoluçãodos sistemas lineares, duas
se destacam no LINPACK: desenrolamento dos laços e reuso dosdados. Observa-se que, fre-
qüentemente, o maior volume de computação de um programa está localizado em menos de 3%
do código fonte. Essa porcentagem do código, também chamadade código crítico, em geral
consiste em um ou alguns poucos laços de repetição imersos, ou seja, em um nível maior de
aninhamento. O desenrolamento do laço consiste em replicarseu conteúdo, fazendo os devidos
ajustes, e aumenta a performance porque há uma diminuição direta dosoverheadsinerentes ao
loop. Nas máquinas com instruções vetoriais, no entanto, essa técnica tem o efeito oposto [26].
Em relação ao reuso de dados, sabe-se que a cada passo do processo de fatoração do LIN-
PACK são feitas operações vetoriais para modificar uma submatriz inteira dos dados. Essa atu-
alização faz com que um bloco dos dados seja lido, atualizadoe escrito novamente na memória
central. O número de operações de ponto flutuante é2/3n3 e o número de referências a dados
é, em ambos os casos (leitura e escrita), de2/3n3. Assim, para cada par adição/multiplicação é
feita a leitura e escrita dos elementos, levando a um baixo reuso dos dados. Mesmo quando as
operações são vetoriais, existe um gargalo significante na movimentação dos dados, o que resulta
em uma performance ruim. Em máquinas vetoriais, isso se traduz em duas operações de vetor
e três referências vetoriais á memória. Nos computadores super-escalares isso resulta em uma
60
grande movimentação e atualização dos dados. Esse contextofaz com que o LINPACK tenha
desempenho reduzido em computadores de alta performance emque o custo do movimento dos
dados é semelhante ao das operações de ponto flutuante. Uma possível solução é reestruturar
os algoritmos de modo que explorem a hierarquia de memória das arquiteturas, o que pode ser
feito, por exemplo, armazenando os dados o maior tempo possível nas memórias de nível mais
próximo do processador, ou seja, aumentando o reuso dos dados [26]. Conforme visto na Seção
5.1, essas otimizações estão presentes nos níveis 2 e 3 da BLAS, adotos por outros pacotes,
como o LAPACK (Seção 5.3.
As soluções dos sistemas lineares no LINPACK benchmark são consideradas corretas quando
possuem a mesma precisão relativa que as apresentadas pelastécnicas padrão como a Eliminação
de Gauss no pacote LINPACK. Essa precisão é dada pela Equação5.1, ondeA ∈ ℜnxn;x, b ∈
ℜn, ǫ é a precisão de máquina para operações de aritmética de pontoflutuante de 64 bits e‖.‖ é
qualquer norma consistente.
‖A~x − ~y‖
‖A‖.‖x‖.n.ǫ= O(1) (5.1)
5.3 LAPACK – Linear Algebra PACKage
O Lapack é uma biblioteca escrita em Fortran 77 que provê subrotinas para resolução de
sistemas de equações lineares, problemas de auto-valores,dentre outros. As fatorações de ma-
trizes utilizadas são LU, Cholesky, QR, SVD Schur e Schur Generalizada. O LAPACK suporta
matrizes densas e do tipo banda, porém, não suporta matrizesesparsas. Os elementos das ma-
trizes podem ser números reais ou complexos e pode-se utilizar precisão simples ou dupla [24].
O LAPACK é o programa estado da arte para resolução de problemas de equações representadas
em matrizes densas e do tipo banda, além de outros tipos de operações da álgebra linear [1].
O objetivo do LAPACK, quando iniciado seu projeto, era fazeras já amplamente utilizadas
bibliotecas EISPACK e LINPACK rodarem eficientemente em computadores paralelos de memória
compartilhada e vetorias. Nessas máquinas, o LINPACK e EISPACK são ineficientes porque
seus padrões de acesso à memória negligenciam a hierarquia de memória dessas máquinas, o
61
que faz com que o custo do acesso aos dados seja bastante alto.No LAPACK, visando contornar
esse problema, os algoritmos foram reorganizados para utilizar as operações de matrizes em
bloco, tais como multiplicação de matrizes, nos laços mais internos. Essas operações em bloco
podem ser otimizadas para cada arquitetura de modo a tirar proveito da arquitetura de memória
e prover uma forma de se atingir alta eficiência nas diversas máquinas modernas [24, 13].
As rotinas do LAPACK são escritas de modo que a computação ocorra, tanto quanto pos-
sível, através de chamadas para rotinas da BLAS. Enquanto o LINPACK e EISPACK são basea-
dos nas operações vetoriais da BLAS (nível 1), o LAPACK foi desenvolvido para explorar o
nível 3, tendo, inclusive, influenciado seu desenvolvimento posteriormente [24]. Essa influência
se deu ao fato de que, somente após a implementação do LAPACK,é que algumas operações
da BLAS como rotinas paracopiar uma matriz(GE_COPY) ecalcular a norma de uma matriz
(GE_NORM), dentre outras, passaram a ser utilizadas com maior freqüência e foram implemen-
tadas como rotinas separadas [1]. Devido à granularidade grossa do nível 3 das operações da
BLAS, seu uso provê alta eficiência em muitos computadores dealta performance, principal-
mente naqueles em que implementações especiais do código são oferecidas pelo fabricante da
máquina [24].
5.4 scaLAPACK - Scalable LAPACK
O scaLAPACK é uma biblioteca de software para fazer computação de matrizes densas e do
tipo banda em computadores do tipo MIMD de memória distribuída, utilizando, para tal, troca
explícita de mensagens e em redes de estações de trabalho quesuportam PVM (Parallel Virtual
Machine) e/ou MPI (Message Passing Interface). É um pacote para ambientes de computação
heterogêneos implementado com paradigma de programação SPMD[25].
Assim como no LAPACK, as rotinas do scaLAPACK são baseadas emalgoritmos que uti-
lizam o paralelismo por decomposição em blocos de movo a diminuir a movimentação dos
dados entre os diferentes níveis da hierarquia de memória. OscaLAPACK mantém suas roti-
nas, quando possível, compatíveis com suas equivalentes noLAPACK. As bases fundamentais
do scaLAPACK são as versões para memória distribuída dos níveis 1, 2 e 3 da BLAS, ou seja,
62
a PBLAS, e um conjunto de subprogramas para comunicação de programas de álgebra linear
(BLACS - Basic Linear Algebra Communication). Nas rotinas do scaLAPACK, toda a comuni-
cação entre os processos ocorre com o uso da PBLAS e da BLACS[25].
5.5 HPL - High Performance Linpack
O HPL [28] é um benchmark para resolução de sistemas densos e aleatórios de equações
lineares. Prove programas que testam o tempo de execução e a precisão das soluções obtidas. O
HPL é na verdade uma evolução do LINPACK, é a implementação domesmo para máquinas de
memória distribuídas. O HPL é o benchmark utilizado na listado top 500 [14] e sua implemen-
tação é baseada nas bibliotecas MPI, para comunicação e na BLAS, para as rotinas de álgebra
linear [28].
63
6 Conclusão
Nesse trabalho foi apresentada a relação da computação paralela com a resolução dos sis-
temas de equações lineares de grande porte. Dentro desse contexto, foram revisados tópicos de
ambas as áreas, além dos métodos numéricos utilizados para resolução dos SELAs. A resolução
eficiente dos sistemas de equações lineares é importantíssima para diversas áreas do conheci-
mento, e, somente através da computação paralela, é que tem-se conseguido atingir os níveis
necessários para o correto entendimento dos problemas por eles modelados.
Observou-se no entanto que uma série de questões carecem de padronização na computação
paralela, como, por exemplo, a implementação de algoritmosportáveis e as medidas de desem-
penho. Não há como comparar arquiteturas diferentes de modogenérico, embora isso ocorra
através, por exemplo do TOP 500, pois, os diferentes benchmarks tendem a beneficiar diferentes
abordagens. Assim, a escolha da arquitetura e dos métodos numéricos utlizados para resolução
dos SELAs deve ser feita de acordo com o problema que se desejamodelar, ou seja, se o prob-
lema resulta em matrizes esparas ou densas, se possuem estruturas específicas na distribuição
dos elementos ou não, dentre outros fatores.
64
Bibliografia
[1] An updated set of basic linear algebra subprograms (blas). ACM Trans. Math. Softw.,
28(2):135–151, 2002.
[2] Automatically tuned linear algebra software (atlas), 2008. http://math-atlas.sourceforge.net
- Acesso em 15 jun. 2008.
[3] Beowulf.org overview, 2008. http://www.beowulf.org/overview - Acesso em 20 mai. 2008.
[4] The message passing interface (mpi) standard, 2008. http://www-unix.mcs.anl.gov/mpi/ -
Acesso em 22 jun. 2008.
[5] Pvm: Parallel virtual machine, 2008. http://www.csm.ornl.gov/pvm/ - Acesso em 22 jun.
2008.
[6] Top500 supercomputing sites, 2008. http://www.top500.org - Acesso em 28 mai. 2008.
[7] George Karypis Ananth Grama, Anshul Gupta.Introduction to Parallel Computing.
Addison-Wesley, Lebanon, Indiana, U.S.A, 2nd edition, 2003.
[8] Dimitri P. Bertsekas and John N. Tsitsiklis.Parallel and Distributed Computation: Nu-
merical Methods. Athena Scientific, 1997.
[9] H.A. Luther Brice Carnahan and James O. Wilkes.Applied Numerical Methods. John
Wiley and Sons, New York, 1969.
[10] Rajkumar Buyya.High Performance Cluster Computing, volume 2. Prentice Hall PTR,
Upper Saddle RIver, New Jersey, 1999.
65
[11] Patrícia Nunes da Silva.Equações Diferenciais Ordinárias, volume 1 ofCálculo Diferen-
cial e Integral. Rio de Janeiro, 1st edition, 2005.
[12] Francisco Chagas da Silva Filho. Modelagem de problemas de engenharia: solução de
equações diferenciais parciais pelo método dos elementos finitos. Revista Tecnologia,
1:134–144, 2005.
[13] J. J. Dongarra, Jeremy Du Croz, Sven Hammarling, and I. S. Duff. A set of level 3 basic
linear algebra subprograms.ACM Trans. Math. Softw., 16(1):1–17, 1990.
[14] Jack Dongarra. Frequently asked questions on the linpack benchmark and top500,
2008. http://www.netlib.org/utk/people/JackDongarra/faq-linpack.html - Acesso em 16
jun. 2008.
[15] Jack Dongarra, Thomas Sterling, Horst Simon, and ErichStrohmaier. High-performance
computing: Clusters, constellations, mpps, and future directions. Computing in Science
and Engg., 7(2):51–59, 2005.
[16] Iain S. Duff and Henk A. van der Vorst. Developments and trends in the parallel solution
of linear systems.Parallel Comput., 25(13–14):1931–1970, 1999.
[17] Dalcidio Moraes Cláudio e Jussara Maria Marins.Cálculo Numérico Computacional. Ed-
itora Atlas, São Paulo, 3rd edition, 2000.
[18] José Alberto Cuminato e Messias Meneguette Junior.Discretização de Equações Diferen-
ciais Parciais: Técnicas de Diferenças Finitas. Universidade de São Paulo e Universidade
Estadual Paulista, Agosto 2002.
[19] Cesar Augusto Fonticielha De Rose e Philippe O. A. Navaux. Arquiteturas paralelas.
Livros Didáticos. Sagra Luzzatto, Porto Alegre, 1st edition, 2003.
[20] Gerson Geraldo H. Cavalheiro e Rafael R. dos Santos.Atualizações em Informática 2007.
SBC, Porto Alegre, 1st edition, 2007. Capítulo de Livro. Capítulo 8: Multiprogramação
leve em arquiteturas multi-core.
66
[21] Michael J. Flynn. Some computer organizations and their effectiveness.IEEE Transactions
on Computers, C-21(9):948–960, September 1972.
[22] Michael J. Flynn and Kevin W. Rudd. Parallel architectures. ACM Computing Surveys,
28(1):67–70, 1996.
[23] National Science Foundation. Blas (basic linear algebra subprograms), 2008.
http://www.netlib.org/blas - Acesso em 18 jun. 2008.
[24] National Science Foundation. Lapack – linear algebra package, 2008.
http://www.netlib.org/lapack - Acesso em 15 jun. 2008.
[25] National Science Foundation. The scalapack project, 2008.
http://www.netlib.org/scalapack/index.html - Acesso em18 jun. 2008.
[26] Antoine Petitet Jack J. Dongarra, Piotr Luszczek. The linpack benchmark: past, present
and future. Concurrency and Computation: Practice and Experience, 15(9):803–820,
2003.
[27] Eric E. Johnson. Completing an mimd multiprocessor taxonomy. SIGARCH Comput.
Archit. News, 16(3):44–47, 1988.
[28] Innovative Computing Laboratory. Hpl - a portable implementation of the
high-performance linpack benchmark for distributed-memory computers, 2008.
http://www.netlib.org/benchmark/hpl/ - Acesso em 23 jun.2008.
[29] Raquel S. Lotti, André Wilson Machado, Ênio Tonani Mazzieiro, and Janes Landre
Júnior. Aplicabilidade científica do método dos elementos finitos. Revista Dental Press de
ortodontia e ortopedia facial, 11(2):35–43, 2006.
[30] Yousef Saad.Iterative Methods for Sparse Linear Systems. SIAM, 2nd edition, January
2000.
[31] David B. Skillicorn. A taxonomy for computer architectures. Computer, 21(11):46–57,
1988.
67
[32] David B. Skillicorn and Domenico Talia. Models and languages for parallel computation.
ACM Comput. Surv., 30(2):123–169, 1998.
[33] Gilbert Strang.Linear Algebra and Its Applications. Thomson Learning, Massachusetts
Institute of Technology, 3rd edition, 1988.
[34] Adenauer Corrêa Yamin. Um estudo das potencialidades elimites na exploração do par-
alelismo. Trabalho individual ii, Universidade Federal doRio Grande do Sul, Porto Alegre,
Agosto 1999. 80p.
68
