Printer Friendly

HPC en simulacion y control a gran escala.

1. Introduccion

Simular y particularmente controlar un fenomeno fisico, tiene un elevado costo computacional, especialmente, cuando el proceso inmerso esta modelado mediante ecuaciones en derivadas parciales (EDPs). Los requerimientos de memoria y tiempo de ejecucion estan ligados estrictamente a la discretizacion del dominio, esto es, a los puntos en los cuales se aproxima la EDP. Esta discretizacion se puede realizar por diferentes metodos como elementos finitos o diferencias finitas. Por esta razon el costo computacional puede crecer, incluso exponencialmente, si se requiere una alta precision para resolver el problema numericamente. Por otro lado, en el proceso de produccion de circuitos integrados de las ultimas decadas, el tamano de los dispositivos se ha reducido y la velocidad de los mismos se ha incrementado en una relacion inversamente proporcional; asi lo afirmo Gordon E. Moore, co-fundador de Intel, el 19 de abril de 1965 en su conocida Ley de Moore (esta ley empirica expresa que, aproximadamente, cada 18 meses se duplica el numero de transistores en un circuito integrado). Es importante senalar que la simulacion numerica ha cumplido, y lo sigue haciendo, un rol fundamental en el diseno de estos dispositivos; ya que permite determinar con alta precision el comportamiento y desempeno del dispositivo.

El sistema dinamico que describe el circuito presenta una gran escala y su resolucion numerica es un reto para la comunidad cientifica, y constituye incluso un problema abierto si la dimension es muy grande y/o si el problema es denso no estructurado. Por la amplia gama de aplicaciones, los sistemas con caracteristicas similares aparecen en general en problemas gobernados por ecuaciones diferenciales parciales, y las numerosas tecnicas que han sido propuestas para tratar estos problemas persiguen reducir su elevado costo computacional [4,5,15,75].

Si el costo computacional esta altamente relacionado con la dimension del sistema, parece razonable y necesario reducir la dimension del mismo para resolverlo en tiempo real. Los metodos que explotan esta idea son conocidos como reduccion de modelos o reduccion de orden. Estas tecnicas han probado ser particularmente eficientes en la simulacion de circuitos electricos [24,39,56,67,74]. En general estos metodos se clasifican en:

* Metodos de proyeccion

* Metodos Moment Matching

* Metodos de truncamiento balanceado

En el presente trabajo nos centraremos en metodos de truncamiento balanceado, particularmente en la resolucion numerica de ecuaciones de Lyapunov que caracterizan esta tecnica. Tambien revisaremos brevemente las ecuaciones matriciales cuadraticas que aparecen en control optimo y cuya resolucion numerica requiere el tratamiento numerico de ecuaciones de Lyapunov y la aplicacion de tecnicas de computacion de alto desempeno. Es importante senalar que existen muchas variantes del metodo de truncamiento balanceado que emplean ecuaciones algebraicas de Riccati [16]; estas ultimas no se revisaran en este trabajo.

El articulo esta organizado de la siguiente manera: en la Seccion 2 se describe el sistema dinamico relacionado con el modelo matematico; en la Seccion 2.1, se introduce la reduccion de modelos y la estructura de cierto tipo de ecuaciones matriciales asociadas a esta tecnica. En la Seccion 2.2 se describen las ecuaciones matriciales que aparecen en el control lineal cuadratico. Los metodos de resolucion para estas ecuaciones matriciales asi como sus algoritmos se ilustran en la Seccion 3, mientras que en la Seccion 4 se describen tecnicas de HPC para implementar eficientemente los algoritmos descritos en la seccion anterior. Finalmente, en la Seccion 5 se presentan algunos ejemplos numericos, con el proposito de comprobar la precision y factibilidad de las tecnicas propuestas en el presente trabajo. Algunas conclusiones se aportan en la Seccion 6.

2. Modelo matematico

Tras linealizar un sistema no lineal o un sistema de mayor orden, se obtiene un sistema lineal no dependiente del tiempo de la forma:

[??](t) = Ax(t) + Bu(t), t > 0,

x(0) = xo,

y(t) = Cx(t) + Du(t), t> 0, (1)

donde A [elemento de] [R.sup.nxn], B [elemento de] [R.sup.nxm], C [elemento de] [R.sup.pxn], y D [elemento de] [R.sup.pxn] son matrices constantes. En el tipo de problemas que nos interesa la dimension del sistema n, en general, varia entre [10.sup.3] - [10.sup.7], dependiendo de la aplicacion especifica que se considere. En cambio, el numero de entradas o controles, representado por u(t), y el numero de salidas del sistema, representado por y(t), a menudo es mucho menor que el orden del sistema, esto es, m,p [mucho menor que] n.

Por ejemplo, si se aplica analisis nodal modificado para generar un sistema de ecuaciones para una red RLC multi-puerto, el sistema tiene la forma siguiente:

[Laplace][??](t) = -Gx(t) + Bu(t) y(t) = Cx(t) (2)

donde

[EXPRESION MATEMATICA IRREPRODUCIBLE EN ASCII].

Los vectores y y u representan las corrientes de los puertos y sus voltajes, respectivamente. Las matrices E, L, G son simetricas y denotan la capacitancia nodal, inductancia y conductividad, respectivamente. M representa la matriz de incidencia asociada con el inductor de corriente. N indica donde se encuentran las fuentes de voltaje, y la matriz C especifica en que punto se miden los puertos de corriente. En este ejemplo, la matriz [Laplace] es no singular, esto es, existe su inversa, por lo cual podemos escribir el sistema (2) en la forma (1). Si la matriz [Laplace] es singular, (2) corresponde a una ecuacion diferencial algebraica, [78].

Existen muchos metodos para resolver numericamente sistemas a gran escala de la forma (1) [15,28]. En la siguiente seccion se introducen los metodos de reduccion de modelos basados en truncamiento balanceado, asi como su relacion con el problema de control lineal cuadratico.

2.1. Reduccion de modelos

El objetivo de la reduccion de modelos (RM) es encontrar un sistema cuyo orden sea mucho menor que el sistema original, y al mismo tiempo este nuevo sistema sea una buena aproximacion del sistema original. Concretamente, se busca para la misma secuencia de entradas, u(t), la salida del sistema, y(t), represente de buena manera lo que ocurre en el sistema original; en otras palabras, se busca un sistema de la forma:

[EXPRESION MATEMATICA IRREPRODUCIBLE EN ASCII], (3)

cuyo orden, r, sea mucho menor que el orden del sistema original, r [mucho menor que] n, y adicionalmente [??] sea una buena aproximacion, de acuerdo a cierto criterio, de la salida del sistema original y.

Entre los metodos de RM, se destacan los de reduccion/substraccion-Guyan, truncamiento modal, aproximaciones de Pade, metodos de subespacio de Krylov, truncamiento balanceado, etc. Una caracteristica de muchos de estos metodos es realizar proyecciones de tipo Galerkin o Petrov-Galerkin del espacio de estado en un subespacio de dimension menor. Para saber mas sobre reduccion de modelos, ver [4,5].

En este trabajo abordaremos los metodos de truncamiento balanceado. En ingenieria uno de los criterios mas utilizados para medir la aproximacion de [??] a y es la denominada funcion de transferencia. Esto se debe a que, bajo ciertas suposiciones, se puede probar que:

[EXPRESION MATEMATICA IRREPRODUCIBLE EN ASCII], (4)

donde

[EXPRESION MATEMATICA IRREPRODUCIBLE EN ASCII]

son las funciones de transferencia de los sistemas (1) y (3), respectivamente, y [[paralelo]x[paralelo].sub.H] denota la norma en el espacio H.[paralelo]x[paralelo] h denota la norma en el espacio H.

Debido a la relacion (4), muchas tecnicas de RM tratan de minimizar [paralelo]G - [??][paralelo][infinito], pues con ello implicitamente estan minimizando [EXPRESION MATEMATICA IRREPRODUCIBLE EN ASCII]; este es precisamente el caso del truncamiento balanceado. Otros criterios de minimizacion consideran modelos de frecuencia pesados y reduccion de controles [85]. El truncamiento balanceado se fundamenta en separar el espacio de transformaciones de estado mediante una aplicacion no singular de la forma:

[EXPRESION MATEMATICA IRREPRODUCIBLE EN ASCII]

donde las matrices del modelo reducido se definen como

[EXPRESION MATEMATICA IRREPRODUCIBLE EN ASCII]

La aplicacion T balancea el Gramian de controlabilidad, Wc, y el Gramian de observabilidad, Wo, que corresponden a soluciones de ecuaciones matriciales de Lyapunov de la forma:

AWc + Wc[A.sup.T] + [BB.sup.T] = 0,

AWo + Wo[A.sup.T] + [CC.sup.T] = 0. (5)

Antes de introducir algunos metodos para resolver (5), en los que se profundiza en las Secciones 3 y 4, es importante senalar que, como Wc y Wo son matrices semidefinidas positivas, se pueden factorizar de la forma Wc = [S.sup.T]S, Wo = RTR. Ademas, si se escogen los factores S, R [elemento de] [R.sup.nxn] triangulares, estos corresponden a los factores de Cholesky de Wc y Wo, respectivamente. Esta propiedad permite definir el truncamiento balanceado usando el producto [SR.sup.T] en lugar del producto de los Gramians. Algunas implementaciones de algoritmos de reduccion de modelos y truncamiento balanceado, se describen en [15,85].

El costo computacional de la RM, via truncamiento balanceado, depende principalmente de la resolucion de las ecuaciones de Lyapunov asociadas a los Gramians de observabilidad y controlabilidad.

En sistemas no estables, la aplicacion del truncamiento balanceado conlleva la necesidad de resolver ecuaciones algebraicas de Riccati, en lugar de ecuaciones de Lyapunov. Las ecuaciones de Riccati aparecen en muchas aplicaciones de ciencia e ingenieria, especialmente en la teoria de control. En la siguiente seccion se describe brevemente el problema de control lineal cuadratico linear-quadratic regulator (LQR) y su relacion con las ecuaciones de Riccati.

2.2. Ecuaciones matriciales en control

Los problemas de optimizacion gobernados por EDPs de tipo parabolico pueden formularse como problemas de Cauchy abstractos. Imponiendo un funcional de costo cuadratico, se obtiene un problema LQR para un sistema de dimension infinita. El control optimo para este problema esta definido via feedback en terminos de la solucion de la ecuacion de Riccati para operadores, ya sea esta de tipo algebraico o diferencial, dependiendo de si el horizonte de integracion que se considera es infinito o finito, respectivamente. Existen varios resultados de convergencia respecto al problema infinito-dimensional y las ecuaciones de Riccati asociadas, ver por ejemplo [65]. Estos resultados nos brindan el marco teorico necesario para resolver el problema numericamente.

Consideremos entonces el problema LQR semidiscretizado:

[EXPRESION MATEMATICA IRREPRODUCIBLE EN ASCII], (6)

con respecto a (1). Por simplicidad tomemos D = 0, Q = I, y R = I. Luego el control optimo esta dado por:

u(t) = -[B.sup.T] [X.sub.[infinito]x](t),

donde [X.sub.[infinito]] es la solucion de la ecuacion algebraica de Riccati

0 = [R.sub.h](X) = [C.sup.T]C + [A.sup.T]X + XA - [XBB.sup.T]X, (7)

o la solucion de la ecuacion diferencial de Riccati

[??] = -[K.sub.h](X) = -[C.sup.T]C - [A.sup.T]X - XA + [XBB.sup.T]X, (8)

si [T.sub.f] = [infinito] y G = 0 o [T.sub.f] < [infinito], respectivamente.

El orden de las ecuaciones matriciales (7) y (8) es cuadratico respecto al orden del sistema, esto es, su solucion implica resolver un sistema de EDOs de orden [n.sup.2]. Por esta razon, en problemas a gran escala, la aplicacion de tecnicas tradicionales no es factible. Si se tiene en cuenta que tipicamente, los coeficientes de (7) y (8) tienen una estructura definida (p.e., dispersion, simetria o rango bajo) es posible definir metodos numericos capaces de explotar eficientemente esta estructura. La relacion entre estas ecuaciones y la ecuacion de Lyapunov, asi como los metodos numericos utilizados se describen en la siguiente seccion.

Observacion 1 En sistemas no estables, la aplicacion de truncamiento balanceado de tipo linear-quadratic gaussian LQG conlleva la necesidad de resolver ecuaciones algebraicas de Riccati en lugar de ecuaciones de Lyapunov, puesto que en este caso particular se deben considerar Gramians de lazo cerrado en el compensador LQG.

3. Algoritmos

Como se ha planteado en la Seccion 2.1, la aplicacion de RM -via truncamiento balanceado a sistemas lineales de gran escala- exige resolver ecuaciones de Lyapunov o ecuaciones de Riccati, dependiendo de si el sistema es estable o no. Estas ultimas ecuaciones aparecen tambien en control optimo, y estan estrictamente ligadas unas a otras, como veremos mas adelante.

La estructura de la ecuacion (7) es compartida por la ecuacion (5) excepto por el termino no lineal [XBB.sup.T]X. Lo mismo ocurre con el lado derecho de la ecuacion (8).

Observemos que la derivada de Frechet (es una generalizacion de la derivada usual a espacios de Banach [61]) del operador iKh en P viene dada por el operador de Lyapunov:

[EXPRESION MATEMATICA IRREPRODUCIBLE EN ASCII];

aplicando el metodo de Newton se tiene que

[EXPRESION MATEMATICA IRREPRODUCIBLE EN ASCII].

Y, por tanto, para una matriz inicial dada, este metodo puede implementarse como se recoge en el Algoritmo 3.1:
Algoritmo 3.1 Metodo de Newton para EAR

Require: [P.sub.0], tal que [A.sub.0] es estable

1: [A.sub.l] [flecha siniestra] [A.sub.h] - [B.sub.h][B.sup.T.sub.h]
[P.sub.l]

2: Resolver la ecuacion de Lyapunov
   [A.sup.T.sub.l] + [N.sub.l][A.sub.l] = -[R.sub.h]([P.sub.l])

3: [P.sub.l+1] [flecha siniestra] [P.sub.l] + [N.sub.l]


En las aplicaciones que consideramos, es posible enunciar que, para intervalos pequenos, la solucion aproximada de [X.sub.k] [aproximadamente igual a] X([t.sub.k]), en general, es un valor estabilizado muy bueno para empezar la iteracion; ver [29]. Por lo tanto, toda [A.sub.l] es estable y la iteracion [P.sub.l] converge a la solucion, [P.sub.*], en orden cuadratico.

El metodo de Newton para ecuaciones algebraicas de Riccati (EAR) puede ser reformulado como una iteracion de un solo paso, si se reescribe de tal manera que la siguiente iteracion sea calculada directamente de la ecuacion de Lyapunov como se muestra en el paso 2 del Algoritmo 3.1, esto es:

[EXPRESION MATEMATICA IRREPRODUCIBLE EN ASCII].

De esta manera, para resolver la EAR, es necesario resolver una ecuacion de Lyapunov en cada iteracion, de la forma:

[F.sup.T] X + XF = -[WW.sup.T], (9)

siendo F una matriz estable en cada iteracion. Usualmente, el metodo de Newton converge muy rapido; por ejemplo, en el problema de transferencia de calor en dos dimensiones descrito en la Seccion 5, la convergencia se alcanza de 3 a 5 iteraciones.

Para resolver ecuaciones de Lyapunov, dispersas a gran escala, usaremos el metodo ADI (Alternating Direction Implicit) que se puede formular como: (ver [87])

[EXPRESION MATEMATICA IRREPRODUCIBLE EN ASCII],

donde p denota el conjugado de p [elemento de] C. Si los parametros [p.sub.j] se escogen apropiadamente, entonces [lim.sub.[flecha diestra][infinito]] [Q.sub.j] = Q converge en orden super-lineal.

Sin embargo, en problemas a gran escala, no es posible aplicar directamente el metodo ADI; incluso no es posible almacenar en memoria las matrices que aparecen en cada paso, pues lo anterior excede facilmente la capacidad de memoria de un computador. Por esta razon, es necesario utilizar tecnicas de rango bajo para su resolucion. Concretamente utilizaremos el metodo de rango bajo Newton-ADI presentado en [25,74].
Algoritmo 3.2 Metodo de rango bajo ADI

Require: F, W y el conjunto de parametros ADI {[p.sub.i], ..., [p.sub.k]}

Ensure: [EXPRESION MATEMATICA IRREPRODUCIBLE EN ASCII]

1: [V.sub.1] = [raiz cuadrada de (-2Re ([p.sub.1]))][([F.sup.T] +
[p.sub.1]I)-1.sup.W]

2: [Z.sub.1] = [V.sub.1]

3: for j = 2, 3, ... do

4: [EXPRESION MATEMATICA IRREPRODUCIBLE EN ASCII]

5: [Z.sub.j] = [[Z.sub.j-1] [V.sub.j]]

6: end for


En el Algoritmo 3.2, Z es una matriz de rango bajo que aproxima X, es decir, Z [elemento de] [R.sup.nxq] con q [mucho menor que] n. La aproximacion de la matriz Z es muy precisa. De hecho, no se requiere un costo computacional adicional para conseguir resultados similares a los que se obtienen con metodos convencionales.

Por otro lado, cabe senalar que todas las matrices [V.sub.j] tienen igual numero de columnas que [EXPRESION MATEMATICA IRREPRODUCIBLE EN ASCII], es decir, en cada iteracion j, es necesario resolver w sistemas de ecuaciones lineales con los mismos coeficientes matriciales [F.sup.T] + [p.sub.j] I. En la implementacion del algoritmo se explota esta propiedad. Si la convergencia del metodo de rango bajo ADI se alcanza con respecto a cierto criterio de parada despues imax pasos, entonces [EXPRESION MATEMATICA IRREPRODUCIBLE EN ASCII]. Cuando n es grande y [n.sub.w] es pequeno, se espera que [EXPRESION MATEMATICA IRREPRODUCIBLE EN ASCII]. En este caso, se ha calculado una aproximacion de rango bajo [EXPRESION MATEMATICA IRREPRODUCIBLE EN ASCII] del factor Z de la solucion X, esto es [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Sin embargo, si [n.sub.w] x [i.sub.max] crece demasiado, el metodo pierde eficiencia. Para controlar el crecimiento de [n.sub.w] x [i.sub.max], pueden aplicarse tecnicas de compresion de columnas que reducen el numero de columnas en [EXPRESION MATEMATICA IRREPRODUCIBLE EN ASCII], sin introducir un aumento del error significativo; ver [31,51].

Si combinamos el Algoritmo 3.2 con el metodo de Newton (Algoritmo 3.1), se puede definir un metodo de rango bajo para resolver ecuaciones algebraicas de Riccati denominado Newton-ADI.

Por otro lado, las ecuaciones diferenciales de Riccati aparecen en reduccion de modelos cuando se consideran sistemas no estables dependientes del tiempo. La resolucion numerica de estas ecuaciones constituyen un reto mayor al de resolver ecuaciones de Lyapunov o ecuaciones algebraicas de Riccati. Mas aun, se espera que las EDRs sean rigidas (stiff), por lo que se deben utilizar metodos que puedan tratar este fenomeno eficientemente. Tambien es necesario explotar la estructura de las matrices de coeficientes, y la estructura de la ecuacion en si, lo que implica aplicar metodos numericos que preserven la estructura de la EDR. Entre las tecnicas, de un paso y paso multiple, para resolver sistemas de ecuaciones diferenciales ordinarias rigidas, los metodos BDF (Backward Differentiation Formulae) y los metodos de tipo Rosenbrock son usados comunmente. Versiones matriciales de estos algoritmos aplicables a EDRs a gran escala, se han propuesto en [29,71]. La idea fundamental en estas tecnicas es realizar una aproximacion de rango bajo a la solucion de la EDR y resolver la ecuacion para los factores de la solucion. Cabe senalar, que el control de paso y orden se realizan tambien directamente sobre los factores.

3.1. Caso no disperso

Cuando las matrices son densas, utilizamos el metodo de la funcion de signo de la matriz (Matrix Sign Function) para resolver las ecuaciones de Lyapunov que aparecen en los problemas de reduccion de modelos y de control. En la literatura existen varios metodos iterativos para calcular la funcion de signo de una matriz. Entre ellos, la iteracion de Newton es una opcion competitiva por su simplicidad, eficiencia, inherente paralelismo y convergencia asintoticamente cuadratica [53].

El Algoritmo 3.3 describe este metodo. El costo computacional es aproximadamente de [2n.sup.3] operaciones aritmeticas de punto flotante por iteracion.
Algoritmo 3.3 Funcion de matriz signo para ecuaciones de Lyapunov

Require: F [elemento de] [R.sup.nxn], W [elemento de] [R.sup.jxn]

Ensure: [??] such that [EXPRESION MATEMATICA IRREPRODUCIBLE EN ASCII] and
[F.sup.T] X + XF = [WW.sup.T]

1: [F.sub.0] := F, [[??].sub.0] := [W.sup.T]

2: k := 0

3: repeat

4: [F.sub.k+1] := 1/[raiz cuadrada de 2] ([F.sub.k]/[c.sub.k] +
[C.sub.k] [F.sub.k.sup.-1]).

Compute the rank-revealing QR (RRQR) decomposition

5: [EXPRESION MATEMATICA IRREPRODUCIBLE EN ASCII]

6: [S.sub.k+1] := [U.sub.s][[PI].sub.s]

7: k := k + 1

8: until [EXPRESION MATEMATICA IRREPRODUCIBLE EN ASCII]


En el Algoritmo 3.3, S representa el factor de la solucion aproximada, es decir, X [aproximadamente igual a] [S.sup.T]S. Es importante tener en cuenta que la convergencia del metodo puede acelerarse mediante algunas tecnicas [53]. En nuestra aproximacion empleamos un escalado definido por el parametro:

[c.sub.k] = [raiz cuadrada de [paralelo][F.sub.k.sup.-1][paralelo][infinito/[paralelo][F.sub.k][paralelo][infinito]]

En la prueba de convergencia, [tau] es la tolerancia para la iteracion, y se determina usualmente como funcion de la dimension del problema y de la precision de la computadora [elemento de].

En la siguiente seccion, se describira una implementacion eficiente de este algoritmo sobre plataformas multinucleo equipadas con varias GPUs. Esta implementacion incorpora numerosas tecnicas de HPC y alcanza un desempeno notable.

4. Computacion de alto desempeno

La alta exigencia computacional de los algoritmos mostrados en la seccion anterior, precisan el uso de tecnicas de computacion y hardware de alto desempeno (HPC).

La utilizacion de este tipo de tecnicas y de hardware presenta como principal inconveniente los altos costos economicos asociados tanto por la adquisicion de hardware y software, como por otros costos derivados; por ejemplo la administracion y mantenimiento del equipamiento o la energia consumida. Sin embargo, en los ultimos anos, han irrumpido nuevas tecnologias de HPC de bajo costo, entre las que se destacan los procesadores graficos, o GPUs. La exigencia ejercida por la industria de los videojuegos sobre esta tecnologia, provoca una constante evolucion de la capacidad de computo de las GPUs a la vez que mantiene su costo en unos niveles moderados. Por todo ello, las GPUs se han posicionado como dispositivos de referencia en computacion de alto desempeno.

4.1. Software de alto desempeno en problemas de algebra lineal

Como quedo explicito en la seccion anterior, la resolucion de los problemas afrontados en este articulo se encuentra estrechamente ligada a la resolucion de problemas de algebra lineal numerica (ALN).

Existen numerosas bibliotecas de alto desempeno que incluyen en su funcionalidad operaciones de ALN. Entre ellas, destaca la especificacion BLAS (Basic Linear Algebra Subprograms), que da soporte a las operaciones basicas de ALN sobre vectores y matrices. Existen numerosas implementaciones de BLAS, la mayor parte de ellas especificas para una plataforma hardware. Es habitual que los propios fabricantes de hardware desarrollen versiones de BLAS especificas para cada uno de sus productos. Con el fin de replicar este caso, se han desarrollado diversos esfuerzos similares, como por ejemplo LAPACK [40] y ScaLAPACK [32], para operaciones mas complejas de ALN.

La utilizacion de estas bibliotecas es fuertemente recomendable por la alta eficiencia que proporcionan sus rutinas, asi como por la simplicidad y legibilidad del codigo o la portabilidad sin perdida de desempeno.

La especificacion LAPACK (Linear Algebra PACKage), incluye soporte para matrices densas, simetricas, banda y triangulares para diferentes operaciones, como la resolucion de sistemas de ecuaciones lineales, la aproximacion de minimos cuadrados y la busqueda de valores propios y de valores singulares. Habitualmente, las implementaciones de LAPACK estan basadas en el uso de rutinas de BLAS, facilitando su portabilidad y adaptacion a diferentes arquitecturas sin afectar negativamente su desempeno.

El aumento en la utilizacion de computadoras paralelas de bajo costo, y en especial, de aquellas con memoria distribuida, motivo la evolucion de LAPACK a una nueva version basada en tecnicas de programacion paralela distribuida; esta nueva biblioteca se denomino ScaLAPACK y la primera version fue publicada en el ano 1995.

Una de las principales caracteristicas de la biblioteca es que su diseno especifica distintos niveles de abstraccion, buscando encapsular diferentes tareas en distintos componentes.

Al igual que LAPACK, ScaLAPACK posee funciones para las operaciones mas comunes en el campo del ALN, tales como la factorizacion LU y la factorizacion QR. Sin embargo, ScaLAPACK esta disenada para su ejecucion en plataformas de memoria distribuida, en las que tanto los datos (matrices y vectores) como los computos se distribuyen entre las diferentes unidades de procesamiento.

Al disenar la biblioteca ScaLAPACK se decidio encapsular las operaciones sobre vectores que se debian realizar de forma distribuida, para lo cual se creo la especificacion de la biblioteca PBLAS (Parallel BLAS) [37]. La especificacion de PBLAS posee las mismas operaciones que dispone BLAS, pero en este caso las operaciones estan disenadas para ser ejecutadas en paralelo sobre plataformas de memoria distribuida mediante el paradigma de pasaje de mensajes (MPI, PVM, etc.) y la invocacion a funciones de la biblioteca de comunicaciones BLACS.

Existen multiples bibliotecas de caracteristicas similares a las presentadas, asi como esfuerzos exitosos para desarrollar soluciones a problemas mas especificos. En el repositorio NetLib [1] se encuentra un buen compendio del trabajo en el area, con una amplia gama de bibliotecas, algoritmos y articulos cientificos en este ambito.

4.2. Aplicacion de GPUs a computacion de proposito general

En los primeros anos del desarrollo de GPUs para computacion de proposito general, el crecimiento de la capacidad computacional del hardware no fue acompanado por un avance en el software de manejo de las GPUs. La programacion de GPUs se realizaba mediante herramientas desarrolladas para el tratamiento de imagenes, como por ejemplo las APIs graficas (OpenGL y DIRECTX). Esto exigia para el desarrollador la necesidad de conocer profundamente los topicos de computacion grafica, para poder, de esta forma, mapear el problema a resolver en conceptos de computacion grafica. Posteriormente, se empezaron a utilizar diferentes lenguajes de tipo ensamblador especificos de cada modelo de GPU, lo que implicaba la existencia de varios lenguajes de desarrollo y una baja portabilidad de las aplicaciones.

Para tratar de paliar esta carencia, se presentaron varios lenguajes de alto nivel, que permitian generar codigos para diferentes modelos de GPUs. No obstante, cada herramienta seguia siendo muy dependiente de la arquitectura de la GPU, el modelo, etc. En el ano 2007, con el proposito de dar respuesta a estas carencias NVIDIA presento CUDA.

Este salto importante en el software se desarrollo en forma conjunta con un cambio radical en la arquitectura de las GPUs (de NVIDIA), ya que, junto con CUDA, NVIDIA presento la arquitectura G80 que, como principal caracteristica, presenta la arquitectura unificada sin distincion entre procesadores de pixeles y de vertices. CUDA se ha convertido en un estandar valido para todas las GPUs desarrolladas por NVIDIA.

CUDA es una arquitectura de software para el computo de proposito general enfocada hacia el calculo masivamente paralelo en tarjetas de NVIDIA. CUDA se compone de una pila de capas de software: un controlador de hardware, una interfaz de programacion de aplicaciones y bibliotecas matematicas de alto nivel, (inicialmente dos, CUFFT y CUBLAS, siendo esta ultima una implementacion de BLAS en GPU). El driver CUDA se encarga de la transferencia de datos entre la GPU y la CPU. Conceptualmente, la arquitectura de CUDA se construye alrededor de una matriz escalable de multiprocesadores (SMs). Un multiprocesador en las GPU actuales (arquitectura G80) consiste de ocho procesadores escalares (SP), asi como de elementos adicionales como la unidad de instrucciones multihilo y un chip de memoria compartida. Los multiprocesadores crean, gestionan y ejecutan hilos concurrentemente sin apenas costo de planificacion. Para lograr explotar de manera optimizada los recursos de la GPU, CUDA provee un nuevo conjunto de instrucciones accesibles a traves de lenguajes de alto nivel como son Fortran y C.

Esta arquitectura permite abstraer la GPU como un conjunto de multiprocesadores compuestos a su vez por un conjunto de procesadores orientados a la ejecucion de hilos. La nocion fundamental en el nuevo modelo de ejecucion es que cada multiprocesador ejecutara en cada uno de sus procesadores el mismo conjunto de instrucciones, pero sobre distintos datos, es decir, paralelismo bajo el paradigma de programacion SPMD. En particular, y debido a la gran cantidad de hilos que manejan, es comun referirse al paralelismo ofrecido por las GPUs como SIMT (del ingles Single Instruction Multiple Threads). Los programas a ejecutar en la GPU se denominan kernels, y se organizan en un grid tridimensional de hilos de ejecucion (coordenadas x, y y z). A su vez, los hilos se agrupan en bloques tridimensionales; estas coordenadas se utilizaran para identificar el conjunto de datos sobre el que se quiere que actue cada hilo. Cada bloque sera ejecutado en un multiprocesador en forma independiente (esto es, no puede haber comunicacion de datos ni sincronizacion entre los bloques) y cada hilo en un procesador de hilo (entre estos hilos si puede haber sincronizacion y comunicacion). Al conjunto de hilos que se ejecutan en el mismo instante de tiempo en un mismo bloque se le denomina warp y su tamano es 32. Organizar la ejecucion de esta manera permite escalar sin necesidad de recompilar el programa.

El acceso a memoria por parte de los hilos es un componente fundamental en el desempeno de un programa sobre CUDA. Se cuenta con una jerarquia de memoria en la cual los accesos deben ser definidos explicitamente (salvo a los registros y la memoria local). La memoria global es comun a todos los bloques de hilos que se ejecutan dentro de una tarjeta grafica; es la memoria de mayor capacidad y los datos almacenados pueden persistir durante la ejecucion de multiples kernels. Es la unica memoria visible desde el host y permite la escritura desde GPU. La memoria compartida es una memoria a la cual solo pueden acceder los hilos que corren dentro un mismo bloque. Su capacidad es menor a la memoria global y su ciclo de vida es de una ejecucion de kernel. Los registros son privados de cada hilo y es la memoria mas rapida pero su tamano es muy limitado y su uso es determinado por el compilador. La memoria local es privada a cada hilo, tiene la misma velocidad que la memoria global y su ciclo de vida es de una ejecucion de kernel. Tambien se dispone de la memoria constante y la memoria de texturas.

Utilizacion de GPUs para la resolucion de problemas de proposito general. La constante evolucion de las GPUs, su precio reducido y su alto desempeno han propiciado que el uso de GPUs en HPC sea en la actualidad una importante area de investigacion que cuenta ya con numerosos y relevantes avances. A continuacion se detallan algunos trabajos relacionados con el topico de este articulo, con interes particular en ALN.

Como se mostro en la Seccion 4.1, en ALN es habitual el uso de la especificacion BLAS para desarrollar aplicaciones. En este sentido, gran parte del esfuerzo inicial para acelerar nucleos de ALN con GPU, hasta la presentacion en 2007 de CUBLAS por parte de NVIDIA, consistio en la implementacion de operaciones incluidas en la especificacion BLAS. En esta categoria destacan trabajos pioneros dedicados a la multiplicacion de matrices de Hall et al. [52] y Larsen et al. [64]; el de Bajaj et al. [8] para resolver expresiones lineales; y el articulo de Kruger et al. [63], primer trabajo en delinear el desarrollo de nucleos al estilo BLAS; donde los vectores son almacenados en memoria de texturas 2D, organizando las matrices completas por diagonales. Tambien afronta el caso en que las matrices son dispersas utilizando la misma idea para matrices de banda y empleando un formato comprimido para matrices dispersas no estructuradas. Implementan los metodos del GC y Gauss-Seidel. Posteriormente, en 2005, los mismos autores extienden el trabajo [62] utilizando DirectX 9.0 y cambian el formato de almacenamiento de las matrices dispersas no estructuradas por una estrategia comprimida a bloques. Otro autor que estudio de forma temprana el desarrollo de nucleos tipo BLAS sobre GPU es Moravansky [72]. Tras la aparicion de la primera version de la biblioteca CUBLAS, muchos trabajos se centraron en presentar mejoras sobre las rutinas desarrolladas por NVIDIA, algunos de los cuales han sido incluidos en posteriores versiones de la biblioteca. En este sentido, destaca el trabajo de Barrachina et al. [10] donde los autores evaluan el desempeno de rutinas de CUBLAS, y posteriormente muestran como conseguir mejoras con base en diferentes estrategias, inteligentes particionamientos, implementaciones hibridas de la multiplicacion de matrices, modificaciones de la estructura del nucleo de resolucion de sistemas lineales triangulares, y en especial la inclusion de estrategias de padding. Tambien Volkov y Demmel [86] presentan estrategias para la aceleracion de la multiplicacion de matrices. Los autores validan sus desarrollos acelerando los metodos de factorizacion LU y QR; su rutina fue incluida en la version 2.1 de CUBLAS. Con base en el trabajo de Volkov y Demmel, L. Chien [36] mostro como alcanzar mejoras (entre un 10 % y un 15 %) sobre la arquitectura Tesla. Posteriormente, Nath et al. [76] presentaron una version optimizada de la multiplicacion de matrices para la arquitectura Fermi; esta propuesta fue incluida en la version 3.2 de CUBLAS. Otro trabajo interesante es el de Fatica [45], que extendio estas ideas para desarrollar una version del benchmark LINPACK que permite incluir GPUs como elementos de computo. Gran parte del avance en lo que respecta a la implementacion de BLAS sobre GPU, se compendia en el Capitulo 4 del libro de Nath et al. [73].

Operaciones relacionadas con LAPACK, como la factorizacion LU, tambien han sido abordadas por diversos autores. En el trabajo de Galoppo et al. [48] se puede encontrar una implementacion pionera de metodos para la resolucion de sistemas lineales densos; en particular se presentan los metodos de Gauss-Jordan, factorizacion LU y factorizacion LU con pivotamiento parcial y completo. Los autores ofrecen un conjunto amplio de pruebas sobre los algoritmos y estudian la influencia de las caracteristicas de las GPUs en el desempeno de los metodos. En el mismo ano, Ino et al. [55] presentaron una implementacion de la factorizacion LU en GPU utilizando estrategias right-looking. Posteriormente, Jung y O'Leary [57,58] presentaron una implementacion de la factorizacion de Cholesky sobre una GPU. Los autores proponen almacenar las matrices (simetricas) en formato comprimido rectangular. Posterior a la presentacion de CUBLAS destaca el trabajo de Ries et al. [79], donde se estudia la inversion de matrices triangulares sobre GPU. Esta operacion es necesaria para computar la matriz inversa mediante la factorizacion LU. Implementaciones eficientes de la factorizacion LU y la factorizacion de Cholesky tambien fueron presentadas por Barrachina et al. [11,9]. En el estudio se incluye la evaluacion de estrategias de computo mixto GPU+CPU, padding y precision mixta con refinamiento iterativo.

Ideas similares son estudiadas en los trabajos de Baboulin et al. [7,83], donde se propone acelerar la factorizacion LU modificando el pivotamiento por un metodo que implica menor costo computacional (aunque sufre de ciertos problemas de estabilidad). Finalmente, los autores extienden los trabajos para acelerar la factorizacion de Cholesky mediante el uso de planificacion dinamica [68]. La inversion de matrices generales con GPUs es estudiada por Ezzatti et al. en [44], mientras que en [43] se extiende el trabajo para utilizar multiples GPUs. Por otro lado, la inversion de matrices simetricas y definidas positivas es cubierto en el articulo de Benner et al. [23].

La arquitectura masivamente paralela de las GPUs es mas adecuada para la resolucion de operaciones de ALN con matrices densas; sin embargo, en diversos trabajos se alcanzan desempenos importantes sobre matrices dispersas. La mayor parte de las investigaciones en esta area se han centrado en la aceleracion de metodos iterativos, ya que la resolucion de estos metodos generalmente se basa en el producto matrix dispersa por vector (spmv) siendo esta operacion mas sencilla de implementar en GPU que las involucradas en los metodos directos. Ademas de los trabajos pioneros con implementaciones del metodo de Jacobi y de Gauss-Seidel de Bajaj et al. [8], y del metodo GC de Bolz et al. [34], Goodnight et al. [50], Hillesland et al. [54], y Kruger et al [62], conviene destacar los trabajos que se comentan a continuacion. Buatois et al. [35] presentan implementaciones del metodo GC incluyendo el uso de un precondicionador de Jacobi. Para resolver el metodo GC implementan las operaciones spmv, dot y SAXPY de BLAS. La implementacion es para operaciones con matrices dispersas, para las que emplean una estrategia de almacenamiento CRS a bloques BCRS. Ademas de los trabajos mencionados, cuyo objetivo era acelerar diferentes metodos iterativos, varios estudios han analizado directa o indirectamente la aceleracion de la funcion spmv. Este es el caso del trabajo de Sengupta et al. [82] donde se presenta una implementacion de scan segmentado sobre GPU. Las matrices dispersas se almacenan con el formato CSR, pero tambien se estudian otras estrategias de almacenamiento. Mas recientemente Bell y Garland [14] presentaron el estudio de la multiplicacion de matriz dispersa por vector sobre CUDA. Evaluan diferentes formatos de almacenamiento para las matrices dispersas e implementan la operacion para los diferentes formatos evaluados. Otros trabajos destacados son el de Gaikwad y Toke [46] que implementan el metodo BiCGTAB y CGS para diferentes formatos de almacenamiento y distintas implementaciones de spmv, la de NVIDIA [14], IBM [13] y la presentada en [35]; y el de Anzt et al. [6] tambien sobre implementaciones de los metodos del CG y GMRES donde se estudian estrategias de precision mixta, reordenamiento y el impacto en el consumo energetico. Para el computo con matrices banda, en el caso de las matrices tridiagonales los trabajos se han centrado en implementaciones del metodo de reduccion ciclica (RC) o sus variantes, en este campo destacamos los trabajos de Zhang et al [88], Maciol et al. [70], Goddeke et al. [49] y Alfaro et al. [3]. En cuanto a la implementacion de metodos directos para matrices dispersas existen pocos trabajos en la literatura. El pionero dentro de tema es el desarrollo de Christen et al. [38], en el que los autores estudian distintas opciones para extender la biblioteca PARDISO [81] utilizando la GPU. La biblioteca utiliza el metodo multifrontal y evalua la aplicacion a nivel de operaciones basicas (mediante invocaciones a nucleos computacionales de BLAS) o a nivel de resolucion de los frentes (mediante la invocacion de nucleos de LAPACK). Mas recientemente, Lucas et al. [69] extienden el trabajo estudiando implementaciones del metodo multifrontal en GPU.

4.3. Aplicacion de tecnicas de HPC para la reduccion de modelos

Uno de los esfuerzos pioneros en el desarrollo de herramientas que incluyen estrategias de HPC para la resolucion de problemas de control, y en particular de reduccion de modelos, es la biblioteca SLICOT (Subroutine Library in Control Theory) (5) [84], que fue desarrollada a finales de la decada de los 90, en el marco del proyecto NICONET (Numerics in Control Network) [2]. La biblioteca esta implementada en Fortran 77 y tambien es posible su utilizacion desde Matlab. Su diseno se basa en el uso de las especificaciones LAPACK y BLAS, por lo cual se puede obtener facilmente una version paralela de la biblioteca sobre memoria compartida utilizando una implementacion multihilo de BLAS. SLICOT es fuertemente utilizada por investigadores del area en la actualidad. Los esfuerzos previos a la introduccion de SLICOT por desarrollar bibliotecas para la resolucion de problemas de control estan revelados en el trabajo de Benner et al. [27].

Posteriormente, se desarrollaron diferentes esfuerzos para incorporar estrategias de paralelismo de memoria distribuida en la biblioteca SLICOT, que llevaron a la realizacion de la biblioteca PLiC (Parallel Library in Control); los principales trabajos relacionados a este proceso se resumen a continuacion.

El trabajo de Blanquer et al. [33] presento algunas propuestas tendientes a paralelizar la resolucion de la ecuacion de Lyapunov proponiendo PSLICOT (Parallel SLICOT). Posteriormente, uno de los principales esfuerzos se vio plasmado en el desarrollo de la primera version de la biblioteca PLiC [30] propiamente dicho. La biblioteca PLiC posee dos componentes fundamentales, PLiCOC y PLiCMR (del ingles Parallel Library in Control: Model Reduction). La primera da soporte a la resolucion de problemas de control, mientras que el segundo componente esta especialmente disenado para la resolucion de problemas de reduccion de modelos. PLiC incluye rutinas para la resolucion de problemas de analisis y diseno de SDLs invariantes en el tiempo y de gran escala sobre arquitecturas de computadoras distribuidas y paralelas. Contiene mas de 30 subrutinas que permiten resolver ecuaciones matriciales lineales y cuadraticas (Lyapunov, Sylvester, Stein y la ecuacion algebraica de Riccati), reduccion de modelos, estabilizacion parcial, y diversas rutinas de soporte. El modelo de paralelismo emula el utilizado por la biblioteca ScaLAPACK. En particular, la resolucion de las operaciones de algebra se hacen invocando a rutinas de ScaLAPACK, siempre que es posible, y las comunicaciones se efectuan mediante invocacion a rutinas de la biblioteca BLACS.

Otra extension a la biblioteca PLiC fue el desarrollo de una version accesible mediante servicios web. Esta propuesta permitio resolver problemas de reduccion de modelos en sistemas remotos a traves de la red, eliminando la necesidad de contar con plataformas de HPC locales para la resolucion de este tipo de problemas. En el trabajo de Benner et al. [26] se puede profundizar en esta tematica. Con el mismo objetivo que el trabajo anterior, de facilitar el uso de la biblioteca PLiC, en el trabajo de Galiano et al. [47] se desarrollo el paquete PyPLiC, en el cual se ofrecio una interfaz de alto nivel en el lenguaje de scripting Python para la biblioteca PLiCMR.

En el trabajo de Remon et al. [77], se aplicaron diferentes tecnicas de paralelismo al metodo LR-ADI para matrices banda. Los autores presentaron dos implementaciones paralelas sobre memoria compartida (arquitecturas SMP) para la factorizacion LU de matrices de banda. En el trabajo se mostro tambien como aplicar las variantes de factorizacion desarrolladas para acelerar los metodos LR-ADI en la resolucion de problemas de reduccion de modelos, consiguiendo porcentajes de mejora superiores al 50%.

En el caso de los metodos de reduccion de modelos basados en subespacios de Krylov, existe gran caudal de trabajos referentes a la aplicacion de tecnicas de HPC [12,80]. Generalmente, dada la estrategia de resolucion de estos metodos, las tecnicas de HPC se aplicaron a la resolucion del sistema, es decir, al metodo de Lanczos o alguna de sus variantes.

Finalmente, en Ezzatti et al. [42], los autores estudiaron la aceleracion del metodo BST sobre matrices simetricas definidas positivas sobre arquitecturas multinucleo.

4.4. Reduccion de modelos y problemas de control con GPUs

Como se ha comentado anteriormente, los problemas de reduccion de modelos y de control estan fuertemente ligados a la resolucion de operaciones de algebra lineal. Dados los excelentes desempenos alcanzados por las GPUs en este tipo de operaciones, en los ultimos anos se han realizado diversos esfuerzos para su aplicacion en la aceleracion de la resolucion de problemas de reduccion de modelos y de control.

El primer esfuerzo tendiente al uso de GPUs para acelerar problemas de reduccion de modelos fue el de Benner et al. [21] sobre el computo de la funcion de signo de una matriz empleando una GPU. En el mismo se propone una estrategia hibrida en la que tanto la CPU como la GPU colaboran en el computo de la funcion de signo de una matriz. Posteriormente, el trabajo fue extendido [18] para la resolucion de ecuaciones de Lyapunov al utilizar una GPU y explotando estrategias de precision mixta para alcanzar niveles de precision altos a un costo computacional bajo. Tambien es interesante el articulo [17], donde se acelera la resolucion del problema de reduccion de modelos generalizado, mediante el metodo BT sobre arquitecturas hibridas GPU-CPU. El trabajo se basa en el uso de implementaciones eficientes propietarias, y particularmente desarrolladas sobre GPU, de las principales operaciones matriciales requeridas por el metodo: la factorizacion LU, la resolucion de sistemas triangulares y la multiplicacion de matrices. Los metodos de error relativo tambien fueron abordados con metodologias de HPC sobre procesadores graficos en [22]. En ese trabajo, los autores emplearon una NVIDIA Tesla C2050 para resolver la ecuacion de Lyapunov y la ecuacion algebraica de Riccati (EAR), utilizando aritmetica en doble precision para implementar una version eficiente del metodo Balanced Stochastic Truncation (BST).

La resolucion de ecuaciones diferenciales de Riccati en plataformas equipadas con una GPU fue abordado en [19], y extendido al caso en el que la plataforma destino estaba equipada con multiples GPUs en [20]. El uso de diversas GPUs, ademas de acelerar los computos, permitio exteder la dimension de los problemas tratados alcanzando problemas de dimension 34.000, trabajando con matrices completas, mientras se alcanzaban valores de aceleracion de hasta 20x.

5. Ejemplos numericos

Modelo PEEC de un inductor espiral [66]. En la Figura 1 se muestra un inductor espiral con una placa de cobre en la superficie. Un sistema dinamico de la forma (1) es obtenido utilizando el software Fasthenry, ver [59,60]. El sistema que se obtiene es:

[EXPRESION MATEMATICA IRREPRODUCIBLE EN ASCII], (10)

donde

[EXPRESION MATEMATICA IRREPRODUCIBLE EN ASCII],

la matriz de inductancia [??] es simetrica y definida positiva, [[??].sup.1/2] representa la raiz cuadrada, esto es, se satisface [EXPRESION MATEMATICA IRREPRODUCIBLE EN ASCII] es la matriz de resistencia la cual es cuadrada y dispersa. Referimos al lector a [59,66,60] para una version detallada del modelo.

Cabe senalar que [??] no se calcula explicitamente, por su elevado costo computacional; en su lugar las operaciones para la matriz [??] se realizan implicitamente.

[FIGURA 1 OMITIR]

Aplicamos reduccion de modelos via truncamiento balanceado explotando la dispersion mediante el Algoritmo 3.2 para resolver las ecuaciones de Lyapunov. El orden original del sistema (10) es n = 1434. El modelo reducido se calculo de manera que satisficiera una tolerancia de [10.sup.-8], para lo que se obtuvo un modelo de orden r = 18. En la Figura 2 se muestra el valor de las funciones de transferencia. La sobreposicion de estos valores constituye un indicador cuantitativo de la eficiencia del modelo reducido para describir la dinamica del modelo original. En efecto, el modelo reducido aproxima al modelo original con 10 digitos de precision. En la Figura 3 se presentan el error absoluto y relativo.

[FIGURA 2 OMITIR]

[FIGURA 3 OMITIR]

Enfriamiento optimo de rieles de tren [28,41]. Este problema aparece en el proceso de produccion de rieles de tren. La idea es reducir la temperatura del riel lo mas rapido posible maximizando las propiedades del material, como maleabilidad, porosidad, etc. El proceso de enfriamiento se realiza rociando fluidos a bajas temperaturas sobre la superficie del material. Se asume que el riel tiene longitud infinita, lo que permite definir el problema en dos dimensiones. Explotando la simetria del dominio, un mallado inicial se presenta en la Figura 4. Si ademas se linealiza el modelo, entonces las ecuaciones que describen el fenomeno tienen la forma

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (11)

donde x([xi], t) representa la temperatura en el tiempo t en el punto [xi], [g.sub.i] incluye las diferencias de temperatura entre el fluido que se rocia y la superficie del material, parametros de intensidad y coeficientes de transferencia de calor.

[FIGURA 4 OMITIR]

Tras realizar la discretizacion se obtiene un sistema como el siguiente:

[EXPRESION MATEMATICA IRREPRODUCIBLE EN ASCII], (12)

donde M, [elemento de] G [R.sup.nxn], M es invertible y [M.sup.-1]N es estable. Luego la forma estandar del sistema es:

A = [M.sup.-1.sub.L] N [M.sup.-1.sub.U], B = [M.sup.-1.sub.L] B, C = C[M.sup.-1.sub.U].

Cabe senalar que A no se calcula explicitamente, por su elevado costo computational. En su lugar las operaciones para la matriz A se realizan implicitamente.

Nuevamente se aplico reduccion de modelos via truncamiento balanceado explotando la dispersion mediante el Algoritmo 3.2 para resolver las ecuaciones de Lyapunov. El orden del sistema original (12) es n = 5177. El modelo reducido se calculo de manera que satisficiera una tolerancia de [10.sup.-6] y [10.sup.-4], se obtuvieron modelos de orden r = 47 y r = 18, respectivamente. En la Figura 5 se muestra el valor de las funciones de transferencia y tambien el error obtenido.

[FIGURA 5 OMITIR]

6. Conclusiones

Con base a lo expuesto es posible aplicar modelos matematicos de gran escala en problemas de simulacion y control, con la ayuda de tecnicas de computacion de alto desempeno.

En particular, en el presente trabajo, se describen algoritmos eficientes de truncamiento balanceado para la resolucion de la ecuacion de Lyapunov, asi como para la resolucion de ecuaciones diferenciales y algebraicas de Riccati. La resolucion de este tipo de ecuaciones constituye el reto computacional mas importante que aparece al tratar problemas de control lineal cuadratico y problemas no lineales mediante tecnicas de control predictivo. La idea fundamental, en todos los casos, es aplicar tecnicas de algebra lineal numerica que exploten la estructura de cada problema particular; de esta manera el costo de resolucion, desde el punto de vista de la cantidad de operaciones aritmeticas necesarias y de la cantidad de memoria requerida, se reducen considerablemente. En la actualidad, las arquitecturas hibridas (compuestas por procesadores multinucleo conectados a uno o varios procesadores graficos) ofrecen una elevada capacidad de computo y al mismo tiempo presentan un costo economico moderado. El uso de los algoritmos citados en este articulo, y de software y tecnicas de computacion de alto desempeno sobre arquitecturas hibridas CPU-GPU, permite la resolucion de problemas de gran escala que anteriormente requerian del uso de plataformas hardware de elevado costo.

Referencias

[1.] Repositorio Netlib. www.netlib.org/. Consultado en octubre (2011)

[2.] Sitio Web oficial de la biblioteca SLICOT www.slicot.org/

[3.] Alfaro, P., Igounet, P, and Ezzatti, P.: Resolucion de matrices tri-diagonales utilizando una tarjeta grafica (GPU) de escritorio. Mecanica Computacional, 30 (2010) 2951-2967

[4.] Antoulas A.C.: Lectures on the approximation of linear dynamical systems. Encyclopedia of Electrical and Electronics Engineering. John Wiley and Sons (1999) 403-422

[5.] Antoulas, A. C., Sorensen, D. C., and Gugercin, S.: A survey of model reduction methods for large-scale systems. Contemporary Mathematics, 280 (2001) 193-219

[6.] Anzt, H., Rocker, B. and Heuveline, V.: Energy efficiency of mixed precision iterative refinement methods using hybrid hardware platforms - An evaluation of different solver and hardware configurations. Computer Science - R & D, 25 (2010) 141-148.

[7.] Baboulin, M., Dongarra, J. and Tomov, S.: Some Issues in Dense Linear Algebra for Multicore and Special Purpose Architectures. Manchester Institute for Mathematical Sciences, University of Manchester, Manchester, UK, jan (2009)

[8.] Bajaj, C., Ihm, I., and Min, J. and Oh, J.: SIMD Optimization of Linear Expressions for Programmable Graphics Hardware. Computer Graphics Forum, 23 (2004) 697714

[9.] Barrachina, S., Castillo, M., Igual, F. D., Mayo, R., Quintana-Orti, E. S.: Solving Dense Linear Systems on Graphics Processors. in Euro-Par '08: Proceedings of the 14th international Euro-Par conference on Parallel Processing, Berlin, Heidelberg, Springer-Verlag, (2008) 739-748

[10.] Barrachina, S., Castillo, M., Igual, F. D., Mayo, R., Quintana-Orti, E. S., QuintanaOrti, G.: Evaluation and Tuning of the Level 3 CUBLAS for Graphics Processors. Departamento de Ingenieria y Ciencia de Computadores, Universidad Jaime I, Campus de Riu Sec, s/n 12.071 - Castellon, Espana, (2008)

[11.] Barrachina, S., Castillo, M., Igual, F. D., Mayo R., Quintana-Orti, E. S., QuintanaOrti, G.: Exploiting the capabilities of modern GPUs for dense matrix computations, Concurrency and Computation: Practice and Experience, 21 (2009) 2457-2477

[12.] Barrett, R., Berry, M., Chan, T. F., Demmel, J., Donato, J., Dongarra, J., Eijkhout, V., Pozo, R., Romine, C., Van der Vorst, H.: Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods. 2nd Edition, SIAM, Philadelphia, PA, (1994)

[13.] Baskaran, M., Bordawekar, R.: Optimizing sparse matrix-vector multiplication on GPUs, IBM Research Report 24704 (2009).

[14.] Bell, N., Garland, M. Implementing sparse matrix-vector multiplication on throughput-oriented processors. Proceedings of the Conference on High Performance Computing Networking, Storage and Analysis, SC '09, New York, NY, USA, ACM, (2009) 18:1-18:11

[15.] Benner, P.: Solving large-scale control problems. IEEE Control Systems Magazine, 14(1) (2004) 44-59

[16.] Benner, P.: System-theoretic methods for model reduction of large-scale systems: Simulation, control, and inverse problems. Proceedings of MathMod 2009, Vienna, February 11-13, 2009, I. Troch and F. Breitenecker, eds., vol. 35 of ARGESIM Reports, (2009) 126-145

[17.] Benner, P., Ezzatti, P., Kressner, D., Quintana-Orti, E. S., Remon, A.: Accelerating model reduction of large linear systems with graphics processors. In Lecture Notes in Computer Science, State of the Art in Scientific and Parallel Computing, Springer, (2010)

[18.] Benner, P., Ezzatti, P., Kressner, D., Quintana-Orti, E. S., Remon, A.: A mixedprecision algorithm for the solution of Lyapunov equations on hybrid CPU- GPU platforms. Parallel Computing, 37 (2011) 439-450

[19.] Benner, P., Ezzatti, P., Mena, H., Quintana-Orti, E. S., Remon, A.: Solving differential Riccati equations on multi-GPU platforms. In 2nd Meeting on Linear Algebra, Matrix Analysis and Applications ALAMA10, (2010)

[20.] Benner, P., Ezzatti, P., Mena, H., Quintana-Orti, E. S., Remon, A.: Solving differential Riccati equations on multi-GPU platforms. In 10th International Conference on Computational and Mathematical Methods in Science and Engineering CMMSE11, (2011) 178-188

[21.] Benner, P., Ezzatti, P., Kressner, D., Quintana-Orti, E. S., Remon, A.: Using hybrid CPU-GPU platforms to accelerate the computation of the matrix sign function. In Euro-Par Workshops, H.-X. Lin, M. Alexander, M. Forsell, A. Knupfer, R. Prodan, L. Sousa, and A. Streit, eds., vol. 6043 of Lecture Notes in Computer Science, Springer, (2009) 132-139

[22.] Benner, P., Ezzatti, P., Kressner, D., Quintana-Orti, E. S., Remon, A.: Accelerating BST methods for model reduction with graphics processors. In Proceedings of the 9th International Conference on Parallel Processing and Applied Mathematics, (2011)

[23.] Benner, P., Ezzatti, P., Kressner, D., Quintana-Orti, E. S., Remon, A.: Hing performance matrix inversion of SPD matrices on graphics processors. In Workshop on Exploitation of Hardware Accelerators WEHA 2011, (2011) 640-646

[24.] Benner, P., Hinze, M., Ter Maten, J.: Model Reduction for Circuit Simulation. Vol. 74 of Lecture Notes in Electrical Engineering, Springer-Verlag, Berlin/Heidelberg, Germany, (2011)

[25.] Benner, P., Li, J.-R., Penzl, T.: Numerical solution of large Lyapunov equations, Riccati equations, and linear-quadratic control problems. Numer. Linear Algebra Appl., 15 (2008) 755-777

[26.] Benner, P., Mayo, R., Quintana-Orti E. S., Quintana-Orti, G.: Enhanced services for remote model reduction of large-scale dense linear systems. In PARA, J. Fagerholm, J. Haataja, J. Jarvinen, M. Lyly, P. Raback, and V. Savolainen, eds., vol. 2367 of Lecture Notes in Computer Science, Springer, (2002) 329-338

[27.] Benner, P., Mehrmann, V., Sima, V., Huffel, S. V., Varga, A.: SLICOT -a subroutine library in systems and control theory. Applied and Computational Control, Signals, and Circuits, Birkhuser, (1997) 499-539

[28.] Benner, P., Mehrmann, V., Sorensen, D.: Dimension Reduction of Large-Scale Systems. Vol. 45 of Lecture Notes in Computational Science and Engineering. Springer-Verlag, Berlin/Heidelberg, Germany, (2005)

[29.] Benner, P., Mena, H.: BDF methods for large-scale differential Riccati equations. In Proc. of Mathematical Theory of Network and Systems, MTNS 2004, B. D. Moor, B. Motmans, J. Willems, P. V. Dooren, and V. Blondel, eds., (2004)

[30.] Benner, P., Quintana-Orti E. S., Quintana-Orti, G.: A portable subroutine library for solving linear control problems on distributed memory computers. In Workshop on Wide Area Networks and High Performance Computing, London, UK, SpringerVerlag, (1999) 61-87

[31.] Bischof, C.H., Quintana-Orti, G.: Computing rank-revealing QR factorizations of dense matrices. ACM Transactions on Mathematical Software, 24(2) (1998) 226-253.

[32.] Blackford, L. S., Choi, J., Cleary, A., Petitet, A., Whaley, R. C., Demmel, J., Dhillon, I., Stanley, K., Dongarra, J., Hammarling, S., Henry, G., Walker, D.: ScaLAPACK: a portable linear algebra library for distributed memory computers - design issues and performance. In Proceedings of the 1996 ACM/IEEE conference on Supercomputing (CDROM), Supercomputing -96, Washington, DCUSA, IEEE Computer Society (1996)

[33.] Blanquer, I., Guerrero, D., Hernandez, V., Quintana-Orti, E. S., Ruiz, P. A.: ParallelSLICOT implementation and documentation standards. Tech. rep., SLICOT Working Note (1998)

[34.] Bolz, J., Farmer, I., Grinspun, E., Schrooder, P.: Sparse matrix solvers on the GPU: conjugate gradients and multigrid. ACM Trans. Graph., 22 (2003) 917-924

[35.] Buatois, L., Caumon, G., Levy, B.: Concurrent number cruncher: An efficient sparse linear solver on the GPU. In High Performance Computation Conference (HPCC), Springer Lecture Notes in Computer Sciences, (2007). Award: Second best student paper.

[36.] Chien, L. S.: Hand Tuned SGEMM on GT200 GPU. Tech. rep., Department of Mathematics, Tsing Hua University, Taiwan, Feb. (2010)

[37.] Choi, J., Dongarra, J., Walker, D.: PB-BLAS: A set of parallel block basic linear algebra subprograms. In Proc. of the 1994 Scalable High Performance Computing Conference, IEEE Computer Society Press, (1994)

[38.] Christen, M., Schenk, O., Burkhart, H.: General-purpose sparse matrix building blocks using the NVIDIA CUDA technology platform. Tech. rep., (2007)

[39.] Cong, J., Shinnerl, J. R., Xie, M., Kong, T., Yuan, X.: Large-scale circuit placement. ACM Trans. Des. Autom. Electron. Syst., 10 (2005) 389-430.

[40.] Demmel, J., Dongarra, J., Croz, J. D., Greenbaum, A., Hammarling, S., Sorensen, D.: Prospectus for the development of a linear algebra library for high-performance computers. Tech. Rep. ANL/MCS-TM-97, 9700 South Cass Avenue, Argonne, IL 60439-4801, USA, (1987)

[41.] Eppler, K., Troltzsch, F.: Discrete and continuous optimal control strategies in the selective cooling of steel profiles., Z. Angew. Math. Mech., 81 (2001) 247-248

[42.] Ezzatti, P., Quintana-Orti, E. S., Remon, A.: Effcient model order reduction of large-scale systems on multi-core platforms. In ICCSA (5), B. Murgante, O. Gervasi, A. Iglesias, D. Taniar, and B. O. Apduhan, eds., vol. 6786 of Lecture Notes in Computer Science, Springer, (2011) 643-653

[43.] Ezzatti, P., Quintana-Orti, E. S., Remon, A.: High performance matrix inversion on a multi-core platform with several GPUs. IEEE Computer Society, (2011) 87-93

[44.] Ezzatti, P., Quintana-Orti, E. S., Remon, A.: Using graphics processors to accelerate the computation of the matrix inverse. The Journal of Supercomputing, online (2011).

[45.] Fatica, M.: Accelerating LINPACK with CUDA on heterogenous clusters. In GPGPU, (2009) 46-51

[46.] Gaikwad, A., Toke, I. M.: Gpu based sparse grid technique for solving multidimensional options pricing pdes. In Proceedings of the 2nd Workshop on High Performance Computational Finance, WHPCF -09, New York, NY, USA, ACM, (2009) 6:1-6:9

[47.] Galiano V., Martin A., Migallon, H. Migallon, V. Penades, J., Quintana-Orti, E.S.: PyPLiC: A high-level interface to the parallel model reduction library PLiCMR. In Proceedings of the Eleventh International Conference on Civil, Structural and Environmental Engineering Computing, B. H. V. Topping, ed., Stirlingshire, United Kingdom, (2007), Civil-Comp Press. paper 62.

[48.] Galoppo, N., Govindaraju, N. K., Henson, M., Manocha, D.: LU-GPU: Effcient algorithms for solving dense linear systems on graphics hardware. In SC 05: Proceedings of the 2005 ACM/IEEE conference on Supercomputing, Washing- ton, DC, USA, IEEE Computer Society, (2005) 3

[49.] Goddeke, D., Strzodka, R.A.: Cyclic reduction tridiagonal solvers on GPUs applied to mixed precision multigrid. IEEE Transactions on Parallel and Distributed Systems, doi: 10.1109/TPDS.2010.61, 22 (2011) 22-32

[50.] Goodnight, N., Woolley, C., Lewin, G., Luebke, D., Humphreys, G.: A multigrid solver for boundary value problems using programmable graphics hardware. In HWWS '03: Proceedings of the ACM SIGGRAPH/EUROGRAPHICS conference on Graphics hardware, Aire-la-Ville, Switzerland, Switzerland, Eurographics Association, (2003) 102-111

[51.] Gugercin, S., Sorensen, D., Antoulas, A.: A modified low-rank Smith method for large-scale Lyapunov equations. Numer. Algorithms, 32(1) (2003) 27-55

[52.] Hall, J., Carr, N., Hart, J.: Cache and bandwidth aware matrix multiplication on the GPU. Tech. rep., UIUCDCS-R-20032328, University of Illinois, (2003)

[53.] Higham, N.: Functions of Matrices: Theory and Computation. SIAM, Philadelphia, USA, (2008)

[54.] Hillesland, K. E., Molinov, S. Grzeszczuk, R.: Nonlinear optimization framework for image-based modeling on programmable graphics hardware. In ACM SIGGRAPH 2005 Courses, SIGGRAPH '05, New York, NY, USA, ACM, (2005)

[55.] Ino, F., Matsui, M., Goda, K., Hagihara, K.: Performance study of LU decomposition on the programmable GPU. In HiPC, (2005) 83-94

[56.] Iordache, M., Dumitriu, L.: Efficient decomposition techniques for symbolic analysis of large-scale analog circuits by state variable method. Analog Integr. Circuits Signal Process., 40 (2004) 235-253

[57.] Jung, J. H., O'leary. D.: Exploiting structure of symmetric or triangular matrices on a GPU. In First Workshop on General Purpose Processing on Graphics Processing Units, Northeastern Univ., Boston, (2007)

[58.] Jung, J. H., O'leary. D.: Implementing an interior point method for linear programs on a CPU-GPU system. Electronic Transactions on Numerical Analysis

[59.] Kamon, M., Tsuk, M., White, J.: Fasthenry: A multipole-accelerated 3-d inductance extraction program. IEEE Transactions on Microwave Theory and Techniques, 42 (1994) 1750-1758

[60.] Kamon, M., Wang, F., White, J.: Generating nearly optimal compact models from krylov-subspace based reduced order models. IEEE Transactions On Circuits and Systems-II: Analog and Digital Signal Processing, 47 (2000) 239-248

[61.] Kolmogorov, A., Fomin, S. V.: Elements of the Theory of Functions and Functional Analysis. Dover Publications, (1999)

[62.] Kruger, J., Schiwietz, T., Kipfer, P., Westermann, R.: Numerical simulations on PC graphics hardware. In ParSim 2004 (Special Session of EuroPVM/MPI 2004, Budapest, Hungary, (2004) 442-450

[63.] Kruger, J., Westermann, R.: Linear algebra operators for GPU implementation of numerical algorithms. ACM Transactions on Graphics, 22 (2003) 908-916

[64.] Larsen, E. S., McAllister, D.: Fast matrix multiplies using graphics hardware. In Proceedings of the 2001 ACM/IEEE conference on Supercomputing (CDROM), Supercomputing '01, New York, NY, USA, (2001), ACM, 55-55

[65.] Lasiecka, I., Triggiani, R.: Control Theory for Partial Differential Equations: Continuous and Approximation Theories I: Abstract Parabolic Systems. Cambridge University Press, Cambridge, UK, (2000)

[66.] Li, J.-R., Kamon, M.: PEEC model of a spiral inductor generated by Fasthenry, in Dimension Reduction of Large-Scale Systems. P. Benner, V. Mehrmann, and D. Sorensen, eds., vol. 45 of Lecture Notes in Computational Science and Engineering, Springer-Verlag, Berlin/Heidelberg, Germany, (2005) 373-377

[67.] Li, J.-R., White, J.: Reduction of large circuit models via low rank approximate gramians. International Journal of Applied Mathematics and Computer Science, 11 (2001) 101-121

[68.] Ltaief, H., Tomov, S., Nath, R., Du, P., Dongarra, J.: A scalable high performant Cholesky factorization for multicore with GPU accelerators. In VECPAR, vol. 6449 of Lecture Notes in Computer Science, Springer, (2010) 93-101

[69.] Lucas, R. F., Wagenbreth, G., Davis, D. M., Grimes, R.: Multifrontal computations on GPUs and their multi-core hosts. In Proceedings of the 9th international conference on High performance computing for computational science, VECPAR'10, Berlin, Heidelberg, Springer-Verlag, (2011) 71-82

[70.] Maciol, P., Banas K.: Testing tesla architecture for scientific computing: the performance of matrix-vector product. vol. 3, (2008)

[71.] Mena, H.: Numerical Solution of Differential Riccati Equations Arising in Optimal Control Problems for Parabolic Partial Differential Equations. PhD thesis, Escuela Politecnica Nacional, Quito, Ecuador, (2007)

[72.] Moravanszky., A., Ag., N.: Dense matrix algebra on the GPU. In Direct3D ShaderX2, Engel W. F., (Ed.). Wordware Publishing, NovodeX AG, (2003) 2

[73.] Nath, R., Tomov, S., Dongarra, J.: BLAS for GPUs. In Scientific Computing with Multicore and Accelerators, J. Kurzak, D. A. Bader, and J. a. Dongarra, eds., CRC Press, Dec. (2010)

[74.] Penzl, T.: Lyapack Users Guide. Tech. Rep. SFB393/00-33, Sonderforschungsbereich 393 Numerische Simulation auf massiv parallelen Rechnern, TU Chemnitz, 09107 Chemnitz, Germany, (2000). Available from http://www.tu-chemnitz.de/sfb393/ sfb00pr.html.

[75.] Penzl, T.: Algorithms for model reduction of large dynamical systems. Linear Algebra Applications, 415 (2006) 322-343. (Reprint of Technical Report SFB393/99 40, TU Chemnitz, (1999)

[76.] Nath, S. T. R., Dongarra, J.: An Improved MAGMA GEMM for Fermi Graphics Processing Units. International Journal in High Performance Computing and Architectures, 24 (2010) 511-515

[77.] Remon, A., Quintana-Orti, E., Quintana-Orti, G.: Parallel solution of band linear systems in model reduction. In Parallel Processing and Applied Mathematics, R. Wyrzykowski, J. Dongarra, K. Karczewski, and J. Wasniewski, eds., vol. 4967 of Lecture Notes in Computer Science, Springer Berlin / Heidelberg, (2008) 678-687

[78.] Riaza R.: Differential-Algebraic Systems. Analytical Aspects and Circuit Applications, World Scientific, (2008)

[79.] Ries, F., De Marco, T., Zivieri, M., Guerrieri, R.: Triangular matrix inversion on graphics processing unit. In Proceedings of the Conference on High Performance Computing Networking, Storage and Analysis, SC '09, New York, NY, USA, ACM, (2009) 9:1-9:10

[80.] Saad, Y.: Iterative Methods for Sparse Linear Systems. Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 2nd ed., (2003)

[81.] Schenk, O. Gartner, K.: Sparse factorization with two level scheduling in pardiso. In PPSC, (2001)

[82.] Sengupta, S., Harris, M., Zhang, Y., Owens, J. D. Scan primitives for GPU computing. In GH '07: Proceedings of the 22nd ACM SIGGRAPH/EUROGRAPHICS symposium on Graphics hardware, Aire-la-Ville, Switzerland, Switzerland, (2007), Eurographics Association, 97-106

[83.] Tomov, S., Nath, R., Ltaief, H., Dongarra, J.: Dense linear algebra solvers for multicore with GPU accelerators. In IPDPS Workshops, IEEE, (2010) 1-8

[84.] Varga, A.: Task II.B.1 - selection of software for controller reduction. SLICOT Working Note 1999-18, The Working Group on Software (WGS), http://www. slicot.org/index.php?site=SLmodredR, (1999)

[85.] Varga, A.: Model reduction software in the SLICOT library. In Applied and Computational Control, Signals, and Circuits, volume 629 of The Kluwer International Series in Engineering and Computer Science, Kluwer Academic Publishers, (2000) 239-282

[86.] Volkov, V, Demmel, J.: Benchmarking GPUs to tune dense linear algebra. In SC '08: Proceedings of the 2008 ACM/IEEE conference on Supercomputing, Piscataway, NJ, USA, (2008), IEEE Press, (2008) 1-11

[87.] Wachspress, E.L.: Iterative solution of the Lyapunov matrix equation. Appl. Math. Letters, 107 (1988) 87-90

[88.] Zhang, Y., Cohen, J., Owens, J. D.: Fast tridiagonal solvers on the GPU. In PPOPP, (2010) 127-136

FECHA DE ENTREGA: 23 DE ENERO DE 2012

FECHA DE APROBACION: 5 DE NOVIEMBRE DE 2012

Peter Benner (1) *, Pablo Ezzatti (2) **, Hermann Mena (3) ***, Enrique S. Quintana-Ortfo, Alfredo Remon (4) *

(1) Max Planck Institute for Dynamics of Complex Technical Systems, Alemania

(2) Centro de Calculo-Instituto de Computacion, Universidad de la Republica, Uruguay

(3) University of Innsbruck, Austria

(4) Departamento de Ingenieria y Ciencia de Computadores, Universidad Jaime I, Espana

(5) Disponible en http://www.win.tue.nl/niconet/NIC2/slicot.html.

* Director del Max Planck Institute for Dynamics of Complex Technical Systems, Sandtorstr. 1, 39106 Magdeburg, Alemania. benner@mpi-magdeburg.mpg.de

** Profesor del Centro de Calculo-Instituto de Computacion, Universidad de la Republica, 11.300-Montevideo, Uruguay. pezzatti@fing.edu.uy

*** University of Innsbruck, Austria. Hermann.Mena@uibk.ac.at

([cruz]) Catedratico del Departamento de Ingenieria y Ciencia de Computadores, Universidad Jaime I, 12.071-Castellon, Espana. quintana@icc.uji.es

([cruz doble]) Investigador del Departamento de Ingenieria y Ciencia de Computadores, Universidad Jaime I, 12.071-Castellon, Espana. remon@icc.uji.es
COPYRIGHT 2013 Politecnico Grancolombiano
No portion of this article can be reproduced without the express written permission from the copyright holder.
Copyright 2013 Gale, Cengage Learning. All rights reserved.

Article Details
Printer friendly Cite/link Email Feedback
Author:Benner, Peter; Ezzatti, Pablo; Mena, Hermann; Quintana-Orti, Enrique S.; Remon, Alfredo
Publication:Elementos
Date:Jun 1, 2013
Words:11687
Previous Article:Editorial.
Next Article:Boosting en el modelo de aprendizaje PAC.
Topics:

Terms of use | Privacy policy | Copyright © 2019 Farlex, Inc. | Feedback | For webmasters