Paralela8
-
Upload
abraham-zamudio -
Category
Education
-
view
981 -
download
1
description
Transcript of Paralela8
1
Programación y Computación paralela (8)
Resolviendo Sist.Ec.Linealesprovenientes de Ec.Dif.Parc. (2)
Glen D. Rodríguez R.
Basado en material de J. Demmel
2
Revisión de este tema
• Se vió la Ec. Poisson
• Y los métodos numéricos para solucionarla
• Método de Jacobi
• Método de Sobre relajación SOR rojo-negro
• Gradiente conjugada
• FFT
• Multigrid
• Comparación de métodos Se
red
uce
n a
mu
ltip
l. d
e ve
cto
r co
n m
atri
z d
isp
ersa
Pre
req
uis
ito
sp
ara
ente
nd
er M
ult
igri
d
3
Ec. Poisson en 1D: ∂∂∂∂2u/∂∂∂∂x2 = f(x)
2 -1
-1 2 -1
-1 2 -1
-1 2 -1
-1 2
T =2-1 -1
Grafo y “estencil”
4
Ec. Poisson 2D
• Similar al caso 1D, pero la matriz T es ahora así:
• 3D es análogo
4 -1 -1
-1 4 -1 -1
-1 4 -1
-1 4 -1 -1
-1 -1 4 -1 -1
-1 -1 4 -1
-1 4 -1
-1 -1 4 -1
-1 -1 4
T =
4
-1
-1
-1
-1
Grafo y “estencil”
5
Algoritmos para Poisson 2D (3D) (N = n2 (n3) vars)
Algoritmo Serial PRAM Memoria #Procs
• Dense LU N3 N N2 N2
• Band LU N2 (N7/3) N N3/2 (N5/3) N (N4/3)
• Jacobi N2 (N5/3) N (N2/3) N N
• Explicit Inv. N2 log N N2 N2
• Conj.Gradients N3/2 (N4/3) N1/2(1/3) *log N N N
• Red/Black SOR N3/2 (N4/3) N1/2 (N1/3) N N
• Sparse LU N3/2 (N2) N1/2 N*log N (N4/3) N
• FFT N*log N log N N N
• Multigrid N log2 N N N
• MINIMO N log N N
Referencia: James Demmel, Applied Numerical Linear Algebra, SIAM, 1997.
6
Motivación del Multigrid
• Recordemos que Jacobi, SOR, CG, y cualquier otro
algoritmo basado en multiplicación de vector por matriz dispersa puede mover información sólo en una celda a
la vez (por iteración)• Mover la información a lo largo de una malla n x n toma por lo
menos n pasos
• Se puede mostrar que disminuir el error por un factor fijo
c<1 toma mínimo unos Ω(log n) pasos límite asintótico mínimo
• Convergencia a un error fijo < 1 toma Ω(log n ) pasos
• Por lo tanto, converger en O(1) pasos (número constante de pasos para cualquier n) implica mover
información a lo largo de la malla o grilla más rápido que una celda por paso
7
Motivación de Multigrid
8
Bases de Multigrid
• Algoritmo básico:• Reemplazar problema en una grilla o malla fina por una
aproximación en una malla más gruesa (o malla burda)
• Resolver aproximadamente el problema en la grilla gruesa, y usar la solución como una “starting guess” para el problema en malla fina, el cual es luego resuelto iterativamente
• Resolver el problema en malla gruesa recursvamente, o sea, usando una aproximación más burda todavía, etc.
• El éxito depende de que la solución en malla burda sea
una buena aproximación a la solución de malla fina
9
Esa misma Gran Idea se usa en otros lados
• Reemplazar el problema fino por una aproximación burda
en forma recursiva.
• Multilevel Graph Partitioning (METIS): • Reemplazar el grafo a ser particionado por un grafo más burdo,
obtenido por el “Maximal Independent Set”
• Dada la partición de la grilla burda, refinar usando Kernighan-Lin
• Barnes-Hut (y el método de Multipolo Rápido) para
computar las fuerzas gravitacionales de n partículas en
tiempo O(n log n):• Aproximar partículas encajas por masa total, centro de gravedad
• Suficientemente bueno para partículas lejanas; para las cercanas, dividir la caja en forma recursiva.
• Todos estos ejemplos dependen de que la aproximación
burda sea suficientemente buena (por lo menos si estamos lejos en la recursión)
10
Multigrid usa “Divide y Vencerás” en 2 formas
• 1ra forma:• Resolver el problema en una malla dada invocando al Multigrid en
una aproximación burda para conseguir una buena “initial guess” a ser refinada después.
• 2da forma:• Piense que el error es la suma de senos de diferente frecuencias.
• La misma idea que con el método de las FFT, pero no en forma explícita.
• Cada llamada al Multgrid es responsable de suprimir los coeficientes de los senos en el 50% de las frecuencias más bajas en el error (gráficos después)
11
Multigrid en 1D a grosso modo
• Considere una grilla 1D de tamaño 2m+1 por simplicidad• Sea P(i) el problema de resolver la ec. discreta de Poisson en una
grilla 2i+1 en 1D. La ec.linear sería T(i) * x(i) = b(i)
• P(m) , P(m-1) , … , P(1) es la secuencia de problemas del más fino al más burdo
12
Multigrid en 1D y 2D a grosso modo
• Considere la grilla 2m+1 en 1D (2m+1 x 2m+1 en 2D) por simplicidad• Sea P(i) el rpoblema de resolver Poisson discreto en una grilla 2i+1 en
1D (2i+1 by 2i+1 en 2D)
• El sistema linear sería T(i) * x(i) = b(i)
• P(m) , P(m-1) , … , P(1) es la secuencia de problemas del más fino al más burdo
13
Operadores Multigrid• Para el problema P(i) :
• b(i) es el vector del lado derecho de la ecuación lineal y
• x(i) es la solución estimada actualmente
• Todos los sgtes. operadores solo promedian valores en puntos aledaños en la grilla (así la info se mueve más rápido en mallas burdas)
• El operador restricción R<i> mapea P(i) a P(i-1)
• Redefine el problema en la grilla fina P(i) a la grilla burda P(i-1)
• Usa muestreo o promedios
• b(i-1)= R<i> (b(i))
• El operador interpolación In<i-1> mapea la solución aprox. x(i-1) a x(i)
• Interpola la solución en la malla burda P(i-1) a la fina P(i)
• x(i) = In<i-1>(x(i-1))
• El operador solución S(i) toma P(i) y mejora la solución x(i)
• Usa Jacobi “ponderado” o SOR en un solo nivel de la malla
• x mejorado (i) = S<i> (b(i), x(i))
Ambos en la grillade tamaño 2i-1
14
Algoritmo Multigrid de ciclo V
Function MGV ( b(i), x(i) )
… Resuelve T(i)*x(i) = b(i) dado un b(i) y un “initial guess” para x(i)
… retorna un x(i) mejorado
if (i = 1)
computar solución exacta de x(1) of P(1) solo 1 incógnita
return x(1)
else
x(i) = S<i> (b(i), x(i)) mejorar la solución
atenuando el error de alta frecuencia
r(i) = T(i)*x(i) - b(i) computar el residual
d(i) = In<i-1> ( MGV( R<i> (r(i)), 0 ) ) resolver T(i)*d(i) = r(i) recursivamente
x(i) = x(i) - d(i) corregir solución en grilla fina
x(i) = S<i> ( b(i), x(i) ) mejorar solución de nuevo
return x(i)
15
¿Por qué es llamado ciclo V?
• Miren el dibujo del orden de las llamadas
• Después de un tiempo, el ciclo V luce así:
16
Complejidad de un Ciclo V
• En una computadora serial• En cada punto del Ciclo V, el cómputo es del orden O(número
de incógnitas)
• Costo del nivel i-esimo es (2i-1)2 = O(4 i)
• Si el nivel más fino es el m, el tiempo total es:
• Σ O(4 i) = O( 4 m) = O(# incógnitas)
• En un sistema paralelo (PRAM)• Con un procesador por cada punto de la grilla y comunicación
con costo despreciable, cada paso del ciclo V tarda un tiempo constante
• Tiempo Total del ciclo V es O(m) = O(log #incógnitas)
m
i=1
17
Full Multigrid (FMG)
• Intuición:
• Mejorar la solución haciendo varios ciclos V
• Evitar los costosos ciclos en grilla fina (alta frecuencia)
• El por qué funciona es para otro curso
Function FMG (b(m), x(m))
… retorna un x(m) mejorado dado un “initial guess”
computa la solución exacta x(1) del P(1)
for i=2 to m
x(i) = MGV ( b(i), In<i-1> (x(i-1) ) )
• En otras palabras:
• Resolver el problema con 1 incógnita
• Dada la solución de un problema burdo P(i-1) , mapear los valores de
la solución como el guess inicial del problema P(i)
• Resolver el problema más fino usando el ciclo V
18
Análisis del costo del Full Multigrid
• Una V para cada llamada al Full MultiGrid• Algunos usan Ws y otras formas
• Tiempo serial: Σ O(4 i) = O( 4 m) = O(# incógn.)
• Tiempo PRAM: Σ O(i) = O( m 2) = O( log2 # incógn.)
m
i=1
m
i=1
19
Complejidad de resolver la Ec.Poisson
• Teorema: error después de una llamada al Full
Multigrid<= c* error antes, donde c < 1/2, independentemente del # incógnitas
• Corolario: se puede hacer que el error se vuelva < que cualquier tolerancia fija arbitraria en n número fijo de
pasos, independentemente del tamaño de la grilla fina
• Esta es la más importante propiedad de convergencia del MultiGrid, que lo diferencia de otros métodos, que
convergen más lento en grillas grandes (finas)
• Complejidad total es proporcional al costo de una
llamada al FMG
20
El operador Solución S<i> - Detalles
• El operador Solución, S<i>, es un Jacobi ponderado
• Considere este problema 1D
• En el nivel i, el Jacobi puro hace:
x(j) := 1/2 (x(j-1) + x(j+1) + b(j) )
• En su lugar, Jacobi ponderado hace:
x(j) := 1/3 (x(j-1) + x(j) + x(j+1) + b(j) )
• En 2D, de forma similar se promedia la variable con sus vecinos.
21
Jacobi ponderado para amortiguar error de alta frecuencia
Error inicial“Áspero”Componentes de alta frec. son fuertesNorma = 1.65
Error después de 1 Jacobi“Suave”Comp. de alta frec. menos fuertesNorma = 1.06
Error después de 2 Jacobi“Más Suave”Débil comp. de alta frec.Norma = .92,
No disminuirá mucho más
22
Multigrid como algoritmo Divide y Vencerás
• Cada nivel en un ciclo V reduce el error en una parte del dominio de la frecuencia
23
El Operador Restricción R<i> - Detalles
• El operador restricción, R<i>, toma
• un problema P(i) con vector del lado derecho b(i) y
• lo mapea en un problema más burdo P(i-1) con vector b(i-1)
• En 1D, se promedia los valores de los vecinos
• xburdo(i) = 1/4 * xfino(i-1) + 1/2 * xfino(i) + 1/4 * xfino(i+1)
• En 2D, se promedia con los 8 vecinos (N,S,E,W,NE,NW,SE,SW)
24
Operador Interpolación In<i-1>: detalles
• El operador Interpolación In<i-1>, toma una función en una malla burda P(i-1) , y produce una función en una malla fina P(i)
• En 1D, se hace interpolación lineal a los vecinos burdos más próximos
• xfino(i) = xburdo(i) si el punto i de la grilla fina está también en la burda, si no:
• xfino(i) = 1/2 * xburdo(izquierda de i) + 1/2 * xburdo(derecha de i)
• En 2D, la interpolación debe promediar los 4 vecinos (NW,SW,NE,SE)
25
Convergencia del Multigrid en 1D
26
Convergencia del Multigrid en 2D
27
Multigrid paralela en 2D
• Multigrid en 2D necesita
computar con los vecinos más próximos
(hasta 8) en cada nivel
de la grilla
• Empieza con una grilla de n=2m+1 x 2m+1 (aquí
m=5)
• Se usa una malla de
procesadores de s x s (aquí s=4)
28
Modelo de performance del Multigrid 2D paralelo (ciclo V)
• Asumir grilla de 2m+1 x 2m+1 puntos, n= 2m-1, N=n2
• Asumir p = 4k procesadores, en una malla de 2k x 2k
• Procesadores comienzan con una sub grilla de 2m-kx 2m-k variables
• Considerar que el Ciclo V empieza en el nivel m
• En los niveles del m al k del ciclo V, cada procesador trabaja.
• En los noveles k-1 al 1, algunos procesadores están ociosos, por que una
malla de 2k-1 x 2k-1 variables (o menos) no puede ocupar cada procesador
• Costo de un nivel en el ciclo V
• If nivel j >= k, then costo =
O(4j-k ) …. Flops, proporcional al número de puntos/procesador
+ O( 1 ) α …. Envía un número constante de mensajes a los procs. vecinos
+ O( 2j-k) β …. Número de words (doubles o floats) enviados
• If nivel j < k, then costo =
O(1) …. Flops, proporcional al número de puntos/procesador
+ O(1) α …. Envía un número constante de mensajes a los procs. vecinos
+ O(1) β .… Número de words enviados
• Sume para todos los niveles 1,2,… y todos los ciclos V en FMG!!!
29
Comparición de Métodos (en O(.))
# Flops # Mensajes # Words enviadas
MG N/p + (log N)2 (N/p)1/2 +
log p * log N log p * log N
FFT N log N / p p1/2 N/p
SOR N3/2 /p N1/2 N/p
• SOR es más lento en todo
• Flops para el MG y FFT dependen de la precisión del MG
• MG comunica menos data (ancho de banda)
• El total de mensajes (latencia) depende …• Este análisis simplista no nos dice si MG o FFT es mejor cuando
α >> β
30
Consideraciones prácticas
• En la práctica, no tenemos que iterar hasta P(1)
• En programa secuencial, las grillas más burdas son
rápidas, pero en un prog.paralelo no.• Considere 1000 puntos por procesador
• En 2D, las vars. a comunicar es 4xsqrt(1000) ~= 128, o 13%
• En 3D, sería 1000-83 ~= 500, o 50%
• Ver Tuminaro y Womble, SIAM J. Sci. Comp., v14, n5, 1993 para el análisis de MG en nCUBE2 de 1024 procs.
• En grilla 64x64 vars., solo 4 por procesador
• eficiencia de 1 ciclo V fue 0.02, y en FMG 0.008
• En grilla 1024x1024
• eficiencias fueron 0.7 (MG ciclo V) y 0.42 (FMG)
• Aunque su eficiencia paralela es peor, FMG es 2.6 veces más
rápido que sólo los ciclos V
• nCUBE tenía comunicación rápida, procesadores lentos
31
Multigrid en un Mesh adaptativo
• Para problemas con un
rango dinámico muy grande, se necesita otro
nivel de refinamiento
• Se construye un mesh
adaptativo y se resuelve el multigrid (típicamente)
en cada nivel
• No se puede costear el uso de mesh fino en todo el
problema
32
Aplicaciones multibloque
• Resolver el sistema de ecuaciones en una unión de
rectángulos• subproblema de AMR
• Ej:
33
Refinando el Mesh Adaptativo (AMR)
• Usar estructuras de datos en el AMR
• El paralelismo usual es asignar grillas en cada nivel a
los procesadores
• Balance de carga es problemático
34
Soporte para AMR
• Dominios en ciertos sistemas (Titanium, etc.) diseñados
para este problema
• Librerías: Kelp, Boxlib, y AMR++
• Primitivas para ayudar con los updates en las fronteras
entre procesadores, etc.
35
Multigrid en un Mesh no estructurado
• Otro enfoque para
manejar trabajos variables es usar un
mesh no estructurado
que se refina más en las áreas de interes
• Rectangular adaptivo o
no estructurado?• Fórmulas más simples en
el rectangular
• Supuestamente más fácil de implementar (arrays) pero las fronteras tienden a dominar el código.
Hasta 39M variables en 960 procesadores,Con 50% eficiencia (Fuente: M. Adams)
36
Multigrid en una Mesh no estructurada
• Se necesita particionar el gráfico para paralelizar
• Recordar: Cuál es el propósito del Multigrid?
• Debe poderse crear grillas burdas (problema “duro”)• No se puede “saltearse los puntos uno si otro no” como antes
• Usar de nuevo “maximal independent sets”
• Como hacer que un grafo burdo aproxime bien al fino???
• Se debe redefinir R<> y In<>• Cómo convertir de grilla burda a fina y viceversa?
• Definir S<>• Cómo se computa la matriz? (fórmulas varían)
• Procesar meshes burdas eficientemente• Sería mejor limitarse a usar pocos procesadores en meshes
burdas?
• Debería cambiarse el “solver” en meshes burdas?
37
Fuente de Mesh no estructurada (elem.finito): Vertebra
Fuente: M. Adams, H. Bayraktar, T. Keaveny, P. Papadopoulos, A. Gupta
38
Multigrid para análisis elástico no lineal de hueso
Test mecánico depropiedades del material
Tomografía por computadora @ resolución 22 µµµµm
Imagen 3dµµµµFE mesh2.5 mm3
elementos 44 µµµµm
Hasta537M variables4088 Procesadores (ASCI White)70% eficiencia en paralelo
Fuente: M. Adams et al