A
proposal for pseudospectra computation on graphic processor units
Zenaida Natividad Castillo Marrero. [1],
Gustavo Adolfo Colmenares Pacheco. [2],
Paulina Elizabeth Valverde Aguirre. [3] &=
amp;
VÃctor Oswaldo Cevallos Vique. [4]
Recibido:
09-03-2021 / Revisado: 16-03-2021 /Aceptado: 06-04-2021/ Publicado: 05-05-2=
021
Keywords:
Pseudospectra,
Eigenvalues, Krylov methods, GPU, NVIDIA.
Resumen.
Palabras claves:
Pseudoespectro, Autovectores, Métodos de Krylov, GPU,
NVIDIA.
Introducción
En los modelos lineales =
la
investigación se basa en la información provista por los autovalores, y p=
ara muchos
problemas relacionados con la ciencia y la ingenierÃa este análisis es
satisfactorio, sobre todo cuando las matrices son normales. Sin embargo, en
algunas áreas de la industria, las matrices que se derivan del modelaje
matemático, no poseen columnas ortogonales, y en caso de que puedan
ortogonalizarse el proceso numérico es costoso o resulta en aproximaciones=
muy
pobres o poco valederas para la toma de decisiones. Adicionalmente, la
no-normalidad influye negativamente en la convergencia de los métodos
iterativos. En estos casos el pseudoespectro es una herramienta útil, que
proporciona una alternativa visual para reforzar el análisis espectral. Ã=
reas
de aplicación exitosa de técnicas con autovalores y pseudoautovalores inc=
luyen
el estudio de la acústica, análisis estructural, mecánica cuántica, mec=
ánica de
fluidos, procesos de revestimiento de superficies, entre otros.
Con la llegada de los
ordenadores de alto rendimiento, que permiten cálculos en paralelo, se abr=
ió la
posibilidad de avanzar con el procesamiento masivo de datos en aplicaciones=
a
gran escala. De esta manera muchos autores se avocaron a la tarea de redefi=
nir
o adaptar sus algoritmos para sacar el mayor provecho al uso de estas
tecnologÃas, entre ellos vale la pena mencionar los trabajos de Bekas et a=
l.
(2001), Mezher & Philippe (2011), y Noschese & Reichel (2015). Estos
algoritmos han sido programados en su mayorÃa a través de una arquitectur=
a de
pase de mensajes como MPI o OpenMP, ver por ejemplo el trabajo de Bekas et =
al.
(2002) con MPI_Matlab o el de Minini et al. (2011) con un hÃbrido MPI-Open=
MP,
ambos implementados en computadores con decenas de miles de procesadores.
Por otra parte, en la ú=
ltima
década, las unidades de procesamiento gráfico o GPU (por sus siglas en in=
glés)
adaptadas a la gran mayorÃa de microcomputadores y computadores personales=
, han
probado tener una capacidad de cómputo que puede ser usada para aligerar
significativamente la carga de la unidad principal de procesamiento o CPU. =
Este
hecho, sumado a que no siempre se dispone de una máquina de alto rendimien=
to,
nos lleva a considerar la implementación de estos métodos bajo un esquema
hÃbrido que involucre cálculos numéricos densos en las GPU y manejo de
entrada-salida a cargo del CPU.
En este trabajo se descr=
ibe
una propuesta para la implementación de un esquema paralelo para el cálcu=
lo del
pseudoespectro de matrices en un modelo CPU-GPU.
Metodologia
Para el desarrollo de la
propuesta se describe el problema de hallar el pseudoespectro de matrices de
gran tamaño y se mencionan algunos métodos que han resultado efectivos en=
su
resolución en computadores de alto rendimiento. En particular se detalla el
funcionamiento de los algoritmos que integran la propuesta y su posible
implementación en un modelo de arquitectura CPU-GPU.
El problema de los autovalore=
s y
el pseudoespectro
El
espectro de una matriz
Ahora bien, si
Otra definición equival=
ente,
muy usada en teorÃa de perturbaciones, expresa al pseudoespectro en térmi=
nos de
los autovalores de una matriz perturbada en un épsilon=
(
Si <=
![if !msEquation]>
Los diferentes métod=
os
para el cálculo del pseudoespectro de una matriz usan estas definiciones u
otras equivalentes. En este trabajo se propone un esquema para el cálculo =
del
pseudoespectro usando esta última definición basada en el valor singular =
mÃnimo
de la matriz
Representación Gráfica del pseudoespectro
Uno de los aportes del pseudoespectro es que el análisis puede
realizarse al observar la gráfica que se genera variando los valores de é=
psilon
(
Si una matriz u operador
La figura 1 muestra curvas de aproximación del pseudoespectro=
de
una matriz hermitiana de orden
Figura
1. Pseudoespectro
de matrices Normales (a) y No-Normales (b)
Fuente:=
Elaboración
propia.
Las
gráficas se presentan en el plano complejo, y los puntos negros representa=
n los
autovalores de la matriz, mientras que los trazos coloreados representan los
lÃmites del pseudoespectro para cada valor de épsilon, el cual se indica =
en la
barra de colores a la derecha. Por ejemplo, la curva más externa de color
amarillo representa una perturbación de la matriz en =
Por otra parte, tenemos los algoritmos de malla, que tienen un
esquema simple y organizado de ejecución. En estos se define una región de
interés en el plano complejo, la cual se discretiza mediante una malla tan
refinada como se desee y por cada punto de la malla se calcula
A continuación, en la figura 3, adaptada de Otero et. al (201=
5),
se presenta un esquema básico del cálculo, basados en algoritmo de malla,=
para
el cálculo del pseudoespectro. La idea general de estos esquemas fue propu=
esta
por Trefethen (1999). Algunos autores han implementado con éxito este tipo=
de
métodos, véase por ejemplo Trefethen & Wright (2001) y las referencia=
s en
Trefethen & Embree (2005).
Figura =
3. Algorit=
mo
propuesto para matrices de gran tamaño
Fuente:
Elaboración propia.
Los algoritmos que proponemos contienen intrÃnsecamente
paralelismo de grano grueso (mÃnima comunicación entre unidades) una vez =
que
están basados en la discretización de la región de interés, y paralelis=
mo fino
(mucha comunicación entre subtareas) en las operaciones matriz-vector que =
se
utilizan para el cálculo de
Proyecciones en espacios de Krylov
Los
métodos basados en Arnoldi se fundamentan en la descomposición Hessenberg=
de la
matriz A o
Pa=
ra
obtener la descomposición de Hessenberg
En este caso, dado un vector inicial
Cuando la matriz
La precisión y la velocidad de convergencia de este tipo de
métodos mejora con el incremento de
Figura 4. Tiempo
de CPU: Iteración m de Arnoldi
Fuente: Elaboración
propia.
Discretización de la región de interés
De
acuerdo a la aplicación, el investigador pudiera estar interesado en la pa=
rte
del espectro donde los autovalores toman un valor en particular, por ejempl=
o,
los de mayor módulo en estudios de estabilidad, o para el análisis de
convergencia de métodos iterativos. Este interés define la región en el =
plano
en la cual se analiza el pseudoespectro.
Una
vez definida la región de interés
Cada punto
Figura 5(a). Â Malla de
Fuente: Elaboración prop=
ia.
Consideremos también que tenemos una arquitectura paralela co=
n
Figura
5(b).
Distribución de la carga en los procesadores
Fuente:
Elaboración propia.
Cada subregión se asocia a un proceso y cada proceso calcula =
Siempre podremos hacer una distribución equitativa del númer=
o de
puntos que atiende cada procesador si el número total de puntos Ad=
icionalmente,
cada cálculo de Unidades de Procesamiento Gráfico (GPUs) En=
la
actualidad los computadores personales y de oficina poseen múltiples núcl=
eos o
cores (procesadores) en la llamada Unidad Central de Procesamiento o CPU,
además la tecnologÃa incorpora, a bajo costo, tarjetas gráficas con una =
Unidad
de Procesamiento Gráfico o GPU que también contiene muchos núcleos diseÃ=
±ados
para el cómputo masivo, por lo que pudiéramos aceptar que la computación
paralela está a nuestro alcance. Los procesadores que se encuentran en la GPU pueden ser
administrados eficientemente desde la unidad central de procesamiento o CPU
para realizar cómputo intensivo. La idea es sustituir los procesos
centralizados por procesos distribuidos entre los procesadores (cores) del =
CPU
y los de las GPU. Para este fin se utilizan arquitecturas vectoriales o
paralelas manejadas por plataformas de enlace que permiten consolidar
resultados entre CPU y GPU. En el año 2006 la compañÃa NVIDIA introduce la Arquitectura=
de
Dispositivos de Cómputo Unificado o CUDA (por sus siglas en inglés), y ha=
ce
posible el manejo efectivo de los núcleos de la GPU a través de un modelo=
de
programación que permite la interacción CPU-GPU desde un programa de usua=
rio.
Detalles de esta tecnologÃa, su uso y programación pueden encontrarse en =
la
guÃa de programación de CUDA en NVIDIA (2021). Con estas herramientas tecnológicas el potencial de la tarjeta
gráfica se extiende a la implementación de algoritmos altamente paraleliz=
ables
que consumen mucho tiempo de CPU; tal es el caso de los algoritmos basados =
en
productos de matrices y de vectores, que aparecen frecuentemente en muchas
aplicaciones del álgebra lineal, como el cálculo del espectro y/o el
pseudoespectro de matrices dispersas y de gran magnitud. Algunos investigadores se han dedicado desde entonces a la
implementación de estos algoritmos, véase por ejemplo el trabajo de Angel=
es et
al. (2011). Â Adicionalmente existen
librerÃas de acceso libre con los códigos más usados. Estas librerÃas p=
oseen
interfaces para brindar flexibilidad en la programación y permitir la ejec=
ución
de instrucciones para la CPU y para la GPU en el mismo programa. Como ejemp=
lo,
véase Thrust, librerÃa desarrollada por Bell y
Hoberock (2011), que permite manejar estructuras de datos en algoritmos en
paralelo, y la librerÃa Cusp que contiene módulos para el manejo eficient=
e de
cálculo matricial y la ejecución de métodos clásicos para la resolució=
n de
sistemas lineales, desarrollada por Bell y Garland (2009). El modelo de arquitectura induce a dejar las tareas con cálcu=
los
secuenciales, como el manejo de entrada y salida de datos, para que se ejec=
uten
en la CPU (usualmente llamada host) y las tareas de cómputo masivo paralel=
izable
que se ejecuten sobre la GPU (denominada device). Actualmente se cuenta con
computadores de escritorio de 8 cores con una tarjeta Nvidia de miles de co=
res
y este gap seguirá ampliándose en el futuro. Cada núcleo de la GPU tiene recursos compartidos, como regist=
ros y
memoria, integrados en el mismo chip, lo cual permite que las tareas (hilos)
que se ejecutan en un mismo núcleo compartan datos sin usar un bus de memo=
ria.
La arquitectura permite al programador usar funciones escritas en lenguajes
clásicos de alto nivel como C y C++ (conocidas como kernels) que luego se
ejecutan en paralelo en la GPU como un conjunto de hilos que organizan dent=
ro
una jerarquÃa de bloques. Para la programación en CUDA se debe considerar los siguientes
elementos: jerarquÃa de bloques de hilos, compartición de memoria y
sincronización. Estos elementos dirigen la programación hacia la división del
problema en tareas que pueden ser ejecutadas en paralelo por bloques de hil=
os
independientes, que a su vez podrÃan agruparse en una malla de procesadore=
s (grid), véase figura 6. La sincronización es manejad=
a por el
programador usando herramientas de CUDA que garantizan la escalabilidad; es
decir, la aplicación de usuario debe soportar el incremento de cores en la=
GPU
sin necesidad de reprogramación. Figura
6. Â Malla de bloques de hilos Fuente:=
Adaptada
de NVIDIA (2021). La propuesta de llevar estos cálculos a un esquema hÃbrido q=
ue
involucre el cómputo masivo en la GPU no solo se sustenta en la posibilida=
d de incrementar
el nivel de paralelismo y disminuir los tiempos de respuesta, sino que adem=
ás
el poder de cómputo de los cores en las actuales GPU es muy superior a los=
de
la CPU, véase por ejemplo el trabajo de Arce et. al (2011). Resultados Una vez analizados cada componente involucrado en el cálculo =
del
pseudoespectro, se resume a continuación las bases de la propuesta: 1) Entrada de datos (CPU): 1.1) Lectura o definición de la matriz 1.2) Lectura o definición de los parámetros de
Arnoldi. 1.3) Distribución de la matriz
2) Proyección de la matriz
2.1) Distribución de la matriz
2.2) Ejecución de bloques de operaciones
matriz-vector para aplicación del método de Arnoldi con reinicio implÃci=
to a la
matriz
2.3) Generación de matriz
3) Discretización o definición de malla.
3.1) Lectura o definición de parámetros de la m=
alla.
(CPU).
3.2) Distribución de regiones en los cores la GP=
U.
(GPU)
4) Cálculo del mÃnimo valor singular de los p=
untos
en cada región (GPU).
3. <=
/span>
4. <=
/span>
4.1) Para cada punto
4.2) EnvÃo de la submatriz correspondiente a los=
Â
5) Visualización de contornos (CPU)
5.1) Graficación de curvas de contornos <=
/span>
Conclusiones
· =
Se han descrito los aspectos fundamentales de una propuesta pa=
ra
el cálculo del pseudoespectro de matrices usando un esquema hÃbrido CPU-G=
PU que
añade un nivel adicional de paralelismo con el fin de acelerar este cálcu=
lo que
generalmente resulta costoso computacionalmente. Adicionalmente se han
expuestos métodos que han logrado tener éxito en todas las implementacion=
es
paralelas del pseudoespectro en máquinas de alto rendimiento. Se espera
implementar estos algoritmos en un futuro inmediato bajo una plataforma CUDA
con MPI usando GPU dispuestas en tarjetas NVIDIA en computadores de gama me=
dia
o baja, a fin de probar el potencial de la propuesta.
·=
=
La propuesta conserva caracterÃsticas generales y cada módulo
puede ser implementado o bien utilizando los métodos recomendados en cada
tarea, u otros métodos de preferencia.
El presente documento estuvo dirigido a:
1. Establecer los fundamentos
teóricos para el cálculo del pseudoespectro; cubriendo detalles sobre su =
uso e
importancia.
2. Justificar el uso de la
herramienta y de los métodos de cálculo.
3. Desarrollar una propuesta
general que contemple las distintas etapas del cálculo del pseudoespectro.=
4. Proponer algoritmos bási=
cos
en cada actividad.
5. Describir las caracterÃs=
ticas
de las unidades de procesamiento gráfico, y su utilidad para hacer cálcul=
os en
paralelo.
6. Señalar detalles de
implementación y consideraciones de rigor en base a la arquitectura propue=
sta.
Referencias bibliográficas
Ãngeles, M., Flores, G., Vidal,= A. (2011).  Implementación CPU-GPU y comparativa de las bibliotecas BLAS-CUBLAS, LAPACK-CULA, Universidad Politécnica de Valencia.
Bekas,
C., Kokiopoulou, E., Kouti=
s,
I., Gallopoulus, E. (2001) Towards the effective
Parallel Computation of Matrix Pseudospectra. IC’s 01: Proceeding of=
the
15th International Conference on Supercomputing, June 2001, 203-226.
Bekas,
C., Kokiopoulou, E., Gallo=
poulos
E., Simoncini, V., Parallel Computation of
Pseudospectra using Transfer Functions on a MATLAB-MPI Cluster Platform.
Lecture Notes in Computer Science. DOI:10.1007/3-540-45825-5_35.
Bell,
N., Garland, M. (2009). Implementing sparse matrix-vector multiplication on
throughput-oriented processors. S=
C’09:
Proceeding of the Conference on High Performance Computing Networking. Stor=
age
and Analysis, New York, NY, USA (ACM), 1-11.
Bell,
N., Hoberock, J. (2011). Thrust: A
Productivity-Oriented Library for CUDA. GPU Computing Gems, Jade Edition. <=
span
class=3DSpellE>Editores Wen-mei W. Hwu.
Castillo,
Z. (2004). A new algorithm for continuation and bifurcation analysis of
large-scale free surface flows. P=
hD
thesis, Rice University, Houston, Texas.
Embree, M., Trefethen. L. N., Pseudospectra Ga= teway (2021). [Online]. Disponible en: https://www.cs.ox.ac.uk20/pseudospectra/index.ht= ml.<= o:p>
Guevara, R. (2012). Paralelización de datos para el cál=
culo
del pseudoespectro. Tesis de Licenciatura.
Universidad Central de Venezuela, Escuela de Computación, Caracas, Venezue=
la.
Golub, G., & Van Loan, C. (1996). Matrix Computations (3rd ed.). The Johns
Hopkins University.
Lui,
S. (1997). Computation of pseudos=
pectra
by continuation. SIAM J. Sci.
Mezher, D., Philippe, B. (2002). Parallel computation of pseudospectra of
large sparse matrices. Parallel Computing, 28, 199-221.
Minini, P., Rosenberg, D., Reddy, R., Pouquet,=
A. (2011). A hybrid MPI-OpenMP scheme
for scalable parallel pseudospectral computatio=
ns for
fluid turbulence. Parallel Computing, 37(6-7), 316-312. Noschese, S. & Reichel, L. (2015). Approximated structured pseudospectra=
. Numerical
Linear Algebra and Applications, 00, 1-15. NVIDIA Corporation, NVIDIA CU=
DA C Programming Guide (Version 11.3.0) (2021). [Online]. Disponible en: <=
span
style=3D'color:black;mso-themecolor:text1;background:white;text-decoration:=
none;
text-underline:none'>https://docs.nvidia.com/cuda/cuda-c-programming-guide<=
/span> Otero, B., Astudillo, R., Castillo, Z. (2015). Un esquema
paralelo para el cálculo del pseudoespectro de
matrices de gran magnitud. Revista Internacional de Métodos Numéricos =
para el
Cálculo y Diseño en IngenierÃa. 31-1, 8-12.<=
o:p> Sorensen D. C. (1997). Implicitly
restarted Arnoldi/Lanczos<=
/span>
methods for large scale eigenvalue calculations, Parallel Numerical
Algorithms, Springer, Dordrecht, 119-165. Trefethen, L. N. & Bau III, D. (1997). N=
umerical
linear algebra. Siam, Philadelphia. Trefethen, L.N. & Embree,
M. (2005). Spectra and Pseudospectra: the behavior of =
nonnormal
Matrices and Operators, Princeton University Press, New Jersey, USA. Trefethen L.N & Wright, T. (2001). Large-scale computation of pseudospect=
ra
using ARPACK and eigs, SIAM J. Sci. Comput. 23 (2001), 591–605. Trefethen, L.N. (1999). =
Computation of pseudospectra. Acta Numerica 8,
1, 247–295. PARA CITAR EL ARTÃCULO INDEXADO. Castillo Marrero, Z. N., Colmenares Pacheco, G. A., Valverde Aguirre=
, P.
E., & Cevallos Vique, V. O. (2021). Una pro=
puesta
para el cálculo del pseudoespectro en unidades=
de
procesamiento gráfico. ConcienciaDigital, 4(2.=
1),
6-20. ht=
tps://doi.org/10.33262/concienciadigital.v4i2.1.1703 El artÃculo que se
publica es de exclusiva responsabilidad de los autores y no necesariamente
reflejan el pensamiento de la Revi=
sta Conciencia
Digital. El
artÃculo queda en propiedad de la revista y, por tanto, su publicación pa=
rcial
y/o total en otro medio tiene que ser autorizado por el director de la Revista Conciencia Digital.
[1] Escuela
Superior Politécnica de Chimborazo, Facultad de Ciencias, Escuela de
EstadÃstica, Grupo CIDED. Riobamba, Ecuador, zenaida.castillo=
@espoch.edu.ec, =
https://orcid.org/0000-0002-4424-8652<=
/span>
[2] Universidad de
Investigación Experimental Yachay, Escuela de Matemáticas y Ciencias
Computacionales, UrcuquÃ, Ecuador, gcolmenares@yach=
aytech.edu.ec, https://orcid.org=
/0000-0003-4789-0859
[3] Escuela Superior
Politécnica de Chimborazo, Facultad de Ciencias, Grupo CIDED. Riobamba,
Ecuador, paulina.valverde@espoch.edu.ec<=
span
lang=3DES style=3D'font-family:"Times New Roman",serif'>, https://orcid.org=
/0000-0003-0458-7083
[4] Escuela Superior
Politécnica de Chimborazo, Facultad de Administración de Empresas, Escuel=
a de
Finanzas, Grupo CIDED, Riobamba, Ecuador, victor.cevallos@=
espoch.edu.ec, https://orcid.org=
/0000-0001-5525-5818
[5] La n=
orma de (ATA
– AAT) proporciona una medida de la no-normalidad.
www.concienciadigital.org
                        Â=
                                     =
                             Vol.
4, N°2.1, p. 6-20, mayo, 20