1 / 41

Operaciones con Matrices Densas

Operaciones con Matrices Densas. Glen D. Rodríguez R. Algoritmos Paralelos Basado en material de J. Demmel. Revisión de BLAS. Bloques básicos de operaciones de alg. Lineal. Las librerías paralelas llaman a las funciones seriales en cada CPU: deben ser veloces!

zahir-davenport
Download Presentation

Operaciones con Matrices Densas

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Operaciones con Matrices Densas Glen D. Rodríguez R. Algoritmos Paralelos Basado en material de J. Demmel

  2. Revisión de BLAS • Bloques básicos de operaciones de alg. Lineal. • Las librerías paralelas llaman a las funciones seriales en • cada CPU: deben ser veloces! • Recordar la eficiencia del algor. q = # flops / # refs mem • A mayor q, más rápido el algoritmo en computadorasreales con jerarquía de memoria • “axpy” (ax plus y): y = a*x + y, con a escalar, x e y vectores

  3. Varias particiones de data en proc. paralelos para matrices 1) Bloques de columnas en 1D 2) Columnas cíclicas en 1D 4) Cambiar columnas por filas en los casos 1,2,3 y salen 3 casos análogos 3) Bloques cíclicos de cols. 1D Generalizar otros 6) Bloques cíclicos 2D 5) Bloques (sub matrices) en 2D

  4. Multiplicación paralela matriz vector • Computar y = y + A*x, con A matriz densa • Partición: • Bloques de filas en 1D • A(i) es el bloque de filas tamaño n por n/p asignado al procesador i, • x(i) e y(i) son los segmentos de x,y asignados al proc.“i” • Algoritmo: • Para todos los procesadores i • Broadcast x(i) • Computar y(i) = y(i) + A(i)*x • Algoritmo use la formula y(i) = y(i) + A(i)*x = y(i) + Sj A(i,j)*x(j) P0 P1 P2 P3 x P0 P1 P2 P3 y

  5. Multip. matriz vector y = y + A*x • Una partición por columnas elimina el broadcast de x • Pero implicaría un mpi_reduce para obtener el valor de y • Una partición de bloques 2D use broadcast y reducción, ambos en un subconjunto de procesadores • sqrt(p) si tenemos una distribución cuadrada de procesadores (2x2, 3x3, 4x4, 5x5, etc.) P0 P1 P2 P3 P0 P1 P2 P3 P4 P5 P6 P7 P8 P9 P10 P11 P12 P13 P14 P15

  6. Multiplicación de matrices en paralelo • Computar C=C+A*B • Usando algoritmo básico: 2*n3 Flops • Variables: • Partición de la data • Topología de las máquinas • Scheduling (calendarización) de la comunicación • Uso de modelos de performance para diseñar el algoritmo • Tiempo comunicación= “latencia” + #words *tiempo-por-word = a + n*b • Eficiencia (en cualquier modelo): • serial time / (p * parallel time) = speedup / p • speedup perfecto (lineal) es “p” efficiency = 1 = 100%

  7. Algoritmo trivial • Si hay suficiente espacio de memoria en cada procesador para guardar todos los elementos de A, B y C; entonces dividir solo el cálculo pero no la data. C A B + x Todos guardan la totalidad de A, B, C. Los valores iniciales se comunican con broadcast Para matrices de tamaño nxn, cada procesador debe usar memoria= tamaño word x 3n2

  8. Mejora al trivial • Se ahorra algo de memoria si se divide la memoria de C y de A, pero no de B. • La división es por bloques de filas 1D • Ver programa ejemplo Proc. 0 Proc. 1 Proc. 2 Proc. 3 … Proc. p-2 Proc. p-1 A Para matrices de tamaño nxn, si hay p procesadores, cada procesador debe usar memoria = tamaño word x (n2 + 2n2/p)

  9. p0 p1 p2 p3 p4 p5 p6 p7 Multiplic. matrices con partición columnas 1D • Asumir matrices de n x n con n divisible por p • A(i) es el bloque de columnas tamaño n por n/p asignado al proc.”i” (lo mismo para B(i) y C(i)) • B(i,j) es un sub-bloque tamaño n/p por n/p de B(i) • desde la fila j*n/p hasta la fila ((j+1)*n/p)-1 • Algoritmo usa la formula C(i) = C(i) + A*B(i) = C(i) + Sj A(j)*B(j,i) Suposición razonable para análisis, no para programar

  10. Partición 1D en Bus ó Anillo • Algoritmo usa la formula C(i) = C(i) + A*B(i) = C(i) + Sj A(j)*B(j,i) • Primero, considerar una máquina que usa bus y sin broadcast: solo un par de procesadores pueden comunicarse a la vez (ej.: ethernet) • Luego, considerar procesadores en un anillo: todos los procesadores se pueden comunicar con sus vecinos inmediatos a la vez.

  11. MatMul: Partición 1D en Bus sin Broadcast Algoritmo simplón: C(myproc) = C(myproc) + A(myproc)*B(myproc,myproc) for i = 0 to p-1 for j = 0 to p-1 except i if (myproc == i) send A(i) al procesador j if (myproc == j) receive A(i) del procesador i C(myproc) = C(myproc) + A(i)*B(i,myproc) barrier Costo del loop interno : computación: 2*n*(n/p)2 = 2*n3/p2 comunicación: a + b*n2 /p

  12. Continua… Costo del loop interno: computación: 2*n*(n/p)2 = 2*n3/p2 comunicación: a + b*n2 /p … aprox. Solo un par de procesadores (i , j) están activos en una iteración, y de ellos, sólo i esta computando => el algoritmo es prácticamente serial Tiempo de ejecución: = (p*(p-1) + 1)*computación + p*(p-1)*comunicación ~= 2*n3 + p2*a + p*n2*b Es peor que el tiempo serial y crece con p. Quién lo usaría si es malo?  ojo, no acelera tiempo de ejecución pero por lo menos aumenta la memoria disponible.

  13. Matmul en columnas 1D y en un anillo • Pares de procesadores se pueden comunicar a la vez Copiar A(myproc) en Tmp C(myproc) = C(myproc) + Tmp*B(myproc , myproc) for j = 1 to p-1 Send Tmp al procesador myproc+1 mod p Receive Tmp del procesador myproc-1 mod p C(myproc) = C(myproc) + Tmp*B( myproc-j mod p , myproc) • Misma idea que en el paralelismo de N-body, usar un buffer de variables temporales • Tiempo del loop interno = 2*(a + b*n2/p) + 2*n*(n/p)2

  14. Matmul en columnas 1D y en un anillo • Tiempo del loop interno= 2*(a + b*n2/p) + 2*n*(n/p)2 • Tiempo Total = 2*n* (n/p)2 + (p-1) * Tiempo loop interno ~ 2*n3/p + 2*p*a + 2*b*n2 • Optimo para partición 1D en Ring o Bus, aún con Broadcast: • Speedup perfecto para aritmética • A(myproc) debe moverse a los demás procesadores, cuesta por lo menos: (p-1)*costo de enviar n*(n/p) words • Eficiencia paralela = 2*n3 / (p * Tiempo Total paralelo) = 1/(1 + a * p2/(2*n3) + b * p/(2*n) ) = 1/ (1 + O(p/n)) • Se acerca a 1 si n/p crece (o si a , b se achican)

  15. MatMul en particiones 2D • Considerar procesadores en malla 2D (física o lógica) • Procs. se pueden comunicar con 4 vecinos inmediatos • Broadcast a lo largo de filas y columnas • Asumir p procesadores forman una malla de sxs p(0,0) p(0,1) p(0,2) p(0,0) p(0,1) p(0,2) p(0,0) p(0,1) p(0,2) = * p(1,0) p(1,1) p(1,2) p(1,0) p(1,1) p(1,2) p(1,0) p(1,1) p(1,2) p(2,0) p(2,1) p(2,2) p(2,0) p(2,1) p(2,2) p(2,0) p(2,1) p(2,2)

  16. Algoritmo de Cannon k … C(i,j) = C(i,j) + S A(i,k)*B(k,j) … se asume que s = sqrt(p) es entero forall i=0 to s-1 … “mover” A desplaz. circular-izq. de la fila i de A en “i” posiciones … tal que A(i,j) es sobreescrito por A(i,(j+i)mod s) forall i=0 to s-1 … “mover” B desplaz.circular-arriba de columna i de B en “i” pos. …tal que that B(i,j) es sobreescrito por B((i+j)mod s), j) for k=0 to s-1 … secuencial forall i=0 to s-1 and j=0 to s-1 … todos procs. en paralelo C(i,j) = C(i,j) + A(i,j)*B(i,j) despl.circular-izq para cada fila de A en 1 despl.circular-arriba para cada columna de B en 1

  17. MatMul con Cannon C(1,2) = A(1,0) * B(0,2) + A(1,1) * B(1,2) + A(1,2) * B(2,2)

  18. Paso inicial para mover Matrices en Cannon A(0,0) A(0,1) A(0,2) B(0,0) B(0,1) B(0,2) • Input inicial por bloques • Después de 1er movimiento, antes de multiplicar A(1,0) A(1,1) A(1,2) B(1,0) B(1,1) B(1,2) A(2,0) A(2,1) A(2,2) B(2,0) B(2,1) B(2,2) A(0,0) A(0,1) A(0,2) B(0,0) B(1,1) B(2,2) A(1,1) A(1,2) A(1,0) B(1,0) B(2,1) B(0,2) A(2,2) A(2,0) A(2,1) B(2,0) B(0,1) B(1,2)

  19. A(0,0) A(0,1) A(0,2) B(0,0) B(1,1) B(2,2) A(1,1) A(1,2) A(1,0) B(1,0) B(2,1) B(0,2) A(2,2) A(2,0) A(2,1) B(2,0) B(0,1) B(1,2) Siguientes movimientos en Cannon • 1er. paso • 2do. paso • 3er. paso A(0,1) A(0,2) A(0,0) B(1,0) B(2,1) B(0,2) A(1,2) A(1,0) A(1,1) B(2,0) B(0,1) B(1,2) A(2,0) A(2,1) A(2,2) B(0,0) B(1,1) B(2,2) A(0,2) A(0,0) A(0,1) B(2,0) B(0,1) B(1,2) B(0,0) B(1,1) B(2,2) A(1,0) A(1,1) A(1,2) A(2,1) A(2,2) A(2,0) B(1,0) B(2,1) B(0,2)

  20. Costo del Algoritmo de Cannon forall i=0 to s-1 … recordar que s = sqrt(p) despl.circ.izq. fila i de A en i … costo = s*(a + b*n2/p) forall i=0 to s-1 desp.circ.arr. columna i de B en i … costo = s*(a + b*n2/p) for k=0 to s-1 forall i=0 to s-1 and j=0 to s-1 C(i,j) = C(i,j) + A(i,j)*B(i,j) … costo = 2*(n/s)3 = 2*n3/p3/2 despl.circ.izq. cada fila de A en 1 … costo = a + b*n2/p desp.circ.arr. cada columna de B en 1 … costo = a + b*n2/p • Tiempo Total = 2*n3/p + 4*s* + 4**n2/s • Eficiencia Paralela = 2*n3 / (p * Tiempo Total paralelo) • = 1/( 1 + a * 2*(s/n)3 + b * 2*(s/n) ) • = 1/(1 + O(sqrt(p)/n)) • Se acerca a 1 si crece n/s = n/sqrt(p) = sqrt(data por procesador) • Mejor que 1D, que tiene eficiencia= 1/(1 + O(p/n))

  21. Pros y Contras de Cannon • Se optimiza la computación local • Difícil de generalizar para: • Valores de p que no son cuadrados perfectos • Matrices A y B no son cuadradas • Dimensiones de A, B que no son divisibles por s=sqrt(p) • A y B no “alineadas” en la forma en que están guardadas en los procesadores • Bloques cíclicos • Consume memoria (copias extra de matrices locales)

  22. Algoritmo SUMMA • SUMMA = Scalable Universal Matrix Multiply • Un poco menos eficiente, pero es más sencillo y más fácil de generalizar • Link: slides de van de Geijn y Watts • www.netlib.org/lapack/lawns/lawn96.ps • Ideas similares aparecieron en muchas ocasiones • Usado en PBLAS = Parallel BLAS • www.netlib.org/lapack/lawns/lawn100.ps

  23. SUMMA B(k,j) k j k Proc(1,1) Proc(2,1) * = C(i,j) i A(i,k) Proc(pr,pc) • i, j representan todas las filas, columnas asignadas a un procesador • k es una sola fila o columna; o un bloque de b filas o columnas • C(i,j) = C(i,j) + Sk A(i,k)*B(k,j) • Asumir una malla de pr por pc procesadores (pr = pc = 4 en el ejemplo de arriba) . El número de procesadores no necesariamente debe ser cuadrado perfecto.

  24. SUMMA B(k,j) k j k * = C(i,j) i A(i,k) For k=0 to n-1 … ó (n/b)-1 donde b es el tamaño de bloque … b= # columnas en A(i,k) y # filas en B(k,j) for all i = 1 to pr… en paralelo dueño de A(i,k) le hace broadcast a todos los procs. de esa fila for all j = 1 to pc… en paralelo dueño de B(k,j) le hace broadcast a todos los procs. de esacol. Receive A(i,k) en Acol Receive B(k,j) en Brow C_myproc = C_myproc + Acol * Brow

  25. Performance de SUMMA • Solo para simplificar análisis, asumir s = sqrt(p) For k=0 to n/b-1 for all i = 1 to s … s = sqrt(p) dueño de A(i,k) le hace bc. a todos los procs. de la fila … tiempo = log s *( a + b * b*n/s), usando árbol for all j = 1 to s dueño de B(k,j) le hace bc. a todos los procs. de la col. … tiempo = log s *( a + b * b*n/s), usando árbol Receive A(i,k) en Acol Receive B(k,j) en Brow C_myproc = C_myproc + Acol * Brow … tiempo= 2*(n/s)2*b • Tiempo Total = 2*n3/p + a * log p * n/b + b * log p * n2 /s

  26. Performance de SUMMA • Tiempo total = 2*n3/p + a * log p * n/b + b * log p * n2 /s • Eficiencia Paralela = • 1/(1 + a * log p * p / (2*b*n2) + b * log p * s/(2*n) ) • ~Mismo término b de Cannon, excepto por el factor log p • log p crece lento, así que no es problema • Término con latencia (a) puede ser mayor, dependiendo de b • Cuando b=1, se tiene a * log p * n • Si b se acerca a n/s, término se achica hacia • a * log p * s (log p veces Cannon) • Storage temporal crece como 2*b*n/s • Se puede cambiar b para disminuir latencia a costa de la necesidad de más memoria

  27. Librería paralela ScaLAPACK

  28. PDGEMM = rutina PBLAS para multipl. matrices Observaciones: Para N fijo, si P aumenta, Mflops aumenta, pero menos que 100% eficiencia Para P fijo, si N aumenta, Mflops (eficiencia) aumenta DGEMM = rutina BLAS para multip. matrices Max. Velocidad para PDGEMM = # Procs * veloc. de DGEMM Observaciones (mismas arriba): Eficiencia por lo menos 48% Para N fijo, si P aumenta, eficiencia cae Para P fijo, si N aumenta, eficiencia aumenta

  29. Notación para el algoritmo de Fox Los índices comienzan en 0 para las sub matrices, pero comienzan en 1 para filas y columnas de procesadores

  30. Algoritmo de Fox (Broadcast, Multiply, y Roll)

  31. Primer paso -- index n=0 en suma de sub-bloque

  32. Segundo paso -- index n=1 en suma de sub-bloque

  33. Segundo paso, continuado

  34. El algoritmo completo para un elemento

  35. MPI: grupos de procesadores y comunicación colectiva • Necesitamos “broadcasts parciales” a lo largo de las filas en Fox (y de las columnas también en SUMMA) • Y rolls (desplazamientos en 1) en las columnas en Fox (parecido en Cannon) • Estas son “comunicaciones colectivas” • “Broadcasts de filas” son broadcasts es sub-grupos especiales de procesadores • Rolls se hacen como una variante del MPI_SENDRECV con “condiciones de frontera ciclicas” • También hay rutinas especiales MPI para definir la malla 2D de procesadores

  36. Broadcast en el caso de la Matriz completa • MatMul usando Fox o Summa con MPI emplea broadcast como comunicación básica • Usaremos esta aplicación para discutir 3 enfoques del broadcast • Trivial • Logarítmico • Pipe (tubo) • Que tiene diferente performance dependiendo en el tamaño del mensaje y la arquitectura del hardware.

  37. Implementación Broadcast trivial y del logarítmico

  38. Broadcast de pipe o tubo • En el caso que el tamaño del mensaje es grande, otras implementaciones son posibles o deseables, ya que será necesario hacer broadcast del mensaje previa sub división en una secuencia de mensajes más chicos. • El broadcast puede establecer un camino (o varios) desde el procesador origen y que visita a los demás procesadores de un grupo. • El mensaje es enviado desde la fuente a lo largo del camino en un “tubo”, donde cada procesador recibe un bloque del mensaje desdensupredecesor y lo envia luego a su sucesor. • La performance de este broadcast es por lo tanto el tiempo para enviar el mensaje al procesador al final del camino más la latencia de inicializar y terminar el “tubo”. • Tiempo = (Tamaño del mensaje +Tamaño de un paquete(√N – 2)) *tcomm • Para granos suficientemente gruesos, el pipe es mejor que el log. • Alta latencia es mala para este enfoque de “tubo” • MPI usa logarítmico, pero se puede “adoptar” el trivial.

  39. Esquema de Operación

  40. Análisis de performance

  41. Resumen de MatMul paralela • 1D • Bus sin broadcast – más lento que serial • Comunicación al vecino más próximo en anillo (o bus con broadcast): Eficiencia = 1/(1 + O(p/n)) • 2D • Cannon • Eficiencia = 1/(1+O(a * ( sqrt(p) /n)3+b* sqrt(p) /n)) • Difícil de generalizar para arbitrarios p, n, bloques cíclicos, alineamientos • SUMMA • Eficiencia = 1/(1 + O(a * log p * p / (b*n2) + b*log p * sqrt(p) /n)) • Más genérico • b chico => menos memoria, menos eficiencia • b grande => más memoria, más eficiencia

More Related