PUNTO DE PARTIDA -...
-
Upload
nguyenthien -
Category
Documents
-
view
218 -
download
1
Transcript of PUNTO DE PARTIDA -...
PROGRAMACIÓN EN EL ENTORNO CUDA EN APLICACIONES DE MECÁNICA COMPUTACIONAL
CÓDIGO DEL PROYECTO
146
VIII. Código Bloque II: GENERACIÓN DE LA Matriz DE
RIGIDEZ A
Este bloque se encarga de la generación de las matrices locales para
cada triángulo, su ensamblado en una matriz única de rigidez y la
simplificación de ésta, eliminando de ella, mediante sumas, las
multiplicidades de los términos compartidos. Esta matriz será almacenada
en formato COO para poder convertirla en un archivo MATRIX MARKET de
manera sencilla y poder manejarla en otros puntos del problema.
1. PUNTO DE PARTIDA:
En Matlab y mediante la herramienta de mallado se generó una trama,
definiendo el dominio, las condiciones de contorno y el grado de
mallado. Este sistema se exportó a archivos mediante el código
“trataMalla” que genera los archivos triangulos.dat, coordenadas.dat,
dirichlet.dat y neumann.dat.
A partir de ahí el objetivo consiste en:
1- Lectura de los archivos que describen el mallado.
2- Cálculo de las matrices locales para cada uno de los elementos del mallado. (Matrices 3x3).
3- El ensamblado de dichas matrices en una matriz dispersa de tamaño Num_nodosxNum_nodos, eliminando duplicidades y valores nulos.
4- El almacenamiento de la matriz generada en formato COO que podría ser volcada en un archivo en formato Matrix Market que sería la entrada
del solver.
METODOLOGÍA
La idea en la que se basaría la solución sería desarrollar un código
usando exclusivamente lenguaje C con las extensiones CUDA. Sin usar
bibliotecas ni códigos previamente desarrollados. La ventaja de este
desarrollo es que el código es una solución optimizada y a medida del
problema que nos ocupa.
Se optó por dividir el problema en dos partes:
1- Generador de matrices que depende de los triángulos.
2- Simplificador de la Matriz a formato COO que depende de los nodos.
Se ha usado la CUDA API C Runtime para CUDA. Se han usado como lenguajes
de programación C y C++. No se han usado bibliotecas estándar ni
"tamplates": CUBLAS, CUSP, THRUST.
PROGRAMACIÓN EN EL ENTORNO CUDA EN APLICACIONES DE MECÁNICA COMPUTACIONAL
CÓDIGO DEL PROYECTO
147
En la figura 1 se ilustra la estrategia usada en este bloque.
Fig 1: Diagrama de flujo de la adaptación objeto de este bloque.
2. MODELO MATEMÁTICO Y SIMPLIFICACIONES
Uno de los elementos fundamentales para la optimización de código es
realizar simplificaciones previas a la codificación. Esto se vuelve
especialmente importante en CUDA cuando tenemos en cuenta que si la
simplificación conlleva la disminución en el número de variables y
operaciones, esto se traduce en una disminución de registros necesarios
para la ejecución y por tanto una mayor velocidad potencial de la
aplicación.
Además un análisis previo del problema algebraico nos permite deducir
características de la matriz final y de las matrices intermedias que
simplificarán los cálculos, disminuyendo como veremos las operaciones
necesarias. Por tanto, a partir de cada triángulo del mallado se
construye una matriz de 3x3 que será la matriz local. Para la
determinación de dichos elementos intervienen las coordenadas de los
tres vértices.
Estos elementos corresponden a combinaciones de los vértices de cada
triángulo tomados de dos en dos.
Para cada triángulo tenemos tres vértices: v1, v2 y v3; que les
corresponden 2 coordenadas de tipo (x,y).
Se forma la matriz de coordenadas que corresponde a los vértices de cada
triángulo: [x1 y1
x2 y2
x3 y3]
BLOQUE DE PARTIDA
EN CPU
ADAPTACIÓN CUDA
Generación de la
Matriz de
Rigidez
Eliminación de
Redundancias y
Ceros en la matriz
ensamblada, en GPU
Traslado de la
información a
GPU.
Volcado a archivo
MM de la matriz
de Rigidez.
Lectura de
archivos
tiengulos.dat y
coordenadas.dat
Generación de las
Matrices Locales y
su ensamblado, en
GPU.
PROGRAMACIÓN EN EL ENTORNO CUDA EN APLICACIONES DE MECÁNICA COMPUTACIONAL
CÓDIGO DEL PROYECTO
148
El jacobiano será:
∣[1 x1 y1
1 x2 y2
1 x3 y3]∣ =(x2*y3)+(x3*y1)+(x1*y2)-(x2*y1)-(x1*y3)-(x3*y2) También se construye la matriz Bkt que es la matriz de transformación
que lleva los vértices del triángulo (0,0) (0,1) (1,0) a la posición de
los vértices V1, V2 y V3:
[x2 y2
x3 y3]- [x1 y1
x1 y1]= [x2− x1 y2− y1
x3− x1 y3− y1] La inversa de Bkt es inmediata por ser una matriz 2x2, inv Bkt:
Bkt− 1=
1
Det * [y3− y1 y1− y2
x1− x3 x2− x1] A partir de ahí construimos la matriz CK definida como: CK= Bkt− 1T
* Bkt− 1
De ahí, CK=
1
Det2 * [y3− y1 x1− x3
y1− y2 x2− x1]* [y3− y1 y1− 23
x1− x3 x2− x1]
CK=
1
Det2 *
2 2
2 2
y3 y1 y1 y2 + x1 x3 x2 x1y3 y1 + x1 x3
y3 y1 y1 y2 + x1 x3 x2 x1 y1 y2 + x2 x1
Se usan las siguientes matrices auxiliares definidas para las matrices
locales:
K XX =
12 [
1 − 1 0
− 1 1 0
0 0 0], KYY =
12 [
1 0 − 1
0 0 0
− 1 0 1 ],
K XY =
12 [
1 0 − 1
− 1 0 1
0 0 0 ] y K XYT
=
12 [
1 − 1 0
0 0 0
− 1 1 0] Para cada triángulo se construye la matriz local:
R=jacob*{CK11* K XX +CK22* KYY +CK12*[ K XY + K XYT
]}
K XY + K XYT
=
12 [
1 0 − 1
− 1 0 1
0 0 0 ]+ 12 [
1 − 1 0
0 0 0
− 1 1 0]= 12 [
2 − 1 − 1
− 1 0 1
− 1 1 0 ]
PROGRAMACIÓN EN EL ENTORNO CUDA EN APLICACIONES DE MECÁNICA COMPUTACIONAL
CÓDIGO DEL PROYECTO
149
R=
12 *
1
Det2 *jacob*
2Ck12 0 0
0 0 0 0 0
0 0 0 0 0
CK11+CK22+ CK11+ CK12 CK22 CK12
CK11+ CK12 CK11+ + + +CK12
CK22 CK12 + +CK12 +CK22+
Observamos que:
R11= R22+ R33+2*R23
R12=-R22+0*R33- R23
R13=0*R22- R33- R23
y que:
R21=R12
R31=R13
R32=R23
Concluimos que las matrices locales son reales, cuadradas y simétricas
con sólo 3 elementos independientes. Por tanto necesitamos obtener sólo
3 elementos que formarán la base de nuestros cálculos: R22, R33 y R23
tal y como están definidos en función de CK. Los 6 restantes son
deducidos de manera inmediata. Estos cálculos permiten simplificar el
código.
Otro elemento esencial es que estas matrices locales una vez ensambladas
forman una matriz de rigidez simétrica real, dispersa, y con diagonal no
nula. Todas estas características son útiles para generarla y
simplificarla.
3. ALGORITMO DE GENERACIÓN
Para la parte de generador, se parte del código MATLAB, y se expresa en
ALGEBRA LINEAL. Tras realizar los cálculos pertinentes para simplificar
las operaciones y permitir en la medida de lo posible la obtención de
resultados inmediatos, se transcribió a lenguaje de programación C.
Rápidamente se observa que los cálculos de mayor exigencia, son los
bucles y entre éstos los bucles anidados (dobles o incluso triples). En
particular el bucle para calcular las matrices locales. Esta parte del
código se transcribió a CUDA mediante inserción de código Host y código
device.
PROGRAMACIÓN EN EL ENTORNO CUDA EN APLICACIONES DE MECÁNICA COMPUTACIONAL
CÓDIGO DEL PROYECTO
150
En la CPU se realiza la lectura de los archivos y el almacenamiento de
los datos en matrices. Por un lado se realiza la lectura del archivo de
coordenadas y se almacena su contenido en la matriz “Cor” de coordenadas
donde la fila “i” corresponde al vértice “V=i+1” y almacena dos valores
x(V) e y(V) que corresponden a sus coordenadas. Por otro lado se realiza
la lectura del archivo de triángulos y se almacena su contenido en la
matriz “Tri” donde la fila “i” corresponde al triángulo “T=i+1” y
almacena tres valores V1, V2 y V3 que corresponden a los vértices que lo
forman.
Es de señalar que no existe ni el triángulo "0", ni tampoco el vértice
“0”. Esto es importante debido al indexado de matrices en C, donde,
recordémoslo, sí existen las posiciones 0.
Estas matrices se copian a la GPU donde se generan las matrices
resultado. Para cada uno de los vértices se leen sus coordenadas
correspondientes y se realizan los cálculos mediante la invocación del
kernel. Finalmente los resultados se trasladarán a la CPU.
La sección del kernel es sin duda la más importante del código porque es
la que implanta el algoritmo de generación de matrices locales y además
lo hace de manera paralela y no secuencial. Cada hilo ejecuta el
algoritmo para un triángulo individual.
Las operaciones son básicamente las expuestas en el desarrollo
algebraico descrito arriba. Se ha buscado por un lado reducir al máximo
los pasos y por otro usar el menor número posible de registros. Esto
permite acelerar la ejecución y aumentar los recursos disponibles para
ejecutar un gran número de hilos simultáneamente.
La lectura de los datos de entrada para cada triángulo se hace casi de
manera independiente leyendo secciones independientes de la matriz de
triángulos. Aunque en teoría estas secciones coinciden o se repiten
entre triángulos (porque los triángulos comparten vértices), en la
práctica no se ha observado una saturación debido a la lectura, ni
tampoco una aceleración cuando se han almacenado las coordenadas en
memoria compartida en lugar de hacerlo en los registros.
Sin embargo la escritura de los resultados se debe hacer de manera
independiente. No se controla la secuencia de escritura de los hilos
independientes por ello si dos o varios hilos escribieran en una misma
posición actualizando el dato almacenado pueden surgir inconsistencias.
Para ello se ha optado por reservar memoria suficiente con posiciones
contiguas para que cada hilo almacene los 9 datos que genera (Que
corresponden a los elementos R de la matriz local) usando 3 arrays que
forman una matriz en el formato COO.
La matriz I almacena la fila que corresponde al vértice tomado como
índice inicial, la matriz J almacena las columnas correspondientes a
todas las combinaciones usando los demás vértices del triángulo para
obtener todas las combinaciones posibles y la matriz Val almacena los
valores que corresponden.
PROGRAMACIÓN EN EL ENTORNO CUDA EN APLICACIONES DE MECÁNICA COMPUTACIONAL
CÓDIGO DEL PROYECTO
151
El índice de almacenamiento y el orden se basan en el índice del
triángulo. De tal modo que al triángulo T le corresponden las posiciones
de 9*T hasta la 9*T+8 todas contiguas en cada uno de los 3 arrays.
El sistema de almacenamiento de resultados que se ha usado se basa en
reservar en cada array un bloque de posiciones contiguas para la
escritura. Como índice se usan los triángulos. Por tanto cada hilo
escribe en una posiciones identificadas e independientemente de los
demás hilos en ejecución o los que puedan ejecutarse o arrancarse
durante su existencia.
La reserva de estos espacios con una dirección de comienzo determinada
por el índice del hilo (en total 9 posiciones en cada matriz para cada
triángulo), evita colisiones en escritura e independiza los accesos de
los hilos.
El orden descrito en este apartado será esencial para realizar las
tareas de simplificación.
En las matrices generadas, sin embargo la posición "i" corresponde al
triángulo “i+1” ya que en C la matrices empiezan en el índice “0”.
Para evitar limitaciones de memoria, se ha optado por volcar el
contenido de las matrices a un archivo provisional: “prov.dat”.
4. ALGORITMO DE SIMPLIFICACIÓN
La matriz generada en el apartado anterior es una matriz que contiene
valores nulos y duplicidades. En realidad su formato no es el de una
matriz COO propiamente dicho. Por ello se hace necesario un paso
adicional que es someterla a un ALGORTIMO SIMPLIFICADOR que elimina los
valores nulos y las DUPLICIDADES.
Las duplicidades Se generan debido a que los vértices son compartidos
por varios triángulos. Además los triángulos contiguos comparten un
lado. Las duplicidades dependen del patrón de mallado, pero generan
pares (I,J) que se repiten aunque puedan tener el mismo valor. Especial
relevancia tienen los pares (I,I) que constituyen la diagonal de la
matriz de ensamblado.
El número de estas duplicidades puede llegar a ser muy importante y de
hecho lo demuestran los datos de salida durante las pruebas.
PROGRAMACIÓN EN EL ENTORNO CUDA EN APLICACIONES DE MECÁNICA COMPUTACIONAL
CÓDIGO DEL PROYECTO
152
Para eliminar estas duplicidades se probaron 3 soluciones:
a. MÉTODO DIRECTO:
La idea es convertir el código CPU siguiente:
for(i=0; i<N; i++){
for(j=i+1; j<N; j++){
if((I[j]==I[i])&&(J[j]==J[i])){
Val[i]+=Val[j]
Val[j]=0
}
}
En código GPU de manera directa:
Si definimos t como índice de hilos, una solución sería ejecutar para
cada hilo el bucle interno
for(j=t+1; j<N; j++){
if((I[j]==I[t])&&(J[j]==J[t])){
T[t]=j
j=N
}
}
Las posiciones se recorren una a una y se comparan.
En la CPU este bucle tarda horas, sin embargo en la GPU la versión
paralela no se ejecuta y la GPU lo interrumpe.
Se debe señalar que al no poder usar el depurador CUDA-gdb por no
disponer de una salida vídeo adicional, se optó por incorporar código
adicional para ayudar a depurar el programa. Esto se hizo mediante la
inserción de código para producir mensajes por pantalla para controlar
los puntos de ejecución (Una especie de orden Break que no interrumpe la
ejecución) y bloques de escritura a archivo para conocer y guardar el
contenido de las variables.
En especial para los kernels se mostró muy útil el uso del temporizador
menor que mide la duración de ejecución del kernel. El temporizador
mayor se usa para medir el tiempo total de ejecución del programa o de
un módulo; incluye por tanto el tiempo CPU y tiempo GPU.
PROGRAMACIÓN EN EL ENTORNO CUDA EN APLICACIONES DE MECÁNICA COMPUTACIONAL
CÓDIGO DEL PROYECTO
153
En condiciones exitosas la duración de un kernel es del orden de decenas
o centenares de milisegundos; en caso de fracaso se detectan valores
absurdos (enormes o infinitamente pequeños) o bien sencillamente "0".
Al no disponer de acceso a los valores de registros de la CUDA en tiempo
de ejecución, sólo se puede especular acerca de la causa. Sin embargo
esta especulación está basada en experiencias parecidas expuestas en los
foros CUDA y NVIDIA.
Un kernel "rebota" cuando los recursos son insuficientes para su
ejecución. Especialmente cuando se agotan las memorias tipo “shared
memory” y registros. La reserva de recursos se hace previamente a la
ejecución por ello simplemente se interrumpe.
Al examinar el código expuesto y sus variaciones se pueden concluir 2
hipótesis:
1- No sólo las escrituras deben ser no competitivas entre hilos;
sino también las lecturas.
2- Recorrer el array en orden implica el almacenamiento de un gran
número de datos. Estos espacios o no se liberan durante la
ejecución del hilo (por ejemplo escribiendo encima de ellos en un
bucle “for”) o lo hacen con lentitud ocasionando el agotamiento de
los registros.
Se debe observar que incluso cuando se modificaron los parámetros de
lanzamiento de kernel bajando al máximo el número de hilos (llegando a 4
hilos por bloque), no se logró que arrancara.
b. MÉTODO DE ORDENACIÓN
Reconsiderando el problema y especialmente la disposición de los datos,
en el archivo “prov.dat”; se observa que los datos están dispuestos
según los triángulos. Esta disposición no es la más idónea para una
matriz donde el orden natural es según nodos. Por ello se buscó algún
método para ordenar los arrays según los nodos. Y una vez ordenados, era
de esperar que las posiciones contiguas contuvieran los valores
repetidos. Por tanto el problema se reduciría a recorrer los arrays en
bloques y no en su conjunto.
Se partiría del primer elemento y se compararía con los k siguientes que
contiene los datos que poseen la misma fila.
En resumen: Se ordenan I, J, y VAL según los nodos fila de I, obteniendo
en las posiciones contiguas el mismo I y distintos J que son los que se
recorrerían buscando repeticiones.
PROGRAMACIÓN EN EL ENTORNO CUDA EN APLICACIONES DE MECÁNICA COMPUTACIONAL
CÓDIGO DEL PROYECTO
154
Para llevar a cabo esta idea se debía usar un código de ordenación
(sorting).
- THRUST dispone de funciones templates (plantillas) para arrays
lineales.
- La biblioteca estándar de C contiene la función qsort.
- La SDK de NVIDIA también contienen implantaciones de algoritmos
Sorting.
En general esas soluciones manejan arrays lineales. Una solución
barajada era convertir cada línea en una cadena; es decir, generar un
array lineal de cadenas (strings) y ordenarlo, porque al ordenador el
array de columnas debemos reorganizar el de valores. Sin embargo esta
solución tenía dos inconvenientes:
1. Conversión adicional de datos con especial cuidado a los
formatos y longitud de cada campo.
2. No aprovechamiento en absoluto de la disposición de datos en la memoria que nos es conocida a priori.
Esta opción sigue siendo válida, pero para un diseño adaptado era
preferible particularizar la solución.
c. MÉTODO DE MAPEADO
Para este método se tuvo en cuenta la disposición de los datos en la
memoria.
Se decidió, por tanto, mapear los datos de los triángulos a vértices.
Recordemos que triangulos.dat y la matriz que de él se genera contienen
para la posiciones “i” que corresponde al triángulo "i+1" los vértices
V1, V2 y V3. En cambio la matriz Vértices se construye haciendo que la
línea "i" corresponda al vértice "i+1" y contenga los triángulos T1, T2,
hasta Tk. Es decir a cada vértice le hace corresponder todos los
triángulos en los que aparece.
Una vez mapeados los triángulos a vértices y usando matrices auxiliares
se ejecuta el código de ordenación en la GPU.
En concreto se usa la matriz V, que contiene las veces que se repite un
vértice en el mallado (en cuantos triángulos entra a formar parte) y la
matriz P, que almacena las posiciones de comienzo de almacenamiento de
todos los elementos que tienen la misma fila.
P[i] es la posición de comienzo de los elementos cuya fila es "i+1". Por
otro lado V[i] nos indica las veces que un nodo se repite. Para
almacenar todos los elementos de la misma fila se necesitan 3xV[i]
posiciones de memoria. Que mediante P se direccionan y se hacen
contiguas.
PROGRAMACIÓN EN EL ENTORNO CUDA EN APLICACIONES DE MECÁNICA COMPUTACIONAL
CÓDIGO DEL PROYECTO
155
No hay colisiones porque los hilos que aquí corresponden a vértices leen
de posiciones determinadas, en concreto, de aquellas en que su índice
coincide con la fila.
Estas posiciones corresponden a 3 posiciones seguidas en las matrices de
resultados según los triángulos. Para escribir también lo hacen en
posiciones determinadas de las matrices de salida previamente fijadas
evitando colisiones.
Se debe señalar que se eligió almacenar los índices de los nodos filas
de un modo lineal en los arrays. Es decir primero todos los elementos
cuyo elemento fila es "1", luego todos los que van en la fila "2" etc.
Por tanto durante todo el proceso se ahorra el manejar al array de
filas.
Este kernel se ha llamado ordenadorSegunI.
Seguidamente se llama un segundo kernel llamado “eliminador” que realiza
comparaciones para cada fila entre los nodos de las columnas. En este
caso se obtiene, como las posiciones son contiguas y conocidas
previamente (porque dependen de las veces que un vértice es compartido
por los triángulos), bucles pequeños de un máximo de k, que son
recorridos cada uno por un hilo independiente. En estos bucles se
realizan comparaciones. Cuando hay una coincidencia se acumula en el
valor inicial y se anula en los siguientes.
En lugar de realizar desde el kernel una llamada a función device, se
prefiere implantar un nuevo kernel. El objetivo es disminuir los
requisitos de ejecución permitiendo que no se agoten en especial los
registros.
Finalmente los datos no nulos obtenidos son volcados a un archivo. Esto
se hace de manera sencilla insertando la condición de no escribir los
valores nulos.
5. PRUEBAS Y OBSERVACIONES
En las pruebas iniciales realizadas en el equipo de transición, se trató
un sistema de mallado con 1277952 triángulos y 640257 nodos. El programa
se ejecutó correctamente con los siguientes datos relevantes:
PRIMERA PARTE: GENERADOR:
La ejecución del programa ha requerido: 11 segundos en total.
La ejecución del programa ha requerido: 298.816 milisegundos en la GPU.
SEGUNDA PARTE: SIMPLIFICADOR:
La ejecución del programa ha requerido: 32.09 segundos en total.
La ejecución del programa ha requerido: 923.906 milisegundos en la GPU.
PROGRAMACIÓN EN EL ENTORNO CUDA EN APLICACIONES DE MECÁNICA COMPUTACIONAL
CÓDIGO DEL PROYECTO
156
Se observó lo siguiente:
- Los códigos generados son bastante rápidos. Recordemos que la
tarjeta usada tiene 4 procesadores con un total de 32 núcleos.
- Los resultados coinciden con los ofrecidos para MATLAB. Esto se
comprueba mediante inspección de los archivos generados por
ambas aplicaciones.
- Los arrays no se pueden inicializar en la GPU pasando el tamaño
como parámetro de la función de llamada. Es decir A[k] no es
válido. En cambio A[6] sí lo es.
- El tamaño de los bloques actúa como des-acelerador de la
velocidad, pero permite acomodar un máximo de hilos.
- Dentro de un mismo código host no se pueden variar las
dimensiones de la malla y de los bloques entre los distintos
kernels a los que se llama.
- Los códigos se desarrollaron de forma compartimentada para
permitir un análisis de cada sección del problema, aunque luego
su fusión no produjo inconvenientes.
- Hay que observar que cualquier bucle "for" excepto los de
lectura y escritura de archivos son susceptibles de ser
trasladados a CUDA. Incluyendo por ejemplo los de mapeado, el
de determinación del máximo orden de repeticiones y los de
inicialización de matrices y arrays (zero padding). La
paralelización de estas secciones supondrá sin duda una mejora.
- No se ha considerado el uso de otro tipo de memorias por
ejemplo texture.
- A la hora de determinar los tamaños de malla y bloques se debe
garantizar que todos los hilos "caben" para que sean
ejecutados, de ahí el redondeo al entero superior y la
condición de “t<número total de tríangulos”; por ejemplo.
- El formato de los datos contenidos en los archivos de entrada
generados por Matlab (triangulos.dat y coordenadas.dat) produjo
problemas al intentar almacenarlos directamente como “int”. Se
deben leer como “float” y reconvertir mediante “cast” a “int”.