1 Sistemas de Ecuaciones LinealesGlen D. Rodríguez R. Algoritmos Paralelos Basado en material de J. Demmel
2 Eliminación Gaussiana y matrices densasAlgebra lineal densa “Gaussian Elimination” (GE) para resolver Ax=b Optimizar GE para máquinas secuenciales con caché Usando multiplicación matriz x matriz (BLAS) LAPACK Partición de la Data en comp.paralelas GE paralelo ScaLAPACK Problemas de Eigenvalue
3 Sca/LAPACK Overview
4 Exitos usando Sca/LAPACKMuy usado: Adoptado por Mathworks, Cray, Fujitsu, HP, IBM, IMSL, NAG, NEC, SGI, … >56M visitas en la web de Netlib (incl. CLAPACK, LAPACK95) Descubrimientos hechos a través de la solución de sistemas de ec. densos Artículo en Nature sobre el universo plano usó ScaLAPACK Otros en Physics Review B también lo usaron 1998 Premio Gordon Bell These are two cover stories that direct results of numerical tools supported by NPACI. SuperLU is a sparse linear system solver, and ScaLAPACK is a dense and band solver for linear systems, least squares problems, eigenvalue problems, and singular value problems. There are other more recent articles in places has Phys Rev B that use ScaLAPACK as well. Cosmic Microwave Background Analysis, BOOMERanG collaboration, MADCAP code (Apr. 27, 2000). ScaLAPACK
5 Usos de Solución de sistemas de ecuaciones linealesResistencia de estructuras (ing. civil) Stress, tensiones en piezas mecánicas Electromagnetismo, electrodinámica Circuitos eléctricos Ing. de tráfico. DSP, óptica (de-convolución) Economía y finanzas.
6 Motivación (1) 3 problemas básicos de Alg. LinealEcuaciones Lineales: Resolver x en Ax=b Mínimos cuadrados: Hallar x que minimiza ||r||2 S ri2 donde r=Ax-b Estadística: Ajustar datos a funciones simples 3a. Eigenvalues: Hallar l , x tal que Ax = l x Análisis de vibraciones (terremotos, circuitos) 3b. Singular Value Decomposition: ATAx=2x Ajuste de data, Obtención de información Muchas variaciones para diferentes estructuras de A Simétrica, definida positiva, en bandas, …
7 Motivación (2) Por qué A densa, en vez de A dispersa?Muchas matrices grandes son dispersas, pero … Alg. Densos son más fáciles de entender Algunas aplicaciones producen matrices densas grandes LINPACK Benchmark (www.top500.org) “Qué tan rápida es tu computadora?” = “Qué tan rápido resuelves Ax=b denso?” Alg. de matriz dispersa grande con frecuencia producen menores (pero aún grandes) problemas densos.
8 Records en la solución de Sist. Densos (2006) Gigaflops Máquina n=100 n= Any n Pico IBM BlueGene/L K K (131K procs) (281 Teraflops) (n=1.8M) NEC SX 8 (8 proc, 2 GHz) (1 proc, 2 GHz) … Efficiency of NEC SX8 at n=1000 = (75140/14960)/8 = .63 Palm Pilot III (1.69 Kiloflops)
9 Records en la solución de Sist. Densos (2011/2) Performance (Gflops) , Problem size on 4 nodes. 4 AMD Athlon K7 500 Mhz (256 Mb) – (2x) 100 Mbs Switched – 2 NICs per node (channel bonding) 8 Duals Intel PIII 550 Mhz (512 Mb) - Myrinet Compaq 64 nodes (4 ev Mhz processors / node) AlphaServer SC GRID 2000 5000 8000 10000 1 x 4 1.28 1.73 1.89 1.95 2 x 2 1.17 1.68 1.88 1.93 4 x 1 0.81 1.43 1.70 1.80 Performance (Gflops), Problem size on 8- and 16-processors grids. GRID 2000 5000 8000 10000 15000 20000 2 x 4 1.76 2.32 2.51 2.58 2.72 2.73 4 x 4 2.27 3.94 4.46 4.68 5.00 5.16 Efficiency of NEC SX8 at n=1000 = (75140/14960)/8 = .63 GRID N 1/2 Nmax Rmax (Gflops) Efficiency 1 x 1 150 6625 1.136 1.000 2 x 2 800 13250 4.360 0.960 4 x 4 2300 26500 17.00 0.935 8 x 8 5700 53000 67.50 0.928 16 x 16 14000 106000 263.6 0.906
10 Eliminación Gaussiana (GE) en Ax=bAñadir múltiplos de cada fila a las demás para volver A en triangular sup. Resolver el sist.triang. resultante Ux = c por sustitución … para cada columna i … “cerear” debajo de la diagonal añadiendo múltiplos de la fila i a las de abajo for i = 1 to n-1 … para cada fila j debajo de la fila i for j = i+1 to n … añadir un multiplo de la fila i a la fila j tmp = A(j,i); for k = i to n A(j,k) = A(j,k) - (tmp/A(i,i)) * A(i,k) . . . … . . . . . . After i=1 After i=2 After i=3 After i=n-1
11 Refinando el algoritmo GE (1)Versión inicial Elimina computación de constante tmp/A(i,i) de loop interno … para cada columna i … “cerear” debajo de la diagonal añadiendo múltiplos de fila i a las de abajo for i = 1 to n-1 … para cada fila j debajo de fila i for j = i+1 to n … añadir un múltiplo de la fila i a la fila j tmp = A(j,i); for k = i to n A(j,k) = A(j,k) - (tmp/A(i,i)) * A(i,k) for i = 1 to n-1 for j = i+1 to n m = A(j,i)/A(i,i) for k = i to n A(j,k) = A(j,k) - m * A(i,k) m
12 Refinando el algoritmo GE (2)Última versión No computa lo que ya se conoce: ceros bajo la diagonal en la columna i for i = 1 to n-1 for j = i+1 to n m = A(j,i)/A(i,i) for k = i to n A(j,k) = A(j,k) - m * A(i,k) for i = 1 to n-1 for j = i+1 to n m = A(j,i)/A(i,i) for k = i+1 to n A(j,k) = A(j,k) - m * A(i,k) m No calcula los ceros
13 Refinando el algoritmo GE (3)Última versión Guarda los factores m bajo la diagonal en el espacio donde irían los ceros para uso posterior for i = 1 to n-1 for j = i+1 to n m = A(j,i)/A(i,i) for k = i+1 to n A(j,k) = A(j,k) - m * A(i,k) for i = 1 to n-1 for j = i+1 to n A(j,i) = A(j,i)/A(i,i) for k = i+1 to n A(j,k) = A(j,k) - A(j,i) * A(i,k) m Guarda m aquí
14 Refinando el algoritmo GE (4)Última versión for i = 1 to n-1 for j = i+1 to n A(j,i) = A(j,i)/A(i,i) for k = i+1 to n A(j,k) = A(j,k) - A(j,i) * A(i,k) Partir el loop for i = 1 to n-1 for j = i+1 to n A(j,i) = A(j,i)/A(i,i) for k = i+1 to n A(j,k) = A(j,k) - A(j,i) * A(i,k) Guardar todos los m aquí antes de modificar el resto de la matriz
15 Refinando el algoritmo GE (5)for i = 1 to n-1 for j = i+1 to n A(j,i) = A(j,i)/A(i,i) for k = i+1 to n A(j,k) = A(j,k) - A(j,i) * A(i,k) Última versión Rápido usando operaciones de matrices (BLAS) for i = 1 to n-1 A(i+1:n,i) = A(i+1:n,i) * ( 1 / A(i,i) ) A(i+1:n,i+1:n) = A(i+1:n , i+1:n ) - A(i+1:n , i) * A(i , i+1:n)
16 Qué computa en realidad la GE?Llamemos a la matriz triangular inferior de multiplicadores m como M, y hagamos L = I+M Llamemos al triangulo superior de la matriz final como U Lema (Factorización LU): Si el algoritmo anterior termina (no hay división por cero) entonces A = L*U Solución de A*x=b usando GE Factorizar A = L*U usando GE (costo = 2/3 n3 flops) Resolver L*y = b para y, usando substitución (costo = n2 flops) Resolver U*x = y para x, usando substitución (costo = n2 flops) Por lo tanto A*x = (L*U)*x = L*(U*x) = L*y = b como se desea for i = 1 to n-1 A(i+1:n,i) = A(i+1:n,i) / A(i,i) A(i+1:n,i+1:n) = A(i+1:n , i+1:n ) - A(i+1:n , i) * A(i , i+1:n)
17 Problemas con el alg. GE básicoQué pasa si algún A(i,i) es cero? O es muy pequeño? Resultado inválido (div.por cero), o “inestable”, se necesita pivotear Se computan operaciones BLAS 1 o BLAS 2, pero sabemos que las BLAS 3 (mult. de matrices) es más rápida (clase anterior…) for i = 1 to n-1 A(i+1:n,i) = A(i+1:n,i) / A(i,i) … BLAS 1 (escala un vector) A(i+1:n,i+1:n) = A(i+1:n , i+1:n ) … BLAS 2 (update en rango 1) - A(i+1:n , i) * A(i , i+1:n) Pico BLAS 3 BLAS 2 BLAS 1
18 Pivotear en GE A = [ 0 1 ] falla por que no se puede dividir entre A(1,1)=0 [ 1 0 ] Pero resolver Ax=b es facilísimo! Cuando elemento diagonal A(i,i) es diminuto (no sólo cero), el algoritmo puede terminar pero dar una respuesta completamente errada Inestabilidad numérica La causa es el error de redondeo a punto flotante Cura: usar Pivot (intercambiar filas de A) tal que A(i,i) sea grande
19 GE con Pivot Parcial (GEPP)Pivoteo Parcial: intercambiar filas tal que A(i,i) es la mayor en su columna for i = 1 to n-1 Hallar y guardar k donde |A(k,i)| = max{i <= j <= n} |A(j,i)| … o sea, el mayor element en el resto de la columna i if |A(k,i)| = 0 exit con advertencia que A es singular, o casi singular elseif k != i intercambiar filas i , k de A end if A(i+1:n,i) = A(i+1:n,i) / A(i,i) … cada cociente cae en [-1,1] A(i+1:n,i+1:n) = A(i+1:n , i+1:n ) - A(i+1:n , i) * A(i , i+1:n) Lema: este algoritmo computa A = P*L*U, donde P es una matriz de permutación. Es numéricamente estable en la práctica Para más detalles ver código LAPACK en
20 Problemas con el alg. GE básicoQué pasa si A(i,i) es cero? O es muy pequeño? Resultado inválido, o “inestable”, así que se debe pivotear Se computan operaciones BLAS 1 o BLAS 2, pero sabemos que las BLAS 3 (mult. de matrices) es más rápida (clase anterior…) for i = 1 to n-1 A(i+1:n,i) = A(i+1:n,i) / A(i,i) … BLAS 1 (escala un vector) A(i+1:n,i+1:n) = A(i+1:n , i+1:n ) … BLAS 2 (update en rango 1) - A(i+1:n , i) * A(i , i+1:n) Pico BLAS 3 BLAS 2 BLAS 1
21 Convirtiendo BLAS2 a BLAS3 en GEPPBloques Usado para optimizar mult.matrices Más difícil en GE por las dependencias de data en GEPP GRAN IDEA: Updates retrasados Salvar varios updates consecutivos BLAS2 a una “matriz relegada” Aplicar muchos updates simultáneamente en una operación BLAS2 La misma idea funciona en otras partes del alg. lineal Quedan preguntas por resolver…. Primer enfoque: necesito escoger un tamaño de bloque b Algoritmo guarda y aplica b updates b debe ser suficientemente pequeño de tal forma que la sub matriz activa consistente en b columnas de A quepa en el cache b debe ser suficientemente grande para que BLAS3 sea rápido
22 GEPP en bloques (www.netlib.org/lapack/single/sgetrf.f)for ib = 1 to n-1 step b … Procesa b columnas end = ib + b … Apunta al final del bloque de b columns aplica versión BLAS2 de GEPP para obtener A(ib:n , ib:end) = P’ * L’ * U’ … LL denota la parte triangular inferior de A(ib:end , ib:end) + I A(ib:end , end+1:n) = LL-1 * A(ib:end , end+1:n) … update siguientes b filas de U A(end+1:n , end+1:n ) = A(end+1:n , end+1:n ) - A(end+1:n , ib:end) * A(ib:end , end+1:n) … aplica updates retrasados con una mult.matrices … de dimensiones w x b , b x w Need to apply row interchanges from A(ib:n, ib:end) to trailing submatrix
23 Eficiencia del GEPP en bloques
24 Vista general de LAPACK y ScaLAPACKLibrería estándar para algebra lineal densa y en bandas Sistemas lineales: A*x=b Problemas de mínimos cuadrados: minx || A*x-b ||2 Problemas de Eigenvalues: Ax = lx, Ax = lBx Descomposición en valores singulares (SVD): A = USVT Alg. reorganizados para usar lo más posible BLAS3 Base de librerías matemáticas en muchos sistemas, Matlab, etc … Aún faltan más mejoras por implementar Proyectos para maestrías, doctorados?
25 Performance de LAPACK (n=1000)Performance de Eigen-values, SVD, etc.
26 Performance de LAPACK (n=100)Eficiencia es mucho menor para matrices más chicas
27 Paralelizando la GE Pasos de la paralelización DescomposiciónDescomposición: identificar suficiente trabajo paralelo, pero no demasiado Asignación: balance de carga entre procesos Coordinación: comunicación y sincronización Mapeo: qué procesadores ejecutan qué tareas? Descomposición En la parte BLAS 2, casi cada flop en el loop interno puede hacerse en paralelo, así que con n2 procs., necesito 3n pasos paralelos Esto es grano muy fino, mejor es llamar a mult.matrices locales Necesario usar mult.matrices paralelo Asignación Que procesadores son responsables de qué submatrices? for i = 1 to n-1 A(i+1:n,i) = A(i+1:n,i) / A(i,i) … BLAS 1 (escala un vector) A(i+1:n,i+1:n) = A(i+1:n , i+1:n ) … BLAS 2 (update rango 1) - A(i+1:n , i) * A(i , i+1:n)
28 Diferentes particiones de data para GE paralelo1 2 3 1 2 3 Mal balance de carga: P0 ocioso luego de n/4 pasos Mejor balance, pero difícil usar BLAS2 o BLAS3 1) Bloque de columnas 1D 2) Columnas 1D cíclicas Se puede negociar balance de carga y performance BLAS2/3 cambiando b, pero la factorización del bloque de columnas es el cuello de botella 1 2 3 1 2 3 Difícil de direccionar b 4) Bloques de diagonales 3) Bloques cíclicos de columnas 1D Another problem with Skewed Layout: BLAS3 updates all nonlocal! Bad load balance in smaller trailing submatrices 1 2 3 1 2 3 Mal balance de carga: P0 ocioso luego de n/2 pasos El mejor! 6) Bloques cíclicos de filas y columnas 2D 5) Bloques de filas y columnas 2D
29 Revisión: GEPP con BLAS 3 (con bloques)for ib = 1 to n-1 step b … Procesa matriz de b columnas end = ib + b … Apunta alfinal del bloque de b columnas aplica ver.BLAS2 de GEPP para onbtener A(ib:n , ib:end) = P’ * L’ * U’ … LL denota parte triangular inferior de A(ib:end , ib:end) + I A(ib:end , end+1:n) = LL-1 * A(ib:end , end+1:n) … update sgtes.b filas de U A(end+1:n , end+1:n ) = A(end+1:n , end+1:n ) - A(end+1:n , ib:end) * A(ib:end , end+1:n) … aplicar updates retrasados con una mult.matrices … de dimensiones w x b, x w BLAS 3
30 Bloque cíclico de Filas y ColumnasProcesadores y bloques de matrices están distribuidos en un array 2d Array de procesadores prow x pcol Bloques de matrices brow x bcol Paralelismo “pcol” veces en cualquier columna, y llamadas a BLAS2 y BLAS3 en matrices de tamaño brow x bcol Cuello de botella serial es menos problemático prow pcol y brow bcol posible, es más deseable bcol brow 1 2 3
31 GE distribuido con Bloques cíclicos 2DTamaño de bloque b en el algoritmo y los tamaños de bloque brow y bcol en la partición satisfacen b=bcol. Regiones sombreadas indican procesadores ocupados con computación o comunicación. No es necesario tener barrera entre pasos sucesivos del algoritmo, ej: pasos 9, 10, y 11 se pueden hacer en “paralelo”
32 Distributed GE with a 2D Block Cyclic LayoutContador de columnas Contador de filas
33 verde = verde -azul * rosa–1 verde = verde -azul * rosa Mult. matrices de
34 Revisión de MatMul paralelaSe quiere problemas grandes por procesador PDGEMM = PBLAS mat.mul. Observaciones: Para N fijo, si P crece, Mflops crece pero menos que 100% de eficiencia Si P es fijo, si N crece, Mflops (eficiencia) crece DGEMM = rutina BLAS para mat.mul. Máxima velocidad para PDGEMM = # Procs * velocidad DGEMM Observaciones: Eficiencia siempre es por lo menos 48% Para N fijo, si P crece, eficiencia disminuye Para P fijo, si N crece, eficiencia crece
35 PDGESV = LU paralelo de ScaLAPACKComo no puede correr más rápido que el loop interno (PDGEMM), hacemos: Eficiencia = Veloc(PDGESV)/Veloc(PDGEMM) Observaciones: Eficiencia muy por encima de 50% para problemas suficientemente grandes Para N fijo, si P crece, eficiencia disminuye (igual que en PDGEMM) Para P fijo, si N crece, eficiencia crece (igual que en PDGEMM) Según la tabla inferior, el coso de resolver Ax=b es aprox. la mitad del mat.mul. para matrices suficientemente grandes. Del contador de flops se deduce que debería ser (2*n3)/(2/3*n3) = 3 veces más rápida, pero la comunicación lo hace algo más lento.
36 Modelos de performance ScaLAPACK (1)Contador de operaciones de ScaLAPACK tf = 1 tm = a tv = b NB = brow=bcol P = prow = pcol Assume brow = bcol = NB, prow = pcol = sqrt(P) t_f = time per flop , t_m = time per message = latency = alpha*t_f, t_v = time per word = 1/BW = beta*t_f Note for GE: bandwidth term same as SUMMA (matmul), latency NB times larger
37 Modelos de performance ScaLAPACK (2)Comparar Predicciones y Mediciones (LU) (Cholesky)
38 Escalabilidad de LAPACK y ScaLAPACKProblemas “unilaterales” son escalables En GE, A factorizada como el producto de 2 matrices A = LU por medio de pre-multiplicar A por una secuencia de matrices más simples Asintóticamente es 100% BLAS3 LU (“Linpack Benchmark”) Cholesky, QR Problemas “Bilaterales” som menos escalables A factorizada como producto de 3 matrices por medio de pre y post multiplicación 50% BLAS2, no llega a 100% BLAS3 Eigenproblemas, SVD descomposición de valor singular Eigenproblemas asimétricos son peores Problemas de banda estrecha son los peores (para llevarlo a BLAS3 o paralelizar) Ecuaciones lineales y eigenproblemas
39 Algoritmos RecursivosTambién usa updates retrasados, pero organizados en forma diferente Puede explotar particiones de data recursiva 3x speedups en mínimos cuadrados para matrices altas y delgadas En teoría, la performance de la jerarquía de memoria es óptima Referencias “Recursive Block Algorithms and Hybrid Data Structures,” Elmroth, Gustavson, Jonsson, Kagstrom, SIAM Review, 2004
40 GE usando un Algoritmo RecursivoF. Gustavson y S. Toledo Algoritmo LU : 1: Partir la matriz en 2 rectangulos (m x n/2) si solo queda 1 columna, escale por el reciproco del pivot y “return” 2: Aplicar el algoritmo LU a la parte izquierda 3: Aplique transformaciones a la parte derecha (solución triangular A12 = L-1A12 y mult. matrices A22=A22 -A21*A12 ) 4: Aplicar algoritmo LU a la parte derecha L A12 A21 A22 La mayor parte del computo es multiplicar matrices de tamaño n/2, n/4, n/8, … Source: Jack Dongarra
41 Factorizaciones recursivasTan precisas como los métodos convencionales Mismo número de operaciones Se forman bloques de tamaño variable en forma automática Sólo operaciones BLAS nivel 1 y 3! Simplicidad de expresión Potencialmente eficientes pero a la vez “indiferente” al caché Pero no se debería llevar la recursión hasta el extremo de 1 columna! La formulación recursiva es sólo un reajuste de algoritmos más viejos del LINPACK
42 Recursive LU Dual-processor LAPACK Recursive LU LAPACK UniprocessorSource: Jack Dongarra
43 Limites de los algortimos recursivosDos clases de descomposiciones de matrices densas “Unilaterales” Secuencia de operaciones sencillas aplicadas en la izquierda de la matriz GE: A = L*U o A = P*L*U GE simétrica: A = L*D*LT Cholesky, sólo matrices simétricas: A = L*LT Descomposición QR para mínimos cuadrados: A = Q*R Puede acercarse a 100% BLAS 3 Susceptible de recursión “Bilaterales” Secuencia de operaciones sencillas aplicadas en ambos lados, alternadamente Algoritmos de Eigenvalores, SVD Por lo menos ~25% BLAS 2 Enfoque recursivo imposible, difícil? Algunos progresos la última decada: SVD (25% vs 50% BLAS2)
44 Algoritmos “Out of Core”Out-of-core significa que la matriz reside en disco; muy grande para entrar completa en la RAM Mucho más difícil evitar efecto de latencia del disco QR más fácil que LU por que no hace falta pivotear para QR Source: Jack Dongarra
45 Resolviendo Sist.Ec.Lineales provenientes de Ec.Dif.Parc.
46 Repaso de “Fuentes de paralelismo”Sistemas de eventos discretos: Ejs: “Game of Life,” simulación de circuitos lógicos. Sistemas de partículas: Ejs: bolas de billar, simulación de dispositivos semiconductores, galaxias. Variables agregadas dependientes de parámetros continuos: EDOs, ejs: simulación de circuitos (Spice), mecánica estructral, cinética química. Variables contínuas dependientes de parámetros continuos: EDPs, ejs: calor, elasticidad, electrostática. Un fenómeno dado puede ser modelado en múltiples niveles. Muchas simulaciones combinan más de una técnica.
47 Repaso: Resolviendo EDPs° Problemas hiperbólicos (ondas): Onda sonora (posición, time) Usan iteración en el tiempo en forma explícita Solución en cada punto depende de los vecinos en Δt anterior Problemas elípticos (estado estable): Potencial electrostático (posición) Todo valor depende de todos lo demás valores Localidad es más difícil de hallar que en los problemas hiperbólicos Problemas parabólicos (dependientes del tiempo): Temperatura (posición, tiempo) Involucra una solución elíptica cada Δt Nos enfocaremos en los elípticos Ejemplo canónico es la ec. Poisson 2u/x2 + 2u/y2 + 2u/z2 = f(x,y,z)
48 Ec. Poisson aparece en muchos lados3D: 2u/x2 + 2u/y2 + 2u/z2 = f(x,y,z) f representa las fuentes; además se necesita condiciones de frontera 2D: 2u/x2 + 2u/y2 = f(x,y) 1D: d2u/dx2 = f(x) Potencial electrostático o gravitacional: Potencial (posición) Flujo de calor: Temperatura (posición, tiempo) Difusión: Concentración (posición, tiempo) Flujo de fluidos: Velocidad, Presión ,Densidad (posición, tiempo) Elasticidad: Tensión, estiramiento (posición, tiempo) Variaciones de Poisson tienen coeficientes variables
49 Relación de ec. Poisson con la Gravedad, ElectrostáticaFuerza en la partícula en (x,y,z) debido a una partícula en (0,0,0) 3D: -(x,y,z)/r3, where r = (x2 +y2 +z2 ) 2D: -(x,y)/r2, where r = (x2 + y2 ) Fuerza es gradiente del potencial V (con signo opuesto) 3D: V = -1/r, V = (V/x, V/y, V/z) 2D: V = log r, V = (V/x, V/y) V satisface la ec. de Poisson (tarea)
50 Ec. Poisson en 1D: 2u/x2 = f(x)Discretizar d2u/dx2 = f(x) en la malla regular ui = u(i*h) Para conseguir [ u i+1 – 2*u i + u i-1 ] / h2 = f(x) Se resuelve como Tu = -h2 * f Para array desconocido u donde Grafo y “estencil” -1 2 -1 T =
51 Ec. Poisson en 2D Similar al caso 1D, pero la matriz T se vuelve:3D es análogo Grafo y “estencil” -1 -1 4 -1 T = -1 -1 -1 -1 6 -1 -1 -1
52 Algoritmos para Poisson 2D (3D) (N = n2 (n3) vars)Algoritmo Serial PRAM Memoria #Procs LU denso N3 N N2 N2 LU bandas N2 (N7/3) N N3/2 (N5/3) N (N4/3) Jacobi N2 (N5/3) N (N2/3) N N Inv.Explicit N2 log N N2 N2 Gradiente Conj. N3/2 (N4/3) N1/2(1/3) *log N N N SOR r/n N3/2 (N4/3) N1/2 (N1/3) N N LU dispersa 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 PRAM es un modelo de paralelismo ideal con costo de comunicación despreciable Referencia: James Demmel, Applied Numerical Linear Algebra, SIAM, 1997. Fix table in earlier lecture for 3D case! Band LU serial complexity: dimension * bandwidth^2 = n^2 * n^2 in 2D and n^3 * (n^2)^2 in 3D Band LU storage: dimension * bandwidth Band LU #procs = bandwidth^2 Jacobi complexity = complexity(SpMV) * #its, #its = cond# = n^2 in 2D and 3D CG complexity = [ complexity(SpMV)+complexity(dotprod) ] * #its, #its = sqrt(cond#) = n in 2D and 3D RB SOR complexity = complexity(SpMV) * #its, #its = sqrt(cond#) = n in 2D and 3D
53 Resumen de Algoritmos en malla n x n grid (n x n x n) (1)El orden de presentación es así (aprox.): Del más lento al más rápido en máquinas secuenciales. Del más general (trabaja en cualquier matriz) al más especializado (trabaj en matrices “parecidas” a T). LU denso: GE; trabaja en cualquier matriz N-by-N. LU en bandas: aprovecha el hecho que T es ≠ de cero solo en las N1/2 (N2/3) diagonales más cercanas a la diagonal principal. Jacobi: básicamente hace mult. matriz vector por T en el loop interno de un algoritmo iterativo. Inverso explícito: Asume que queremos resolver muchos problemas con la misma T, así que se precalcula y guarda inv(T) “gratis”, y se multiplica por ella (pero todavía caro!). Gradiente conjugada: usa mult. matriz vector, como Jacobi, pero aprovecha propiedades matemáticas de T que Jacobi no aprovecha. SOR (successive over-relaxation) rojo/negro: Variación de Jacobi que explota otras propiedades de T. Usado en multigrid. LU disperso: GE explotando la estructura particular de ceros de T. FFT (fast Fourier transform): Trabaja sólo en matrices muy similares a T. Multigrid: trabaja en matrices similares a T, que vienen de ec. elípticas. MINIMO: estadística del mejor Serial (tiempo para tener la respuesta en output); del mejor paralelo (tiempo para combinar los N outputs). Jacobi – converges if matrix diagonally dominant (our problem is nearly so) CG – converges if matrix symmetric positive definite (SPD) Red Black SOR – converges if matrix is SPD, analysis of convergence rate depends on more structure (RB ordering)
54 Resumen de Algoritmos en malla n x n grid (n x n x n) (2)Dimensión = N = n2 en 2D (=n3 en 3D) Ancho de banda = BW = n en 2D (=n2 en 3D) Número de Condición= k = n2 en 2D ó 3D #its = número de iteraciones SpMV = Sparse-matrix-vector-multiply LU densa: costo = N3 (serial) ó N (paralelo, con N2 procs) LU en bandas: costo = N * BW2 (serial) ó N (paralelo con BW2 procs) Jacobi: costo = costo(SpMV) * #its, #its = k. Gradiente conjugada: costo = (costo(SpMV) + costo(prod vectorial))* #its, #its = sqrt(k) SOR roja/negra: costo = costo(SpMV) * #its, #its = sqrt(k) Algoritmos Iterativos necesitan diversos supuestos para análisis de convergencia Jacobi – converges if matrix diagonally dominant (our problem is nearly so) CG – converges if matrix symmetric positive definite (SPD) Red Black SOR – converges if matrix is SPD, analysis of convergence rate depends on more structure (RB ordering)
55 Comentarios en mallas en la prácticaMallas regulares 1D, 2D, 3D Importantes como bloques básicos para fromar mallas más complicadas Se discutirán primero Mallas en la prácticas son irregular con frecuencia Mallas compuestas, consistentes de vrias mallas regulares “torcidas” unidas en los bordes Mallas no estructuradas, Con puntos de malla y coneccones arbitrarias Mallas adaptativas, cambian en resolución (tamaño, finura) a lo largo del proceso de solución para computar más donde se necesita más precisión
56 Malla compuesta en una estructura mecánica
57 Malla no estructurada: Perfil de ala 2D (NASA)
58 Refinamiento de malla adaptativa (AMR)Malla adaptativa alrededor de una explosión Fuente: John Bell y Phil Colella en LBL
59 Malla irregular : tubo que se estrecha (multigrid)
60 Método de Jacobi Para derivar este método, escribir Poisson como:u(i,j) = (u(i-1,j) + u(i+1,j) + u(i,j-1) + u(i,j+1) + b(i,j))/4 Hacer u(i,j,m) la aproximación de u(i,j) después de m pasos u(i,j,m+1) = (u(i-1,j,m) + u(i+1,j,m) + u(i,j-1,m) + u(i,j+1,m) + b(i,j)) / 4 O sea, u(i,j,m+1) es un promedio ponderado de vecinos Motivación: u(i,j,m+1) escogido para satisfacer exactamente la ecuación en el punto (i,j) Pasos para converger proporcional al tamaño del problema, N=n2 Ver Lección 24 de Por lo tanto, complejidad serial es de O(N2)
61 Convergencia de los métodos de vecinosJacobi involucra la computación de los vecinos más próximos en una malla nxn (N = n2) Así que toma O(n) = O(sqrt(N)) iteraciones para que la información se propague Ej., considerar un vector b = 0, excepto en el centro =1 La solución exacta es más o menos así: Aún en el mejor caso, cualquier computación de vecinos más próximos toma n/2 pasos para propagarse en una malla nxn
62 Convergencia de los métodos de vecinosBest 5 step solution is not a particular method, but uses the true solution 5 steps out from center, and zero elsewhere.
63 Paralelizando el método de JacobiSe reduce a una mult. matriz dispersa T (o casi) x vector U(m+1) = (T/4 - I) * U(m) + B/4 Cada valor de U(m+1) puede ser calculado por separado mantener 2 copias para las iteraciones m y m+1 Requiere que los valores en los límites se comuniquen Si cada procesador posee n2/p elementos para update Cantidad de data comunicada, n/p por vecino, es relativamente poco si n>>p
64 Optimización de Localidad en JacobiUpdate de vecinos más próximos en Jacobi tiene: Buena localidad espacial (en mallas regulares): fila / columna Mala localidad temporal: pocos flops por punto Ej: en mallas 2D: 4 sumas y 1 mult. vs. 5 lecturas y 1 escritura Para proc. Paralelo o secuencial, se puede reducir ops. memoria por un aumento de cómputos Idea: para cada submalla, hacer varias iteraciones de Jacobi Tamaño de submalla se achica con cada iteración, así que comenzar con la mayor Usado en uniprocesadores: Rescheduling para Localidad, Michelle Mills Strout: www-unix.mcs.anl.gov/~mstrout Algoritmos multigrid eficientes en uso de cache, Sid Chatterjee Y comput.paralelas, incluyendo “ajustes” de mallas Ian Foster et al, SC2001
65 Nodos fantasmas redundantes en JacobiResumen de optimización de jerarquía de memoria Puede ser usado en mallas no estructuradas Tamaño de región fantasma (y cómputo redundante) depende en velocidad de red o memoria vs. computación Para computar verde Copiar amarillo Computar azul
66 Sobre-Relajación Sucesiva (SOR)Similar a Jacobi: u(i,j,m+1) se calcula como una combinación lineal de los vecinos Coefic. numéricos y orden de updates son diferentes Basado en 2 mejoras respecto a Jacobi Usa los “valores más recientes” disponible para u, por que estos valores posiblemente son más precisos Gauss Seidel Actualizar el valor de u(m+1) “más agresivamente” a cada paso Primero, note que cuando evaluamos secuencialmente u(i,j,m+1) = (u(i-1,j,m) + u(i+1,j,m) … algunos de los valores para m+1 ya están disponibles u(i,j,m+1) = (u(i-1,j,latest) + u(i+1,j,latest) … donde “latest” es ya sea m ó m+1
67 Gauss-Seidel Actualizar de izq.-a-der. en orden por filas, es el algoritmo de Gauss-Seidel for i = 1 to n for j = 1 to n u(i,j,m+1) = (u(i-1,j,m+1) + u(i+1,j,m) + u(i,j-1,m+1) + u(i,j+1,m) + b(i,j)) / 4 No puede paralelizarse, debido a las dependencias, así que mejor se usa el orden “rojo-negro” forall puntos negros u(i,j) u(i,j,m+1) = (u(i-1,j,m) + … forall puntos rojos u(i,j) u(i,j,m+1) = (u(i-1,j,m+1) + … Para un gráfico general, usar coloreo de grafos Se puede usar conjuntos independ. máximos repetidos para colorear Grafo(T) es bipartito => 2 colores (rojo y negro) Nodos para cada color pueden ser actualizados simultáneamente Sigue siendo multip. vector matriz dispersa, usando submatrices
68 Sobre-Relajación Sucesiva (SOR)Este Gauss-Seidel rojo-negro converge dos veces más rápido que Jacobi, pero hay el doble de pasos paralelos, así que es lo mismo en la práctica. Para motivar la siguiente mejora, que el paso básico sea: u(i,j,m+1) = u(i,j,m) + corrección(i,j,m) Si la “corrección” es una buena dirección a donde moverse, entonces uno debería moverse aún más lejos en esa diección por un factor w>1 u(i,j,m+1) = u(i,j,m) + w * correction(i,j,m) Se le llama “successive overrelaxation” (SOR) Paraleliza como Jacobi (Multip. vector matriz dispersa…) Se puede probar que w = 2/(1+sin(p/(n+1)) ) da mejor convergencia Número de pasos para converger = complejidad paralela = O(n), en vez del O(n2) de Jacobi Complejidad Serial O(n3) = O(N3/2), en vez del O(n4) = O(N2) de Jacobi
69 Gradiente Conjugada (CG) para resolver A*x = bEste método puede usarse cuando la matriz A es simétrica, o sea, A = AT positiva definida, o sea cualquiera de estas 3 definiciones equivalentes: Todos los eigenvalores son positivos xT * A * x > 0 para todo vector x diferente de 0 Tiene una factorización Cholesky A = L*LT El algoritmo mantiene 3 vectores x = la solución aproximada, mejorada tras cada iteración r = el residual, r = A*x - b p = dirección de búsqueda, o la gradiente conjugada Costo de una iteración Mult. vector x matriz dispersa A (costo más importante) 3 productos vectoriales, 3 saxpy (escala*vector + vector) Converge en O(n) = O(N1/2) pasos, como SOR Complejidad serial = O(N3/2) Complejidad paralela = O(N1/2 log N), factor log N del producto vectorial
70 Gradiente Conjugada (CG) para resolver A*x = bAproximación inicial x (“initial guess”,adivinada) r = b – A*x, j=1 Repetir rho = rT*r … producto vectorial If j=1, p = r, else beta = rho/viejo_rho, p = r + beta*p, endif … saxpy q = A*p … MatMul dispersa alpha = rho / pT * q … producto vectorial x = x + alpha * p … saxpy r = r – alpha * q … saxpy viejo_rho = rho; j=j+1 Hasta que rho sea lo suficientemente pequeño
71 Cómo escoger el algoitmo iterativo para Ax=bUsar sólo mult.matriz-vector, prod.vect., saxpy, … Ver Simétrica? N S AT disponible? Definida? N N S S Almacenam. es “caro”? trata Minres o CG Eigenvalores conocidos? trata QMR N S N S trata CG o Chebyshev trata CGS o Bi-CGSTAB trata CG trata GMRES
72 Sumario de Jacobi, SOR y CGTodos (Jacobi, SOR, y CG) hacen mult. Matriz dispersa x vector Para Poisson, esto significa comunicación entre vecinos más próximos en una malla n-by-n grid Toma n = N1/2 pasos para que la información viaje a lo largo de una malla n-por-n Como la solución en un lado de la malla depende de la data en el otro lado de la malla, se necesitan métodos que muevan la información más rápidos FFT Multigrid
73 Resolviendo la Ec.Poisson con FFTMotivación: expresar la solución continua como serie de Fourier u(x,y) = Si Sk uik sen(p ix) sen(p ky) uik es el coeficiente de Fourier de u(x,y) La ecuación de Poisson 2u/x2 + 2u/y2 = b se convierte en: Si Sk (-pi2 - pk2) uik sen(p ix) sen(p ky) = Si Sk bik sen(p ix) sen(p ky) donde bik son los coefic. de Fourier de b(x,y) Unicidad: como la serie de Fourier es única, uik = bik / (-pi2 - pk2) Algoritmo continuo (Algoritmo discreto) Computa coefic.Fourier bik del lado derecho de la ecuación Aplicar FFT 2D a los valores b(i,k) en la malla Computa coefic.Fourier uik de la solución Divida cada transformada b(i,k) entre la función(i,k) Computa la solución u(x,y) usando los coef.Fourier Aplicar la FFT inversa 2D a los valores de b(i,k) divididos
74 FFT Serial Hacer i=sqrt(-1) y que los índices de matrices y vectores empiezen en 0. La Transformada discreta de Fourier de un vector v de m elementos es: F*v Donde F es la matriz m*m definida como: F[j,k] = v (j*k) Donde v es: v = e (2pi/m) = cos(2p/m) + i*sen(2p/m) v es un número complejo cuyo m-esima potencia vm =1 y por ello es una raiz m-esima de la unidad Ej., for m = 4: v = i, v2 = -1, v3 = -i, v4 = 1,
75 Usando la FFT 1D para filtrarSeñal = sen(7t) + .5 sen(5t) en 128 puntos Ruido = número aleatorio no mayor que 0.75 Filtrar cereando los componentes FFT < 0.25
76 Usando la FFT 2D para compresión de imágenesImagen = matriz 200x320 pixels (valores) Comprimido manteniendo el 2.5% más significativo de los componente del FFT, descarta el resto. Similar idea usada en formato jpeg (transf. del coseno) Jpeg breaks an image into 8x8 pixel blocks, does a DCT, quantizes the coefficients (round them so they can be Stored in very few bits), and stored these bits
77 Transformadas relacionadasLa mayor parte de aplicaciones necesitan multiplicar por F y por la inversa de F. Multiplicando por F y inversa(F) es esencialmente lo mismo (inversa(F) es la conjugada compleja de F dividida por n) Para resolver la ec. Poisson y varias otras, usaremos las variaciones de la FFT La transformada del seno – parte imaginaria de F La transformada del coseno – parte real de F Algoritmos similares, así que nos enfocaremos en la “forward FFT”
78 Algoritmo serial para la FFTComputar la FFT de un vector v con m elementos, F*v (F*v)[j] = S F(j,k) * v(k) = S v (j*k) * v(k) = S (v j)k * v(k) = V(v j) Donde V es el polinomio V(x) = S xk * v(k) m-1 k = 0 m-1 k = 0 m-1 k = 0 m-1 k = 0
79 Motivación del MultigridRecordemos 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 W(log n) pasos límite asintótico mínimo Convergencia a un error fijo < 1 toma W(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
80 Motivación de Multigrid
81 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 recursivamente, 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
82 Esa misma Gran Idea se usa en otros ladosReemplazar 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)
83 Multigrid usa “Divide y Vencerás” en 2 formas1ra 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)
84 Multigrid en 1D a grosso modoConsidere 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
85 Multigrid en 1D y 2D a grosso modoConsidere la grilla 2m+1 en 1D (2m+1 x 2m+1 en 2D) por simplicidad Sea P(i) el problema 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
86 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 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 (b(i)) El operador interpolación In
87 Algoritmo Multigrid de ciclo VFunction 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 (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
88 ¿Por qué es llamado ciclo V?Miren el dibujo del orden de las llamadas Después de un tiempo, el ciclo V luce así:
89 Complejidad de un Ciclo VEn 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: S 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
90 Full Multigrid (FMG) Intuición: Function FMG (b(m), x(m))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
91 Análisis del costo del Full MultigridUna V para cada llamada al Full MultiGrid Algunos usan Ws y otras formas Tiempo serial: S O(4 i) = O( 4 m) = O(# incógn.) Tiempo PRAM: S O(i) = O( m 2) = O( log2 # incógn.) m i=1 m i=1
92 Complejidad de resolver la Ec.PoissonTeorema: 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
93 El operador Solución S - DetallesEl operador Solución, S, 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.
94 Jacobi ponderado para amortiguar error de alta frecuenciaError inicial “Áspero” Componentes de alta frec. son fuertes Norma = 1.65 Error después de 1 Jacobi “Suave” Comp. de alta frec. menos fuertes Norma = 1.06 Error después de 2 Jacobi “Más Suave” Débil comp. de alta frec. Norma = .92, No disminuirá mucho más
95 Multigrid como algoritmo Divide y VencerásCada nivel en un ciclo V reduce el error en una parte del dominio de la frecuencia
96 El Operador Restricción R - DetallesEl operador restricción, R, 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) /2 * xfino(i) /4 * xfino(i+1) En 2D, se promedia con los 8 vecinos (N,S,E,W,NE,NW,SE,SW)
97 Operador Interpolación In
98 Convergencia del Multigrid en 1D
99 Convergencia del Multigrid en 2D
100 Multigrid paralela en 2DMultigrid 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)
101 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-k x 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 ) a …. Envía un número constante de mensajes a los procs. vecinos + O( 2j-k) b …. 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) a …. Envía un número constante de mensajes a los procs. vecinos + O(1) b … Número de words enviados Sume para todos los niveles 1,2,… y todos los ciclos V en FMG!!!
102 Comparición de Métodos (en O(.))# Flops # Mensajes # Words enviadas MG N/p (log N) (N/p)1/2 + log p * log N log p * log N FFT N log N / p p1/ N/p SOR N3/2 /p N1/ 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 a >> b
103 Consideraciones prácticasEn 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 ~= 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
104 Multigrid en un Mesh adaptativoPara 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
105 Aplicaciones multibloqueResolver el sistema de ecuaciones en una unión de rectángulos subproblema de AMR Ej:
106 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
107 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.
108 Multigrid en un Mesh no estructuradoOtro 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)
109 Multigrid en una Mesh no estructuradaSe 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?
110 Fuente de Mesh no estructurada (elem.finito): VertebraFuente: M. Adams, H. Bayraktar, T. Keaveny, P. Papadopoulos, A. Gupta
111 Multigrid para análisis elástico no lineal de huesoFuente: M. Adams et al Test mecánico de propiedades del material Imagen 3d mFE mesh 2.5 mm3 elementos 44 mm Tomografía por resolución 22 mm Hasta 537M variables 4088 Procesadores (ASCI White) 70% eficiencia en paralelo