Implementación paralela de algoritmos de control y optimización en CUDA

Equation Chapter 1 Section 1 Trabajo Fin de Grado Grado en Ingeniería en Tecnologías Industriales Implementación paralela de algoritmos de control y optimización en CUDA Autor: Eduardo Mayoral Briz Tutor: Daniel Rodríguez Ramírez Dpto. Ingeniería de Sistemas y Automática Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2018 iii Trabajo Fin de Grado Grado en Ingeniería en Tecnologías Industriales Implementación paralela de algoritmos de control y optimización en CUDA Autor: Eduardo Mayoral Briz Tutor: Daniel Rodríguez Ramírez Profesor titular de Universidad Dpto. de Ingeniería de Sistemas y Automática Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2018 v Trabajo Fin de Grado: Implementación paralela de algoritmos de control y optimización en CUDA Autor: Eduardo Mayoral Briz Tutor: Daniel Rodríguez Ramírez El tribunal nombrado para juzgar el Proyecto arriba indicado, compuesto por los siguientes miembros: Presidente: Vocales: Secretario: Acuerdan otorgarle la calificación de: Sevilla, 2018 El Secretario del Tribunal vii A mi familia A mis maestros ix Agradecimientos Me gustaría agradecer a toda mi familia el apoyo y cariño recibido durante estos años, sin el cual no podría haber llegado a escribir estas líneas. Gracias a todos mis amigos y compañeros de estudio, que han estado en los malos y en los buenos momentos durante estos duros años. También agradezco la ayuda a mi tutor, Daniel, quien me ha orientado durante el trabajo y me ha abierto las puertas a nuevos conocimientos. Eduardo Mayoral Briz Sevilla, 2018 xvii ÍNDICE DE FIGURAS Figura 1-1. Arquitectura de CPU y GPU. 1 Figura 1-2. Rejilla de bloques de hilos bidimensional. 2 Figura 1-3. Tipos de memoria en el device. 4 Figura 1-4. Ejecución de un programa CUDA C. 5 Figura 1-5. Distribución de bloques por SMs. 6 Figura 2-1. Resultados tras ejecutar DeviceQuery. 7 Figura 2-2. Configuración del proyecto: librerías. 8 Figura 2-3. Configuración del proyecto: compilación separada. 9 Figura 2-4. Configuración del proyecto: capacidad computacional. 10 Figura 3-1. Influencia del tamaño. 18 Figura 3-2. Influencia del número de hilos. 19 Figura 3-3. Influencia del número de bloques. 20 Figura 3-4. Resultado tras aplicar nvprof al kernel de factorización. 20 Figura 3-5. Ocupación. 21 Figura 3-6. Comparativa CPU vs GPU. 22 Figura 3-7. Comparativa CPU vs GPU para número de matrices bajo. 22 Figura 4-1. Paralelismo dinámico en CUDA. 26 1 1 DESCRIPCIÓN DE CUDA 1.1 Descripción general CUDA (Compute Unified Device Architecture) es una arquitectura de cálculo paralelo desarrollada por NVIDIA, con el fin de sacar el máximo provecho de las GPUs (Graphics Processing Unit) en labores más allá de las tareas gráficas que siempre han desempeñado. Desde su creación, ha proporcionado un incremento en el rendimiento de diferentes aplicaciones enfocadas a vídeo, biología, química computacional, machine learning o dinámica de fluidos, entre otros. Gracias a CUDA, en las tareas anteriormente mencionadas, y en muchas otras, se ha optado por utilizar una GPU como alternativa a la tradicional CPU, gracias a su gran capacidad de programación en paralelo. Mientras que en una CPU actual encontramos un numero de núcleos en torno a 4-8, en las GPUs más modernas podemos llegar a ver hasta 4000 núcleos. Estos, aunque tienen una capacidad computacional mucho menor, se benefician de su trabajo en paralelo para superar a la CPU en determinados tipos de tarea. El lenguaje de programación puede ser C/C++, Fortran, o incluso algunas directivas de lenguaje abierto como OpenACC. Dispondremos de dos partes a la hora de programar, el llamado host (anfitrión, asociado a la CPU) y el device (dispositivo, relacionado con la GPU). 1.2 Modelo de programación 1.2.1 Kernels Los kernels toman el papel de una función de C, pero son ejecutados por la GPU. Por lo general, se ejecutan desde el host, aunque también puede hacerse desde el device mediante el uso del paralelismo dinámico. Se ejecutan un determinado número de veces en paralelo, a diferencia de una función en la CPU, que solo se ejecuta una vez. Definimos un kernel anteponiéndole la declaración __global__, y lo llamamos especificando los parámetros de ejecución, el número de bloques e hilos que se utilizarán en su ejecución, entre <<<…>>>, de la forma <<< número de bloques, número de hilos >>>. Dichos parámetros pueden ser de tipo int o dim3. A continuación, se propone un kernel de ejemplo, donde se realiza la suma de dos vectores en paralelo: Ejemplo. Suma de dos vectores en paralelo. // Definicion del kernel __global__ void VecAdd(float* A, float* B, float* C) { int i = threadIdx.x; C[i] = A[i] + B[i]; } int main() { Figura 1-1. Arquitectura de CPU y GPU. Descripción de CUDA 2 ... // Invocacion del kernel con N hilos VecAdd << <1, N >> >(A, B, C); ... } 1.2.2 Jerarquía de hilos Cada uno de los kernels mencionados anteriormente se ejecutarán en tantos hilos como se desee. Es tarea del programador establecer el número de hilos, así como la organización y estructura de estos a la hora de lanzar el kernel. Los hilos (threads) forman parte de bloques (blocks), mientras que un grupo de estos constituye una rejilla (grid). Un hilo puede ser definido con un índice de hasta 3 dimensiones, al igual que cada uno de los bloques que contienen a estos hilos. Este índice se representa con threadIdx, y lo acompañamos de la dimensión correspondiente usando threadIdx.x (por ejemplo, para la dimensión x). De igual forma indexamos los bloques, con la variable blockIdx. Se pueden definir tantos índices como se deseen, y el kernel trabajará en todos y cada uno de los hilos indicados. Sin embargo, existen limitaciones en el número de hilos por bloque, así como en el número de bloques por rejilla, dependiendo de los recursos de la GPU que se trate. Por lo general, se tienen como máximo 1024 hilos por bloque. Aun así, existe la posibilidad de lanzar varios bloques de igual cantidad de hilos para alcanzar un número de hilos en ejecución mayor al límite mencionado anteriormente. A la hora de programar, podemos conocer este límite de hilos por bloque con la variable blockDim. Al igual que las anteriores, irá seguida de su dimensión (blockDim.x), ya que en cada dimensión los bloques tendrán diferentes limitaciones. Un kernel de ejemplo para la suma de dos matrices, usando dos índices, uno para filas y otro para columnas, es el siguiente: Figura 1-2. Rejilla de bloques de hilos bidimensional. 3 3 Implementación paralela de algoritmos de control y optimización en CUDA Ejemplo. Suma de dos matrices en paralelo con dos índices. // Definicion del kernel __global__ void MatAdd(float A[N][N], float B[N][N], float C[N][N]) { int i = threadIdx.x; int j = threadIdx.y; C[i][j] = A[i][j] + B[i][j]; } int main() { ... // Invocacion del kernel con un bloque de N*N hilos int numBlocks = 1; dim3 threadsPerBlock(N, N); MatAdd << <numBlocks, threadsPerBlock >> >(A, B, C); ... } Cada uno de los hilos se ejecutará de forma independiente, aunque pueden compartir datos de memoria si están en el mismo bloque. Esto es gracias a una memoria compartida, shared memory, de una latencia muy baja, que coordina los accesos de cada hilo a la memoria. Debemos hacer uso de la llamada __syncthreads(), que actúa como barrera para que los hilos se coordinen, y no se de la operación por finalizada hasta que acaben todos los hilos. 1.2.3 Jerarquía de memoria En CUDA, los hilos pueden hacer accesos a diferentes tipos de memoria. Cada hilo tiene una memoria local privada asociada, mientras que cada bloque de hilos tiene una memoria compartida que es visible ante todos los hilos que componen el bloque. Todos los hilos de este bloque tendrán acceso a dicha memoria compartida. Además, tienen acceso a los espacios de memoria de sólo lectura constant y texture. También tienen acceso a la memoria global de la GPU, la memoria asociada al device. Sin embargo, aunque está a disposición de cualquier hilo, los accesos son más lentos que si se usa la memoria compartida. El uso de dicha memoria shared puede darnos en ocasiones un rendimiento mayor que la global, aunque añadirá cierta complejidad al problema, llegando a casos en los que su uso es inviable, debido a su pequeño espacio. Para alojar espacio en la memoria global, disponemos de funciones como cudaMalloc() o cudaFree(), análogas a las ya conocidas malloc() o free() en C/C++ estándar. Si queremos transferir datos desde el host al device, debemos usar alguna función del tipo cudaMemcpy(). Esto es algo fundamental, ya que no podemos acceder directamente al contenido de la memoria global del device sin antes transferirla al host. Es interesante mencionar que algunas funciones como cudaMemcpy() son síncronas, es decir, no vuelven para seguir ejecutando instrucciones hasta que esta finalice. En muchas ocasiones, estaremos desaprovechando tiempo si las siguientes instrucciones en la cola de ejecución no tienen relación con esta y están paradas. Para ello se usan funciones asíncronas, como en este caso cudaMemcpyAsync(). Esta función hace la transferencia de memoria que se le indique, mientras ya ha dado paso a la siguiente instrucción. Sin embargo, el uso de estas funciones asíncronas tiene algunos inconvenientes. Por un lado, el uso de cualquier dato que dependa de una función que aún no ha terminado puede desencadenar errores. Por otro, el programa puede volverse más complejo, ya que en muchos casos es necesario el uso de objetos auxiliares, como los streams. Descripción de CUDA 4 1.2.4 Programación heterogénea El modelo de programación de CUDA asume que los hilos se ejecutan en un dispositivo físicamente separado (el llamado device). Los kernels lanzados serán ejecutados en la GPU, mientras que el resto de código C se ejecuta en la CPU. El programa escrito es el que gestiona los accesos a los diferentes tipos de memoria, así como las transferencias de datos, mediante el llamado CUDA Runtime (Figura 1-4). El programa escrito en C se ejecutará de forma secuencial. Así mismo, las órdenes correspondientes a la CPU, también lo harán de forma secuencial. Sin embargo, a la hora de ejecutar una función asociada a la GPU, como puede ser el lanzamiento de un kernel, está se ejecutará en paralelo y dará paso a la siguiente instrucción. Existen algunas funciones que, aunque se ejecutan en el device, son asíncronas. Es responsabilidad del programador evitar que las funciones que se ejecuten asíncronamente no usen datos aun no disponibles de funciones anteriores. De lo contrario, se podrían producir problemas y la forma de actuar del programa no estaría definida. 1.2.5 Capacidad computacional La capacidad computacional de un dispositivo, o compute capability, identifica a cada GPU, definiendo las características soportadas por cada una. Es usada por las aplicaciones en tiempo de ejecución para determinar que características de hardware e instrucciones están disponibles en la GPU que se trate. Se representa por un número X.Y, donde X refleja la arquitectura de la GPU: 7 para arquitecturas Volta, 6 para Pascal, 5 para Maxwell, 3 para Kepler, 2 para Fermi y 1 para Tesla. No debe ser confundido con la versión de CUDA. En ocasiones, muchas actualizaciones del Toolkit introducidas por NVIDIA tienen limitaciones de hardware. Mediante la capacidad computacional, establecen umbrales en los que crean grupos de GPUs que podrán o no disfrutar de dichas actualizaciones. Un claro ejemplo es el paralelismo dinámico, disponible únicamente a partir de la capacidad computacional 3.5. Figura 1-3. Tipos de memoria en el device. 5 5 Implementación paralela de algoritmos de control y optimización en CUDA 1.3 Arquitectura Una GPU NVIDIA está construida en torno a unos procesadores multihilo, llamados streaming multiprocessors (SMs). Cuando un programa en CUDA invoca un kernel, los bloques de un grid (red, rejilla) son enumerados y distribuidos a lo largo de los multiprocesadores con capacidad de ejecución disponible. Los hilos de un bloque de hilos se ejecutan concurrentemente en un multiprocesador, así como múltiples bloques de hilos son ejecutados en un multiprocesador. Cuando estos bloques terminan, se lanzan nuevos bloques en los multiprocesadores vacantes (Figura 1-5). Un multiprocesador está diseñado para ejecutar cientos de hilos concurrentemente. Para gestionar tal cantidad de hilos, se emplea una arquitectura llamada SIMT (Single-Instruction, Multiple-Thread). Las diferentes instrucciones son encoladas y, a diferencia de las CPUs, no se puede predecir en qué orden serán ejecutadas. Cada multiprocesador crea, ejecuta y gestiona los hilos en grupos de 32, llamados warps (urdimbres). Los hilos que componen un warp comienzan a la vez, en la misma dirección del programa, aunque cada hilo tiene su propio registro y contador de instrucciones, por lo que son libres para ejecutarse independientemente. Cuando un multiprocesador recibe uno o más bloques de hilos para ser ejecutados, los reparte en warps y hace una esquematización de la ejecución. La forma en que los bloques se reparten en warps es siempre la misma: cada warp contiene hilos consecutivos y de identificador creciente, comenzando por el hilo 0. Cada warp ejecuta una instrucción al mismo tiempo, por lo que la máxima eficiencia se alcanza cuando los 32 hilos del warp concuerdan en la orden de ejecución. En particular, cada multiprocesador tiene un conjunto de registros de 32 bits distribuidos a lo largo de los Figura 1-4. Ejecución de un programa CUDA C. Descripción de CUDA 6 warps, y una memoria compartida (shared memory) distribuida a lo largo de los bloques. El número de bloques e hilos que pueden alojarse y ser procesados a la vez en un multiprocesador para un kernel determinado, depende de la cantidad de registros y memoria compartida usados por el kernel y de la cantidad de registros y memoria compartida disponibles en el multiprocesador. También hay un límite de warps y bloques de hilos por cada multiprocesador. Todos estos límites, así como la cantidad de registros y memoria compartida, viene dado por la capacidad computacional del dispositivo. Si en un momento dado no existieran suficientes registros o memoria compartida disponible para procesar al menos un bloque, el kernel fallaría al ser lanzado. La ocupación es una relación establecida entre los registros y memoria disponibles, y los utilizados. Teóricamente, una alta ocupación hará que los multiprocesadores estén siempre trabajando, lo que ocultará mucho mejor la latencia, y aportará resultados más eficientes. El Toolkit de CUDA nos provee de una calculadora de ocupación teórica, asi como de otras herramientas para calcular la ocupación real, que serán tratadas más adelante. Figura 1-5. Distribución de bloques por SMs. 7 2 CONFIGURACIÓN DE CUDA EN WINDOWS 2.1 Instalación de CUDA CUDA está disponible para distintas plataformas, aunque nos centraremos en su instalación en Windows, sistema operativo donde se ha usado para este proyecto. Para usar CUDA en Windows debemos contar con los siguientes recursos básicos:  Una GPU que soporte CUDA. Casi cualquiera de los últimos años no debería presentar problemas de compatibilidad, más allá de las capacidades computacionales de cada una.  Un sistema operativo Windows compatible. Esto es variable según la versión de CUDA en cuestión. La versión más reciente, CUDA 9.2, podría ejecutarse sin problemas en Windows 10, Windows 8 y Windows 7.  Una versión de Microsoft Visual Studio compatible.  El Toolkit de NVIDIA CUDA. Se puede descargar de forma gratuita desde la web de desarrolladores de NVIDIA, tanto el más reciente como versiones anteriores. Una vez descargado el Toolkit, se procede a ejecutar el instalador. Este se encargará de instalar CUDA, actualizar drivers, librerías y códigos de ejemplo, entre otros. Una vez finalizado, debemos asegurarnos de que se ha instalado correctamente. Una forma de hacerlo es abrir una ventana de comandos, e introducir nvcc -V, comando que comprobará la versión que hemos instalado. Usaremos nvcc siempre que queramos usar el compilador de CUDA. En caso de obtener un error, puede deberse a una mala instalación o a una mala configuración del path de CUDA en las variables de entorno del sistema. Figura 2-1. Resultados tras ejecutar DeviceQuery. Paralelización de algoritmos: factorización de Cholesky 14 Con crea_term_indep, generamos los términos independientes de las matrices. Como es una función auxiliar, igual que la anterior, se generan en función de los coeficientes de la matriz factorizada (por eso se hace tras la factorización, simplemente por ahorrar coste computacional). La forma de generarlos es la siguiente: recorre los coeficientes de cada fila, y los va sumando. De esta forma, genera una serie de términos independientes tales que todas las soluciones del sistema serán 1. Esto facilita mucho la revisión y la comprobación de la factorización. Por último, con la función comprueba_soluciones comprobaremos que se ha factorizado y resuelto con éxito. Esta función es únicamente de uso auxiliar, y solo funcionará con éxito si previamente se ha usado crea_term_indep. Si es así, revisa que todas las soluciones valgan 1, y en ese caso imprime un mensaje diciendo que se ha factorizado con éxito. En caso de haber algún número distinto de 1, sale y avisa de que ha habido fallos. 3.3 Factorización en la CPU Como ya se mencionó antes, las matrices se almacenarán en punteros. Cada puntero hará referencia a todas las submatrices de cada tipo. Por tanto, basta con utilizar un índice para movernos dentro de cada matriz “grande” y localizar la submatriz deseada. Para medir los tiempos de ejecución, utilizaremos distintos métodos para CPU y GPU. En este caso, haremos uso de las funciones QueryPerformanceCounter y performancecounter_diff para establecer puntos de inicio y final, y calcular el tiempo transcurrido entre ambos puntos, respectivamente. Código. Medición de tiempos en rutinas CPU. double performancecounter_diff(LARGE_INTEGER *a, LARGE_INTEGER *b) { LARGE_INTEGER freq; QueryPerformanceFrequency(&freq); return (double)(a->QuadPart - b->QuadPart) / (double)freq.QuadPart; } void main(){ QueryPerformanceCounter(&iniTime); // primera llamada cholesky_cpu(h_m, h_Bfact, tam); // función que se ejecutará QueryPerformanceCounter(&endTime); // segunda llamada printf("Tiempo en crear %d submatrices de dimensiones %d x %d: %.16g milisegundos\n", num, tam, tam, 1000 * performancecounter_diff(&endTime, &iniTime)); } La función principal encargada de la factorización es cholesky_cpu. Los parámetros que recibe la función son dos matrices (*m y *Bfact) que a su vez contienen un total de num submatrices. En *m se guardan las matrices originales, mientras que en *Bfact se guardan las matrices tras la factorización. También recibe dim, que es el tamaño de cada submatriz. Esta función se encarga de realizar la factorización cada submatriz en bucle, en el que el índice idx_sub se encarga de recorrer todas las submatrices, desde 0 hasta num. Los pasos que sigue son poner a 0 los términos por encima de la diagonal en Bfact y ejecutar el algoritmo de la factorización. Por otro lado, la función resuelve_sist_cpu se encarga de resolver cada sistema de ecuaciones. Agrupa todas las soluciones en una matriz de dimensiones num_submats * filas_submats, es decir, devuelve un puntero a una matriz sols que contiene todas las soluciones para cada incógnita. Es una función prácticamente auxiliar, ya que sirve para comprobar las soluciones y verificar que la factorización se ha realizado con éxito. 15 Implementación paralela de algoritmos de control y optimización en CUDA 3.4 Factorización en la GPU Para la GPU se ha desarrollado un programa diferente al anterior, que nos da opción a usar múltiples GPUs. En cualquier caso, el programa es capaz de detectar el número de GPUs disponibles, para no ejecutarse si se le pide que lo haga usando más de las disponibles. Es necesario definir un tipo de estructura que abarque las variables con independencia para cada GPU, como la que aparece en el código a continuación. Código. Estructura y órdenes para ejecución en varias GPU. typedef struct{ int tams; //número de elementos de cada matriz int num_submats; //número de submatrices en esta GPU int fila_submats; //longitud de la fila en esta GPU int num_bloq; int num_hilos; float *sols; float *h_m; float *h_Bfact; float *d_m; float *d_Bfact; cudaStream_t stream; } TGPUplan; int detected_GPUs; // buscamos el numero de GPUs disponibles cudaError_t error = (cudaGetDeviceCount(&detected_GPUs)); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0);} if (NUM_GPUs > detected_GPUs){ printf("El numero de GPUs escogido (%d) es mayor al numero de GPUs instaladas (%d)\nSolo se trabajara con %d GPUs\nSaliendo...\n", NUM_GPUs, detected_GPUs, detected_GPUs); exit(EXIT_FAILURE);} TGPUplan plan[NUM_GPUs]; NUM_GPUs es un macro que es definido por el usuario, donde se indica el número de GPUs que se desea usar. Previamente a la creación de la estructura en cuestión, se verifica que el usuario no haya introducido un número mayor de GPUs al existente, mediante la función cudaGetDeviceCount. 3.4.1 Uso de múltiples GPUs En algunas ocasiones, podemos encontrarnos con un PC con más de una GPU instalada. Mediante ciertas funciones de CUDA seremos capaces de exprimir esta cualidad, haciendo que ambas trabajen a la vez, con la correspondiente mejora de rendimiento que en la mayoría de los casos supone. Para desempeñar el cálculo en varios dispositivos, es interesante recalcar el uso de la función cudaSetDevice(int disp), que nos permitirá elegir la GPU donde queremos que se desarrollen las funciones que se mencionen a continuación de esta. El parámetro disp será un entero que representa el número de la GPU, comenzando en 0. En general, a lo largo del programa se usará un bucle del estilo for (i = 0; i < NUM_GPUs; i++), dentro del cual escribiremos el código que queramos que se desarrolle en cada GPU. La única diferencia en cada bucle es que al comienzo habrá que llamar a cudaSetDevice(i), para seleccionar la GPU donde queremos ejecutar las órdenes. El uso de varios dispositivos, aparte de darnos un extra en rendimiento, nos añadirá un tanto de complejidad. Paralelización de algoritmos: factorización de Cholesky 16 Por ejemplo, a la hora de lanzar un kernel, tendremos que especificar al sistema en qué dispositivo queremos que se lleve a cabo la ejecución. Para ello, tendremos que hacer uso de streams, los cuales nos permitirán añadir un nivel más de concurrencia. Gracias a esto, asociando cada GPU a un stream diferente, haremos que cada una haga su trabajo por separado y no interfieran entre ellas. 3.4.2 Medición de tiempos Se utilizan funciones de la librería de CUDA que tienen un funcionamiento similar al mencionado en el apartado de la CPU. Estas funciones servirán siempre que queramos medir tiempos en CUDA. Para comenzar, debemos crear dos eventos, uno para comenzar la cuenta y otra para finalizar. Esto lo haremos con la función cudaEventCreate. Para definir los puntos o eventos de comienzo-final, utilizamos la función cudaEventRecord, que darán comienzo y final a la cuenta cuando se llamen. Se introducirá como argumento el evento que se desee, y previamente creado con cudaEventCreate. Estas funciones requieren de un stream, al igual que para el cálculo en múltiples GPUs. Sin embargo, en esta ocasión vale con introducir un 0 como argumento, usando así el stream por defecto. Tras llamar dos veces a cudaEventRecord, es necesario sincronizar los eventos con la función cudaEventSynchronize. Por último, con la función cudaEventElapsedTime almacenamos en la variable ms el valor del tiempo medido en milisegundos. 3.4.3 Transferencias de memoria Al ser el kernel ejecutado en la GPU, necesitamos transferir los datos desde el host al device, y viceversa. Para esto, aunque se suele usar la función cudaMemcpy, se usará la función cudaMemcpyAsync. La diferencia es que esta última es asíncrona respecto al host, y por tanto volverá antes de finalizar la copia, lo que implica un ahorro de tiempo. Otra diferencia es que es necesario crear un stream al que asociar la función. Este debe ser creado en cada estructura, ya que habrá un stream diferente para cada llamada de cudaMemcpyAsync y, por tanto, para cada estructura. Es por esto por lo que cada estructura, y por tanto cada GPU, estará asociada a un stream diferente. Código. Funciones asociadas a la transferencia de memoria síncrona y asíncrona. error = cudaMemcpy(h_m, d_m, tam*fila * sizeof(float), cudaMemcpyDeviceToHost); if (error != cudaSuccess) { … } error = cudaMemcpyAsync(plan[disp].d_m, plan[disp].h_m, plan[disp].tams * sizeof(float), cudaMemcpyHostToDevice, plan[disp].stream); if (error != cudaSuccess) { … } Es conveniente también mencionar que, para evitar problemas tras el uso de la función cudaMemcpyAsync, debemos utilizar la función cudaStreamSynchronize. Esta función establece un punto de interrupción donde se esperará que todos los streams hayan acabado con las tareas asociadas a ellos. Gracias a esto evitaremos que se trabaje con zonas de memoria que aún no ha sido transferidas por cudaMemcpyAsync y se produzcan problemas indeseados y difíciles de localizar. 3.4.4 Kernels usados 3.4.4.1 Kernel de la factorización: cholesky_gpu Tiene un trabajo similar a la función cholesky_cpu, pero con las particularidades de un kernel. La primera diferencia es que un kernel no trabaja en un bucle secuencial con un índice, sino que lo hace en paralelo en un determinado número de hilos y bloques elegidos por el usuario. El índice en este caso es int idx_sub = blockIdx.x * ThreadsPerBlock + threadIdx.x, y cada valor corresponderá a una submatriz. Para no causar problemas por desborde, todo el procesamiento que sigue a 17 Implementación paralela de algoritmos de control y optimización en CUDA continuación se realizará dentro de una sentencia if (idx_sub < num_submats) {…}. Los pasos que sigue son los mismos que en la función cholesky_cpu. La única diferencia es que se utilizan los valores num_submats y fila_submats, los cuales son diferentes a num y fila. Esto es debido a que en cada GPU puede haber un numero diferente de submatrices, en caso de que el total de submatrices no fuera múltiplo del número de GPUs utilizadas. Por ello, el kernel recibe como parámetros dichos valores independientes para cada estructura. Otra particularidad del kernel es la necesidad de usar la función __syncthreads() , que establece un punto de interrupción hasta que lleguen todos los warps, para evitar problemas indeseados de sincronización entre hilos. El lanzamiento del kernel se realiza de la siguiente forma: cholesky_gpu <<< BlocksPerGrid, ThreadsPerBlock, 0, plan[disp].stream >>>(plan[disp].d_m, plan[disp].d_Bfact, plan[disp].num_submats, plan[disp].fila_submats , tam). Entre comillas se definen las características de lanzamiento: número de bloques, de hilos, y stream asociado al lanzamiento. El stream no sería necesario si se hubiera usado la función cudaMemcpy en vez de cudaMemcpyAsync. Entre paréntesis se incluyen los argumentos, como una función normal. 3.4.4.2 Kernel de la resolución de sistemas: resuelve_sist_gpu Resuelve los sistemas en paralelo en la GPU, como alternativa a la CPU. El prototipo es igual, si obviamos las diferencias que existen entre una llamada a una función y el lanzamiento de un kernel: resuelve_sist_gpu. El indexado es igual al kernel de la factorización cholesky_gpu, ya que estamos haciendo un tratamiento a cada una de las submatrices que componen el conjunto de submatrices. La única diferencia está en el cuerpo, el cual sí que es prácticamente igual al de la función resuelve_sist_cpu. El kernel se lanza con la siguiente sentencia: resuelve_sist_gpu <<< BlocksPerGrid, ThreadsPerBlock, 0, plan[disp].stream >>>(plan[disp].d_Bfact, plan[disp].d_sols, plan[disp].num_submats, plan[disp].fila_submats, tam). Como es lógico, para lanzar el kernel con éxito habrá que haber reservado y transferido la memoria necesaria, tal y como se explica en apartados anteriores. 3.4.5 Análisis de las variables Como se ha mencionado en puntos anteriores, a la hora de ejecutar nuestro programa en la GPU existen determinadas variables. Es evidente que cualquier modificación de ellas hará que obtengamos tiempos de factorización diferente. Por tanto, se hará un breve análisis para intentar estudiar la influencia de cada uno de los parámetros más destacables en nuestro kernel de factorización. Para estos análisis se ha utilizado una CPU Intel Core i7 a 2,5GHz, y diferentes dispositivos gráficos. Por un lado, se ha usado una GPU NVIDIA GeForce GTX 870M, con 3GB de VRAM, 7 SMs (multiprocesadores) y 1344 CUDA Cores. Por otro lado, se han hecho análisis con 2 dispositivos NVIDIA GeForce GTX TITAN, de 6GB de VRAM, 15 SMs y 2880 CUDA Cores cada uno. Como cada dispositivo de este tipo presenta 2 GPUs, en algunos casos se ha llegado a hacer un uso a la vez de hasta 4 GPUs. 3.4.5.1 Influencia del tamaño de cada submatriz Partimos de un conjunto de submatrices cuadradas cuya dimensión es una de las variables de nuestro problema. Por el estilo de nuestro kernel, en el cual cada hilo hace un procesamiento de factorización individual de cada submatriz, podemos esperar que aumentar el tamaño de cada una de las submatrices implique un aumento del tiempo de factorización en una determinada proporción. El resultado de la comparación en una GTX 870M se puede apreciar en la Figura 3-1. Paralelización de algoritmos: factorización de Cholesky 18 Tal y como se esperaba, encontramos un comportamiento casi proporcional, en el que aumenta el tiempo de factorización considerablemente. Esto es debido a que cada uno de los hilos que procesa una submatriz tiene una capacidad de computación muy baja, si lo comparamos con una CPU. Es decir, hacer que un hilo tenga mucha carga computacional nos puede llevar a tiempos demasiado altos, por lo que lo ideal para tamaños más altos sería repartir la carga de un kernel a varios hilos. 3.4.5.2 Influencia del número de hilos de lanzamiento Como se mencionó en el apartado de la arquitectura de CUDA, la información es procesada por los multiprocesadores en grupos de 32 hilos, llamados warps. Esto es algo clave si queremos observar como influye la variación del número de hilos en los tiempos obtenidos. Para comenzar, hemos obtenido una serie de tiempos en función del número de hilos, quedando de la siguiente forma (Tabla 3-1): Tabla 3-1. Influencia del numero de hilos en la factorización. Hilos Tiempo por matriz (µs) 20 7.69 32 4.81 40 4.92 50 3.31 64 3.23 100 4.0 128 2.89 200 3.88 256 2.81 400 3.78 512 2.63 1000 3.33 1024 3.24 Lo que nos decía la teoría de la arquitectura se cumple con creces. Vemos que, si vamos alternando el número 49,21 33,1 29,49 28,74 26,71 33,14 678,43 703,3 667,4 659,1 646,7 691,8 0 100 200 300 400 500 600 700 800 32 / 320 64 / 160 128 / 80 256 / 40 512 / 20 1024 / 10 Tiempo (ms) Hilos / Bloque Influencia del tamaño Tamaño 20 Tamaño 50 Figura 3-1. Influencia del tamaño. 19 Implementación paralela de algoritmos de control y optimización en CUDA de hilos, conseguimos tiempos mejores para aquellos números múltiplos de 32, el tamaño de un warp. Como ya se dijo, esto es debido a que los multiprocesadores trabajan en grupos de 32 hilos, y el hecho de hacerlo con un número de hilos que no hagan warps completos hará que algunos SMs se ejecuten de forma ineficiente. Por ello, a la hora de ejecutar nuestro kernel, o alguno que tenga una estructura similar y desarrolle tareas diferentes, la elección óptima estará entre los valores marcados en rojo en la Figura 3-2. Sin embargo, esta tendencia de mejor eficiencia a mayor número de hilos no se puede generalizar. En otro problema, y con otra GPU, podemos encontrar mejores resultados con 1024 hilos. 3.4.5.3 Influencia del número de bloques Cada bloque se irá asignando a un SM, que lo procesará. Teóricamente, mientras más warps maneje cada SM, tendrá más facilidad para ocultar la latencia, y con ello conseguir tiempos mejores. Esto es debido a que, si un warp determinado provoca un bloqueo, el SM tiene que elegir otro. Por tanto, mientras más haya, tendrá más elecciones posibles, y con ello estará menos tiempo bloqueado. Intentaremos hacer un análisis viendo si utilizar un número de bloques distinto de un múltiplo de 7 puede influir. Para ello, fijaremos el número de hilos utilizados en 1024, e iremos variando el número de bloques en el lanzamiento, para ver las variaciones que sufre (Figura 3-3). Vemos que, al llegar a 7 bloques, obtenemos un tiempo mejor que si pasamos a 8. Esto es debido a que cuando elegimos 7, estamos ocupando todos los SMs a la vez. Sin embargo, al pasar a 8, habrá un instante en el que solo un SM esté trabajando, y el kernel no finaliza hasta que este acabe. Por ello encontramos ahí un incremento de tiempo. Esta teoría no se mantiene si nos vamos a los siguientes múltiplos de 7. Aún así, puede tener una explicación bastante sencilla, y es el orden en que se van asginando los bloques a los SMs. Que se asignen 14 bloques, no implica que esto vaya a ser en dos “tandas” de 7. Cuando observamos que en 15 se mantienen los tiempos obtenidos en 14, corroboramos esto. Y es que, si ha quedado algún SM vacante, se le ha asignado el bloque sobrante. Es decir, se han conseguido procesar los 15 bloques en sólo 2 “tandas”. Esto toma más sentido si vemos la tendencia de la curva obtenida hacia un número infinito de bloques: a medida que hay más bloques, hay más posibilidades de ocupar SMs vacantes, y con ello mejorar los tiempos. Esto conlleva una reducción de la latencia, que se mencionaba en el apartado de la arquitectura. Esta curva, obtenida de la misma forma, pero en un alto número de bloques, debe ser prácticamente horizontal. 2 3 4 5 6 7 8 20 32 40 50 64 100 128 200 256 400 512 1000 1024 Hilos Tiempo (µs) Influencia de hilos Figura 3-2. Influencia del número de hilos. Paralelización de algoritmos: factorización de Cholesky 20 3.4.5.4 Ocupación Como ya sabemos, cada dispositivo tiene una serie de límites físicos a la hora de almacenar registros y memoria en cada SM. A su vez, según los parámetros que establezcamos nosotros, entre ellos bloques e hilos, conseguiremos que cada SM almacene un determinado número de registros y memoria compartida. La relación porcentual entre los registros y memoria usadas en cada momento, respecto a los límites establecidos por el dispositivo, representa la ocupación. Teóricamente, una ocupación alta es algo positivo, ya que ocultamos bien la latencia y mejoran los tiempos. Por el contrario, una baja nos debe dar peores resultados. Sin embargo, en la práctica esto no siempre es así, y muchas veces los tiempos óptimos se consiguen con una ocupación media-baja. Para el estudio de la ocupación en nuestro kernel hemos usado dos herramientas. Estas nos permiten calcular la ocupación teórica y la ocupación real según el número de hilos usados en nuestro kernel. La teórica se calcula mediante una hoja Excel proporcionada por NVIDIA, descargable desde la web oficial, llamada CUDA_Ocuppancy_calculator.xls. Por otro lado, para la real, utilizamos la herramienta nvprof en la línea de comandos del símbolo del sistema de Windows. Esta herramienta es un profiler (perfilador), que nos proporciona información sobre nuestro archivo, como por ejemplo sobre la ejecución del kernel, transferencias de memoria o configuración de la memoria. En concreto, podemos elegir que comprobar añadiendo el flag –-metrics, y los parámetros 0,05 0,055 0,06 0,065 0,07 0,075 0,08 0,085 0,09 0,095 0,1 0 Bloques 714 21 28 35 Tiempo (ms) Bloques Influencia del número de bloques Figura 3-3. Influencia del número de bloques. Figura 3-4. Resultado tras aplicar nvprof al kernel de factorización. 21 Implementación paralela de algoritmos de control y optimización en CUDA deseados. El resultado de aplicar esta herramienta al kernel se puede consultar en la Figura 3-4. Reuniendo la información necesaria de ambas ocupaciones para la ejecución del programa con diferentes números de hilos, se obtiene la gráfica que podemos ver en la siguiente figura (Figura 3-5): Como podemos ver, no siempre una alta ocupación implica mejores tiempos, ni una baja ocupación malos. Es algo bastante complejo como generalizarlo con un simple caso, ya que es algo muy variable según el tipo de problema y dispositivos que estemos usando. 3.5 Comparativa entre CPU y GPU Vistas las características del algoritmo tanto en la CPU como en la GPU, y tras haber profundizado un poco en el funcionamiento del kernel, procedamos a la comparación de ambos métodos. En la mayoría de los casos, elaborar un algoritmo en código C estándar es mucho más sencillo que en CUDA. Muchas menos líneas de código, menos complejidad y más soporte hacen que la programación en la CPU sea nuestra primera opción. Cuando nos decantamos por el uso de CUDA y la programación en GPU, es porque esto nos va a dar una optimización y rendimiento mayor que la CPU. Sin embargo, pueden existir casos en los que, tras horas y horas de programación en CUDA, obtengamos resultados incluso peores que con un simple programa en C. En este punto se va a hacer una comparativa entre ambos métodos, y se intentará explicar cuando es conveniente usar cada uno, así como los beneficios y contras que nos aportan. Cuando hacemos una ejecución de ambas formas, intentando hacer que los códigos y parámetros de ejecución sean lo más parecidos posibles, obtenemos los resultados esperados: la factorización en la GPU da resultados mucho mejores a partir de un determinado umbral. En este caso, para submatrices de tamaño 20, a partir de 1.000 submatrices comienza a compensar el uso de la GPU frente a la CPU (Figura 3-6). Este umbral será variable en función del tipo de kernel que tratemos. Es posible incluso encontrar funciones que sean totalmente inviables en CPU o GPU, ya que del modo contrario siempre da mejor resultado. 0,055 0,06 0,065 0,07 0,075 0,08 0,085 0,09 10% 20% 30% 40% 50% 60% 70% 80% 64 128 192 256 320 384 448 512 576 640 Hilos Tiempo/mat (ms) Ocupación (%) Ocupación y tiempo de factorización según número de hilos Ocupacion EXCEL (%)(teórico) Ocupacion NVPROF (%)(real) Tiempo/mat (ms) Figura 3-5. Ocupación. Paralelización de algoritmos: factorización de Cholesky 22 En nuestro caso, este límite será distinto según el número o tamaño de las matrices, como es lógico, y será la aplicación la que marque qué es mejor usar, si CPU o GPU. Por ejemplo, para factorizar 50 matrices, compensaría usar CPU, ya que los costes de overhead en la GPU superan a los costes de cálculo de la CPU. Estos costes de overhead abarcan el lanzamiento del kernel, transferencias de memoria, y creación y destrucción de ciertos objetos que no serían necesarios en una CPU, entre otras cosas. Esto se ve de forma clara si hacemos zoom en la zona baja de la gráfica anterior (Figura 3-7). Podremos observar como es el comportamiento de la factorización para un número de matrices bajo, y como en este caso sí que no compensaría desarrollar y usar el kernel, ya que hay un nivel bajo de paralelización, y un código secuencial en CPU sería incluso más efectivo. 0 100 200 300 400 500 600 700 800 900 100 1.000 10.240 20.480 40.960 100.000 Tiempo (ms) Número de matrices CPU vs. GPU CPU GPU GTX 870M GPU TITAN Figura 3-6. Comparativa CPU vs GPU. 0 5 10 15 20 25 30 35 40 45 50 10 100 1.000 10.000 Tiempo (ms) Número de matrices CPU vs. GPU (num. mat. bajo) CPU GPU GTX 870M GPU TITAN Figura 3-7. Comparativa CPU vs GPU para número de matrices bajo. 23 Implementación paralela de algoritmos de control y optimización en CUDA En esta última figura, vemos que cuando trabajamos con pocas matrices tenemos varios escenarios diferentes. Para un número muy bajo (<100), la opción de la CPU es la que obtiene los mejores tiempos. Para tratar menos de 1000 matrices aproximadamente, y suponiendo que nuestra GPU fuera la GTX 870M, la opción de la CPU seguiría siendo más viable. Sin embargo, para un número mayor a partir de ahí, la opción de la GPU es la que predomina, sea de gama media o alta. La CPU toma tiempos del orden de 20 veces más, por lo que sería difícil encontrar algún caso en este contexto en el que mereciese la pena usarla. Bagging en ajuste por mínimos cuadrados 30 30 31 Implementación paralela de algoritmos de control y optimización en CUDA 5 CONCLUSIONES Y LÍNEAS DE CONTINUACIÓN 5.1 Conclusiones Gracias a CUDA, hoy en día tenemos las herramientas necesarias para llevar a cabo la paralelización del algoritmo que queramos. Es evidente que el uso de estas va a ser mayor con el paso de los años, sobre todo en ramas de conocimientos aún por descubrir. La capacidad de hacer una computación en paralelo supone algo nuevo con respecto a lo que conocíamos hasta hace unos años. Sin embargo, por sencilla que sea la cuestión, siempre estará presente la duda de cuál será la forma más eficiente de lograr nuestros objetivos. Esto es algo que, a nivel básico, se ha intentado describir a lo largo del presente trabajo. Y es la propia aplicación que tengamos entre manos, la que definirá si es conveniente usar una CPU, una GPU, o ambas a la vez. Aquí se ha trabajado con una determinada función, que bien puede ser secuencial o incluso única, a la que se le da cierto nivel de paralelismo, como puede ser la factorización de Cholesky para la resolución de sistemas de ecuaciones. A partir de ahí, se exponen diferentes casos y situaciones en los que descubrimos que todos los métodos usados son útiles. Lo único necesario es saber en qué momento y de qué forma usarlos. Para la factorización de Cholesky hemos encontrado que la CPU nos puede ayudar para factorizar un número bajo de matrices. Si quisiéramos hacerlo a miles de ellas, tendríamos un grave problema en cuanto a tiempos de ejecución. Esto se soluciona con la programación en paralelo mediante CUDA, llevando a cabo una operación para miles en el mismo tiempo que una CPU convencional se toma en hacerlo en apenas 100. Por otro lado, también hemos descubierto la gran influencia que tiene conocer la arquitectura interior de CUDA. Desarrollar un programa en paralelo requiere de habilidad y conocimiento para distribuir los recursos, y esto resulta imposible si no conocemos como se comportan estos. Además, hemos verificado como las suposiciones teóricas en cuanto a la estructura de dichos recursos son ciertas para la aplicación que le hemos dado. Hemos complementado todo este estudio con el máximo uso posible de diferentes alternativas, como pueden ser las librerías de NVIDIA (cuBLAS, cuRAND), a la vez que también hay partes donde se desarrollan estas funciones con un código propio. Lo mismo sucede con el paralelismo dinámico, el cual se ha comprobado no ser lo suficientemente eficiente para las labores en las que lo necesitamos. Sin embargo, se ha hecho uso de él ya que se considera que, si bien en esta aplicación no es lo óptimo, en muchas otras sí que lo puede ser. Aporta ventajas que no pueden ser alcanzadas si no es con su uso. También se ha desarrollado un caso particular, el cual encaja a la perfección con la forma de los kernel que se han desarrollado en este trabajo. El bagging para la identificación por mínimos cuadrados nos exige dar un nivel de paralelismo más para conseguir los objetivos, que influye de manera positiva en términos estadísticos. En este caso, no queda comprobado que el uso del bagging sea algo necesario para la estimación por mínimos cuadrados de una recta. Esto es debido a que dicho estimador trabaja de forma óptima en el caso lineal, por lo que el bagging no mejora con creces los resultados que se obtendrían con un solo gran conjunto. Sin embargo, esto no quiere decir que el bagging no sea útil para lo que se propone. Será necesario verificar si esto nos proporciona mejoras en sistemas de grados mayores, más allá de una estimación lineal. Seguramente ahí, donde el uso del estimador de mínimos cuadrados ya no está tan optimizado, el bagging nos aporte sus ventajas estadísticas. 5.2 Líneas de continuación El carácter de este trabajo deja muy bien definidas las líneas de continuación. Por un lado, es necesario llevar todo el conocimiento de la arquitectura CUDA a un nivel más, desarrollando funciones en las que se busque Conclusiones y líneas de continuación 32 32 mucha más optimización y eficiencia, más allá de sentar las bases de este lenguaje de programación y dar a conocer los aspectos básicos. Con ello, también se debería intentar incorporar una concurrencia GPU-CPU, en la que podamos observar que ocurre si combinamos las ventajas de ambos métodos. Esto mismo también es aplicable al bagging. Sería oportuno implementar mejoras en el código, que permitan llevar la estimación por mínimos cuadrados a sistemas más complejos, más allá de lo lineal. A su vez, también se podría trabajar más en el nivel estructural de este, haciendo que se reparta mejor la computación, y no quede toda la carga de cada conjunto en hilos independientes, ya que estos tienen muy poca capacidad. Todo lo mencionado anteriormente no sería posible si no se hubieran sentado las bases aquí escritas, ya que las conclusiones obtenidas pueden ser el origen de otros proyectos más complejos. 33 Implementación paralela de algoritmos de control y optimización en CUDA ANEXOS A. Código de la factorización de Cholesky en CPU/GPU #include <stdio.h> #include <stdlib.h> #include <math.h> #include <cuda_runtime.h> #include <device_launch_parameters.h> #include <time.h> #include <Windows.h> #include "cublas_v2.h" #define BlocksPerGrid 5 // valor maximo: 65535 por dimension #define ThreadsPerBlock 30 // valor maximo: 1024 (x,y) / 64 (z) #define tam 200 // tamaño de cada submatriz #define num (BlocksPerGrid * ThreadsPerBlock) // ---PARA GPU--- numero de submatrices. se genera como combinación de hilos y bloques para evitar problemas de indexado al lanzar el kernel //#define num 40960 // ---PARA CPU--- numero de submatrices. #define fila (tam * num + num) // longitud de una fila en una matriz grande (que contiene muchas submatrices) double performancecounter_diff(LARGE_INTEGER *a, LARGE_INTEGER *b) { /* Funcion para medir tiempos de ejecucion durante la factorizacion en la CPU. */ LARGE_INTEGER freq; QueryPerformanceFrequency(&freq); return (double)(a->QuadPart - b->QuadPart) / (double)freq.QuadPart; } int * crea_mat_sdp(int *randmat, int *randmattransp, int *sdpmat, int dim) { int i, j, k; float sum = 0; // Generación de matriz triangular inferior aleatoria 'randmat' // NOTA: términos de la diagonal >0 para evitar indeterminaciones for (i = 0; i < dim; i++) for (j = 0; j < dim; j++) if (i >= j) *(randmat + i*dim + j) = 1 + (rand() % 5); else *(randmat + i*dim + j) = 0; // Traspuesta de 'randmat' for (i = 0; i < dim; i++) for (j = 0; j < dim; j++) *(randmattransp + i*dim + j) = *(randmat + j*dim + i); // Matriz producto (Matriz definida positiva) for (i = 0; i < dim; i++) { for (j = 0; j < dim; j++) { for (k = 0; k < dim; k++) { sum += (*(randmat + dim*i + k)) * (*(randmattransp + k*dim + j)); } *(sdpmat + i*dim + j) = sum; Anexos 34 34 sum = 0; } } return sdpmat; } void crea_term_indep(float *h_Bfact, int *term, int dim) { int i, j, k; for (k = 0; k < num; k++) { for (i = 0; i < tam; i++) { *(term + i) = 0; for (j = 0; j <= i; j++) { *(term + i) = *(term + i) + *(h_Bfact + i*fila + k*dim + k + j); } *(h_Bfact + i*fila + k*dim + k + tam) = *(term + i); } } } float * resuelve_sistema(float *h_Bfact, float *sols, int dim) { // resolucion de los sistemas de ecuaciones // conjunto de soluciones de tamaño tam*num int i, j, k; float sum; for (k = 0; k < num; k++) { for (i = 0; i < dim; i++) { sum = 0; for (j = 0; j < i; j++) { sum = sum + (*(sols + j*num + k) * *(h_Bfact + i*fila + k*dim + k + j)); } *(sols + i*num + k) = (*(h_Bfact + i*fila + k*dim + k + dim) - sum) / *(h_Bfact + i*fila + k*dim + k + i); } } return sols; } void comprueba_soluciones(float *sols) { // Comprueba que todas las soluciones del sistema valen 1. // Este comprobador solo es válido si se usa la funcion crea_term_indep para... // ... generar los terminos independientes de las matrices. int i, j, aux = 0; for (i = 0; i < tam; i++) { for (j = 0; j < num; j++) { if (*(sols + i*num + j) != 1) { printf("ERROR: Error al resolver el sistema\n"); aux = 0; break; } else { aux = 1; } 35 Implementación paralela de algoritmos de control y optimización en CUDA } } if (aux == 1) printf("Exito al resolver el sistema!\n"); } __global__ void cholesky_gpu(float *m, float *Bfact, int dim) { // -- m es una matriz de tamaño 'tam*tam*num', que contiene a 'num' submatrices // -- Bfact es una matriz de tamaño 'tam*tam*num', que contiene a 'num' submatrices factorizadas // -- dim es la dimensión de cada submatriz (definida cte. por 'tam') int i = 0, j = 0, k = 0; double s = 0, t = 0; // idx_sub es el número de la submatriz que queremos factorizar [0, (num-1)] // es el elemento para indexar la paralelización int idx_sub = blockIdx.x * ThreadsPerBlock + threadIdx.x; if (idx_sub < num) { // para poner a 0 los términos por encima de la diagonal for (i = 0; i < dim; i++) { for (j = 0; j < dim + 1; j++) { if (j > i) *(Bfact + i*fila + idx_sub*dim + j + idx_sub) = 0; if (j == dim) *(Bfact + i*fila + idx_sub*dim + j + idx_sub) = *(m + i*fila + idx_sub*dim + j + idx_sub);}} __syncthreads(); // algoritmo de la factorización for (i = 0; i < dim; i++) { for (j = 0; j < (i + 1); j++){ if (i == j) { s = 0; for (k = 0; k < j; k++) { s += (*(Bfact + i*fila + idx_sub*dim + k + idx_sub)) * (*(Bfact + i*fila + idx_sub*dim + k + idx_sub)); } *(Bfact + i*fila + idx_sub*dim + j + idx_sub) = sqrt(*(m + i*fila + idx_sub*dim + j + idx_sub) - s); } if (i != j) { t = 0; for (k = 0; k < j; k++) { t += (*(Bfact + i*fila + idx_sub*dim + k + idx_sub)) * (*(Bfact + j*fila + idx_sub*dim + k + idx_sub)); } *(Bfact + i*fila + idx_sub*dim + j + idx_sub) = (*(m + i*fila + idx_sub*dim + j + idx_sub) - t) / *(Bfact + j*fila + idx_sub*dim + j + idx_sub); } } } __syncthreads(); } Anexos 36 36 } void cholesky_cpu(float *m, float *Bfact, int dim) { int i = 0, j = 0, k = 0, idx_sub; double s = 0, t = 0; // idx_sub es el número de la submatriz que queremos factorizar [0, (num-1)] for (idx_sub = 0; idx_sub < num; idx_sub++) { // para poner a 0 los términos por encima de la diagonal for (i = 0; i < dim; i++) { for (j = 0; j < dim + 1; j++) { if (j > i) *(Bfact + i*fila + idx_sub*dim + j + idx_sub) = 0; if (j == dim) *(Bfact + i*fila + idx_sub*dim + j + idx_sub) = *(m + i*fila + idx_sub*dim + j + idx_sub); } } // algoritmo de la factorización for (i = 0; i < dim; i++) { for (j = 0; j < (i + 1); j++) { if (i == j) { s = 0; for (k = 0; k < j; k++) { s += (*(Bfact + i*fila + idx_sub*dim + k + idx_sub)) * (*(Bfact + i*fila + idx_sub*dim + k + idx_sub)); } *(Bfact + i*fila + idx_sub*dim + j + idx_sub) = sqrt(*(m + i*fila + idx_sub*dim + j + idx_sub) - s); } if (i != j) { t = 0; for (k = 0; k < j; k++) { t += (*(Bfact + i*fila + idx_sub*dim + k + idx_sub)) * (*(Bfact + j*fila + idx_sub*dim + k + idx_sub)); } *(Bfact + i*fila + idx_sub*dim + j + idx_sub) = (*(m + i*fila + idx_sub*dim + j + idx_sub) - t) / *(Bfact + j*fila + idx_sub*dim + j + idx_sub); } } } } } void main(void) { int i, j, k; //clock_t t_ini, t_fin; // tipos para cálculo de tiempos en CPU //LARGE_INTEGER iniTime, endTime; // tipos para cálculo de tiempos en CPU //double t_crear, t_fact; cudaError_t error; cudaEvent_t start, stop; // tipos para cálculo de tiempos en GPU float tiempo; 37 Implementación paralela de algoritmos de control y optimización en CUDA // declaracion y reserva de memoria para matrices en 'host' y 'device' float *h_m = (float *)malloc(tam*fila * sizeof(float)), *h_Bfact = (float *)malloc(tam*fila * sizeof(float)), *sols = (float *)malloc(tam*num * sizeof(float)); float *d_m, *d_Bfact; error = cudaMalloc((void **)&d_m, tam*fila * sizeof(float)); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } error = cudaMalloc((void **)&d_Bfact, tam*fila * sizeof(float)); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } int *A = (int *)malloc(tam*tam * sizeof(int)), *B = (int *)malloc(tam * sizeof(int)), *randmat = (int *)malloc(sizeof(int)*tam*tam), *randmattransp = (int *)malloc(sizeof(int)*tam*tam); if (h_m == NULL || h_Bfact == NULL || A == NULL || B == NULL || sols == NULL || randmat == NULL || randmattransp == NULL) { printf("ERROR: imposible reservar memoria (h_m, h_Bfact, A, B, sols)\n"); exit(EXIT_FAILURE); } // generacion de matrices apropiadas para la factorizacion printf("\nGenerando %d matrices de dimensiones %d x %d...\n", num, tam, tam); for (k = 0; k < num; k++) { srand(2 * k + time(NULL)); // Cambio de semilla para generar matrices siempre diferentes A = crea_mat_sdp(randmat, randmattransp, A, tam); for (i = 0; i < tam; i++) { for (j = 0; j < (tam + 1); j++) { if (j == tam) { *(h_m + i*fila + k*tam + k + j) = 0; } else { *(h_m + i*fila + k*tam + k + j) = *(A + i*tam + j); } } } } // copia de datos desde 'host' a 'device' error = cudaMemcpy(d_m, h_m, tam*fila * sizeof(float), cudaMemcpyHostToDevice); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } Anexos 38 38 error = cudaMemcpy(d_Bfact, h_Bfact, tam*fila * sizeof(float), cudaMemcpyHostToDevice); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } // lanzamiento del kernel printf("\nFactorizando %d matrices de dimensiones %d x %d...\n", num, tam, tam); error = cudaEventCreate(&start); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } error = cudaEventCreate(&stop); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } error = cudaEventRecord(start, 0); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } cholesky_gpu << < BlocksPerGrid, ThreadsPerBlock >> >(d_m, d_Bfact, tam); // QueryPerformanceCounter(&iniTime); //cholesky_cpu(h_m, h_Bfact, tam); //QueryPerformanceCounter(&endTime); //printf("Tiempo en crear %d submatrices de dimensiones %d x %d: %.16g milisegundos\n", num, tam, tam, 1000*performancecounter_diff(&endTime, &iniTime)); error = cudaEventRecord(stop, 0); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } error = cudaEventSynchronize(start); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } error = cudaEventSynchronize(stop); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } 39 Implementación paralela de algoritmos de control y optimización en CUDA error = cudaEventElapsedTime(&tiempo, start, stop); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } printf("Tiempo en factorizar %d submatrices de dimensiones %d x %d: %3.5f milisegundos\n\n", num, tam, tam, tiempo); error = cudaDeviceSynchronize(); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } // copia de datos de 'device' a 'host' error = cudaMemcpy(h_m, d_m, tam*fila * sizeof(float), cudaMemcpyDeviceToHost); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } error = cudaMemcpy(h_Bfact, d_Bfact, tam*fila * sizeof(float), cudaMemcpyDeviceToHost); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } // creacion de los terminos independientes // se generan como combinacion lineal de los coeficientes crea_term_indep(h_Bfact, B, tam); // resolucion de los sistemas sols = resuelve_sistema(h_Bfact, sols, tam); // comprobacion de los resultados comprueba_soluciones(sols); // liberacion de memoria reservada free(A); free(B); free(h_m); free(h_Bfact); free(sols); free(randmat); free(randmattransp); error = cudaFree(d_m); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } error = cudaFree(d_Bfact); if (error != cudaSuccess) { Anexos 46 46 error = cudaMalloc((void**)&plan[disp].d_sols, tam*plan[disp].num_submats * sizeof(float)); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } error = cudaMallocHost((void**)&plan[disp].h_m, plan[disp].tams * sizeof(float)); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } error = cudaMallocHost((void**)&plan[disp].h_Bfact, plan[disp].tams * sizeof(float)); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } error = cudaMallocHost((void**)&plan[disp].h_sols, tam*plan[disp].num_submats * sizeof(float)); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } } // generacion de matrices apropiadas para la factorizacion printf("\nGenerando %d matrices de dimensiones %d x %d en %d GPUs...\n", num, tam, tam, NUM_GPUs); for (int disp = 0; disp < NUM_GPUs; disp++){ srand((disp+4) * 2 + time(NULL)); for (k = 0; k < plan[disp].num_submats; k++) { srand((disp+k) * (k+2) + time(NULL)); // Cambio de semilla para generar matrices siempre diferentes A = crea_mat_sdp(randmat, randmattransp, A, tam); for (i = 0; i < tam; i++) { for (j = 0; j < (tam + 1); j++) { if (j == tam) { *(plan[disp].h_m + i*plan[disp].fila_submats + k*tam + k + j) = 0; } else { *(plan[disp].h_m + i*plan[disp].fila_submats + k*tam + k + j) = *(A + i*tam + j); } } } } } // copia de datos desde 'host' a 'device' for (int disp = 0; disp < NUM_GPUs; disp++){ error = cudaSetDevice(disp); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); 47 Implementación paralela de algoritmos de control y optimización en CUDA exit(0); } error = cudaMemcpyAsync(plan[disp].d_m, plan[disp].h_m, plan[disp].tams * sizeof(float), cudaMemcpyHostToDevice, plan[disp].stream); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } error = cudaMemcpyAsync(plan[disp].d_Bfact, plan[disp].h_Bfact, plan[disp].tams * sizeof(float), cudaMemcpyHostToDevice, plan[disp].stream); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } } printf("Factorizando %d matrices de dimensiones %d x %d...\n", num, tam, tam); error = cudaEventCreate(&start); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } error = cudaEventCreate(&stop); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } error = cudaEventRecord(start, 0); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } // lanzamiento del kernel for (int disp = 0; disp < NUM_GPUs; disp++){ error = cudaSetDevice(disp); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } cholesky_gpu << < BlocksPerGrid, ThreadsPerBlock, 0, plan[disp].stream >> >(plan[disp].d_m, plan[disp].d_Bfact, plan[disp].num_submats, plan[disp].fila_submats , tam); } error = cudaEventRecord(stop, 0); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } Anexos 48 48 error = cudaEventSynchronize(start); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } error = cudaEventSynchronize(stop); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } error = cudaEventElapsedTime(&tiempo, start, stop); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } printf("\nTiempo en factorizar %d submatrices de dimensiones %d x %d: %3.5f milisegundos\n", num, tam, tam, tiempo); error = cudaDeviceSynchronize(); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } printf("\n"); // copia de datos de 'device' a 'host' for (int disp = 0; disp < NUM_GPUs; disp++){ error = cudaSetDevice(disp); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } error = cudaMemcpyAsync(plan[disp].h_m, plan[disp].d_m, plan[disp].tams * sizeof(float), cudaMemcpyDeviceToHost, plan[disp].stream); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } error = cudaMemcpyAsync(plan[disp].h_Bfact, plan[disp].d_Bfact, plan[disp].tams * sizeof(float), cudaMemcpyDeviceToHost, plan[disp].stream); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } } // liberacion de memoria reservada para la factorizacion for (int disp = 0; disp < NUM_GPUs; disp++){ error = cudaSetDevice(disp); 49 Implementación paralela de algoritmos de control y optimización en CUDA if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } error = cudaStreamSynchronize(plan[disp].stream); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } error = cudaFree(plan[disp].d_m); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } } // creacion de los terminos independientes // se generan como combinacion lineal de los coeficientes printf("Generando terminos independientes en CPU...\n"); for (int disp = 0; disp < NUM_GPUs; disp++){ crea_term_indep(plan[disp].h_Bfact, B, plan[disp].num_submats, plan[disp].fila_submats, tam); } // resolucion de los sistemas en GPU // copia de datos desde 'host' a 'device' for (int disp = 0; disp < NUM_GPUs; disp++){ error = cudaSetDevice(disp); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } error = cudaMemcpyAsync(plan[disp].d_sols, plan[disp].h_sols, tam*plan[disp].num_submats * sizeof(float), cudaMemcpyHostToDevice, plan[disp].stream); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } error = cudaMemcpyAsync(plan[disp].d_Bfact, plan[disp].h_Bfact, plan[disp].tams * sizeof(float), cudaMemcpyHostToDevice, plan[disp].stream); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } } printf("Resolviendo sistemas en la GPU...\n\n"); // lanzamiento del kernel for (int disp = 0; disp < NUM_GPUs; disp++){ error = cudaSetDevice(disp); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); Anexos 50 50 } error = cudaStreamSynchronize(plan[disp].stream); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } resuelve_sist_gpu << < BlocksPerGrid, ThreadsPerBlock, 0, plan[disp].stream >> >(plan[disp].d_Bfact, plan[disp].d_sols, plan[disp].num_submats, plan[disp].fila_submats, tam); } // copia de datos de 'device' a 'host' for (int disp = 0; disp < NUM_GPUs; disp++){ error = cudaSetDevice(disp); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } error = cudaMemcpyAsync(plan[disp].h_sols, plan[disp].d_sols, tam*plan[disp].num_submats * sizeof(float), cudaMemcpyDeviceToHost, plan[disp].stream); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } error = cudaMemcpyAsync(plan[disp].h_Bfact, plan[disp].d_Bfact, plan[disp].tams * sizeof(float), cudaMemcpyDeviceToHost, plan[disp].stream); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } } // liberacion de memoria reservada for (int disp = 0; disp < NUM_GPUs; disp++){ error = cudaSetDevice(disp); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } error = cudaStreamSynchronize(plan[disp].stream); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } error = cudaFree(plan[disp].d_Bfact); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } error = cudaFree(plan[disp].d_sols); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); 51 Implementación paralela de algoritmos de control y optimización en CUDA } error = cudaStreamDestroy(plan[disp].stream); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } } // comprobacion de soluciones aux = 0; for (int disp = 0; disp < NUM_GPUs; disp++){ comprueba_soluciones(plan[disp].h_sols, plan[disp].num_submats, disp); } // liberación del resto de memoria alojada free(A); free(B); free(randmat); free(randmattransp); for (int disp = 0; disp < NUM_GPUs; disp++){ error = cudaSetDevice(disp); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } error = cudaFreeHost(plan[disp].h_m); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } error = cudaFreeHost(plan[disp].h_Bfact); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } error = cudaFreeHost(plan[disp].h_sols); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } error = cudaDeviceReset(); if (error != cudaSuccess) { printf("Cuda failure %s:%d: '%s'\n", __FILE__, __LINE__, cudaGetErrorString(error)); exit(0); } } printf("\n"); } C. Código del bagging: mincuad.cu #include <stdio.h> Anexos 52 52 #include <stdlib.h> #include <math.h> #include <cuda_runtime.h> #include <cublas_v2.h> #include <curand.h> #include <curand_kernel.h> #include <cuda.h> #include <device_launch_parameters.h> #include "subconjuntos.cuh" //#define ThreadsPerBlock 32 const int ThreadsPerBlock = 1024; #define CUDA_CALL(x) do { if((x)!=cudaSuccess) { \ printf("Error at %s:%d\n",__FILE__,__LINE__);\ return EXIT_FAILURE;}} while(0) #define CURAND_CALL(x) do { if((x)!=CURAND_STATUS_SUCCESS) { \ printf("Error at %s:%d\n",__FILE__,__LINE__);\ return EXIT_FAILURE;}} while(0) #define PERR(call) \ if (call) {\ fprintf(stderr, "%s:%d Error [%s] on "#call"\n", __FILE__, __LINE__,\ cudaGetErrorString(cudaGetLastError()));\ exit(1);\ } #define ERRCHECK \ if (cudaPeekAtLastError()) { \ fprintf(stderr, "%s:%d Error [%s]\n", __FILE__, __LINE__,\ cudaGetErrorString(cudaGetLastError()));\ exit(1);\ } __global__ void inv_kernel(int nvect, float *a_ig, float *c_og, int n) { int idx = blockIdx.x*ThreadsPerBlock + threadIdx.x; float *a_i, *c_o; if (idx < nvect){ a_i = a_ig + idx*n*n; c_o = c_og + idx*n*n; int *p = (int *)malloc(n * sizeof(int)); int *info = (int *)malloc(sizeof(int)); int batch; cublasHandle_t hdl; cublasStatus_t status = cublasCreate_v2(&hdl); //printf("handle %d n = %d\n", status, n); info[0] = 0; batch = 1; float **a = (float **)malloc(sizeof(float *)); *a = a_i; const float **aconst = (const float **)a; float **c = (float **)malloc(sizeof(float *)); *c = c_o; status = cublasSgetrfBatched(hdl, n, a, n, p, info, batch); __syncthreads(); if (status!=0) printf("inv_kernel idx %d, status: rf %d info %d\n", idx, status, info[0]); //Depuracion status = cublasSgetriBatched(hdl, n, aconst, n, p, c, n, info, batch); __syncthreads(); 53 Implementación paralela de algoritmos de control y optimización en CUDA if (status != 0) printf("inv_kernel idx %d, status: ri %d info %d\n", idx, status, info[0]); //Depuracion cublasDestroy_v2(hdl); //printf("exito\n"); //Depuracion free(p); free(info); free(a); free(c); } } __global__ void mult_kernel(int nvect, float *Ag, int Afilas, int Acols, float *Bg, int Bcols, float *Cg){ /* Kernel que realiza la multiplicacion matricial INV * M' = RES mediante CUBLAS NOTA: CUBLAS usa como dimension principal las COLUMNAS (al contrario que C estandar, que usa filas) Para usar la misma notacion que en C, basta con intercambiar M' e INV en los argumentos de la llamada cublasSgemm */ int idx = blockIdx.x*ThreadsPerBlock + threadIdx.x; float *A, *B, *C; int tam_A = Afilas*Acols; int tam_B = Acols*Bcols; int tam_C = Afilas*Bcols; if (idx < nvect){ /* asignacion de vector corresp. a cada hilo */ A = Ag + idx*tam_A; B = Bg + idx*tam_B; C = Cg + idx*tam_C; const float alpha = 1.0f; const float beta = 0.0f; cublasHandle_t hdl; cublasStatus_t status = cublasCreate_v2(&hdl); //printf("handle %d\n", status); status = cublasSgemm_v2(hdl, CUBLAS_OP_N, CUBLAS_OP_N, Bcols, Afilas, Acols, &alpha, B, Bcols, A, Acols, &beta, C, Bcols); __syncthreads(); if (status != 0) printf("mult_kernel idx %d, status: %d\n", idx, status); cublasDestroy_v2(hdl); } } __global__ void mult_trans_kernel(int nvect, float *Mg, int Mfilas, int Mcols, float *Cg){ /* Kernel que realiza la multiplicacion matricial M' * M = C mediante CUBLAS NOTA: CUBLAS usa como dimension principal las COLUMNAS (al contrario que C estandar, que usa filas) Para usar la misma notacion que en C, basta con intercambiar M' y M en los argumentos de la llamada cublasSgemm */ int idx = blockIdx.x*ThreadsPerBlock + threadIdx.x; float *M, *C; int tam_M = Mfilas*Mcols; int tam_C = Mcols*Mcols; if (idx < nvect){ /* asignacion de vector corresp. a cada hilo */ M = Mg + idx*tam_M; C = Cg + idx*tam_C; const float alpha = 1.0f; const float beta = 0.0f; cublasHandle_t hdl; cublasStatus_t status = cublasCreate_v2(&hdl); //printf("handle %d\n", status); status = cublasSgemm_v2(hdl, CUBLAS_OP_N, CUBLAS_OP_T, Mcols, Mcols, Mfilas, &alpha, M, Mcols, M, Mcols, &beta, C, Mcols); Anexos 54 54 /* Parametros: - Supongo la matriz M almacenada en filas. Como CUBLAS usa el almacenamiento en columnas, nuestra operacion M' * M se la pasaremos a CUBLAS como M * M'. Por ello marcamos el flag CUBLAS_OP_T en el segundo argumento */ __syncthreads(); if (status != 0) printf("mult_trans_kernel idx %d, status: %d\n", idx, status); cublasDestroy_v2(hdl); } } __global__ void mult_trans_kernel_2(int nvect, float *Ag, int Afilas, int Acols, float *Bg, int Bcols, float *Cg){ /* Kernel que realiza la multiplicacion matricial INV * M' = RES mediante CUBLAS NOTA: CUBLAS usa como dimension principal las COLUMNAS (al contrario que C estandar, que usa filas) Para usar la misma notacion que en C, basta con intercambiar M' e INV en los argumentos de la llamada cublasSgemm */ int idx = blockIdx.x*ThreadsPerBlock + threadIdx.x; float *A, *B, *C; int tam_A = Afilas*Acols; int tam_B = Acols*Bcols; int tam_C = Afilas*Bcols; if (idx < nvect){ /* asignacion de vector corresp. a cada hilo */ A = Ag + idx*tam_A; B = Bg + idx*tam_B; C = Cg + idx*tam_C; const float alpha = 1.0f; const float beta = 0.0f; cublasHandle_t hdl; cublasStatus_t status = cublasCreate_v2(&hdl); //printf("handle %d\n", status); status = cublasSgemm_v2(hdl, CUBLAS_OP_T, CUBLAS_OP_N, Bcols, Acols, Afilas, &alpha, B, Afilas, A, Acols, &beta, C, Bcols); __syncthreads(); if (status != 0) printf("mult_trans_kernel_2 idx %d, status: %d\n", idx, status); cublasDestroy_v2(hdl); } } int crea_subconj_y_mat(int d, int m, int num_sub, int tam_sub, int datomax, float *d_M, float *d_Y) { /* Necesaria libreria "subconjuntos.cuh" */ /* Variables del problema */ float *d_datos_x; // Puntero a los datos (coord x)(GPU) float *d_datos_y; // Puntero a los datos (coord y)(GPU) int *d_prueba_indice; // Puntero a los subconjuntos (indice) (GPU) float *d_prueba_dato; // Puntero a los subconjuntos (dato) (GPU) //int datomin = 10; // Dato minimo para evitar redondeos a 0 al generar los puntos X curandState *d_states_ygen; 55 Implementación paralela de algoritmos de control y optimización en CUDA curandState *d_states_subconj; // states para cuRAND. Se crean tantos como vectores haya /* Reserva de memoria */ CUDA_CALL(cudaMalloc((void **)&d_datos_x, d*sizeof(float))); CUDA_CALL(cudaMalloc((void **)&d_datos_y, d*sizeof(float))); CUDA_CALL(cudaMalloc((void **)&d_prueba_indice, num_sub*tam_sub*m*sizeof(int))); CUDA_CALL(cudaMalloc((void **)&d_prueba_dato, num_sub * 2 * tam_sub*m*sizeof(float))); CUDA_CALL(cudaMalloc((void **)&d_states_subconj, m*sizeof(curandState))); CUDA_CALL(cudaMalloc((void **)&d_states_ygen, d*sizeof(curandState))); //printf("Trabajando en la GPU:\n\tSe generaran %d vectores de %d subconjuntos\n\tValor maximo de cada dato: %d\n\tTotal de puntos generados: %d\n\n", m, num_sub, datomax, d); /* Se genera un puntero de datos de m vectores y n datos por vector Se utilizan ordenes desde el host (host API) */ gen_datos_x(d, d_datos_x); /* Se escalan estos datos entre [0, datomax] */ escala_datos << < ((d) + ThreadsPerBlock - 1) / ThreadsPerBlock, ThreadsPerBlock >> > (datomax, d, d_datos_x); cudaDeviceSynchronize(); ERRCHECK; /* Se generan las coordenadas Y para cada dato */ gen_datos_y << < ((d) + ThreadsPerBlock - 1) / ThreadsPerBlock, ThreadsPerBlock >> > (d, d_datos_x, d_datos_y, d_states_ygen); cudaDeviceSynchronize(); ERRCHECK; /* Se copia al host y se imprimen (también en archivo CSV) */ //copia_datos(d, d_datos_x, d_datos_y); /* Se crean los subconjuntos */ crea_subconj << < (m + ThreadsPerBlock - 1) / ThreadsPerBlock, ThreadsPerBlock >> > (m, d, tam_sub, num_sub, d_prueba_indice, d_prueba_dato, d_datos_x, d_datos_y, d_states_subconj); cudaDeviceSynchronize(); ERRCHECK; /* Se copian al host y se imprimen */ //copia_subconj(m, num_sub, tam_sub, d_prueba_indice, d_prueba_dato); /* Inicializacion de matrices M e Y */ crea_mat << < (m + ThreadsPerBlock - 1) / ThreadsPerBlock, ThreadsPerBlock >> > (m, num_sub, tam_sub, d_prueba_dato, d_M, d_Y); cudaDeviceSynchronize(); ERRCHECK; /* Copia M e Y al host e imprime por pantalla */ //copia_mat(m, num_sub, tam_sub, d_M, d_Y); CUDA_CALL(cudaFree(d_datos_x)); CUDA_CALL(cudaFree(d_datos_y)); CUDA_CALL(cudaFree(d_prueba_indice)); CUDA_CALL(cudaFree(d_prueba_dato)); CUDA_CALL(cudaFree(d_states_subconj)); CUDA_CALL(cudaFree(d_states_ygen)); Anexos 62 62 CUDA_CALL(cudaMemcpy(h_prueba_indice, prueba_indice, total_datos * sizeof(int), cudaMemcpyDeviceToHost)); CUDA_CALL(cudaMemcpy(h_prueba_dato, prueba_dato, total_datos * sizeof(float), cudaMemcpyDeviceToHost)); /* Show result */ printf("Subconjuntos obtenidos de forma aleatoria:\n"); for (i = 0; i < num_vect; i++){ printf("\nVector numero %d:\n\n", i); for (j = 0; j < num_subconj; j++){ printf("\tSubconj num %d: ", j); for (k = 0; k < tam_subconj; k = k + 2){ printf("X[%d]: %7.4f; ", *(h_prueba_indice + i*num_subconj*tam_subconj + j*tam_subconj + k), *(h_prueba_dato + i*num_subconj*tam_subconj + j*tam_subconj + k)); printf("Y[%d]: %7.4f; ", *(h_prueba_indice + i*num_subconj*tam_subconj + j*tam_subconj + k), *(h_prueba_dato + i*num_subconj*tam_subconj + j*tam_subconj + k + 1)); printf("\n"); } printf("\n"); } } free(h_prueba_indice); free(h_prueba_dato); } int copia_mat(int num_vect, int num_subconj, int tam_subconj, float *M, float *Y){ /* Copia las matrices M e Y al host y las imprime */ /* Reserva de memoria auxiliar para el host */ float *h_M, *h_Y; h_M = (float *)malloc(2 * num_subconj*num_vect * sizeof(float)); if (h_M == NULL){ printf("Error al reservar memoria para h_M\nError en %s:%d\n", __FILE__, __LINE__); return EXIT_FAILURE; } h_Y = (float *)malloc(num_subconj*num_vect * sizeof(float)); if (h_Y == NULL){ printf("Error al reservar memoria para h_Y\nError en %s:%d\n", __FILE__, __LINE__); return EXIT_FAILURE; } /* Copia de memoria al host */ CUDA_CALL(cudaMemcpy(h_M, M, 2 * num_subconj*num_vect * sizeof(float), cudaMemcpyDeviceToHost)); CUDA_CALL(cudaMemcpy(h_Y, Y, num_subconj*num_vect * sizeof(float), cudaMemcpyDeviceToHost)); /* Imprime el resultado */ for (int i = 0; i < num_vect; i++){ printf("Matriz M[%d]\n", i); for (int j = 0; j < num_subconj; j++){ for (int k = 0; k < 2; k++){ printf("%7.4f ", *(h_M + i*tam_subconj*num_subconj + 2 * j + k)); } printf("\n"); } printf("\n"); } 63 Implementación paralela de algoritmos de control y optimización en CUDA for (int i = 0; i < num_vect; i++){ printf("Matriz Y[%d]\n", i); for (int j = 0; j < num_subconj; j++){ printf("%7.4f\n", *(h_Y + i*num_subconj + j)); } printf("\n"); } free(h_M); free(h_Y); } __global__ void escala_datos(int dato_max, int ndatos, float *d_dato){ /* Reescala los datos aleatorios a numeros en el conjunto [0, ndatos] y asigna valores a estos datos entre [0, dato_max] */ int idx = blockIdx.x * ThreadsPerBlock + threadIdx.x; if (idx < ndatos){ *(d_dato + idx) = *(d_dato + idx) * dato_max; } __syncthreads(); } __global__ void gen_datos_y(int d, float *x, float *y, curandState *state){ /* Funcion para "obtener" y en funcion de x */ /* Introducimos ruido en el termino independiente */ int idx = blockIdx.x * ThreadsPerBlock + threadIdx.x; unsigned long seed = 1; float ruido; float escala_ruido = 200000; curand_init(seed, idx, 0, &state[idx]); curandState localState = state[idx]; /* El ruido será un flotante entre [ 0, 2*escala_ruido ] */ /* Se resta escala_ruido para dejarlo entre [ -escala_ruido , escala_ruido ] */ ruido = (curand_uniform(&localState) * 2 * escala_ruido) - escala_ruido; if (idx < d){ /* funcion de la recta que forman los puntos (x,y) */ *(y + idx) = (/*2 **/ *(x + idx)) + ruido; } __syncthreads(); state[idx] = localState; } __global__ void crea_subconj(int num_vect, int ndatos, int tam_subconj, int num_subconj, int *prueba_indice, float *prueba_dato, float *d_dato_x, float *d_dato_y, curandState *state){ /* [device API] ---> se trabaja en la GPU Crea 'num_subconj' conjuntos de tamaño 'tam_subconj'. [[Se intentara que los subconjuntos se identifiquen por direccion de memoria. Asi evitamos repetir y ocupar memoria al copiar el dato]] */ int idx = blockIdx.x*ThreadsPerBlock + threadIdx.x; unsigned long seed = 1; unsigned int x; curand_init(seed, idx, 0, &state[idx]); curandState localState = state[idx]; /* Cada hebra procesa un vector */ if (idx < num_vect){ Anexos 64 64 /* Creacion de los subconjuntos */ for (int i = 0; i < num_subconj; i++){ for (int j = 0; j < tam_subconj; j = j + 2){ /* Se asigna el flotante obtenido de curand_uniform al entero x, redondeando hacia el entero mas bajo. curand devuelve un flotante en [0,1] y lo escalamos con tam_vect para elegir una posicion del vector de datos */ x = (curand_uniform(&localState) * ndatos); /* La siguiente asignacion es clave. Es donde podría asignarse la direccion de memoria en vez del dato (?) Las dos asignaciones a continuacion necesitan reservar en memoria dichos punteros (querriamos evitarlo) */ /* Se almacena el indice de cada valor */ *(prueba_indice + idx*num_subconj*tam_subconj + i*tam_subconj + j) = x; /* Se almacena el propio dato */ /* Dato X */ *(prueba_dato + idx*num_subconj*tam_subconj + i*tam_subconj + j) = *(d_dato_x + x); /* Dato Y */ *(prueba_dato + idx*num_subconj*tam_subconj + i*tam_subconj + j + 1) = *(d_dato_y + x); /* Tenemos todos los puntos (x,y) guardados en prueba_dato */ } } } __syncthreads(); state[idx] = localState; } __global__ void crea_mat(int num_vect, int num_subconj, int tam_subconj, float *prueba_dato, float *M, float *Y){ /* Inicializacion de las matrices M e Y */ int idx = blockIdx.x*ThreadsPerBlock + threadIdx.x; if (idx < num_vect){ for (int i = 0; i < num_subconj; i++){ for (int j = 0; j < tam_subconj; j = j + 2){ /* -------------------------- MATRICES ----------------------- --- */ /* Se inicializan las matrices M e Y para el calculo de parametros */ *(M + idx*tam_subconj*num_subconj + i*tam_subconj + j) = 1; *(M + idx*tam_subconj*num_subconj + i*tam_subconj + j + 1) = *(prueba_dato + idx*num_subconj*tam_subconj + i*tam_subconj + j); *(Y + idx*num_subconj + i) = *(prueba_dato + idx*num_subconj*tam_subconj + i*tam_subconj + j + 1); /* ----------------------------------------------------------- --- */ } } } __syncthreads(); } 65 Implementación paralela de algoritmos de control y optimización en CUDA E. Código del bagging: subconjuntos.cuh #include <device_launch_parameters.h> #ifndef SUBCONJUNTOS_CUH #define SUBCONJUNTOS_CUH int gen_datos_x(int d, float *d_dato); int copia_datos(int d, float *d_datos_x, float *d_datos_y); int copia_subconj(int num_vect, int num_subconj, int tam_subconj, int *prueba_indice, float *prueba_dato); int copia_mat(int num_vect, int num_subconj, int tam_subconj, float *M, float *Y); __global__ void gen_datos_y(int d, float *x, float *y, curandState *state); __global__ void escala_datos(int dato_max, int ndatos, float *d_dato); __global__ void crea_subconj(int num_vect, int ndatos, int tam_subconj, int num_subconj, int *prueba_indice, float *prueba_dato, float *d_dato_x, float *d_dato_y, curandState *state); __global__ void crea_mat(int num_vect, int num_subconj, int tam_subconj, float *prueba_dato, float *M, float *Y); #endif Anexos 66 66 67 Implementación paralela de algoritmos de control y optimización en CUDA REFERENCIAS [1] Sanders, J. y Kandrot, E. CUDA by example, An Introduction to General-Purpose GPU Programming. 2010. [2] James, G., Witten, D., Hastie, T., Tibshirani, R. An Introduction to Statistical Learning. 2013. [3] Hastie, T., Tibshirani, R., Friedman, J. The Elements of Statistical Learning. 2009. [4] Zeller, C. CUDA C/C++ Basics, Supercomputing 2011 Tutorial. 2011. [5] NVIDIA Corporation. CUDA Toolkit Documentation. CUDA C Programming Guide. [En línea] Consultado el 14 de febrero de 2018. Disponible en: https://docs.nvidia.com/cuda/cuda-c-programming- guide/ [6] Stack Overflow. nvidia simpleDevLibCUBLAS error. [En línea] Consultado el 5 de junio de 2018. Disponible en: https://stackoverflow.com/questions/37653654/nvidia-simpledevlibcublas-error [7] Stack Overflow. CUDA 5.0 separate compilation and linking option in Visual Studio. [En línea] Consultado el 12 de junio de 2018. Disponible en: https://stackoverflow.com/questions/15192787/cuda-5- 0-separate-compilation-linking-option-and-visual-studio-2010 [8] NVIDIA Corporation. CUDA Toolkit Documentation. cuBLAS. [En línea] Consultado el 20 de mayo de 2018. Disponible en: https://docs.nvidia.com/cuda/cublas/ [9] NVIDIA Corporation. CUDA Toolkit Documentation. cuRAND. [En línea] Consultado el 13 de mayo de 2018. Disponible en: https://docs.nvidia.com/cuda/curand/ [10] Rosetta Code. Cholesky decomposition. [En línea] Consultado el 26 de marzo de 2018. Disponible en: https://rosettacode.org/wiki/Cholesky_decomposition [11] Terraludus. Consejos C++: midiendo tiempos. [En línea] Consultado el 1 de abril de 2018. Disponible en: http://terraludus.blogspot.com/2012/06/consejos-c-ii-midiendo-tiempos.html [12] Stack Overflow. Different occupancy between calculator and nvprof. [En línea] Consultado el 2 de mayo de 2018. Disponible en: https://stackoverflow.com/questions/23469180/different-occupancy-between- calculator-and-nvprof [13] Stack Overflow. Re-seeding cuRAND host API. [En línea] Consultado el 1 de junio de 2018. Disponible en: https://stackoverflow.com/questions/42072377/re-seeding-curand-host-api Referencias 68 68 69 Implementación paralela de algoritmos de control y optimización en CUDA GLOSARIO CUDA: Compute Unified Device Architecture. Tecnología de cálculo paralelo desarrollada por NVIDIA y destinada a aprovechar el perfil de las tarjetas gráficas en tareas comunes 1 GPU: Graphics Processing Unit 1 CPU: Central Processing Unit. Responsible de cálculos, control y supervisión de otras partes en un ordenador1 C/C++, Fortran: Lenguajes de programación de propósito general 1 OpenACC: Estándar de programación de informática paralela 1 Host: En programación paralela, anfitrión, lado asociado a la CPU 1 Device: Al contrario que host, dispositivo, lado asociado a la GPU 1 Kernel: Función que se ejecuta en el dispositivo, en un determinado número de bloques e hilos 1 Thread: Hilo de ejecución en CUDA 2 Block: Bloque formado por hilos en CUDA 2 Grid: Conjunto que abarca varios bloques CUDA 2 SM: Streaming Multiprocessor. Cada uno de los multiprocesadores que componen una GPU 5 SIMT: Single-Instruction, Multiple-Thread. Arquitectura empleada en los multiprocesadores de una GPU 6 Warp: Conjunto de hilos que son procesados a la vez en un multiprocesador de una GPU con una sola instrucción 6 Host API: Hace referencia al conjunto de funciones que son ejecutadas en el lado del host 10 Device API: Hace referencia al conjunto de funciones que son ejecutadas en el lado del device 11