Paralelización mediante procesadores gráficos (GPU) del cálculo de la dinámica de electrones en plasma confinado magnéticamente

Máster Universitario en Computación Paralela y Distribuida Departamento de Sistemas Informáticos y Computación Paralelización mediante procesadores gráficos (GPU) del cálculo de la dinámica de electrones en plasma confinado magnéticamente TRABAJO FINAL DE MÁSTER Autor: Tomás Navarro Cosme Director: José Enrique Román Moltó Septiembre de 2015 . Abstract CUDA parallelization of a program, originally written in Fortran with MPI, for the calculation of the dynamic of electrons, taking particular attention to the formation of runaway electrons in a mangentically confined plasma, under the effect of an electric field, using the Langevin’s approximate calculation method. The final implementation will retain the most basic part source code written in Fortran, a language with a large base of use in the field of Physics, to show its compatibility with the newest NVIDIA’s Graphical Processors Units via CUDA. The implementation is progressive, showing successive technical improvements that can be applied to such codes, analyzing the obtained speedups. Resumen Paralelización en CUDA de un programa, originalmente escrito en Fortran usando MPI, para el cálculo de la dinámica de los electrones, con especial atención a la aparición de electrones fugitivos (runaway), dentro de un plasma confinado magnéticamente y que se mueven por efecto de un campo eléctrico, utilizando el método de cálculo aproximado de Langevin. En la paralelización se mantendrá la parte básica de Fortran, lenguaje con una gran base de uso en el campo de la Física, para mostrar su compatibilidad de uso con los aceleradores gráficos de NVIDIA usando CUDA. La implementación es progresiva, mostrando las sucesivas técnicas y mejoras que pueden aplicarse a este tipo de códigos, analizando las aceleraciones de cálculo obtenidas. Palabras clave: CUDA, Fortran, GPU, GPGPU, computación numérica, código científico, HPC, plasma, Langevin, electrones runaway i ÍNDICE ABREVIADO Índice abreviado · iii Índice general · v Índice de figuras · ix Índice de tablas · ix Índice de listados de código · x 1 Introducción · 1 2 Conocimientos previos · 3 3 Ecuaciones de Langevin para describir un plasma de partículas y la generación de electrones runaway · 13 4 Implementación de la paralelización con CUDA · 17 5 Resultados · 43 6 Trabajo futuro · 71 7 Conclusiones · 73 Glosario · 75 Siglas · 77 Bibliografía · 79 iii ÍNDICE GENERAL Índice abreviado iii Índice general v Índice de figuras ix Índice de tablas ix Índice de listados de código x 1 Introducción 1 2 Conocimientos previos 3 2.1 La simulación computacional en la Física . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.2 Características de los códigos científicos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.3 GPGPU .................................................. 7 2.4 CUDA ................................................... 9 2.4.1 Lenguaje CUDA C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.4.2 Tipos de código en CUDA C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.4.3 Compilación de aplicaciones CUDA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 3 Ecuaciones de Langevin para describir un plasma de partículas y la generación de electrones runaway 13 3.1 El plasma y su aplicación práctica: reacción controlada de fusión . . . . . . . . . . . . . . . . 13 3.2 Confinamiento artificial del plasma . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.3 Física de partículas: la ecuación de Langevin . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 4 Implementación de la paralelización con CUDA 17 4.1 Selección del lenguaje auxiliar de programación . . . . . . . . . . . . . . . . . . . . . . . . . . 17 4.2 Utilización de wrappers en Cpara ser invocados por Fortran . . . . . . . . . . . . . . . . . . 18 4.2.1 Respeto de la nomenclatura (name mangling) de los símbolos según compilador . 19 4.2.2 Esquema general de llamadas a CUDA desde Fortran a través de wrappers C. . . . 21 4.3 Diferencias programáticas importantes entre Fortran yCyCUDA-C . . . . . . . . . . . . . . 23 4.4 Implementación en CUDA del código de Langevin originalmente en Fortran . . . . . . . . 24 4.4.1 Aplanamiento de subrutinas por inclusión en un kernel aglutinador . . . . . . . . . 24 4.4.2 Minimización de la transferencia de datos entre CPU yGPU . . . . . . . . . . . . . . 25 4.4.3 Utilización adecuada de los generadores de números aleatorios de CUDA . . . . . . 26 4.4.4 Mejoras en el código: utilización de funciones estándar frente a métodos numéricos pesados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 v 4.4.5 Algoritmos de reducción en la GPU . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 4.4.5.1 Reducción global iterativa . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 4.4.5.2 Reducción mediante el uso de operaciones atómicas . . . . . . . . . . . . 29 4.4.6 Mejora en el algoritmo de clasificación de velocidades . . . . . . . . . . . . . . . . . . 30 4.4.7 Detección de errores (bugs) en el código . . . . . . . . . . . . . . . . . . . . . . . . . . 31 4.5 Desacople de computación y comunicación usando OpenMP . . . . . . . . . . . . . . . . . . 33 4.5.1 Desacople usando OpenMP con hilo en espera activa . . . . . . . . . . . . . . . . . . 34 4.5.2 Implementación final con cerrojos de OpenMP (omp_locks) . . . . . . . . . . . . . 36 4.5.2.1 Consideraciones con la adición de OpenMP en el código Fortran . . . . . 40 4.6 Valoración cualitativa del proceso de paralelización con CUDA . . . . . . . . . . . . . . . . . 40 5 Resultados 43 5.1 Pruebas realizadas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 5.1.1 Primera etapa: optimización del kernel principal según los parámetros propios de su invocación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 5.1.2 Segunda etapa: medición de prestaciones de las versiones finales en un entorno HPC real . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 5.2 Entorno de ejecución . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 5.2.1 Primera etapa: optimización de los parámetros de invocación del kernel principal . 46 5.2.2 Segunda etapa: entorno HPC real, MinoTauro . . . . . . . . . . . . . . . . . . . . . . . 47 5.3 Ejecución de las pruebas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 5.3.1 Parámetros óptimos de invocación al kernel principal advance . . . . . . . . . . . . 47 5.3.1.1 Breve descripción de la prueba . . . . . . . . . . . . . . . . . . . . . . . . . . 47 5.3.1.2 Resultados obtenidos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 5.3.1.3 Fenómenos observados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 5.3.1.4 Conclusión de la prueba . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 5.3.2 Aceleración en un entorno de ejecución de altas prestaciones . . . . . . . . . . . . . 51 5.3.2.1 Caracterización del comportamiento del programa original . . . . . . . . 52 5.3.2.2 Caracterización básica de la implementación consistente en la paralelización con CUDA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 5.3.2.2.1 Prueba de corrección . . . . . . . . . . . . . . . . . . . . . . . . . . 55 5.3.2.2.2 Aceleración y eficiencia respecto del código en CPU . . . . . . . 55 5.3.2.3 Caracterización básica de la implementación con CUDA yOpenMP . . . . 57 5.3.2.4 Pruebas masivas de escalabilidad de las implementaciones con CUDA . . 58 5.3.2.4.1 Resultados obtenidos . . . . . . . . . . . . . . . . . . . . . . . . . . 58 5.3.2.4.2 Escalabilidad de las implementaciones . . . . . . . . . . . . . . . 61 5.3.2.4.3 Comparación directa de las implementaciones . . . . . . . . . . 64 5.3.2.4.4 Cálculo de la aceleración en una carga grande . . . . . . . . . . . 66 5.4 Valoración . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 5.4.1 Valoración numérica cuantitativa . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 5.4.2 Valoración cualitativa de los resultados . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 5.4.2.1 Facilidad de programación . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 5.4.2.2 Necesidad de la reformulación de los algoritmos . . . . . . . . . . . . . . . 68 5.4.2.3 Depuración semántica de las implementaciones históricas . . . . . . . . . 69 5.4.2.4 Selección de códigos a paralelizar . . . . . . . . . . . . . . . . . . . . . . . . 69 6 Trabajo futuro 71 vi 7 Conclusiones 73 Glosario 75 Siglas 77 Bibliografía 79 vii Paralelización mediante GPU de la dinámica de electrones en plasma confinado magnéticamente Finalmente se expondrán otras tecnologías de paralelización, algunas de las cuales también han sido utilizadas e incluso incorporadas. Así por ejemplo tenemos MPI, una biblioteca de programación basada en el paradigma de memoria distribuida, mediante la utilización de paso de mensajes; y también OpenMP, que está especialmente pensado para memoria compartida. Últimamente han surgido otras posibilidades que intentan facilitar la labor de la paralelización abstrayéndose un poco más del hardware existente por debajo, de entre los cuales destacaremos el caso de OpenACC. Por otra parte, la parte física del problema, por su especial relevancia y entidad, será tratada en el siguiente Capítulo 3 “Ecuaciones de Langevin para describir un plasma de partículas y la generación de electrones runaway”, aunque sin entrar en excesivo detalle científico, sino más descriptiva de su significado, así como de las principales características a tener en cuenta a la hora de paralelizar el algoritmo de simulación. 2.1 La simulación computacional en la Física En el ámbito de la física de particulas es muy habitual el uso de la simulación computacional. Esta utilización no se realiza para su estudio directo, sino como método de comprobación de las de hipótesis y nuevos desarrollos teóricos. En este campo de la física de partículas elementales, el método de observación está limitado por la máxima precisión alcanzable, determinada por el principio de indeterminación de Heisenberg, y ese límite físico acota también la extracción de conocimiento a partir de la experimentación. A esto se une la complejidad de las actuales teorías, que deben compaginar tanto fenómenos bien conocidos y teóricamente bien desarrollados, como la gravitación, la dinámica y el electromagnetismo, con los efectos cuánticos y la aplicación de la teoría de la relatividad, ya que es habitual trabajar con partículas que se mueven a altas velocidades (esto es, son de alta energía). Por todo ello, cuando se trabaja a nivel de partículas tales como electrones, neutrones, protones y partículas alfa, las ecuaciones teóricas no pueden resolverse analíticamente. El procedimiento habitual es adaptarlas al caso particular en estudio, realizando las oportunas simplificaciones e ndo que no afecten al resultado a obtener, y finalmente, con el conjunto de fórmulas simplificadas, y por lo tanto aproximadas, realizar intensos cálculos de simulación para comprobar si se pueden reproducir los resultados obtenidos experimental o empíricamente. En caso de que así sea, ese conjunto de fórmulas aproximadas puede someterse a un nuevo proceso de simulación que se extienda a casos más complejos que difícil o muy onerósamente pueden contrastarse experimentalmente. De esta forma se obtiene información valiosa que permite acelerar el proceso tecnológico y de ingeniería para la construcción de nuevas instalaciones de experimentación. Este es el caso de la física de plasma, en donde partículas elementales de los átomos más ligeros se hallan confinadas por el efecto de fuertes campos electromagnéticos a altas temperaturas, inonizadas, libres entre sí, y produciéndose choques ocasionales de alta energía, que, en condiciones de tiempo, temperatura y densidad suficiente, pueden provocar el inicio del proceso de la fusión nuclear. Desde el punto de vista informático, los procesos de simulación se caracterizan por ser computacionalmente intensos, dado que se debe tratar un número lo más grande posible de partículas, y porque generalmente interesa la dinámica del sistema, esto es, cómo evoluciona desde un determinado estado inicial. Esto implica generalmente que se debe calcular el estado inicial, después suponer un incremento de tiempo pequeño, para intentar linealizar lo más posible la dinámica del sistema, esto es, que los errores por las aproximaciones realizadas no se propaguen en exceso, pero no tan pequeño como para suponer una inmensa cantidad de incrementos “diferenciales” del tiempo, y con ello, una gran cantidad de nuevos estados a calcular hasta alcanzar el estado final buscado (generalmente, ver si alcanza un estado estacionario, o analizar y comprender cómo se comporta si no. Por lo tanto, las ecuaciones son complejas, con muchos cálculos por cada estado, y muchos estados a computar. Esto puede complicarse en función de la complejidad del desarrollo teórico utilizado. De esta forma, interesa obtener unas ecuaciones que calculen el estado 4 Capítulo 2. Conocimientos previos 2.2. Características de los códigos científicos de cada partícula de forma independiente al comportamiento individual de cada una del resto, puesto que si existe interdependencia entre las partículas, el orden del número de operaciones a computar por cada estado pasa de ser de O(N), siendo Nel número de partículas, a orden O(N·m), siendo mel número de particulas vecinas con las que interactúa, o incluso O(N2) si todas interactúan con todas en cada estado. Esta alta necesidad de poder de cálculo implica que, si se utiliza un único procesador, el tiempo de cálculo sea excesivamente elevado, lo que conlleva problemas por: incremento del retraso entre diferentes simulaciones, con la consiguiente paralización de la validación de los desarrollos teóricos aproximativos realizados; mayor riesgo de finalización abrupta de la simulación, tanto por fenómenos externos (cortes de luz) como internos (sobrecalentamiento del nodo, error de software,...). La solución que se ha venido aplicando en los últimos 20 años ha sido la paralelización de los procesos, bien sea mediante OpenMP (en sistemas multiprocesadores de memoria compartida), o mediante MPI (paso de mensajes entre nodos con memoria distribuida). Todo ello es posible utilizando el lenguaje Fortran, que tiene su principal nicho de uso dentro del terreno de cómputo científico, puesto que: está pensado para transcribir las fórmulas matemáticas; permite la programación estructurada, generalmente idónea para implementar los habituales algoritmos de cálculo científico; permite un uso muy simplificado de las operaciones vectoriales y matriciales, muy habituales en ciencia; está orientado al cómputo, con compiladores muy optimizados al efecto; dispone de librerías específicas de cálculo numérico, algebraico y científico muy optimizadas, robustas y contrastadas; tiene una muy buena implementación y soporte en los dos históricamente principales paradigmas de programación paralela, OpenMP yMPI, dado que sus estándares se desarrollan e implementan teniendo muy en cuenta las particularidades y ventajas de Fortran. Este último aspecto es fundamental para continuar la herencia de muchos programas y bibliotecas de cálculo matemático y científico desde los años 60 hasta la actualidad, lo cual tiene aspectos muy positivos (programas contrastados, eficientes y robustos), pero que sutilmente esconde algunos inconvenientes que conviene tener en cuenta y analizar, como se reveló duante la fase de transcripción y adaptación al nuevo paradigma SIMD de CUDA, tal como se verá en este trabajo. 2.2 Características de los códigos científicos Por código científico entenderemos en este trabajo aquel programa o biblioteca que implemente directamente un determinado algoritmo de computación científica. Este algoritmo de computación científica será generalmente un conjunto de ecuaciones interdependientes que permiten resolver un determinado problema físico. Generalmente, este conjunto de fórmulas proviene del desarrollo teórico en una determinada área científica, y más específicamente, nos centraremos en aquellas que sirven para la simulación de sistemas físicos. Estos problemas tienen la característica común de que son lo suficientemente complejos como para que el conjunto de fórmulas derivadas estrictamente de la teoría sean demasiado complejas para su resolución analítica. Esta imposibilidad de resolución analítica impide que se pueda analizar el comportamiento de estos sistemas en diferentes situaciones prácticas. Al objeto de poder estudiarlos, se hace necesario realizar determinadas suposiciones que permiten aligerar la complejidad del conjunto de ecuaciones, mediante la inclusión de aproximaciones razonables para el caso concreto en estudio. Tras estas aproximaciones, es muy frecuente que se obtenga otro conjunto de ecuaciones más simple, pero aún así irresolubles o de difícil interpretación. Para resolverlas se recurren a métodos numéricos mediante la implementación de programas informáticos específicos para ese problema. Estos códigos científicos se ejecutan en múltiples instancias, puesto que es frecuente que estas funciones tengan términos estocásticos, esto es, que para capturar la variabilidad de las condiciones externas o de la complejidad interna innata del sistema, y ante la imposibilidad o gran dificultad de expresión teórica estricita, en el cuerpo teórico se introducen términos 5 Paralelización mediante GPU de la dinámica de electrones en plasma confinado magnéticamente de variabilidad aleatoria (denominados términos estocásticos), que permiten dotar al sistema de un cierto comportamiento y adaptación pseudoaleatoria necesario. Estas múltiples ejecuciones se denominan simulaciones, y permiten extraer nuevo conocimiento del sistema físico en estudio del análisis de sus resultados. Históricamente, estos códigos científicos han sido implementados por los propios científicos teóricos encargados de desarrollar las fórmulas básicas y sus derivaciones aproximadas, y presentan las siguientes características: • Lenguaje de programación: preferencia por lenguajes como Fortran o entornos como MatLab, sobre todo en comparación con C, con sintaxis más sintéticas en las instrucciones de cálculo matemático, especialmente al permitir trabajar con vectores y matrices sin necesidad de generar bucles que los recorran. • Estructuración: pensada para la implementación rápida de algoritmos matemáticos. Los módulos permiten definir bibliotecas de funciones, y la evitación de bucles explícitos para operaciones matriciales permite una mayor claridad del código, permitiendo al programador a centrarse en el algoritmo y no en el código. • Estructura: Los códigos cientíticos son principalmente de computación numérica. Por lo tanto suelen estar claramente estructurados en entrada o lectura de datos, preparación del cálculo, ejecución del cálculo (parte computacionalmente más intensa) y recopilación de resultados para mostrar o guardar. Por este hecho utilizan bibliotecas para la entrada y salida de datos (tal como silo, con formatos especiales como HDF5, y precisan de un lenguaje y bibliotecas muy completas y eficientes computacionalmente. Está fuertemente basado en biblitecas que resuelven los problemas computacionales más generales. • Independiencia entre cálculo y visualización: En relación con lo anterior, el código científico está pensado para ser revisado y reescrito, y por ello se hace independiente de procesos posteriores como el de la visualización. Los programas generan ficheros en formatos estandarizados para su posterior tratamiento por otros códigos, o por paquetes de visualización de datos, como VisIt, o el propio MatLab • Legibilidad: Pensados para ir actualizándolos conforme evolucione la teoría que genera el conjunto de ecuaciones a resolver. Por lo tanto, interesa que el lenguaje de programación evite artefactos propios del lenguaje que enmascare o dificulte la expresión y comprensión del algoritmo implementado. • Fuerte influencia de los algoritmos: lo cual puede ser un problema. Muchas veces, la transcripción directa de la fórmula a código, hace que se incurra en la utilización de algoritmos computacionalmente ineficientes cuando son implementados por un científico no experto y que no considere el coste computacional de su implementación. • Necesidad de tipos específicos de variables: ejemplo claro es la necesidad de trabajar con números complejos, y que la bibliotecas matemáticas y computacionales también lo hagan. También el reconocimiento de números especiales (i, e) y del correcto tratamiento de valores extremos (valores mínimos y máximos predefinidos, generación y trabajo con variables que hayan excedido el intervalo definido (NaN, Not a Number), para lo cual es muy interesante que se sigan el estándar IEEE 754 (Institute of Electrical and Electronics Engineers Standard for Floating-Point Arithmetic). 6 Capítulo 2. Conocimientos previos 2.3. GPGPU 2.3 GPGPU GPGPU [5], General-Purpose computing on Graphics Processing Units, consiste en la realización de cómputos de propósito general utilizando para ello procesadores gráficos, GPU, en lugar de los habituales procesadores genéricos, CPU. Debido a la economía de escala y al exitoso mercado de tarjetas gráficas impulsado por el sector de los juegos en el PC y consolas, los procesadores gráficos que se utilizan en esas tarjetas gráficas han evolucionado desde ser unos coprocesadores que simplemente representaban caracteres o gráficos de mapa de bits, a ser complejos procesadores capaces de realizar no tan simples cálculos matemáticos y vectoriales sobre cada píxel para poder dotar al conjunto de la imagen de un mucho mayor realismo, a la vez que proporcionar un alto refresco en la imagen para suavizar las transiciones entre cada instantánea [11]. De esta forma, el poder computacional de las GPU se ha multiplicado mucho más en estos últimos años que el de los procesadores de propósito general, CPU. Una tarjeta gráfica de alta gama puede tener perfectamente más de 2800 procesadores individuales todos ellos trabajando de forma (casi) independiente en paralelo, y con capacidad de comunicarse a través de rápidas memorias DDR5, incluso de varios GB de capacidad, ubicadas en la misma tarjeta [9]. No obstante, su conjunto de instrucciones es mucho más específico, orientado a la realización de cálculos gráficos, y se resiente de la falta de instrucciones optimizadas muy comunes en la programación estándar no gráfica (saltos y comparación, continua necesidad de trasvasar gran cantidad de información de la CPU a la GPU y viceversa, por ejemplo). Tampoco tienen la universal flexibilidad de las CPU pues tiene diferentes zonas de memoria, cada una con usos específicos y no intercambiables, o la coordinación entre procesadores dentro de un hilo de ejecución. Otra característica que penalizaba mucho las anteriores tarjetas gráficas era la inexistencia de instrucciones específicas para el flujo de control, toda vez que desde el punto de vista meramente gráfico, los programas siempre han sido preferentemente secuenciales. En la actualidad, los procesadores gráficos ya incorporan instrucciones para el control de flujo, pero las ramificaciones (if-then-else) son muy costosas, sobre todo si en cada hilo el resultado es diferente, lo que provoca graves penalizaciones en tiempo, para conseguir una adecuada sincronización en la ejecución del programa, por lo que conviene tenerlo muy en cuenta para obtener buenos resultados [6]. Estas restricciones se han suavizado con la llegada de las nuevas arquitecturas Kepler yMaxwell, pero aún así son aspectos de funcionamiento que conviene seguir teniendo en consideración para la optimización de la ejecución del código en las GPU. Aún así, a mediados de los 2000 empezaron los experimentos de realizar cómputos matemáticos pesados usando órdenes puramente gráficas, principalmente basadas en el estándar gráfico OpenGL. A pesar de la dificultad de programación ofrecida por esta vía, sí que se comprobó empíricamente la gran potencialidad que tenían las GPU para el cómputo numérico, y fruto de ello, NVIDIA lanzó en 2007 su lenguaje de programación propietario denominado CUDA C. Desde entonces prácticamente NVIDIA ha ido proporcionando la capacidad de ejecución de programas CUDA a la gran mayoría de sus tarjetas gráficas. Esta característica está especificada mediante la denominada “capacidad de computación CUDA” (en inglés, CUDA capabilities), que consiste en un número que permite conocer qué características y qué comportamientos tiene la tarjeta gráfica para la ejecución de código CUDA, sobre todo teniendo en cuenta que el lenguaje y las prestaciones ha evolucionado mucho desde 2007 hasta la fecha. Además, la numeración también tiene cierta correlación con las diferentes arquitecturas hardware que han aparecido: Tesla,Fermi,Kepler y la más reciente Maxwell. Posteriormente, Apple, liderando un conjunto de otras compañías (AMD,IBM,Intel yNVIDIA), impulsó el lanzamiento de otro lenguaje para GPGPU, denominado OpenCL [13], que al contrario que CUDA es de especificaciones abiertas. Actualmente el encargado de la especificación es el Grupo Khronos [8], quien también es el responsable de la especificación de OpenGL [14]. 7 Paralelización mediante GPU de la dinámica de electrones en plasma confinado magnéticamente Por ello, mientras CUDA únicamente está implementado en los procesadores gráficos de NVIDIA por ser un lenguaje propietario y totalmente cerrado, OpenCL puede utilizarse no sólo en los procesadores gráficos de diferentes fabricantes, como ATI [2], Intel [7] y la propia NVIDIA [12], sino que incluso hay CPUs que son capaces de ejecutar código de OpenCL (IBM Power 77X eIBM BladeCenters,Intel Ion y la tercera generación de Intel Cores [7] y los recientes procesadores de AMD, incluidos los AMD Fusion [1]). Samsung también dispone de una implementación de OpenCL que puede ejecutarse sobre procesadores Cell BE, ARM y sobre DSP [15]. A pesar del mayor número de fabricantes que proporcionan hardware compatible con OpenCL,CUDA es un lenguaje con una mayor evolución, y su compilador está mucho más optimizado para el hardware al que va destinado, dado que está específicamente diseñado para utilizar todo el potencial de los procesadores gráficos de NVIDIA.OpenCL es mucho más genérico, e incluso, como ya se ha comentado, puede ser ejecutado por CPU, por lo que tiene un mayor grado de abstracción sobre el hardware y por tanto, no permite emplear instrucciones o recursos específicos de ciertos procesadores gráficos. En la Tabla 2.1 se compara el desarrollo de cada uno de estos lenguajes. CUDA OpenCL Versión del Toolkit Fecha Especificación 1.0 julio 2007 1.1 noviembre 2007 junio 2008 1.0 2.0 agosto 2008 2.1 diciembre 2008 2.2 mayo 2009 2.3 julio 2009 3.0 marzo 2010 3.1 mayo 2010 junio 2010 1.1 3.2 noviembre 2010 4.0 mayo 2011 noviembre 2011 1.2 4.1 enero 2012 4.2 abril 2012 5.0 octubre 2012 5.5 julio 2013 noviembre 2013 2.0 6.0 abril 2014 6.5 agosto 2014 enero 2015 2.1 provisional 7.0 marzo 2015 7.5 septiembre 2015 Tabla 2.1: Comparación de la evolución cronológica de CUDA yOpenCL No obstante, es necesario recalcar que la computación a través de GPU puede en algún momento sobrevalorarse. Existen campos en los que las actuales CPU son tan capaces como las más potentes GPU, y es por ello que siempre es necesario valorar primero si el problema a abordar, o cada una de las partes en que se pueda dividir, debe serlo a través de computación genérica o de computación basada en GPU. Así, 8 Capítulo 2. Conocimientos previos 2.4. CUDA en campos tales como las matrices dispersas, los cálculos en que en cada procesador pueden divergir en cuanto a su flujo de control, o en los que existe una gran dependencia entre datos en diferentes ubicaciones espaciales, la computación genérica puede ser más eficiente, tal como se señala en [16]. 2.4 CUDA CUDA son las siglas de Compute Unified Device Architecture. Ha sido propuesta, implementada y desarrollada exclusivamente y de forma cerrada por NVIDIA. Se trata de una arquitectura hardware para procesadores gráficos que permite que éstos puedan realizar cómputos de carácter general (no gráficos) de una forma sencilla y directa, sin necesidad de invocar directamente las primitivas gráficas originales basadas en OpenGL oDirectX. Todo ello sin disminuir su propia capacidad intrínseca de ejecución de tareas gráficas basadas en cualquiera de estos estándares gráficos. 2.4.1 Lenguaje CUDA C Este concepto arquitectónico es aprovechado en la práctica mediante la definición de un nuevo lenguaje que es capaz de abstraer la complejidad del multiparalelismo de los actuales procesadores gráficos manteniendo una estructura y nomenclatura muy similar al lenguaje C, universalmente utilizado para la programación genérica usando CPUs. Este nuevo lenguaje se denomina CUDA C, y en lo sucesivo, en este trabajo, se abreviará simplemente por CUDA, el mismo nombre que la arquitectura. Se trata de un lenguaje de programación de especificaciones cerradas, que desarrolla en exclusiva NVIDIA, que permite la utilización de sus procesadores gráficos para computación genérica, no exclusivamente gráfica. Los programas realizados con CUDA solamente puede ejecutarse en aquellos procesadores gráficos que tengan «capacidad CUDA» (CUDA capability). Es importante diferenciar entre la versión del Toolkit, que representa el API de programación en CUDA con una sintaxis basada en el lenguaje C, y la capacidad de computación CUDA. Esta última se trata de una característica del hardware, esto es, de la GPU, y representa qué subconjunto de instrucciones es capaz de ejecutar directamente cada procesador gráfico. Existe una evidente interrelación, por cuanto en cada versión del Toolkit se menciona expresamente qué capacidad CUDA es necesaria para cada una de las funciones de librería que se proporciona, al objeto de utilizar optimamente el lenguaje en función del hardware que se vaya a utilizar. En [3] se listan todas las GPU que tienen capacidad de computación CUDA indicando expresamente su nivel de capacidad, que actualmente llega hasta la 5.2 para las GPU GeForce más avanzadas (serie GTX 970 y 980). 2.4.2 Tipos de código en CUDA C La programación en CUDA se caracteriza por contener dos tipos bien diferenciados de bloques de código: el código del anfitrión (host code) y el código del dispositivo (device code). • Anfitrión (host): Por anfitrión se entiende el entorno hardware de ejecución genérico proporcionado por las CPU convencionales. Su código es puro C/C++ y no precisa nada de los procesadores gráficos para su ejecución, pero obviamente, no se aprovecha el potencial de ellos. • Dispositivo (device): Por dispositivo se entiende el entorno hardware de ejecución propio de los procesadores gráficos, GPU. Su código, como ya se ha comentado, tiene una estructura y sintaxis muy similar al C, pero extendido para permitir programar fácilmente las nuevas capacidades que posee la arquitectura CUDA. 9 Paralelización mediante GPU de la dinámica de electrones en plasma confinado magnéticamente CUDA proporciona un compilador propio, denominado nvcc, que está especializado en la compilación de código de dispositivo, pero que también compila la parte del anfitrión y efectúa el enlace entre ambos códigos en un único programa ejecutable. El código del dispositivo se caracteriza por estar basado en kernels, que tienen la apariencia de funciones estándar C, pero con una nomenclatura de llamada especial, y que instruyen al procesador genérico de la CPU a enviar el código y los datos necesarios a la GPU, donde dicho kernel efectivamente se ejecuta. Los resultados de esos kernels pueden devolverse directamente a la CPU, para lo cual es necesario hacer una nueva transferencia de datos, en este caso entre la GPU y la CPU, o incluso mantenerse dentro de la memoria de la GPU para reutilizar en posteriores ejecuciones de nuevos kernels. Consiguientemente, la CPU siempre es la que lleva el control de ejecución del programa global, pero va delegando en la GPU aquellas partes del programa que permitan realmente obtener una gran ventaja de la paralelización masiva que ofrecen las GPU. 2.4.3 Compilación de aplicaciones CUDA La compilación de las aplicaciones CUDA se realiza fundamentalmente con el driver compilador nvcc.nvcc puede realizar muchas funciones, siguiendo caminos de compilación diferentes según las necesidades, tal como se muestra en la Figura 2.1. Una de las principales consiste en separar el código del anfitrión (host, que está en C) del dispositivo (device, para el que se usan extensiones específicas CUDA C). Una vez separados los códigos, lanza los compiladores específicos de cada tipo de código por separado para después ensamblarlos en un único ejectuable, generalmente. La parte del anfitrión se delega en un compilador de uso genérico, como pueden ser gcc oicc. La compilación del código de la tarjeta gráfica es mucho más crítica, y sólo puede realizarse con las herramientas específicas proporcionadas por NVIDIA, dado el carácter cerrado de la especificación hardware de sus tarjetas, a pesar de ser pública la especificación del código intermedio, formato denominado PTX (Parallel Thread eXecution), para el cual NVIDIA provee el compilador ptxas.ptxas puede ser invocado directamente por nvcc, pero también puede hacerse en tiempo real por la propia aplicación, realizando una compilación JIT (Just In Time) del código intermedio a código máquina específico de la arquitectura y capacidades CUDA de la tarjeta concreta en la que se va a ejectura el kernel. La compilación del código puede realizarse de tres formas básicas, hacia tres formatos diferentes, a saber: •.cubin: El formato cubin es microcódigo nativo específico de una tarjeta concreta NVIDIA con capacidad CUDA. Por lo tanto, es directamente ejecutable por ella. Los códigos binarios cubin en general no pueden intercambiarse entre diferentes arquitecturas o diferentes capacidades CUDA (CUDA capabilities. •.ptx: Es la representación intermedia del código a ejecutar en la GPU. Éste puede ser compilado a su vez y convertirse en código binario cubin, arriba descrito. El formato ptx es neutro, en el sentido de que es independiente de la arquitectura y de las capacidades computacionales de la tarjeta NVIDIA en la que finalmente se ejecutará. Está diseñado para ser resistente a cambios futuros («future-proof ») por innovaciones tecnológicas hardware y software. Puede ser compilado offline por el compilador driver nvcc, o compilado automáticamente online JIT según se necesite. •.fatbin: Es un formato especial que contiene en su seno diferentes compilaciones o precompilaciones, bien sea por contener tanto código cubin como código ptx, bien por tener código dirigido a diferentes arquitecturas o con diferentes capacidades CUDA, así como todo ello a la vez. 10 Capítulo 2. Conocimientos previos 2.4. CUDA SOFTWARE ARCHITECTURE 58 nvcc and PTX PTX (“Parallel Thread eXecution”) is the intermediate representation of compiled GPU code that can be compiled into native GPU microcode. It is the mechanism that enables CUDA applications to be “future-proof” against instruction set innovations by NVIDIA—as long as the PTX for a given CUDA kernel is available, the CUDA driver can translate it into microcode for whichever GPU the application happens to be running on (even if the GPU was not available when the code was written). PTX can be compiled into GPU microcode both “offline” and “online.” Offline compilation refers to building software that will be executed by some computer .cu file (Mixed CPU/GPU) nvcc Host-only Code GPU Code .ptx, .fatbin Host Executable (with embedded GPU code) CUDA Runtime (e.g., libcudart.so.4) Host Compiler .cu file (GPU only) nvcc GPU Code .ptx, .cubin CUDA Driver Host code (with embedded GPU code) CUDA Driver (e.g., libcuda.so) CPU Source Code Host Executable Offline components of application build process CUDA Runtime Driver API Figure 3.2 nvcc workflows. Figura 2.1: Diagrama de los flujos de compilación más habituales de nvcc Es evidente que los cubin son más rápidos en cargarse en la GPU, puesto que se envían directamente en cuanto debe ejecutarse dicho kernel. Por su parte, los ptx son más portables, pues pueden compilarse a cualquier arquitectura o capacidad igual o superior a su formato, pero incurren en la considerable sobrecarga de la compilación en línea justo en el momento en que debe lanzarse a la GPU. Los fatbin son ficheros más voluminosos, pero con mayor compatibilidad, según lo que se incluya, sobre todo si se incluyen ptx de diferentes arquitecturas y capacidades. Nótese (y esto es crucial para el presente trabajo) que la compilación específica del ptx se realiza justo en el momento en que se va a lanzar el kernel, no durante la carga en memoria de la aplicación al inicio de su ejecución. De igual forma que la selección del cubin concreto, en el caso de que ya estén incluidos varios y no sea necesario compilar. Es interesante indicar que los kernels CUDA pueden estar incrustados en el propio fichero ejecutable ELF (o EXE, en Windows), o también pueden cargarse durante la ejecución desde otros ficheros, en cualquiera de los tres formatos arriba indicados. De hecho, cuando el kernel está incluido en el programa principal, está en formato cadena (string literals), y lo que se hace es leerlos y pasarlos a la GPU mediante unas instrucciones específicas previas, que se encargan también de lanzar su ejecución dentro de la tarjeta gráfica. 11 Paralelización mediante GPU de la dinámica de electrones en plasma confinado magnéticamente Por último, la correcta selección del tipo de tarjeta gráfica para el cual se quiere compilar (directamente a código binario .cubin u obtener código intermedio PTX) es esencial. Esta selección se verá de vital importancia a la hora de la migración, dado que en todo caso, la aplicación que se lanza deberá, como requisito básico, poder ejecutarse tanto en la GPU inicial como en la GPU a la que se migra. 12 CA P Í T U L O 3 ECUACIONES DE LANGEVIN PARA DESCRIBIR UN PLASMA DE PARTÍCULAS Y LA GENERACIÓN DE ELECTRONES runaway En este capítulo resumen brevemente la base física nuleótica que soporta el conjunto de ecuaciones que son resueltas numérica a través del código original de Langevin realizado en Fortran, y que permiten la simulación del comportamiento de los electrones de un plasma a alta temperatura confinado magnéticamente, y eventualmente, bajo la acción de un campo eléctrico perpendicular. Primeramente se introducirá el concepto de plasma, y su importancia económica por el impacto que podría llegar a suponer la obtención de una reacción de fusión nuclear con balance positivo de energía, que resultaría en una fuente relativamente limpia de energía a partir de elementos muy comunes en nuestro entorno, tales como el Helio, el Deuterio, el Tritio y el Litio. Posteriormente, se repesará muy por encima, la base teórica de la física de partículas que explica el comportamiento de los plasmas, haciendo especial mención a una forma de abordarlas denominada ecuaciones de Langevin. Este método de Langevin permite, mediante unas aproximaciones razonables, obtener unas ecuaciones relativamente sencillas, sobre las que es posible aplicar métodos numéricos y con ello conseguir la simulación del comportamiento del plasma. Este trabajo se centrará específicamente en la simulación de un tipo de partículas muy iteresantes, los electrones fugitivos ( runaway). 3.1 El plasma y su aplicación práctica: reacción controlada de fusión El plasma es un estado de la materia en el cual ésta está sometida a altas temperaturas (generalmente millones de grados Kelvin), a resultas de lo cual, sus partículas están totalmente ionizadas. En este estado su comportamiento es similar al de un gas, por la gran libertad de movimiento de sus partículas. Esa gran temperatura y gran libertad, implica una alta energía cinética, esto es, unas altas velocidades dentro del medio. Las condiciones necesarias para la generación y mantenimiento de la materia en estado de plasma se dan comúnmente en la naturaleza, pero no en nuestro inmediato entorno, sino en el interior de las estrellas, por ejemplo. 13 Paralelización mediante GPU de la dinámica de electrones en plasma confinado magnéticamente Para su utilización en C(básicamente, dentro de las funciones wrapper), se crean unos ficheros de cabecera (.h), que mediante directivas de compilador #define, se realiza una traducción del nombre complejo en que aparecen estas variables globales en los ficheros objeto compilados con Fortran. Para ser genérico, conviene incluir directivas compiladores de compilador que sean capaces de detectar los posibles compiladores que se van a usar, y según eso, utilizar la traducción adecuada. A continuación se expone un ejemplo de cómo difernciar entre compiladores para direccionar correctamente el símbolo de una variable global exportado desde código Fortran: Listado 4.1: Compilación condicional mediante detección del compilador 1#if defined(__ICC) || defined(__INTEL_COMPILER) 2/*Intel ICC/ICPC. ------------------------------------------ */ 3#define L_TWO_NOISES __langevin_params_mp_l_two_noises 4#define L_DE_LOR __langevin_params_mp_l_de_lor 5... 6#elif (defined(__GNUC__) || defined(__GNUG__)) && !(defined(__clang__) || defined( __INTEL_COMPILER)) 7/*GNU GCC/G++. --------------------------------------------- */ 8#define L_TWO_NOISES __langevin_params_MOD_l_two_noises 9#define L_DE_LOR __langevin_params_MOD_l_de_lor 10 ... 11 #endif En http://nadeausoftware.com/articles/2012/10/c_c_tip_how_detect_compiler_name_ and_version_using_compiler_predefined_macros se muestra cómo detectar más familias de compilador, y con ello poder ampliar este ejemplo a una mayor variedad de compiladores de C. Diferencias en los prototipos de función entre Fortran yC A pesar de que Fortran no está totalmente estandarizado, especialmente en lo referente al name mangling, los compiladores Fortran gfortran de GNU eifort de Intel comparten las siguientes características particulares en la creación de los símbolos de las subrutinas y funciones: • Las funciones y subrutinas en el fichero objeto están siempre con minúsculas • Además, se les añade a su nombre el carácter de subrayado (también conocido por “guión bajo”) 0_0. Resumen comparativo de las diferentes nomenclaturas según lenguaje y familia de compilador Como referencia y resumen final, en la Tabla 4.1 se ejemplifican algunos casos: Tipo de símbolo Tipo Nombre declarado Símbolo creado por Símbolo creado por nm en Fortran gfortran (GNU)ifort (Intel) Función propia T ADVANCE advance_advance_ Función externa (indefinida) U RHS_COEFFS_rhs_coeffs_rhs_coeffs_ Variable propia inicializada D MYDT __advance_MOD_mydt __advance_mp_mydt Variable propia no inicializada B MY_ITER __advance_MOD_my_iter __advance_mp_my_iter Variable externa (indefinida) U DT __orb_params_MOD_dt __orb_params_mp_dt Tabla 4.1: Ejemplos del name mangling de los compiladores Fortran usados 20 Capítulo 4. Implementación de la paralelización con CUDA Sección 4.2 En la tabla anterior, se muestran los símbolos que se crean en el fichero objeto producto de la compilación de un código fuente de un supuesto fichero ADVANCE.F que define la subrutina ADVANCE(). Se supone que esta subrutina realiza una llamada a la función RHS_COEFFS que está declarada y definida en un fichero externo. En la subrutina ADVANCE se definen dos variables, MYDT yMY_ITER, las cuales son globales para todo el programa Fortran, una de ellas inicializada y otra no, y utiliza la variable externa DT declarada en el fichero ORB_PARAMS. Más información sobre el significado y tipos declarados por la utilidad GNU nm se puede obtener en http://www.freebsd.org/cgi/man.cgi?query=nm&apropos=0&sektion=0&manpath=FreeBSD+9. 0-RELEASE&arch=default&format=html 4.2.2 Esquema general de llamadas a CUDA desde Fortran a través de wrappers C El método elegido para adaptar el programa en Fortran 90 a CUDA consiste en la creación de una función wrapper en C, invocada de forma normal desde Fortran, que recoge los parámetros y accede a las variables globales declaradas en Fortran necesarias para pasar al kernel. A partir de ese momento, todo el código puede hacerse en Co en C++, dado que se compilará con el compilador nvcc proporcionado por CUDA. Se deberán tener en cuenta las cuestiones de nomenclatura y de paso de parámetros ya expuestas en el anterior apartado. En general el proceso será: • La función wrapper C, deberá definirse siempre con los calificadores extern "C" void, dado que una subrutina en Fortran nunca devuelve valor. El extern "C" permite la compatibilidad con la nomeclatura de funciones de Cdesde C++, puesto que nvccen realidad es un envoltorio sofisticado del compilador C++ presente. • En Cla función deberá terminar siempre con el carácter de subrayado adicional (también conocido como guión bajo) ’_’, y con todos sus caracteres en minúsculas. Todos sus parámetros serán punteros, y es recomendable que aquellos que no deban ser modificados por esa función o sus derivadas (parámetro sólo de entrada), lleven el calificador const. • Desde Fortran se invocará a la función por el nombre definido en C, con cualquier combinación de mayúsculas y minúsculas, pero obligatoriamente sin el último carácter de subrayado. Además, debe invocarse como se hace a una subrutina de Fortran (con CALL <nombre_subrutina(lista_parámetros)), y no como se hace a una función de C. • En la función wrapper Cse gestionarán los parámetros de entrada, y sobre todo, se recomienda allí utilizar macros para renombrar las variables globales de Fortran, por ser tremendamente farragosas de escribir tal cual desde C. A partir de allí, con esos sinónimos en Cse operará normalmente, siempre teniendo en cuenta que se trata de punteros. Extremar la precaución con los parámetros que sean matrices multidimensionales por la cantidad de indirecciones que serán necesarias para su correcta referencia al valor, o a la dirección (puntero donde se encuentre el valor real). • En el wrapper Cse gestionarán las transferencias de datos que sean necesarias entre la CPU y la GPU. Se recomienda, por eficiencia, evitar en todo lo posible la creación de nuevas zonas de memoria para cambiar el orden de acceso de las matrices. • Desde el wrapper Cse invocará ya de forma normal a los kernels, siguiendo las particularidades del lenguaje CUDA-C, en especial, de la invocación a dichos kernels. 21 Paralelización mediante GPU de la dinámica de electrones en plasma confinado magnéticamente • Dentro de los kernels, donde se utiliza el lenguaje CUDA-C, se deben tener en cuenta, aparte de todo lo característico de ese lenguaje extensión de C, la procedencia de las zonas de memoria multidimensionales (arrays). Si son en Cdeben accederse de la forma habitual (orden por filas), pero si provienen (fueron definidas) en Fortran, entonces el acceso debe realizarse teniendo en cuenta que el almacenamiento es por columnas, esto es, el segundo elemento en la zona de memoria no es el elemento A(1,2), sino el elemento A(2,1). Es por ello que se recomienda el uso de macros en Cpara direccionar siempre con ellas las zonas de memoria definidas en Fortran. A continuación se muestran las macros utilizadas en este proyecto para realizar tal acceso, en función de si se pretende acceder al valor, o a la dirección de memoria donde está almacenado dicho valor. Estas macros deberán definirse para cada tipo de array multidimensional. Las que se presentan son para matrices bidimensionales de números flotantes de doble precisión (**double odouble[][]). Listado 4.2: Macros en C para el acceso a matrices 2D definidas en Fortran 1// Convert linear vector representation to 2D matrix in C-like file/column order 2// Use: MATRIX2D_DBL_PTR(pointer_to_matrix_start, file, column, file_width) 3#define MATRIX2D_DBL_PTR(m,i,j,N) (((m))+((j)*(N)+(i))) // A pointer to value 4#define MATRIX2D_DBL_VAL(m,i,j,N) (*(((m))+((j)*(N)+(i)))) // The value itsef Siguiendo estas indicaciones, supongamos una función denominada ADVANCE, que está declarada en el código original en Fortran, y que al realizar gran computación numérica la vamos a sustituir por un kernel. Sea PARAM1 sea un parámetro de entrada, y PARAM2 otro de salida. Su definición en Fortran será: Listado 4.3: Definición de una subrutina común en Fortran 90 1SUBROUTINE ADVANCE(PARAM1, PARAM2) 2INTEGER,INTENT(IN) :: PARAM1 3REAL*8, INTENT(OUT):: PARAM2 4... 5END SUBROUTINE ADVANCE Para sustituir esa función por un wrapper en C, su declaración y contenido básico será la siguiente: Listado 4.4: Definición de la función wrapper en C 1extern "C" void kernel_advance_wrapper_(const int *param1, double *param2) { 2... 3// Copy data from CPU to GPU, if apply 4kernel_advance<<<GridSize,BlockSize(param1, ...); 5// Copy data from GPU to CPU, if apply 6... 7return;// It’s void. So we have to return nothing! 8} Finalmente, esta función wrapper en Cse invoca desde Fortran de la siguiente forma: Listado 4.5: Ejemplo de llamada a un wrapper Cdesde Fortran 1CALL KERNEL_ADVANCE_WRAPPER(PARAM1, PARAM2) 22 Capítulo 4. Implementación de la paralelización con CUDA Sección 4.3 4.3 Diferencias programáticas importantes entre Fortran yCyCUDA-C En las secciones anteriores ya se han comentado aspectos diferenciadores entre CyFortran, pero que tenían más que ver con la sintaxis. En esta sección se hará especial hincapié en las principales diferencias semánticas y programáticas que deben considerarse en la paralelización con CUDA-C de programas Fortran, que son las siguientes: Paso de parámetros: Por valor frente a por referencia En Cel paso de parámetros es por valor siempre (para pasar por referencia, en realidad se hace un pase por valor, pero del puntero a la zona de memoria referenciada), en Fortran el paso es siempre por referencia. En un apartado posterior se analizarán las consecuencias de este hecho. De esta forma, en Fortran no se puede incluir una constante como parámetro, y por otra, desde el wrapper en C, al recibir un parámetro siempre será un puntero, y como tal, estará expuesto a ser modificado o eliminado desde C. Es por ello que se recomienda que sea declarado como puntero constante (const) en el prototipo en C. Variables globales Fortran dentro de una subrutina llamada por diferentes hilos OpenMP Las subrutinas Fortran, cuando son llamadas desde diferentes hilos de ejecución, no comparten las variables globales definidas en Fortran. Están disponibles, pero cada hilo tiene su propia versión local, sin que se copie el valor global cuando ésta tenía al ser llamada. Para solucionar este problema, se propone el uso de parámetros adicionales en el prototipo de la subrutina Fortran. Así, por ejemplo, si se precisa tener acceso a la variable global thread del programa Fortran, se debe añadir un nuevo parámetro, con el modificador INTENT adecuado ("IN", o "INOUT"), para que en la subrutina se utilice esta nueva variable local con el valor que tenía su correspondiente variable global en el programa principal en el momento en que fue llamada. De lo contrario se inicializa a cero. Este problema afecta también a las matrices y variables definidas como ALLOCATABLE. Diferente orden de almacenamiento en memoria de matrices multidimensionales Ya ha sido comentado en el apartado anterior, pero es tan importante que merece un nuevo recordatorio. En Cel acceso a las zonas de memoria multidimensionales se define y realiza por filas primero, por columnas después (si es multidimensional, de izquierda a derecha en el orden de las dimensiones tal cual se escriben). En Fortran, para matrices 2D, es por columnas primero y después por filas. Si son multidimensionales, de derecha a izquierda en el orden que aparecen las dimensiones o índices. Diferente inicio del índice de vectores y matrices En C, el índice de inicio siempre es 0, tanto en filas como en columnas. No es lo mismo en Fortran, en donde el índice del primer elemento normalmente empieza por uno (véase la excepción en la siguiente característica). Eso es muy importante a la hora de reutilizar el código en Fortran como base para su traducción a C: los bucles e índices deben revisarse convenientemente para que accedan correctamente al elemento correspondiente en cada momento. Matrices multidimensionales con índices que no empiezan en cero En Fortran es posible definir matrices en las que una (o varias) dimensiones no empiecen en 1 como normalmente se hace. Puede definirse que el índice de una determinada dimensión tenga un intervalo cualquiera, como por ejemplo, entre -3 y +6, en lugar de 1 a 10. En Clos índices siempre serán, en ese ejemplo, entre 0 y 9, por lo que a la hora de migrar el código de Fortran aCdeberán contem23 Paralelización mediante GPU de la dinámica de electrones en plasma confinado magnéticamente plarse estas eventualidades para que el acceso sea siempre al elemento correcto (aparte de que se podría, además, estar accediendo a zona de memoria fuera de la matriz, con riesgo de corrupción de la misma). Compatibilidad entre los diferentes tipos numéricos Es muy importante tener en cuenta realmente cuál es el espacio justo de almacenamiento de las variables en Csegún su declaración en Fortran, donde se pueden definir, por ejemplo, números reales de X decimales e Y cifras enteras. Esto da lugar a tipos de variables que deben analizarse para asegurar que en Cse le asigna un tipo compatible en signo, en capacidad (4, 8, 16 Bytes) y en precisión. Especial atención a los tipos inexistentes en C, como los LOGICAL de Fortran. Otras características menores No merece la pena entrar en detalle de la multitud de pequeños detalles de divergencia entre estos lenguajes. Como mero ejemplo, en este proyecto en algún momento se pensó en la conveniencia de utilizar un enumerate de C, pero se observó que en Fortran no existen las denominadas “constantes nombradas”. Por lo tanto, se tiene que trabajar con variables normales directamente, con valores enteros. Por ello se precisa un experto conocedor de ambos lenguajes (Fortran,CyCUDA-C) para realizar una correcta y óptima migración o adaptación del código histórico en Fortran aCUDA. Como corolario, es tremendamente importante tener presente en todo momento de dónde proceden las matrices, puesto que según hayan sido definidas en Co en Fortran, su acceso y tipología son diferentes. Además, es muy importante recordar que según en qué lenguaje estemos programando, el acceso a las matrices y su definición, son diferentes. 4.4 Implementación en CUDA del código de Langevin originalmente en Fortran En esta sección se van a comentar brevemente las principales características de la implementación realizada resultante de la migración parcial del código de Langevin en Fortran aCUDA. En una sección posterior se comentará la modificación que a ésta primera implementación en CUDA debe realizarse para obtener una mejora de rendimiento por el desacoplamiento entre computación e intercambio de mensajes y escritura en disco. 4.4.1 Aplanamiento de subrutinas por inclusión en un kernel aglutinador En CUDA la invocación a los kernels, y los cambios de contextos internos entre los diferentes hilos de la GPU son menos costosos que en los lenguajes tradicionales de CyFortran. Por otro lado, es muy importante el número de registros y otros recursos, como las unidades lógicas de coma flotante (simple y doble, pues son independientes) no se saturen, porque eso disminuye la cantidad de hilos que pueden estar en ejecución simultáneamente en el mismo núcleo interno (bien sea en estado activo, con actualización constante del contador de programa, bien en estado pasivo, por estar a la espera de la finalización de alguna operación larga, como las aritméticas de doble precisión). Consiguientemente, es muy factible en CUDA obtener buenos resultados incluso aunque se utilicen muchos kernels simples, siempre que éstos precisen pocos recursos, pero el número total de recursos queden bien saturados, de forma lo más homogénea posible, por la invocación y ejecución en cada núcleo del 24 Capítulo 4. Implementación de la paralelización con CUDA Sección 4.4 número de bloques óptimo (o un poco más, nunca menos), como se demostrará en el siguiente Capítulo 5 “Valoración cualitativa del proceso de paralelización con CUDA”. A pesar de ello, y como demuestra este trabajo, es posible también contar con un buen rendimiento ocupacional aunque el kernel sea más bien complejo, si la utilización de los recursos es homogéneo. Tal es el caso que se produce en el kernel_advance, que es el principal responsable del tiempo de ejecución global del programa. Este kernel se encarga de, para cada partícula, calcular la nueva velocidad vectorial en función de la velocidad anterior, del valor de los campos magnético y eléctrico, y de un componente estocástico (fundamental al ser unas ecuaciones del tipo Langevin) que simula el efecto colisional del resto de partículas (electrones e iones) a su alrededor. En el código original, tal como se observa en la Figura 5.1, la subrutina ADVANCE es responsable del 5,45% del tiempo de computación, pero ésta llama exclusivamente a otras subrutinas de ayuda. Estas subrutinas de ayuda sirven para calcular algunos parámetros concretos de la fórmula de Langevin, que es la que se computa en la subrutina padre ADVANCE, tales como el efecto de la colisión con iones (RHS_ION) y el ángulo de la velocidad resultante tras impactar con iones o electrones, (PHI_AND_PSI, llamada a través de la otra subrutina RHS_COEFFS). También se utilizan subrutinas auxiliares, como GET_GAUSSIAN_NOISE, encargada de realimentar el plasma con un electrón estadísticamente similar a la condición inicial por cada electrón fugitivo (runaway) que se escapa del confinamiento magnético. Esta subrutina, GET_GAUSSIAN_NOISE, a su vez, llama a una subrutina de biblioteca, RAN2, que es la que genera los números aleatorios necesarios para el comportamiento estocástico de la ecuación de Langevin. Por lo tanto, en realidad, la subrutina ADVANCE, y las dependientes de ella, es responsable del 98,5% del tiempo de ejecución del código original. Esta cadena de llamadas de subrutinas no puede realizarse de forma equivalente con kernels. No en el sentido de que un kernel llame a otro kernel. Aunque eso es posible a partir de las tarjetas Kepler, esas invocaciones de un kernel por otro, está permitido sólo con un nivel de anidamiento de uno, e implica una nueva creación de muchos más bloques con sus respectivos hilos. Está pensado para permitir, hasta cierto punto, la programación dinámica, un paradigma de programación que no es aplicable a nuestro programa. Lo que se hace, es invocar a un único kernel, el cual puede llamar a otras funciones en CUDA-C que lo que hacen es simplificar la expansión del código completo en una única función kernel. Pero en realidad, el kernel es único, lo que pasa es que éste llama a funciones, compiladas en el código propio de la GPU (PTX). Por lo tanto, en el código presentado, el kernel_advance_wrapper_incluye otras funciones CUDA-C, rhs_coeffs yrhs_ion, que son la migración de sus correspondientes subrutinas de Fortran. Se observa que “desaparecen” las subrutinas GET_GAUSSIAN_NOISE yRAN2. Eso es debido a que las tarjetas NVIDIA, a través de CUDA tienen sus propios (y variados) generadores de números aleatorios, muy probablemente más eficaces y seguros que cualquier otro método que se quisiera implementar directamente con CUDA-C en la GPU. 4.4.2 Minimización de la transferencia de datos entre CPU yGPU Las implementaciones presentadas han minimizado las transferencias de datos en ambos sentidos entre la CPU y la GPU. Únicamente se transmiten los siguientes, que son inevitables: • Vectores de velocidades iniciales cuando éstos son leídos de un fichero de reinicio (ficheros con el nombre langevin.restart.????): Evidentemente, pues la GPU no puede acceder directamente al disco duro. Nótese que si el programa debe generar las velocidades iniciales, según los parámetros indicados en el fichero de configuración (distribución uniforme o gaussiana), éstas se generan ya 25 Paralelización mediante GPU de la dinámica de electrones en plasma confinado magnéticamente directamente en la GPU, sin intervención de la CPU, quien desconcerá hasta la finalización de la ejecución el valor concreto de la velocidad de cada partícula del plasma. • Vectores de velocidades finales: Es necesario pasarlos de la GPU a la CPU al final de la computación, por la misma razón: la GPU no puede guardar esos datos en ficheros langevin.restart.???? en disco duro. • Valores estadísticos intermedios: La propia GPU es la encargada de catalogar la velocidad de cada partícula en su correspondiente casillero para finalmente enviar a la CPU sólo los datos estadísticos ya elaborados que periódicamente, si así se ha indicado, deben grabarse en disco o mostrarse por consola. Estas estadísticas pueden ser de diferentes valores y características, pero siempre son computadas por la GPU y transmitidos los resultados estadísticos elaborados a la CPU. • La contabilización del número de electrones fugitivos generados: Que en realidad es redundante con lo anterior, puesto que es otra estadística más. Con ello se consigue eliminar la necesidad de intercambiarse la información sobre las velocidades de las partículas entre la GPU y la CPU. La versión que definitivamente eliminó esa necesidad proporcionó una aceleración interna entre versiones de prácticamente el 100%. Esto es, el intercambio en cada iteración temporal de los vectores de velocidad suponía un coste semejante al coste computacional de los kernels. 4.4.3 Utilización adecuada de los generadores de números aleatorios de CUDA CUDA, mediante la biblioteca estándar CURAND, proporciona un muy variado y eficiente conjunto de generadores de números aleatorios. Tal es su variedad, que se pueden seleccionar diferentes métodos de obtención de números, que se basan en diferentes métodos, con ventajas y desventajas respecto de la agrupación, seguridad, distribución y velocidad de generación de los números aleatorios. En este trabajo no se ha entrado a valorar las diferentes técnicas base, y se han seleccionado los que parecían más usualmente usados en la bibliografía, en tutoriales y otros programas. Además, las funciones que aporta, permiten la generación directamente de números con distribución uniforme (la más usualmente implementada de forma estándar en los demás lenguajes), pero también de números con distribución normal (también conocida como distribución gaussiana). Es más, incluso es posible generar los números aleatorios por pares, con una mejora sobre un 25% respecto a generar dos números aleatorios de forma individual. Ambas características han sido utilizadas en estas implementaciones, permitiendo que fuera la propia GPU la encargada de sintetizar los términos estocásticos directamente dentro del kernel. Otro aspecto de rendimiento muy importante a tener en cuenta es que los números aleatorios deben ser inicializados. Esta inicialización no debe realizarse en cada invocación del kernel, sino que si las variables de los generadores aleatorios son guardados en la zona de almacenamiento global de la GPU, se asegura su permanencia durante toda la ejecución de la aplicación, aunque el kernel concreto que los haya inicializado o usado por última vez finalice. De esta forma, es posible tener los generadores ya inicializados para cada partícula, e invocar la petición de nuevos números aleatorios aprovechando las semillas existentes de una invocación a otra del kernel. Cuando este fenónemo se descubrió, se obtuvo una aceleración relativa del 35% en el tiempo global de ejecución del programa global, entre las dos versiones implicadas. Esto significa que aunque los números aleatorios pueden ser rápida y eficientemente generados en la GPU, no es así su proceso de inicialización. Es por ello que es conveniente inicializarlos una única vez durante toda la ejecución de la aplicación, y reusarlos siempre que sea posible. 26 Capítulo 4. Implementación de la paralelización con CUDA Sección 4.4 Otra mejora fue la de asegurarse que estos generadores están en memoria compartida en el momento de la ejecución del kernel. Esto es, en general están en memoria global, pero cuando son necesarios se invocan varias veces. Para asegurarnos de que únicamente se bajan y suben de la memoria global a la más cercana a los núcleos CUDA una vez (al inicio y al final de cada invocación al kernel) y no depender de que estén o no todavía presentes en la cache, se hace una copia de ellas en memoria shared y después, al terminar el kernel, se sube el estado del generador aleatorio de nuevo a la memoria global. Esto es posible porque tenemos un generador por cada una de las partículas. No obstante, hay que señalar un problema que se produjo con los números aleatorios: en determinadas simulaciones, es necesario contar con números aleatorios distribuidos uniformemente para unas cosas, pero otros distribuidos de forma normal (gaussiana) para otras. En esos casos, no puede utilizarse la misma variable generadora, siendo necesario crear generadores diferentes para distribuciones diferentes. Esto implica un sobrecoste, pero se realiza una única vez durante el programa. Como idea de optimización se deja la posibilidad de probar diferentes generadores, para ver si su inicialización y constante generación es más rápida en unos que en otros. Nótese que nuestro programa, en principio, no requiere de especiales condiciones de extrema seguridad en cuanto a la calidad de los números aleatorios generados. Otra idea de optimización consiste en seguir los consejos que en diversos foros aparecen, consistentes en que se inicializa sólo una variable generadora de números aleatorios, o un número reducido de ellos, por ejemplo, una por cada bloque. Posteriormente, cuando se solicitan números aleatorios, en lugar de pedir el siguiente, se pide el segundo, tercero o cuarto, etc. según el número de hilo. Dicen que si se hace avanzar la cuenta, es posible que no existan colisiones de números generados. No obstante, hay críticas a estos métodos por cuanto que su utilización implica la pérdida de la seguridad teórica de los algoritmos base de la implementación de estos números aleatorios. Por lo tanto, los críticos indican que los números obtenidos pueden estar sesgados, o tener ciertas agrupaciones locales indeseadas. Es un tema abierto a la experimentación, tanto para ver si acelera realmente las actuales implemenaciones, y también sobre si los resultados finales se muestran estadísticamente semejantes o diferentes a los obtenidos con la actual implementación, más conservadora y posiblemente más lenta, pero segura. 4.4.4 Mejoras en el código: utilización de funciones estándar frente a métodos numéricos pesados En el código original, la subrutina advance debe implementar una determinada integral. En el código original, la implementación es tal cual, esto es, se realiza una integración numérica usando el método trapezoidal con unos 200 pasos. Esto computacionalmente es costoso. Si se analiza la integral, y como en algunos artíclos puramente científicos donde se desarrollan estas ecuaciones de Langevin para plasma, esa integral ha dado lugar a la definición como tal de una función ampliamente conocida: la función error. La sorpresa está en que esta función es tan frecuente, que prácticamente todos los lenguajes computacionales (esto es, que mínimamente sirvan o se usen para cálculos matemáticos), tienen prácticamente de serie, alguna función propia o en una biblioteca numérica estándar, implementada esta función. Tal como incluso se ha comprobado experimentalmente en este proyecto, la llamada a esta función es mucho más precisa, exacta y rápida que el método de integración implementado en Fortran dentro del código original. Es más, por ejemplo, en C, el error de la función está asegurada en ser como máximo dos unidades en la última posición representativa del número flotante de doble precisión resultante, siendo de una como máximo en la mayor parte del intervalo razonable y normal de uso. Es evidente que doscientas sumas y multiplicaciones para un cálculo numérico aproximado va a tener una exactitud de varios órdenes mucho peor. 27 Paralelización mediante GPU de la dinámica de electrones en plasma confinado magnéticamente La función error en Cy en CUDA-C tiene el mismo nombre, y se trata de erf(), la cual se ha utilizado en el kernel_advance, con mejora de prestaciones. Esta mejora ha sido del 1–2% en el momento de realizar el cambio. 4.4.5 Algoritmos de reducción en la GPU El secreto evidente de la buena aceleración que se obtiene en los resultados de las implementaciones presentadas en este trabajo, es la absoluta independencia de los cálculos entre cada una de las partículas simuladas. Por ello, no existe dependencia de datos en la GPU y prácticamente se asegura que todos los accesos a memoria son coalescentes, propiedad que mejora las prestaciones computacionales de la tarjeta gráfica, al mejorar el trasiego de información interno entre la memoria y los núcleos de la GPU. No obstante, existe un par de sitios en el algoritmo, en donde la GPU se ve forzada a realizar una reducción. Esto es, debe obtener un resultado a partir de los valores calculados por el resto de hilos y bloques, esto es, de todas las partículas: • Cálculos estadísticos para obtener estadísticas sobre el valor de las componentes de la velocidad de los electrones (perpendicular y paralela al campo magnético): caracterizado porque se deben siempre tener en consideración el valor de todas las partículas (por lo tanto, de todos los hilos participantes en el kernel). • Recuento del número de electrones que en cada iteración pasan a ser electrones fugitivos (runaway), por exceder la velocidad crítica de confinamiento: se caracteriza porque, en condiciones normales del plasma, se trata de una baja proporción de electrones (hilos) los que deben computarse. Estas reducciones se pueden abordar al menos desde dos perspectivas diferentes: • Mediante algoritmos propios de reducción, generalmente con coste computacional de orden O(log2(N)): Son los indicados para cuando el porcentaje de hilos (partículas) a participar en la reducción es elevada respecto al total. • Mediante operaciones atómicas que cada partícula realiza sobre una variable global común: Son los indicados cuando el porcentaje de hilos (partículas) a participar es despreciable frente al conjunto total de partículas. La elección del método es muy importante por su impacto en el tiempo de computación. Esta decisión es un claro ejemplo de las dificultades que suelen tener los programadores que no son expertos en ingeniería informática, dado que frecuentemente se obvia la realización de un estudio previo (o aunque sea posterior, experimental) de costes del algoritmo a implementar. La elección para cada uno de los dos casos es diferente, y se basa en un estudio de costes: Estadística de velocidades Cada hilo se ubica en su contenedor según sus componentes de velocidad perpendicular y paralela, pero hay que ver cuántas partículas hay en cada contenedor. Por lo tanto es una reducción en las que siempre participan todas las partículas. Es por ello que se requiere un algoritmo de reducción basado en reducción binaria en árbol, que tiene un coste de orden O(log2(N))- Conteo de electrones runaway En condiciones normales de plasma, el número de electrones fugitivos es muy pequeño. Si fuera grande, el tokamak estaría descontrolado y la reacción nuclear de fusión se pararía en pocos milisegundos. Por lo tanto, no tiene sentido estar simulando en esas condiciones. Consiguientemente, interesa el algoritmo basado en operaciones atómicas sobre una misma variable, 28 Capítulo 4. Implementación de la paralelización con CUDA Sección 4.4 pues en el caso extremo de que coincidan todas las operaciones de adición en el mismo momento, éstas se serializarían, y el coste sería de ordenlineal, O(log2(n)), con nsiendo el número de colisiones. Analicemos cada una de este tipo de reducción por separado. 4.4.5.1 Reducción global iterativa Esta reducción tiene sus triquiñuelas cuando se tiene que realizar en GPU, sobre todo cuando la capacidad computacional CUDA de la tarjeta usada es más baja. Esto es así, porque éstas manejan mucho peor los accesos no coalescentes a la memoria. Para resolver esta eventualidad, es importante planificar estos accesos para que sean lo más coalescentes posibles, y siempre bajo el prisma de que será necesario realizar una entrada recursiva (en realidad, iterativa) de orden O(log2(N)), siendo Nel número de partículas, para su final reducción. En el SDK de CUDA aparece un método terriblemente eficiente que utiliza unas macros de compilador para obtener la máxima potencia de la tarjeta gráfica, independientemente de su capacidad computacional, al optimizar mucho el acceso coalescente a la memoria. No obstante, es un método muy rebuscado como para poder reutilizarlo. Es por ello la implementación se ha realizado utilizando la técnica del árbol binario, con coste O(log2(m)), siendo mel número de hilos en cada bloque (32, o sea, 6 iteraciones), para posteriormente, hacer otra pasada y volver a totalizar esas sumas parciales (total de bloques dividido por 32), una por cada bloque lanzado. Como ejemplo, se puede acudir al código fuente de la implementación, al fichero kernel_calculate_statistics.cu, concretamente a los kernels kernel_calculate_statistics_ratios ykernel_calculate_statistics_ratios2D, que son demasiado extensos para su inclusión en esta memoria. 4.4.5.2 Reducción mediante el uso de operaciones atómicas CUDA provee de operaciones atómicas sobre variables globales. El único problema con las operaciones atómicas en la GPU es que hasta el momento (hasta CUDA versión 7.0, capacidades computacionales hasta 5.2), no está permitido realizar una reducción sobre variables de longitud de 64 bits. Esto es, no se puede hacer un atomicAdd sobre un entero largo (long int) o sobre un número en coma flotante de doble precisión (double). Como en nuestro caso, el conteo se realiza con enteros normales, de 4 octetos, es perfectamente viable, y recomendable desde el punto de vista de teoría de coste computacional. Simplemente consiste en seleccionar sólo aquellos hilos (partículas) que deban contabilizarse, y con ello, invocar directamente a la función atomicAdd. Veamos como ejemplo, la forma en que se utiliza en el kernel que recuenta los electrones con velocidad mayor que la crítica, considerados por ello en la simulación como electrones runaway. El Código 4.6 es un ejemplo de su utilización. 29 Paralelización mediante GPU de la dinámica de electrones en plasma confinado magnéticamente El problema con todo este esquema es que las detenciones en espera del semáforo, es que se trata de esperas activas, en un bucle vacío, simplemente comprobando constantemente la condición del cambio de estado. Esto es ineficiente energéticamente, pues obliga a que la CPU esté siempre al 100% de utilización. Además, según el planificador, puede que le esté robando ciclos preciosos al otro hilo que sí que está haciendo tareas provechosas (computando, o transmitiendo o escribiendo en disco). Se propusieron soluciones alternativas que realmente no parecían llegar a solucionar totalmente el problema: • Uso de señales: La idea es crear los servicios de la variable condición utilizando los servicios del sistema operativo. Se desconoce cómo usar, lanza y gestionar señales del sistema operativo en Fortran. Por ello, se descarta inmediatamente. • Dormir el hilo en espera, para que no esté en espera activa, y despertarlo pasado cierto tiempo. El problema está de nuevo en que Fortran no provee de funciones de detención inferiores a un segundo. Se llegó a crear un wrapper para invocar a una nueva función en Cque permitía realizar esperas desde un microsegundos a varios milisegundos. En las pruebas realizadas, no obstante, se comprobó que la granularidad efectiva mínima era de unos 250 µs (0,25 ms), y que este valor era variable según el entorno de ejecución. Una espera de 0,5 ms implica 2.000 esperas así cada segundo. Los tiempos de ejecución indicaban que en ese tiempo, incluso el hilo 0 computacional habría realizado algunas iteraciones. Esto implica que si el hilo computacional, el teóricamente más lento, invocara a esta llamada, se perdería un mínimo de 5 iteraciones, lo cual parecía poco razonable, puesto que a lo mejor el hilo de comunicaciones ya había terminado antes. Por estas primeras pruebas, también se descartó esta opción. Descartadas estas opciones, se pensó en utilizar los cerrojos de OpenMP (omp_locks). 4.5.2 Implementación final con cerrojos de OpenMP (omp_locks) La versión r27 final recoge las ideas y la experiencia de la anterior r26, pero la aplica a partir de la r24. Añade a esta versión dos cerrojos, con los cuales se gestiona el acceso a las iteraciones en que se realiza el cómputo estadístico, y las inmediatamente siguientes a ella. Se trata pues de una gestión de grado más grueso, puesto que ya no se pretende gestionar el caso de cada una de las estadísticas por separado, sino que simplemente se actúa sobre toda la iteración, lo cual incluso tiene cierto sentido computacional por la existencia de dependencia de datos entre las diferentes estadísticas de una misma iteración. Los requisitos son los mismos: se crean permanentemente zonas de memoria auxiliar (buffers) para almacenar los resultados estadísticos de una iteración, de tal forma que en la siguiente los datos originales puedan ser sobreescritos sin problema. De paso, toda la gestión de memoria dinámica se pasa a estática para prácticamente todas las variables, dado que se considera también muy costoso el estar constantemente reservando, inicializando y liberando memoria en cada iteración estadística. Los cerrojos en OpenMP permiten bloquear en una espera no activa al bloque que intenta tomar el cerrojo pero éste ya está tomado por otro hilo. Sólo cuando el otro hilo lo libera, el sistema operativo avisa el hilo bloqueado para que éste pueda coger el cerrojo y seguir su camino. Se trata pues de una gestión propia de OpenMP, y por lo tanto optimizada para la arquitectura subyacente y teóricamente más eficiente que el intento realizado con la versión r26 anteriormente explicada. 36 Capítulo 4. Implementación de la paralelización con CUDA Sección 4.5 Para el caso que nos ocupa, la solución planteada ha sido la utilización de dos cerrojos diferentes, así como dos variables auxiliares de seguimiento de quién es su propietario (es importante, dado que está prohibido liberar un cerrojo que no se posea, ni volver a coger un cerrojo que ya se posee, en ambos casos, la ejecución termina abruptamente con error). La idea básica es que inicialmente el hilo 0 posee ambos cerrojos, y siempre tendrá al menos uno de ellos. Nunca se le estará permitido liberar los dos cerrojos. Por su parte, el hilo 1 podrá tener o ninguno, o sólo uno de los dos cerrojos. Podrán haber cerrojos libres, sin ningún propietario, que siguiendo las reglas anteriores, y otras por definir, podrán ser tomados por cualquiera de los dos hilos, si se le permite. Dependiendo de qué cerrojo o cerrojos poseea cada hilo, podrá avanzar en su tarea, o deberá quedarse bloqueado en espera de un cambio de situación, esto es, coger uno de los otros cerrojos. Cada vez que se coge el cerrojo, el nuevo propietario anota en la variable auxiliar correspondiente que él es el propietario de ese cerrojo. Cada vez que se libera, esa variable se pone a -1, indicando que no tiene propietario. Así, cualquier hilo puede en cualquier momento saber quién es el propietario del cerrojo, o si éste está libre. Los hilos 0 y 1 tienen asignadas las mismas funciones que las descritas en el caso anterior de r26. El esquema detallado de funcionamiento es el siguiente: 1. La ejecución empieza con el hilo 0 poseyendo ambos cerrojos, que denominaremos abreviadamente WORK yWAIT. 2. Al inicio de la iteración de generación de estadísticas, el hilo 0 comprueba por las variables auxiliares si es el propietario de WORK. Si no lo es, intenta cogerlo, quedando bloqueado si no lo consigue (llamada a omp_set_lock()). 3. Si ya era propietario de WORK, entonces intenta coger WAIT, de tal forma que si falla, también se queda bloqueado hasta obtenerlo. 4. Sólo cuando el hilo 0 posee ambos cerrojos, continua la ejecución de la iteración de estadísticas. Esto es, precisa tener los dos cerrojos para elaborar estadísticas. 5. No obstante lo anterior, cuando ha obtenido los dos cerrojos, y va a continuar la ejecución, lo primero que hace es liberar el cerrojo de WAIT. 6. Eventualmente el hilo 0 termina el ciclo de generación de las estadísticas, que las guarda en las variables auxiliares al efecto, y llega a la siguiente iteración, la de estadísticas+1. 7. Cuando el hilo 0 llega al inicio de la iteración estadísticas+1 (o termina el bucle, que es equivalente), el hilo 0 intenta coger el WAIT. Si no puede cogerlo, queda bloqueado hasta conseguirlo. Si lo consigue, o cuando, tras estar bloqueado, lo consigue, es porque el hilo 1 ya lo ha tomado y lo ha devuelto, y por lo tanto, sirve para que el hilo 0 obtenga el reconocimiento del hilo 1 de que éste sabe que el hilo 0 está preparando los nuevos datos. 8. El hilo 0, tras obtener de nuevo el WAIT, libera el WORK, con lo cual ha respetado en todo momento tener al menos uno de los dos cerrojos en su posesión. 9. Pasemos al hilo 1. Cuando éste llega a la iteración de generación de estadísticas, opera de forma diferente al hilo 0. Primero mira si él es el poseedor de alguno de los dos cerrojos, WORK oWAIT, si es así, se ha producido un error inesperado en el protocolo, y directamente aborta la ejecución, terminando anómalamente el programa. 37 Paralelización mediante GPU de la dinámica de electrones en plasma confinado magnéticamente 10. Si sigue vivo (no tenía ninguno de los dos cerrojos), el hilo 1 intenta obtener el WAIT, o se queda bloqueado en espera de él. 11. Al conseguir el WAIT se entera de que el hilo 1 ya ha llegado también al hilo de generación de estadísticas, así que lo devuelve inmediatamente como reconocimento de haberse enterado. 12. Inmediatamente, intenta coger el cerrojo de WORK, dado que el hilo 0 debe estar en esos momentos generando las estadísticas, y no lo liberará hasta que el hilo 0 no termine esa iteración de trabajo estadístico, y al inicio de la iteración de estadísticas+1, recoja el WAIT para saber que el hilo 1 ya está preparado, y devuelva el WORK como señal al hilo 1 de que los datos ya están generados y en sitio seguro para poder ser gestionados por el hilo 1. 13. El hilo 1, cuando toma el WORK, se desbloquea, y continua con su iteración, que será la de estadísticas acabas de generar por el hilo 0. 14. Cuando el hilo 1 termina su iteración de intercambio y resumen de estadísticas, y almacenamiento en disco de ellas, si procedía, llega al inicio de la iteración de estadísticas+1. 15. Al inicio de la iteración estadísticas+1, el hilo 1 devuelve el cerrojo de WORK como señal de que ya ha terminado su trabajo. 16. Consecuentemente, el hilo 1 sólo ha tenido como máximo un único cerrojo en su posesión, el WAIT para reconocer la situación de trabajo del hilo 0, o el WORK cuando está gestionando los datos estadísticos generados por el hilo 0. Cuando está iterando sin hacer nada hasta la siguiente iteración de estadísticas, o bloqueado al inicio de ésta, no posee ningún cerrojo. Con este protocolo se ha implementado la versión r27, que como se verá en el Capítulo 5 “Resultados”, es más rápida que la versión de partida r24 en aproximadamente un 3–4% en el peor de los casos, excepto en un caso puntual que se explicará en su momento, y que pudiera tener relación con el entorno de ejecución. A continuación se incluye el código al inicio de la iteración principal del bucle, donde se gestionan estos cerrojos. El Código 4.9 está ampliamente comentado para seguir el protocolo anterior. Listado 4.9: Gestión de cerrojos para permitir el correcto solapamiento de cálculo y comunicación 1DO jk = 1, niter - 1 ! Main loop of temporal iterations 2IF(l_production) THEN ! Check if the statistics are active 3if (thread==0) then ! Thread 0 does the computational tasks 4if (nprod*int(jk/nprod) == jk .or. jk==niter-1 .or. 5& ncontrol*int(jk/ncontrol) == jk) then 6! Iteration when the special statistical calculations must be done 7! Master (thread==0) cannot start if thread==1 hast not finished 8! messaging and storing the previous stats 9! The thread==0 (master) has to have both locks to continue 10 ! It always have one at least, don’t need to check if have none 11 ! Then, when he has got both, he inmediatly release lock_wait 12 if (lock_work_owner .ne. thread) then 13 call omp_set_lock(lock_work) 38 Capítulo 4. Implementación de la paralelización con CUDA Sección 4.5 14 lock_work_owner = thread 15 iter_work = jk 16 elseif (lock_wait_owner .ne. thread) then 17 call omp_set_lock(lock_wait) 18 lock_wait_owner = thread 19 iter_wait = jk 20 endif 21 call omp_unset_lock(lock_wait) 22 lock_wait_owner = -1 23 elseif ((jk .ne. 1) .and. (nprod*int(jk/nprod)==jk-1 .or. 24 & ncontrol*int(jk/ncontrol)==jk-1) ) then 25 ! Iteration after the special stats have been done 26 ! The master (thread==0) has only the lwork and waits for lwait 27 ! When it gets lwat, then it releases lwork so it can be taken 28 ! by thread==1 who will message and/or store that info 29 call omp_set_lock(lock_wait) 30 lock_wait_owner = thread 31 iter_wait = jk 32 call omp_unset_lock(lock_work) 33 lock_work_owner = -1 34 endif ! end of checking iteration number for thread==0 35 else ! It’s thread 1, which does the MPI and I/O part 36 if (nprod*int(jk/nprod) == jk .or. 37 & ncontrol*int(jk/ncontrol) == jk) then 38 ! Iteration when the special statistical calculations must be done 39 ! Thread==1 cannot start if master (thread==0) hast not yet finished 40 ! this very same iteration, because needs to be sure the data is OK. 41 ! Thread==1 have no lock, we check that, and release them if so . 42 ! Then thread==1 tries to get lock_wait before continue the mes &IO work. 43 ! When it gets it, releases it inmediatly and then tries to get 44 ! the lock_work 45 if (lock_work_owner==thread .or. 46 & lock_wait_owner==thread) then 47 print *,’thread ’,thread,’have one lock at it=’,jk, 48 &’That is not allowed. STOPING!!!’ 49 STOP ! No allowed, something went wrong. Abnormal termination 50 endif 51 call omp_set_lock(lock_wait) 52 lock_wait_owner = thread 53 iter_wait = jk 54 call omp_unset_lock(lock_wait) 55 lock_wait_owner = -1 39 Paralelización mediante GPU de la dinámica de electrones en plasma confinado magnéticamente 56 call omp_set_lock(lock_work) 57 lock_work_owner = thread 58 iter_work = jk 59 elseif ((jk .ne. 1) .and. (nprod*int(jk/nprod)==jk-1 .or. 60 & (ncontrol*int(jk/ncontrol)==jk-1 .or. jk==niter-1))) 61 &then 62 ! Iteration after the special stats have been done 63 ! The thread==1 has the lock_work, and releases it 64 ! to let know the master it has already messaged and 65 ! stored previous stats 66 call omp_unset_lock(lock_work) 67 lock_work_owner = -1 68 endif ! end of checking iteration number for thread==1 69 endif ! if (thread==0) then 70 ENDIF ! IF(l_production) THEN 71 ! End of lock management, never make this region OMP CRITICAL or will deadlock 4.5.2.1 Consideraciones con la adición de OpenMP en el código Fortran Al utilizar hilos de OpenMP se pueden definir qué variables van a ser locales y cuáles van a ser compartidas entre ambos hilos participantes de la paralelización de la región dada. No obstante, hay un caso especial y es cuando desde Fortran se invoca a una subrutina. En ese momento, las variables que hasta ahora eran globales, desaparecen como tales, y pasan realmente a ser locales. Es más, se reinicializan a su valor por defecto (generalmente cero), sin ni siquiera conservar el valor que la variable global tenía en el momento de invocarse la subrutina. Ese comportamiento no es así si la creación de los hilos es dentro de la subrutina. Para solventar ese problema se hace necesario, pues, pasar como parámetros en la invocación a la subrutina de todas aquellas variables que el hilo precise conocer para la ejecución correcta de la subrutina invocada. Esto hace que sea necesaria la modificación del prototipo de casi todas las subrutinas invocadas desde el programa principal, pues es allí donde se crean los hilos. 4.6 Valoración cualitativa del proceso de paralelización con CUDA La realización de la migración paso a paso ha permitido evaluar en cada momento cuáles han sido los principales saltos cualitativos y cuantitativos obtenidos durante el proceso, así como tener que enfrentarse a decisiones para solventar aquellos problemas que han ido apareciendo. Tras este proces, se presentan de forma resumida las principales características que se consideran interesantes: • La migración a CUDA ha resultado ser, en lo que se refiere a las funciones individuales de cómputo, relativamente sencilla. • Dicha migración sí que ha presentado mayores problemas a la hora de definir la estrategia de migración. • También ha resultado ser problemático analizar el algoritmo en el código implementado, puesto que en ocasiones éste no era el más evidente ni el más óptimo. 40 Capítulo 4. Implementación de la paralelización con CUDA Sección 4.6 • Se han tenido que resolver problemas de programación interesantes en CUDA-C, tales como la utilización eficiente de los generadores de números aleatorios, la utilización y programación de diferentes tipos de reducción, y la optimización en la invocación a los kernels. • El método de los wrappers, junto con los ficheros de cabecera para acceder más fácilmente a las variables globales declaradas en Fortran ha resultado eficiente. No obstante, aporta un grado de complejidad que puede representar un escalón importante para programadores no familiarizados en C. • El análisis del código, junto con la comprensión de la parte física de las ecuaciones de Langevin con ellas implementada, ha permitido la detección de dos errores de programación que se han corregido en la implementación propuesta, así como de la propuesta e implementación de ella de algoritmos y uso de funciones más eficientes que los métodos numéricos utilizados en el código original. • La inclusión de OpenMP en el código ya paralelizado, para obtener una aceleración aún mayor, ha resultado ser mucho más complejo que los resultados con él obtenidos. Evidentemente esto depende del caso concreto, pero el grado de complejidad, no solo en programación, sino también en la necesidad de idear algoritmos correctos libres de bloqueo, hace que esta última parte no sea recomendable intentarla si no se tiene a una persona muy experta en computación paralela y distribuida en el equipo. Se incluye la distribuida, porque el funcionamiento asíncrono de los hilos plantea problemas alejados de la computación paralela, y mucho más propios de la distribuida, con su necesidad de estudio de corrección y viveza de los algoritmos que se propongan implementar. • Desde el punto de vista cuantitativo, la migración es un total éxito, dado que es correcta en sus resultados, y aporta una aceleración en torno a 700 respecto al código original ejecutándose en un núcleo de CPU. Económicamente, si las simulaciones son muchas, implica un importante ahorro en el tiempo de uso del cluster de producción utilizado, y con ello, también económico. 41 Paralelización mediante GPU de la dinámica de electrones en plasma confinado magnéticamente (Página intencionadamente en blanco) 42 CA P Í T U L O 5 RESULTADOS En este capítulo se presentan los resultados obtenidos en diferentes entornos con la ejecución del programa implementado, o de algunas de sus versiones previas, si ya era suficiente para sacar conclusiones válidas. Las pruebas se han realizado en varios entornos de ejecución diferentes, por lo que será interesante presentar una breve reseña de los mismos con antelación a la presentación de dichos resultados. Finalmente, se efectuará un análisis de los resultados, que permitirá valorar la implementación, en sus dos aspectos cuantitativos más importantes a valorar: • su corrección, esto es, si produce los resultados que de él se espera, • su eficiencia, y por lo tanto, qué rendimiento se extrae de su utilización, de forma comparada con otras implementaciones existentes. No obstante, y recordando que la paralelización del código en concreto utilizado es una caso de prueba y demostración, se debe valorar cualitativamente el proceso de paralelización utilizado, en comparación con otras posibles alternativas. Todo ello en vistas a poder informar a la parte científica de los requisitos que debe contar el personal que deba acometer en el futuro esta labor, personal que actualmente, en su mayor parte, está formada por la propia comunidad científico teórica del campo en cuestión, la física de partículas. 5.1 Pruebas realizadas Las pruebas se han realizado en dos etapas diferentes: • inicio de la implementación, para la optimización del tamaño de la malla y estimar ya desde las primeras etapas del trabajo la aceleración mínima obtenible, • con las implementaciones finales terminadas de las dos versiones del código, la puramente MPI con CUDA, y la que se añade OpenMP para desligar las operaciones computacionales de los kernels de las de comunicación y escritura en disco. 43 Paralelización mediante GPU de la dinámica de electrones en plasma confinado magnéticamente A continuación se detallan ambas etapas, y las razones y objetivos de las mismas: 5.1.1 Primera etapa: optimización del kernel principal según los parámetros propios de su invocación En la primera etapa el trabajo se focalizó en comprobar si el proyecto era viable, o mejor dicho, si el código inicialmente escogido era paralelizable con CUDA y susceptible de tener una buena aceleración. Para ello se identificó el mínimo trozo de código que era responsable de una gran proporción del tiempo total de ejecución. En la Figura 5.1 se representa gráficamente el diagrama de llamadas obtenido con el profiling del programa original. En él se detectó que la subrutina advance, junto con aquellas a las cuáles ésta llamaba (rhs_coeffs,phi_and_psi,rhs_ion yget_gaussian_noise) era la consumidora del 85% del tiempo de ejecución, al cual hay que añadir prácticamente otro 13% por la parte correspondiente de las subrutinas estáticas de biblioteca hpsort y sobre todo de ran2, que son principalmente utilizadas por las subrutinas dependientes de advance. 84.54% 999000× 50.91% 1998000× 18.17% 1998000× 10.00% 1998000× 85.45% 1× 39.09% 1998000× MAIN__ 85.45% (0.91%) 1× advance_ 84.54% (5.45%) 999000× rhs_coeffs_ 50.91% (11.82%) 1998000× get_gaussian_noise_ 18.18% (18.18%) 1999000× rhs_ion_ 10.00% (10.00%) 1998000× main 85.45% (0.00%) phi_and_psi_ 39.09% (39.09%) 1998000× ran2_ 13.64% (13.64%) hpsort2_ 0.91% (0.91%) Figura 5.1: Diagrama de llamadas del programa original obtenido del perfilado Es por ello que la primera etapa de implementación consistió en paralelizar mediante kernels de CUDA la subrutina advance y todas aquellas que son llamadas por ésta. La ejecución del código preliminar así obtenido mostró que efectivamente esta parte del código era la que más afectaba al tiempo de ejecución, al obtenerse ya un speedup en torno a 30–50, según el entorno de ejecución. Tras optimizar dicho código, se lanzó una prueba de optimización de los parámetros de llamada a los 44 Capítulo 5. Resultados 5.1. Pruebas realizadas kernels. Adicionalmente a los típicos parámetros de entrada/salida de toda función C/C++, los kernels de CUDA se lanzan indicando entre dos y cuatro parámetros más, que se incluyen en la llamada con una notación especial que precisamente caracteriza su código de llamada: Kernel_Name<<< GridSize, BlockSize, SMEMSize, Stream >>> (arguments,....); Los valores relevantes para nuestro caso son BlockSize yGridSize, que definen respectivamente la cantidad de bloques y cómo éstos se estructuran al lanzar la ejecución del kernel en la GPU. Aunque estos parámetros son de un tipo especial de vector de tres dimensiones (x,y,z) definido por el propio CUDA, en nuestro código no es necesaria tal complejidad y se utilizará únicamente las dimensiones xeyen el del tamaño del bloque y xen el tamaño de la malla. La forma en que se invocan los kernels en la GPU es importante por las siguiente razones: • el número de bloques junto con el tamaño de cada bloque determinará el número de ejecuciones individuales que se realizará del código. Debe ser el número suficiente para que acoja, en nuestro caso, todas y cada una de las partículas a simular, puesto que cada hilo de ejecución se encargará de una de ellas en exclusiva, • el número de bloques (GridSize) tiene que ser suficientemente elevado como para permitir la máxima utilización de los multiprocesadores de la tarjeta gráfica, • el tamaño del bloque (BlockSize) determinará un factor muy importante para la eficiencia de la ejecución paralela en la tarjeta gráfica, que es el grado de ocupación de los registros, y con ello, el grado de paralelismo solapado existente en la ejecución. Téngase en cuenta que hay mucho cálculo de doble precisión, por lo que las unidades de coma flotante estarán constantemente en utilización, y conviene utilizar el tiempo en que ese warp está inactivo en espera del resultado, en otros cálculos utilizando registros y otras unidades lógicas del core. La prueba consiste en variar estos parámetros, teniendo siempre en cuenta que existe un límite máximo según la arquitectura de número de hilos que pueden ejecutarse concurrentemente en un bloque. Para ello, se varían las dimensiones xeyde BlockSize, cambiándose el valor de GridSize para que se cumpla que se lancen los mínimos bloques necesarios para que cada iteración (partícula) sea ejecutada una y al menos una vez por un único hilo. Estas prebas se realizan en diversos entornos reducidos a una única máquina, dado que en este caso, el resultado no depende de la masividad del cómputo, sino de las características hardware de la tarjeta gráfica utilizada, estudiando para ello tres diferentes generaciones de arquitectura CUDA:Fermi,Kepler y Maxwell (primera generación). 5.1.2 Segunda etapa: medición de prestaciones de las versiones finales en un entorno HPC real Una vez obtenidas las dimensiones óptimas para el lanzamiento del principal kernel, el trabajo se centra en paralelizar el resto de funciones que están dentro del bucle principal del programa. De esta forma se obtiene un programa final que realiza prácticamente todos los cálculos sobre cada partícula utilizando CUDA. Esto permite evitar traspasos continuos de memoria entre la CPU y la GPU en cada iteración temporal. Sólo se realizan aquellas transferencias de memoria necesarias, siempre de GPU aCPU, para enviar los datos mínimos necesarios para que la CPU pueda guardar los datos estadísticos y de progreso de la simulación en disco (o presentarlos por consola). 45 Paralelización mediante GPU de la dinámica de electrones en plasma confinado magnéticamente También se obtendrán conclusiones interesantes sobre aspectos inicialmente no planificados, como por ejemplo, influencia del lanzamiento de un número no divisible por 2 del tamaño de los bloques, detección de posibles problemas con la utilización de OpenMPI con InfiniBand. Las pruebas se realizan en el entorno de ejecución de MinoTauro. Se utiliza un número variable de nodos, así como un número variable de GPU. También se varía el número de GPU por nodo que se utiliza en cada prueba, para ver si la mayor localidad de los procesos en menos nodos incrementa la eficiencia de la ejecución de la simulación. Se seleccionan diferentes tamaños de carga, variando tanto el número de iteraciones temporales, como el número de partículas que son consideradas en cada iteración. También se cambia el tipo de simulación (con generación o no de resultados estadísticos intermedios, con fuerte carga de intercambio de mensajes y escritura en disco). Y finalmente se utilizan tres programas: el código original en Fortran+MPI; la nueva implementación sustituyendo la parte computacional con CUDA, y algo de C para facilitar la invocación a la API; y la que mejora ésta última con la adición de OpenMP para desacoplar el cómputo de la comunicación, y así solapar lo más posibles estas dos actividades. Todas estas pruebas se hacen de una forma ortogonal (ortogonalidad no completa, pero sí suficiente para extraer conclusiones adecuadas). De esta forma se puede estudiar la escalabilidad del programa bajo diferentes puntos de vista: espacial, temporal y localidad, así como evaluar si la inclusión de CUDA primero y de OpenMP después compensan respecto de su respectivo tiempo de desarrollo. A continuación se procede a describir, presentar resultados, analizar los resultados, y obtener conclusiones sobre cada una de las pruebas realizadas. 5.3.2.1 Caracterización del comportamiento del programa original Para poder evaluar la nueva implementación, se necesita hacerlo por comparación con el código original de partida. En esta prueba servirá para caracterizar el comportamiento del programa original, ya que la nueva implementación deberá mejorar los aspectos más débiles del original, sin que los aspectos fuertes de éste sufran merma. La caracterización consistirá en analizar su escalabilidad espacial y temporal. Breve descripción de las pruebas Consiste en la ejecución del código original en el entorno de MinoTauro mediante las siguientes pruebas: 1. Con una carga pequeña constante (10.000 iteraciones temporales de 20.000 partículas), variar el número de procesos que se ejecutan en un único nodo. Con esto se comprobará la escalabilidad en un entorno puro de memoria compartida. 2. Con la misma carga pequeña constante anterior, ejecutar siempre con 12 procesos, pero con un número variable de nodos, variando por lo tanto el número de procesos por nodo. Así se podrá comparar la velocidad de intercomunicación en memoria compartida con la existente entre diferentes nodos (memoria distribuida). 3. Con una carga grande constante (500.000 iteraciones temporales de 200.000 partículas), variar el número de nodos, y por lo tanto de procesos, haciendo siempre que cada nodo trabaje al 100% de su capacidad, esto es, con 2x6 =12 procesos por nodo. La idea es comprobar la escalabilidad del programa en un entorno mixto de memoria compartida y distribuida a gran escala. Los tamaños son escogidos con dos ideas: 52 Capítulo 5. Resultados 5.3. Ejecución de las pruebas • La carga pequeña es para que las ejecuciones no sean excesivamente largas, sobre todo cuando se ejecutan en un único nodo. • La carga mayor está relacionada con la memoria máxima que una única GPU del cluster puede soportar, y el número de iteraciones es uno ya suficientemente elevado como para que el fenómeno del plateau (meseta) en la formación de los electrones fugitivos sea claro, se vea claramente estabilizado, en estado estacionario. Resultados obtenidos Las pruebas 1 y 2 se realizan en un único nodo, para ver la escalabilidad en memoria compartida. En la Tabla 5.4 se muestran los resultados obtenidos. Procesos Nodos Procesos Tiempo de Aceleración Eficiencia por nodo ejecución 1 1 1 2045,806 s 1,000 1,000 2 1 2 1021,929 s 2,002 1,001 3 1 3 683,256 s 2,994 0,998 4 1 4 511,999 s 3,996 0,999 6 1 6 343,132 s 5,962 0,994 9 1 9 228,449 s 8,955 0,995 12 1 12 170,748 s 11,981 0,998 12 2 6 171,560 s 11,925 0,994 12 3 4 173,361 s 11,801 0,983 12 4 3 171,696 s 11,915 0,993 12 6 2 172,575 s 11,855 0,988 12 12 1 170,895 s 11,971 0,998 Tabla 5.4: Aceleración y eficiencia del código original de CPU (memoria compartida y distribuida) En la Gráfica 5.2, en representación doble logarítmica, se representan los datos obtenidos de la Prueba 3, con gran cantidad de procesos y una carga grande. Se observa la disminución del tiempo de ejecución conforme se incrementa el número de nodos utilizados en la computacion de forma muy proporcional al número de procesos participantes. Fenómenos observados En la pruebas 1, sobre un único nodo, esto es, sólo con memoria compartida, tal como se ve en la Tabla 5.4 en las filas sombreadas de amarillo, la escalabilidad es prácticamente lineal en todo el intervalo. Lo mismo es cierto para cuando se cambia progresivamente a memoria distribuida (filas sombreadas de azul), donde igualmente la eficiencia es prácticamente la unidad en todos los casos, nunca bajando del 98%. En la prueba 3, mezcla de memoria compartida y distribuida, Figuras 5.2 y 5.3, se observa que la línea es casi recta, con una pendiente de prácticamente la unidad, lo que adelanta que la aceleración será bastante lineal. Ese hecho se comprueba en la Gráfica 5.3 que visualiza que tanto la aceleración como la eficiencia se aproximan a la idealidad (eficiencia muy cercana a 1 incluso con 64 nodos, que son 768 procesos, y la línea muy coincidente con la ideal). Aún así, con 64 nodos (768 procesos MPI) la eficiencia sigue siendo muy alta, prácticamente del 95%. 53 Paralelización mediante GPU de la dinámica de electrones en plasma confinado magnéticamente 1; 66988 4; 16910 16; 4382 32; 2142 64; 1105 1 10 100 1000 10000 100000 1 10 100 Número de nodos (cada nodo ejecuta 12 procesos en sendos núcleos) Tiempo (s) Tiempo ejecución total Figura 5.2: Escalabilidad temporal del código original de CPU (500.000 partículas y 200.000 iteraciones) 16; 15.287 32; 31.276 64; 60.615 4; 3.962 1; 1.000 0.977 0.947 0.990 1.000 0.955 0 10 20 30 40 50 60 70 0 8 16 24 32 40 48 56 64 72 Número de nodos (cada nodo ejecuta 12 procesos en sendos núcleos) Aceleración (-) 0.880 0.900 0.920 0.940 0.960 0.980 1.000 1.020 Eficiencia (-) Aceleración Escalabilidad lineal perfecta Eficiencia Figura 5.3: Aceleración y eficiencia del código original de CPU (500.000 partículas y 200.000 iteraciones) Conclusiones de las pruebas A raíz de estas tres pruebas, se obtienen las siguientes conclusiones respecto al programa original en Fortran paralelizado con MPI ejecutándose sólo en CPU y también sobre el entorno de ejecución de MinoTauro. • El programa original escala muy bien en memoria compartida, pero al ser tan lento, no es suficiente para simulaciones que permitan extraer conclusiones físicas interesantes. • El programa original escala muy bien en un cluster HPC, lo que hasta ahora ha permitido realizar simulaciones suficiente para extraer conclusiones. • La red de comunicaciones del cluster MinoTauro es muy rápida. No hay diferencia sustancial entre el uso de memoria compartida y la memoria distribuida. No obstante, se observa que para poder alcanzar tamaños de simulación interesantes, se necesitangrandes recursos computacionales durante un tiempo nada despreciable. Como se verá, la paralelización con CUDA permitirá reducir drásticamente esta relación entre tiempo necesario y recursos hardware utilizados, lo que implicará o un importante ahorro de explotación, y la posibilidad de lanzar simulaciones mucho más grandes que permitan obtener resultados más conclusivos, o disminuir la cantidad de aproximaciones realizadas para alcanzar un conjunto de ecuaciones más complicadas de computar, pero más fieles respecto del comportamiento esperado del plasma de partículas simulado. 5.3.2.2 Caracterización básica de la implementación consistente en la paralelización con CUDA A continuación se va a proceder a estudiar el comportamiento de la versión final del programa en el que, basándose en el código original en Fortran+MPI ejecutable en CPU, se realiza una paralelización de toda la parte de computación numérica y del cálculo de las estadísticas intermedias mediante CUDA. De esta forma, Fortran queda sólo como soporte básico para la lectura y escritura de ficheros y la invocación a la librería MPI para el intercambio de datos entre los diferentes procesos, así como también del flujo de control básico del programa, incluyendo la gestión de las iteraciones temporales, y de cuándo corresponde realizar las estadísticas. La implementación ha sido tal que tan sólo es necesaria la transferencia de los datos de velocidades de las partículas al principio y al final del programa. Durante todas las iteraciones temporales, la transferencia de datos se realiza únicamente de la GPU a la CPU solamente en aquellas iteraciones en las que hay que 54 Capítulo 5. Resultados 5.3. Ejecución de las pruebas recopilar y guardar en disco las estadísticas intermedias solicitadas, y sólo de los mínimos datos necesarios. De esa forma se ha conseguido disminuir mucho la transferencia de memoria entre la GPU y la CPU, lo que permite incrementar el porcentaje de código que el programa está realizando cómputos, y por lo tanto, por la Ley de Amdahl, permitirá un mucho mejor aprovechamiento del poder de la paralelización, como se observará a través de los resultados que a continuación se presentan y discuten. Al igual que en la sección anterior, primero se detallarán las pruebas realizadas, después se presentarán los datos cuantitativos resumidos obtenidos de dichas pruebas, que permitirán observar las principales características de la ejecución, y con ello, caracterizar la implementación, a partir de lo cual se extraerán conclusiones sobre ella. 5.3.2.2.1 Prueba de corrección Durante todo el proceso de implementación se ha ido comprobando que los resultados de cada parte que se implementaba no cambiaban los resultados intermedios estadísticos que genera el programa. En realidad, los mismos resultados no pueden generarse, dado que se cambió el generador de números aleatorios. Sin embargo, forzando a que unos pocos números sean los mismos se comprobó que las velocidades y su clasificación en contenedores fueran iguales a los del código original. De hecho, durante esta revisión se localizaron varios posibles bugs en el código original que fueron puestos en conocimiento de sus responsables. A partir de ese momento, los números generados aleatoriamente se fijaron para la GPU hasta la versión final, y aunque al procesar con diferente número de GPU esos valores vuelven a cambiar, sí que se observa que la distribución de los valores en cada columna en los ficheros generada es estadísticamente similar a los del código original. Además, visualizando los ficheros .silo con VisIt, se observa que los mapas de velocidades y su distribución son similares a los que se muestran en [4]. Es por todo ello que se considera que las implementaciones presentadas son correctas en cuanto a la fidelidad de la resolución de las ecuaciones físicas que se usan para simular el plasma confinado magnéticamente. 5.3.2.2.2 Aceleración y eficiencia respecto del código en CPU Breve descripción de las prueba La prueba consiste en ejecutar la misma carga que se utilizó en la prueba de la caracterización de la CPU, cuyos datos se han mostrado gráficamente en la Figuras 5.2 y 5.3, esto es, 200.000 iteraciones temporales sobre 500.000 partículas. Se realizan varias simulaciones variando dos parámetros: • El número de procesos, con cada proceso utilizando una GPU. Se varía desde 1 nodo con 1 GPU (caso base) hasta 128 procesos, utilizando 128 GPU. • El número de GPU que se utilizan en cada nodo. Las mismas pruebas se realizarán utilizando una GPU por nodo, o dos. Esto último permitirá establecer si hay alguna diferencia importante entre la ejecución dentro de un mismo nodo o con un mayor número de nodos, y por lo tanto, con una mayor utilización de la red de interconexión. Resultados obtenidos 55 Paralelización mediante GPU de la dinámica de electrones en plasma confinado magnéticamente En la Tabla 5.5 se muestran los datos de esta prueba. Los experimentos se han lanzado por duplicado, siendo muy coincidentes en todos los casos parejos, por lo que en la tabla se muestra el valor promedio medido. Procesos Nodos GPU nodo Tiempo de ejecución Aceleración Eficiencia Aceleración Aceleración Total GPU I/O MPI Otros GPU/CPU relativizada 1 1 1 1181,871 1177,111 4,559 0,057 0,144 1,000 1,000 680 680 2 2 1 600,646 597,715 2,596 0,386 −0,051 1,968 0,984 1338 669 4 4 1 317,959 315,880 1,631 0,525 −0,077 3,717 0,929 2528 632 8 8 1 167,233 165,699 1,141 0,515 −0,122 7,067 0,883 4807 601 16 16 1 93,997 92,619 0,901 0,583 −0,106 12,573 0,786 8552 534 32 32 1 55,610 54,437 0,792 0,565 −0,184 21,253 0,664 14455 452 64 64 1 44,437 38,298 0,716 5,525 −0,101 26,597 0,416 18090 283 2 1 2 601,306 598,022 2,780 0,512 −0,007 1,966 0,983 1337 668 4 2 2 317,888 315,954 1,625 0,370 −0,060 3,718 0,929 2529 632 8 4 2 167,683 166,181 1,151 0,431 −0,079 7,048 0,881 4794 599 16 8 2 93,997 92,619 0,901 0,583 −0,106 12,573 0,786 8552 534 32 16 2 55,610 54,437 0,792 0,565 −0,184 21,253 0,664 14455 452 64 32 2 44,437 38,298 0,716 5,525 −0,101 26,597 0,416 18090 283 128 64 2 26,550 25,318 0,656 1,086 −0,509 44,515 0,348 30277 237 Tabla 5.5: Aceleración y eficiencia de la implementación Fortran+MPI+CUDA con misma carga que Tabla 5.4 Para el cálculo de la aceleración sobre la CPU, y dado que el código original escala linealmente de forma perfecta, se ha tomado el tiempo total de ejecución en un nodo completo, corrigiéndose el hecho de que ese tiempo se corresponde con 12 procesos. Estos datos se visualizan en las Gráficas 5.4. 0 5000 10000 15000 20000 25000 30000 35000 40000 0 16 32 48 64 80 96 112 128 Número de procesos (gpu) Aceleración respecto 1 CPU (-) 0 100 200 300 400 500 600 700 800 Aceleración relativizada (-) Aceleración 2GPU/nodo respecto CPU Aceleración 1GPU/nodo respecto CPU Aceleración 2GPU/nodo relativizada Aceleración 1GPU/nodo relativizada (a) Comparación con el código original en CPU 0 20 40 60 80 100 120 140 0 16 32 48 64 80 96 112 128 Número de procesos (gpu) Aceleración entre GPU (-) 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 Eficiencia (-) Acel. 2GPU/n CPU Acel. 1GPU/n sobre CPU Escalabilidad lineal perfecta Acel. 2GPU/n relat. Acel. 1GPU/n relat. Eficiencia unidad (b) Comperación entre las propias versiones en GPU Figura 5.4: Aceleración y eficiencia de implementación Fortran+MPI+CUDA con misma carga que Tabla 5.4 Fenómenos observados En la Tabla 5.5 se observan los siguientes fenómenos: • La aceleración de una única GPU frente a un único núcleo de CPU alcanza el elevado valor de 680. Esto quiere decir que una única GPU realiza el mismo trabajo prácticamente que 32 nodos biprocesador como los existentes en MinoTauro. 56 Capítulo 5. Resultados 5.3. Ejecución de las pruebas • Conforme se utilizan más GPU la aceleración incrementa, sin notarse que llegue un momento de desaceleración. • Conforme se utilizan más GPU la aceleración relativa decrementa. Esto es, añadir el doble de nodos con GPU al mismo trabajo no disminuye a la mitad el tiempo global de ejecución. Esto es, la eficiencia se observa que baja, hasta llegar a un 35% al usar 128 nodos. Este fenómeno es el esperable, dado que el mantenerse constante trabajo total, al incrementarse el número de GPU, disminuye su carga individual, por lo que en general, al disminuir la carga computacional por nodo, es esperable cierto decremento de prestaciones. En todo caso, en la práctica, es una situación que no suele darse, puesto que el lanzamiento se suele ya realizar sobre los mínimos procesos necesarios para evitar el coste de comunicaciones y aprovechar la localidad de la carga. • Al igual que en pruebas anteriores, la misma carga repartida entre el doble de nodos que trabajan con la mitad de sus recursos (sólo 1 procesador CPU y sólo 1 tarjeta gráfica) tiene prácticamente el mismo tiempo de ejecución que si se ejecutan de forma más concentrada en la mitad de nodos a tope de capacidad (2 procesadores CPU con 2 GPU). • Sorprendentemente, cuando se utilizan 64 nodos en MinoTauro, en todas las ejecuciones realizadas en esta prueba, el tiempo parcial empleado por la aplicación en comunicaciones MPI se incrementa de manera súbita de 0,5 segundos a 5,5 segundos. Cuando se utilizan 128 nodos, vuelve a baja a un valor justificable entre 0,5 y 1,0 segundos. • Los tiempos parciales de escritura (I/O) decrementan conforme aumenta el número de nodos. Esto es así por las características del sistema de ficheros distribuido GPFS que usa MinoTauro. Conclusiones de las pruebas De las observaciones anteriores se concluye que: • La implementación en GPU del método aproximativo de Langevin es muy eficiente comparada con la original en CPU+MPI, con una aceleración de 680. • Esta primera implementación con CUDA no escala bien cuando el número de GPU incrementa más allá de 16, aunque sigue acelerándose el cómputo. Este hecho es el que impulsa la siguiente versión añadiendo OpenMP para desacoplar más el cálculo de la comunicación y la escritura en disco. • Se confirma que la red de interconexión MinoTauro es excelente. No obstante, hay algunas combinaciones de número de nodos, quizá ligada a la topología resultante y probablemente también a algún problema entre OpenMPI eInfiniBand, que hace que el tiempo de comunicación por MPI se dispare, con la consiguiente pérdida de rendimiento global. 5.3.2.3 Caracterización básica de la implementación con CUDA yOpenMP Consiste en las mismas pruebas y metodología que en el anterior apartado 5.3.2.2. Los resultados son muy similares a la versión anterior, sólo que un 3–5% más rápidos, según la carga. La aceleración máxima obtenida en este caso es de 695 sobre el código original en CPU. Por mor de la brevedad, se omite la tabla y las gráficas de resultados, al no aportar mayor información, y ser redundantes con las de la primera implementación por su semejanza. 57 Paralelización mediante GPU de la dinámica de electrones en plasma confinado magnéticamente 5.3.2.4 Pruebas masivas de escalabilidad de las implementaciones con CUDA Tras observar el funcionamiento básico de la nueva implementación, y siendo ésta correcta y obteniendo resultados de aceleración muy importantes respecto al código original en CPU, se procede a realizar varias pruebas masivas a nivel de producción de los dos códigos presentados, al que sólo se le ha incorporado CUDA (referenciada como programa o versión r24, con Fortran+MPI+CUDA) y al que añade a la anterior OpenMP para gestionar dos hilos en cada proceso, uno encargado de la computación numérica en la GPU principalmente y otro responsable de las comunicaciones y la escritura en disco de los resultados (referenciada como r27). Breve descripción de las pruebas Las pruebas se realizan en el entorno de ejecución de MinoTauro. Se realizan diferentes lanzamientos variando la carga, para lo cual básicamente se mantiene constante el número de iteraciones temporales, fijándolas en 200.000 iteraciones (dado que llega un momento en que se alcanza el estado estacionario y no tiene sentido físico alguno continuar la simulación), y variar el número de partículas, desde muy pocas a muchas. Es importante indicar que existe un límite para el número de partículas que se pueden lanzar en una única GPU. La tarjeta gráfica NVIDIA M2090, que dispone de 6GB de memoria RAM, y con la actual implementación, sólo permite un máximo de unas 825.000 partículas. Para trabajar con un poco de margen, se escoge lanzar simulaciones con 500.000 u 800.000 partículas por cada GPU. Presentan las siguientes características, lanzadas de forma ortogonal entre todas ellas: • Carga de entre 1 millón y 50 millones de partículas, siempre con 200.000 iteraciones temporales. • Se varía el número de GPU, desde el mínimo número posible por la capacidad en memoria de las tarjetas gráficas, que coincide con unas 825.000 partículas/GPU, hasta el máximo disponible en el cluster por política de utilización, que es de un máximo de 128 GPU. • Se prueban las dos versiones r24 (con CUDA) y r27 (con CUDA más OpenMP). Estas pruebas se componen de forma ortogonal, barriendo todas las posibilidades combinativas, pero por mor de la brevedad, sólo se mostrarán las más relevantes, sobre todo teniendo en cuenta que se comprueba que alguna de las variables inicialmente prevista no tiene en la práctica, en este entorno de ejecución, ninguna influencia estadística. 5.3.2.4.1 Resultados obtenidos Una importante observación tras las pruebas, y concordante con lo visto con anterioridad, es que da prácticamente lo mismo lanzar Xprocesos en Xnodos (utilizando sólo 1 GPU por nodo), que lanzar esos mismos Xprocesos en X 2nodos (utilizando 2 GPU por nodo). Es por ello que en adelante sólo se mostrarán los datos correspondientes a la utilización de 2 GPU por nodo, excepto en el caso de las pruebas con 128 nodos, que por política del cluster, sólo podían ser lanzadas en 64 nodos. En las Tablas 5.6 y 5.7 se muestran los resultados de los tiempos de ejecución de las diferentes pruebas. Las pruebas se realizan por duplicado, repitiéndose si se observa algún valor extraño, que generalmente se presenta como un valor extrañamente alto en el tiempo parcial de las comunicaciones. En las tablas se muestran el valor promedio de las ejecuciones realizadas, eliminando los valores extremos extraños, que suelen ser en torno al 10–15% de las ejecuciones, especialmente las que se ejecutan con mayor número de GPU() o nodos. 58 Capítulo 5. Resultados 5.3. Ejecución de las pruebas Promedio de Tiempo Partículas Procesos Parte 1,000,000 2,000,000 4,000,000 8,000,000 16,000,000 32,000,000 50,000,000 2 GPU 1177.366 I/O 4.628 MPI 0.570 Otros 0.024 Total 1182.587 4 GPU 597.950 1177.426 I/O 2.663 4.643 MPI 0.364 0.397 Otros -0.065 -0.072 Total 600.912 1182.394 8 GPU 315.947 597.928 1177.343 I/O 1.621 2.784 4.596 MPI 0.550 0.488 0.802 Otros -0.091 -0.087 -0.006 Total 318.027 601.112 1182.736 16 GPU 165.991 316.038 597.802 1177.352 I/O 1.126 1.621 2.618 4.636 MPI 0.447 0.503 0.460 0.470 Otros -0.141 -0.089 -0.105 -0.076 Total 167.422 318.072 600.775 1182.382 32 GPU 92.590 166.030 315.896 597.833 1177.577 I/O 0.848 1.133 1.634 2.621 4.598 MPI 0.751 0.536 0.557 0.629 0.570 Otros -0.195 -0.149 -0.165 -0.149 -0.125 Total 93.994 167.550 317.922 600.929 1182.620 64 GPU 54.031 92.574 166.073 315.856 597.692 1177.630 1861.650 I/O 0.763 0.906 1.139 1.633 2.619 4.569 6.751 MPI 0.773 0.889 0.719 1.238 9.418 0.716 8.898 Otros -0.299 -0.467 -0.223 -0.291 -0.253 -0.243 -0.169 Total 55.267 93.902 167.707 318.436 609.475 1182.672 1877.130 128 GPU 54.140 92.571 166.113 315.868 597.822 954.424 I/O 0.788 0.932 1.209 1.660 2.627 3.681 MPI 11.897 23.777 49.865 101.112 204.925 355.542 Otros -0.588 -0.535 -0.618 -0.691 -0.734 -0.515 Total 66.237 116.745 216.569 417.949 804.639 1313.133 Tabla 5.6: Tiempos de ejecución de las pruebas de producción para la versión r24 (CUDA+MPI) Para poder observar mejor el comportamiento, para cada una de las implementaciones y tipo de simulación, se presentan también las siguientes gráficas: Figura 5.5: Tiempos de ejecución en función del número de GPU, para cada una de las cargas (número de partículas), en escala doble logarítmica. Permitirá visualizar cómo escala el programa, puesto que lo ideal sería obtener una recta de pendiente la unidad. Figura 5.6: Tiempos de ejecución en función del número de partículas de la carga, según el número de GPU utilizado. Figura 5.7: Tiempos parciales de la parte de comunicaciones MPI en función del número de GPU y del número de partículas simuladas. Permite comprobar si la pérdida de escalabilidad se puede achacar principalmente a esta parte del programa, por lo que es interesante observar cómo varía también entre las diferentes implementaciones y tipos de simulación realizadas. 59 Paralelización mediante GPU de la dinámica de electrones en plasma confinado magnéticamente Promedio de Tiempo Partículas Procesos Parte 1,000,000 2,000,000 4,000,000 8,000,000 16,000,000 32,000,000 50,000,000 2 GPU 1127.802 I/O 3.932 MPI 1.416 Otros 5.009 Total 1138.158 4 GPU 570.129 1127.719 I/O 2.413 3.817 MPI 0.573 0.884 Otros 5.479 5.422 Total 578.594 1137.842 8 GPU 299.392 569.961 1127.621 I/O 1.759 2.446 3.700 MPI 0.576 0.626 2.696 Otros 5.536 5.528 4.884 Total 307.263 578.560 1138.901 16 GPU 155.517 299.533 570.070 1127.730 I/O 1.243 1.671 2.334 3.776 MPI 0.471 0.396 0.499 0.552 Otros 5.533 5.711 5.700 5.987 Total 162.763 307.310 578.602 1137.905 32 GPU 85.486 155.836 299.634 570.313 1127.866 I/O 0.945 1.263 1.675 2.397 3.748 MPI 0.750 0.546 0.702 0.540 3.493 Otros 5.460 5.475 5.831 5.578 5.938 Total 92.641 163.119 307.842 578.828 1141.045 64 GPU 50.852 86.174 156.630 300.415 571.260 1128.100 I/O 0.939 0.977 1.281 1.796 2.399 3.830 MPI 0.733 0.628 0.619 0.588 3.734 9.993 Otros 4.704 5.360 5.354 5.510 5.554 5.661 Total 57.228 93.138 163.884 308.309 582.947 1147.586 128 GPU 52.316 87.130 157.008 302.132 657.749 1141.044 I/O 1.017 1.146 1.478 2.077 2.578 3.297 MPI 5.944 19.578 45.918 111.002 365.214 801.091 Otros 6.205 7.089 7.253 -8.189 -241.537 -670.985 Total 65.481 114.941 211.656 407.021 784.003 1274.447 Tabla 5.7: Tiempos de ejecución de las pruebas de producción para la versión r27 (CUDA+MPI) Figura 5.8 a) y b): Eficiencia temporal relativa del programa en función del número de GPU, para cada carga. Se ha denominado relativa porque no se utiliza como referencia el código original Fortran en CPU, sino que se compara con la ejecución que utiliza menos recursos. En este caso, se trata se compara con la ejecución en un único nodo con dos GPU, dado que la carga menor no puede ser ejecutada por una única GPU por falta de memoria. Figura 5.8 c) y d): Eficiencia espacial relativa del programa, esto es, cómo se comportan los programas con el incremento de la carga, representada por el número de partículas, pues el número de iteraciones se mantiene constante en 200,000 en todas estas pruebas. Figura 5.9 Comparación directa entre los tiempos de ejecución de las dos versiones con GPU presentadas (r24 y r27). Es importante explicar por qué aparece el valor negativo en la parte de computación “Otros”, que es 60 Capítulo 5. Resultados 5.3. Ejecución de las pruebas aquella que no está contabilizada ni como de cálculo en la GPU y su gestión, ni de transmisiones MPI, ni de escritura en disco duro (“I/O”). Aparece sobre todo en la versión con OpenMP, y especialmente cuando más procesos están ejecutando, y en las cargas mayores. Esto es debido a que hay dos hilos dentro de cada proceso, y el hilo que se encarga de la transmisión MPI contabiliza los tiempos de comunicación y también la de escritura en disco. Por su parte, el hilo que ejecuta los kernels cronometra el tiempo de GPU utilizado. El hecho de que la suma de la computación numérica más la de comunicaciones y escritura en disco es mayor al tiempo total de ejecución real de la aplicación demuestra que existe asincronía entre ambos conjuntos de operaciones. Así se justifica también que el tiempo de ejecución global sea menor que en la versión r24, sin OpenMP, dado que allí están seralizadas la computación, la computación y la escritura en disco, de tal forma que cuando aumenta el tamaño del problema, o el número de nodos participando de la computación, aumenta la cantidad de datos a transmitir y a escribir en disco. 5.3.2.4.2 Escalabilidad de las implementaciones Para estudiar la escalabilidad de las implementaciones presentadas nos serviremos de las Figuras 5.5 y 5.6, donde se muestran los resultados obtenidos enfocando la representación a un análisis espacial o temporal, respectivamente. Número de partículas en la simulación 10 100 1000 10000 1 10 100 1000 Número de procesos (gpu) Tiempo de ejecución (s) 1,000,000 2,000,000 4,000,000 8,000,000 16,000,000 32,000,000 50,000,000 (a) Código r24 (CUDA+MPI) Número de partículas en la simulación 10 100 1000 10000 1 10 100 1000 Número de procesos (gpu) Tiempo de ejecución (s) 1,000,000 2,000,000 4,000,000 8,000,000 16,000,000 32,000,000 50,000,000 (b) Código r27 (CUDA+MPI+OpenMP) Figura 5.5: Escalabilidad espacial de las implementaciones realizadas De las Tablas 5.6 y 5.7 también se observa un súbito incremento del tiempo empleado por las comunicaciones con MPI cuando se llega a usar gran parte de nodos del cluster MinoTauro. Las Gráficas 5.7 muestran mejor ese fenómeno. De todos estos datos se observa que: • El código con OpenMP se ejecuta en unos tiempos inferiores a la versión con sólo CUDA yMPI. • Parece existir algún problema con la versión r27 con OpenMP cuando se trabaja con 64 nodos (128 procesos) y con más de 30.000.000 de partículas. No se ha analizado la causa, pues el fenómeno de pérdida de escalabilidad en esas condiciones no se reproduce en la versión r24. Una hipótesis podría ser que cuando en ese cluster se ejecutan simulaciones con la mitad de los nodos globales de 61 Paralelización mediante GPU de la dinámica de electrones en plasma confinado magnéticamente de la paralelización de otros códigos científicos históricos, se procede a valorar cualitativamente la experiencia del proceso de implementación, las decisiones que se han tomado en dicha implementación y otras alternativas que se descartaron al inicio del trabajo. Como toda valoración cualitativa, ésta es subjetiva, por lo que a continuación se incluye se expresa como opinión formada a través de todo el proceso de realización del presente trabajo. Sin embargo, se considera muy interesante su expresión explícita en esta memoria porque es la parte quizá más interesante para el trabajo futuro, más que la mera perfección de las implementaciones concretas presentadas. A continuación se presentan diferentes características que se consideran de interés. 5.4.2.1 Facilidad de programación CUDA está basado en CyC++, y éstos ya son menos intuitivos de usar que Fortran para los científicos sin base específica de programación. CUDA representa otro paso más allá en la programación, pues se cambia el paradigma de programación. Por ello, la paralelización de código histórico en Fortran requiere: • que una parte del grupo de científicos se especialice seriamente en programación paralela, • que la tarea de programación sea encargada a un experto ingeniero informático. Esta segunda opción sería la preferible, pero siempre bajo el asesoramiento de algún científico, para poder clarificar las fórmulas, algoritmos, y en definitiva, dar a entender al programador qué se está haciendo. Por la razón anterior, CUDA es un lenguaje complicado de usar, aunque exprime a fondo las tarjetas gráficas de NVIDIA, muy potentes computacionalmente. Es interesante estar atento a otras alternativas, concretamente con el avance de OpenMP, versión 4, y sobre todo, de OpenACC. Ambos permiten el offload (envío de tareas de computación) a los aceleradores (principalmente, pero no exclusivamente, GPU, pero también a otros como Intel Xeon Phi o tarjetas gráficas de AMD) de una forma más intuitiva que CUDA. En este último trimestre de año se espera el lanzamiento de versiones ya finales de compiladores libres, como gcc, con soporte pleno a estas nuevas tecnologías. Se propone estar atentos a este hecho, y comparar la paralelización usando las GPU a través de OpenMP 4.0 y OpenACC. 5.4.2.2 Necesidad de la reformulación de los algoritmos El gran leitmotiv de Fortran es su propio originario acrónimo: FORmula TRANslation system, hecho por el cual es tan ampliamente utilizado en el mundo científico. Pero también es uno de sus principales problemas. Al ser tan sencillo transcribir las fórmulas, los métodos matemáticos y los algoritmos en ellas externamente visualizados, es muy fácil caer en la programación directa de la fórmula, sin pensar en aspectos importantes desde el punto de vista de ingeniería del software: • tener en cuenta los costes espaciales y temporales a la hora de implementar algoritmos (por ejemplo, clasificación en canastas, que puede hacerse con coste 1 por partícula, total O(N), siendo innecesario hacer previamente una ordenación total de los valores, con coste O(N·log(N)), y después recorrer otro bucle con coste O(N)), • implementar eficientemente según la arquitectura hardware y software disponible en el momento, (por ejemplo, estar actualizado en cuanto al avance de la tecnología de la paralelización, tales como OpenMP,OpenACC yCUDA), 68 Capítulo 5. Resultados 5.4. Valoración • amplio conocimiento de las funciones y biblotecas disponibles en el entorno (ejemplo, utilización de la función error, erf en Fortran,CyCUDA), en lugar de calcularla cada vez mediante métodos numéricos de integración (trapecio o Simpson). • actualización del programa base a las nuevas versiones de Fortran (por ejemplo, actualizar los programas desde Fortran 77 o Fortran 90 a Fortran 95 (bindings con C) o incluso Fortran 2003 (interoperabilidad con Cy por lo tanto, más fácil trabajar con CUDA, y objetos, incluso aunque sean como meras estructuras de datos, para ocultación de variables actualmente declaradas globalmente) 5.4.2.3 Depuración semántica de las implementaciones históricas Durante el trancurso del trabajo se detectaron lo que podían ser dos errores de programación que podrían influir de manera muy importante en los resultados que éste genera. A la fecha de hoy, el físico responsable de los códigos no ha confirmado ni negado esos informes de error enviados. Los códigos históricos utilizando tecnologías de paralelización basadas en clusters de CPU tienen cierta limitación en cuanto al tamaño y, por lo tanto, validez, de las simulaciones que con ellos se pueden realizar. Ésa es la razón que fundamenta este trabajo, de hecho, el estudio de viabilidad de su conversión para la utilización de tarjetas gráficas de computación. Dado ese limitado poder de simulación, es posible que en los códigos históricos, aunque estén en constante evolución y revisión, persistan errores de base, difíciles de detectar en código y también en simulaciones a pequeña escala, pero que en simulaciones mayores, provoquen resultados mucho más divergentes con la realidad que lo razonablemente esperado. Esto llevaría a falsas conclusiones y generación de dudas sobre el propio cuerpo teórico científico en que se basan. Es por ello que antes de cada paralelización, se revise completamente el código original, o incluso en mucho casos, que la implementación se realice desde cero para evitar de raíz este problema. Esto tiene la ventaja de que proporcionaría un código mucho más adaptado a las actuales tecnologías de programación paralela. 5.4.2.4 Selección de códigos a paralelizar Finalmente, y es una obviedad, dadas todas las dificultades arriba destacadas, es siempre conveniente realizar un adecuado estudio de selección de qué códigos históricos deben ser paralelizados con GPU. La programación de GPU para GPGPU es delicada, laboriosa y difícil de depurar, pero su tecnología (hardware y software), potencial y viabilidad económica, no paran de crecer. Aún así, no todos los problemas aprovechan esta nueva tecnología. Es una labor crucial tener a algún experto que asesore sobre la viabilidad de una implementación con CUDA u otras tecnologías emergentes (OpenACC,OpenMP 4), y de esta forma, dirigir los esfuerzos en paralelizar primero aquellas aplicaciones más interesantes pero que sean viables, tanto desde de el punto de vista tecnológico (rendimiento obtenido), como temporal (coste económico de la implementación). 69 Paralelización mediante GPU de la dinámica de electrones en plasma confinado magnéticamente (Página intencionadamente en blanco) 70 CA P Í T U L O 6 TRABAJO FUTURO La paralelización realizada de CUDA, con y sin OpenMP, sobre un código científico basado en Fortran con MPI, ha permitido comprobar que la aplicación del modelo SIMT a cierto tipo de códigos permite obtener unas prestaciones extraordinarias comparativamente con el código funcionando sobre CPU. También ha servido para detectar dificultades que deberán ser superadas para la aplicación de este procedimiento de paralelización a otros códigos científicos. Las líneas de futuro a partir de este trabajo se dividen en tres categorías: 1. Ampliación del trabajo sobre las implementaciones en CUDA presentadas, con objeto de mejorar su estabilidad, rendimiento y aplicabilidad: • Depuración del programa, mediante detección de posibles bugs, limpieza de código y otras tareas básicas para que el código quede en un estado de mejor mantenibilidad. • Mejorar la ubicación de los cerrojos de OMP, que permitan, por ejemplo, disminuir la granularidad, y con ello, intentar extraer todavía más prestaciones basadas en un mayor solapamiento de la computación con las comunicaciones y la escritura en disco o consola. • Mejorar programáticamente lasimplementaciones actuales, por ejemplo, evitando la copia, dentro de la GPU de los vectores de velocidad nuevos sobre los viejos, mediante el sistema de rotar el puntero para realizar dichos cambios con copia cero, o investigando sobre la mejora de prestaciones en el uso de los generadores de números aleatorios. • Investigar por qué razón no se permiten ejecuciones de simulaciones de más de 825.000 partículas por cada GPU, cuando éstas parecen disponer de suficiente capacidad de almacenamiento (6GB). Posiblemente a través de una mejor gestión de los diferentes tipos de memoria presentes en la tarjeta gráfica. • Reimplementar el código utilizando otras tecnologías de paralelización emergentes que también se aprovechan de los aceleradores gráficos, pero permiten también hacerlo sobre otros aceleradores, como sobre las tarjetas gráficas de AMD, o sobre todo, del acelerador Intel Xeon Phi. Ejemplos de estas tecnologías emergentes son OpenMP 4 y OpenACC, de los cuales se espera una implementación final libre en este último trimestre gracias a GNU. 71 Paralelización mediante GPU de la dinámica de electrones en plasma confinado magnéticamente 2. Aplicación de la paralelización a otros códigos científicos • Migrar otros códigos similares a CUDA, u otras tecnologías emergentes, para aprovechar el extremo poder computacional de las GPU. • Modificar las ecuaciones base de Langevin obtenidas inicialmente, a través de la progresiva eliminación de algunas de las aproximaciones efectuadas para la deducción del conjunto de ecuaciones implementadas en el presente código. Esto debería permite obtener simulaciones más cercanas a la realidad. • A partir de esta experiencia, y viendo que es posible, añadir complejidad al cuerpo teórico inicial, por ejemplo, mediante la consideración de la relatividad dentro del conjunto de ecuaciones, para de esta forma acercarse más si cabe al comportamiento que tendrá el plasma cuando esté realmente en reacción de fusión controlada de forma continuada, bien sea en un reactor de tipo tokamak o del tipo stellarator. 3. Utilización a fondo de la nueva capacidad de simulación en este campo, y estudiar nuevas condiciones, variando los parámetros de entrada, realizando las oportunas mínimas variaciones en el código si es necesario, de tal forma que se convierta en una herramienta más para la comprensión del comportamiento real del plasma de partículas confinado magnéticamente. 72 CA P Í T U L O 7 CONCLUSIONES En el presente Trabajo Fin de Máster se ha estudiado la viabilidad y metodología para poder migrar código histórico científico, generalmente en lenguaje Fortran, al nuevo paradigma de programación mediante aceleradores externos, concretamente, mediante tarjetas gráficas de NVIDIA a través de su tecnología propia privativa CUDA. Para ello se ha utilizado como caso de estudio un código en Fortran ya paralelizado mediante MPI, que permite la simulación de un plasma de partículas a alta temperatura confinado magnéticamente, que es la física básica en la que se basan los actuales intentos de fusión en caliente en reactores de tipo tokamak. Los resultados cuantitativos son plenamente satisfactorios: se han obtenido aceleraciones de hasta 700 en comparación con el código ejecutándose en un único núcleo de CPU, todo ello, obteniendo resultados correctos en la simulación y generando correctamente los ficheros intermedios estadísticos pertinentes. Esto representa grandes ventajas evidentes: permite la ejecución de simulaciones más grandes en un tiempo razonable, y disminuye el coste de explotación (número de horas) de utilización del cluster de computación de altas prestaciones. Asímismo, se ha redactado una memoria que puede servir de guía sobre los problemas, dificultades, técnicas, planteamientos y decisiones que podrían surgir durante la realización a gran escala de la migración de código científico a las nuevas tecnologías de GPGPU. A pesar de ello, la viabilidad de extender este caso de estudio a otros proyectos más ambiciosos no debe sólo ser informado por esta valoración muy positiva desde el punto de vista cuantitativo. Como ya se ha indicado al final de los Capítulos 4 y 5, hay otros aspectos cualitativos que sean tanto o más importantes a la hora de tomar la decisión: • La migración es técnicamente viable y se han demostrado cómo realizarla. No obstante, no es trivial. El grado de complejidad del código resultado es elevado, sobre todo por la mezcla de diferentes lenguajes y tecnologías. Es por ello que deben concluirse se precisa una persona cona la suficiente formación en computación paralela y en ingeniería del software para poder enfrentarse con garantías y exitosamente al reto. • A pesar de lo anterior, es muy importante destacar la importancia que tiene el significado físico de las fórmulas y algoritmos que implementan, así como el conocimiento del fenómeno físico estudiado, y 73 Paralelización mediante GPU de la dinámica de electrones en plasma confinado magnéticamente con ello, los resultados que a priori pueden obtenerse de la simulación. Es por ello que en el proceso de migración también debe estar presente un experto en la materia científica de que se trate, al objeto de guiar, comprobar y revisar que la implementación se ajuste a la teoría, esto es, efectúe una revisión semántica de los códigos. • Con relativamente poco esfuerzo, se pueden conseguir resultados notables. Mejorarlos es también posible, pero muchas veces requiere de una inversión en tiempo y tecnología más elevado. Pero todo depende del caso concreto. • Hay un aspecto que no se ha podido analizar, y es la existencia de compiladores específicos privativos de pago que permiten el uso de CUDA directamente sobre un dialecto o variación del lenguaje Fortran. Cabe la posibilidad de que su uso sea más sencillo que la solución propuesta, y que los requisitos del personal a cargo del proyecto sean también diferentes. • No todo código científico histórico puede ser fácilmente paralelizable con CUDA, o incluso puede que sea contraproducente hacerlo. Se requiere un buen conocimiento de CUDA y del programa a migrar, para poder determinar si el nuevo paradigma de paralelización masiva de la GPGPU es adecuado al proyecto en cuestión. Como ejemplos, los métodos recursivos deben pasarse a iterativos. Tampoco se tolera la programación dinámica, aunque es poco común en el mundo científico. Los algoritmos con gran dependencia entre datos pueden presentar problemas de acceso a memoria que disminuyan mucho el rendimiento de las GPU. • Es por ello que es necesaria una adecuada selección del código o algoritmo que va a ser migrado o implementado con CUDA, para extraer de esta tecnología todo lo que puede aportar. 74 GLOSARIO buffer Espacio de memoria intermedia, que puede ser software o hardware. bug En español, simplemente, error, aunque a veces se traduce literalmente como bicho, o gusano. Error, generalmente sutil, en el código que provoca el comportamiento anómalo del programa, incluso provocando su finalización abrupta. Especialmente peligrosos cuando sus manifestaciones son erráticas, dependiendo de alguna combinación determinada de condiciones para revelarse. Muy difíciles de depurar en computación paralela y distribuida. cache Memoria rápida intermedia de pequeña capacidad para acelerar el acceso a una memoria principal de gran capacidad más lenta. cluster Conjunto de ordenadores interconectados por una red de altas prestaciones para la ejecución conjunta de aplicaciones paralelas. computación grid Entorno de ejecución HPC formado por varios clusters localizados geográficamente en ubicaciones muy dispersas y administrados por diferentes organizaciones que, no obstante, comparten sus recursos entre sí para lograr una mayor ocupación y utilización de sus instalaciones, ofreciendo una interfaz de lanzamiento común. Las ejecuciones deben ser planificadas y coordinadas, para lo cual existen gestores de colas especializados que pueden funcionar en base a créditos. CUDA kernel Función cuyo código que se ejecuta en una GPU, y que por lo tanto de manera masivamente paralelizada. daemon Apicación servidora que se inicia usualmente en el arranque del sistema operativo y que se encarga de gestionar un determinado servicio para el sistema operativo o las aplicaciones. Aunque normalmente se suele traducir al español por «demonio», en este trabajo se opta por el término «duende», término que parece más correcto, pues demonio tiene connotaciones malignas que no cabe esperar de un servicio legítmo arrancado por el sistema operativo, mientras que el término «duende» evoca más al efecto «mágico» de aquello que se ocupa de hacer algo por nosotros sin casi percibir esa ayuda. disrupción Fenómeno brusco en un plasma en confinamiento magnético por el cual pierde sus condiciones de estabilidad y que provoca que en muy pocos milisegundos desaparezcan las condiciones necesarias para la fusión nuclear, apagándose el reactor como consecuencia. HPC High Performance Computing, computación de alto rendimiento. Consiste en la ejecución de programas reales en entornos masivos de producción dotados de la última tecnología especialmente diseñados y gestionados con el fin de obtener las máximas prestaciones computacionales.. mapear Tecnicismo utilizado para describir la acción de referenciar un objeto o zona de memoria desde otra ubicación, generalmente con el objeto de virtualizar o facilitar el acceso directo al objeto original. 75 Paralelización mediante GPU de la dinámica de electrones en plasma confinado magnéticamente middleware Pieza de software que ofrece servicios a otras aplicaciones más allá de los que ofrece el sistema operativo, tales como permitir la mútua interacción o el intercambio de información. Es comúnmente utilizado en sistemas distribuidos, en donde un único sistema operativo no puede centralizar y gestionar toda la intercomunicación. Otro uso es facilitar al programador el acceso a diferentes recursos de entrada/salida, ya sean locales o remotos, y de tipo software o hardware. name mangling Nomenclatura particular utilizada por cada compilador para cada lenguaje para la denominación de sus símbolos en los ficheros objeto .o. plasma Estado de la materia, similar al gas, en el que por el efecto de las altísimas temperaturas, la materia se haya completamente ionizada, dada la alta velocidad térmica de las partículas. profiling Investigación del comportamiento de un programa a través de la información recopilada desde el análisis dinámico del mismo. Existen programas específicos para tales tareas, siendo el más común gprof, dentro del software libre en Linux. NVIDIA pone a disposición dentro de su Toolkit, la utilidad nvprof. Aunque el término inglés está ampliamente difundido, puede traducirse por “perfilaje” o por “análisis del rendimiento”.. runtime Biblioteca de funciones que dan soporte para la ejecución de aplicaciones que utilicen un determinado lenguaje, en nuestro caso CUDA C. speedup Aceleración. Es una métrica de la mejora relativa en prestaciones del rendimiento para la ejecución de una misma tarea. En el ámbito de la computación paralela, hace referencia al cociente entre el tiempo de ejecución de la tarea secuencial y el tiempo de ejecución del programa paralelizado. stellarator Dispositivo donde se confina magnéticamente un plasma de partículas a alta temperatura con la finalidad de obtener la fusión nuclear controlada. Tiene la particularidad de que el campo magnético es tal que las órbitas de las partículas contenidas tiene una forma helicolidal. tokamak Dispositivo similar al stellarator, pero en el que el campo magnético confina al plasma en una conjunto de órbitas de forma toroidal. warp En español, urdimbre. Representa la agrupación de 32 ó 64 hilos de ejecución (según la arquitectura CUDA) que tienen un ámbito de ejecución común en una tarjeta gráfica de NVIDIA. Comparten recursos esenciales y escasos como registros y memoria cache, pero sobre todo, tienen un mismo flujo de ejecución, de tal forma que si éste es divergente para cada hilo, los hilos que no estén en esa bifurcación permanecen parados, y luego al revés, de tal forma que se acumulan los tiempos. Esto es así porque cada urdimbre sólo tiene activa una instrucción, la cual será o no ejecutada por el hilo que la componga. La urdimbre es, en la industria textil, el hilo longitudinal que se mantienen en tensión en un marco o tela, y que sirve de guía y estructura sobre la que se entretejen los hilos. wrapper Función envoltorio, que ofrece una interfaz conocida para acceder a otra función o secuencia de funciones. Su uso puede estar motivado por diferentes causas, tales como facilitar el acceso a otras funciones, simplificando la forma de invocarlas, o de esconder el funcionamiento interno de un software por temas de seguridad o de protección de la propiedad intelectual. 76 SIGLAS API Application Programming Interface, Interfaz de programación de aplicaciones. ARPACK ARnoldi PACKage, Paquete de (método de) Arnoldi. BLAS Basic Linear Algebra Subprograms, Subprogramas de álgebra lineal básica. CPU Central Processing Unit, Unidad central de procesamiento. CUDA Compute Unified Device Architecture, Arquitectura unificada de dispositivos de computación. FORTRAN Anteriormente FORmula TRAslation system, Sistema de traducción de fórmulas. En la actualidad ya no se interpreta de esta manera, se considera una palabra propia. GPFS General Parallel File System, Sistema de ficheros paralelo general (desarrollado por IBM. GPGPU General-Purpose computing on Graphics Processing Units, Computación de propósito general en unidades de procesamiento gráfico. GPU Graphics Processing Unit, Unidad de procesamiento gráfico. ITER International Thermonuclear Experimental Reactor, Reactor experimental termonuclear internacional. LAPACK Linear Algebra PACKage, Paquete de álgebra lineal. Lis Library of Iterative Solvers, Biblioteca de resolvedores iterativos. MPI Message Passing Interface, Interfaz de intercambio de mensajes. OpenCL Open Computing Language, Lenguaje abierto de computación. OpenMP Open Multi-Processing, Multiprocesamiento abierto. PETSc Portable, Extensible Toolkit for Scientific Computation, Conjunto portable y extensible de herramientas para la computación científica. ScaLAPACK Scalable Linear Algebra PACKage, Paquete de álgera lineal escalable. SIMT Single Instruction Multiple Threads, Instrucción única, múltiples hilos. Modelo de ejecución en paralelo, utilizada por las GPU de NVIDIA con capacidad CUDA, entre otras, en la que múltiples hilos ejecutan concurrentemente la misma instrucción sobre diferentes datos. SLEPc Scalable Library for Eigenvalue Problem Computations, Biblioteca escalable para computaciones sobre el problema de los valores propios o característicos. 77