Source: http://docplayer.es/52921025-Universidad-eafit-universidad-eafit-issn-version-impresa-x-colombia.html
Timestamp: 2017-12-15 07:22:37+00:00

Document:
Universidad Eafit Universidad Eafit ISSN (Versión impresa): X COLOMBIA - PDF
Universidad Eafit Universidad Eafit ISSN (Versión impresa): X COLOMBIA
Download "Universidad Eafit Universidad Eafit ISSN (Versión impresa): X COLOMBIA"
Juan Luis Santos de la Fuente
1 Universidad Eafit Universidad Eafit ISSN (Versión impresa): X COLOMBIA 2001 Carlos Enrique Mejía / Tomás Restrepo / Christian Trefftz LAPACK UNA COLECCIÓN DE RUTINAS PARA RESOLVER PROBLEMAS DE ÁLGEBRA LINEAL NUMÉRICA Universidad Eafit, julio-septiembre, número 123 Universidad Eafit Medellìn, Colombia pp Red de Revistas Científicas de América Latina y el Caribe, España y Portugal Universidad Autónoma del Estado de México
2 73 LAPACK Una colección de rutinas para resolver problemas de Álgebra Lineal Numérica Carlos Enrique Mejía Tomás Restrepo Christian Trefftz 1. INTRODUCCIÓN El modelamiento matemático ha sido siempre fundamental para las ciencias y las ingenierías pero ha cobrado mucha más importancia desde que se empezaron a usar computadores en gran escala. LAPACK (Linear Algebra Package) es una colección de subrutinas escritas en FORTRAN 77 para resolver los problemas matemáticos más comunes que surgen a partir del modelamiento y que se enmarcan en el campo del álgebra lineal numérica. El desarrollo de LAPACK es una tarea que se empezó hacia 1987 y aun no termina. Lo hacen universidades y laboratorios de investigación de Estados Unidos y Europa, un equipo de investigadores de primera línea apoyados financieramente por varias instituciones entre las que se destaca la Fundación Nacional de Ciencias de Estados Unidos (National Science Foundation, NSF.) Las subrutinas de LAPACK se basan en llamadas a unas subrutinas más sencillas que se conocen por la sigla BLAS (Basic Linear Algebra Subprograms.) Hasta el momento se dispone de subprogramas BLAS de tres niveles: Nivel 1, publicado en 1979, que se encarga de operaciones de vector con vector. Ejemplo: Subrutina SAXPY, que sirve para sumar un múltiplo de un vector con otro vector. Nivel 2, publicado en 1988, que se ocupa de operaciones de matriz con vector. Ejemplo: Subrutina SGEMV, que sirve para sumar el producto de un múltiplo de una matriz por un vector, con un múltiplo de otro vector. Nivel 3, publicado en 1990, que se dedica a operaciones entre matrices. Ejemplo: Subrutina SGEMM, que sirve para sumar un múltiplo de una matriz con un múltiplo del producto de otras dos matrices. La construcción de LAPACK con base en los subprogramas BLAS genera un alto nivel de estandarización, proporciona gran rapidez de cálculo y hace más fácil hacer un seguimiento cuando se presentan dificultades. La utilización de LAPACK es universal, no sólo directamente, sino también como motor de cálculo en software con interfase gráfica como OCTAVE, desarrollado en la Universidad de Wisconsin. LAPACK incluye rutinas para resolver sistemas de ecuaciones lineales, sistemas de ecuaciones lineales por mínimos cuadrados, problemas de valores propios y problemas de valores singulares. Se trata pues, de un software de primera calidad, que no requiere de pago alguno CARLOS ENRIQUE MEJÍA SALAZAR. Matemático, Universidad Nacional, Medellín, M.S. Matemáticas, Universidad Nacional, Medellín, Ph.D. Matemáticas, University of Cincinnati, Profesor asociado, Universidad Nacional en Medellín desde TOMÁS RESTREPO. Estudiante de Ingeniería de Sistemas, Universidad EAFIT. CHRISTIAN TREFFTZ G. Ph.D. en Ciencias de la Computación de Michigan State University. Profesor Asistente Grand Valley State University en Allendale Michigan. REVISTA Universidad EAFIT No. 123 julio.agosto.septiembre 2001
3 74 para ser utilizado y que sirve para el tratamiento numérico de los problemas que se encuentran por igual investigadores, profesores y estudiantes, en una gran diversidad de disciplinas científicas y técnicas. Las preguntas que surgen son: Por qué motivo LAPACK es tan desconocido en Colombia? Por qué razón es frecuente que en Colombia se recurra a software de dudosa calidad o a software de calidad comparable a la de LAPACK pero de alto precio? Para evitar que sea necesario seguir haciendo estas preguntas, se propuso divulgar, para una audiencia bastante general, la existencia de LAPACK, la forma de obtenerlo y de instalarlo y una revisión somera de su contenido y su alcance. Este artículo es un resumen de un documento que pretende lograr la divulgación descrita. Se dirige a personas acostumbradas a usar un computador, que conocen rudimentos de FORTRAN y que están familiarizadas con las nociones básicas del álgebra lineal. La construcción de LAPACK con base en los subprogramas BLAS genera un alto nivel de estandarización, proporciona gran rapidez de cálculo y hace más fácil hacer un seguimiento cuando se presentan dificultades. La utilización de LAPACK es universal, no sólo directamente, sino también como motor de cálculo en software con interfase gráfica como OCTAVE, desarrollado en la Universidad de Wisconsin. El contenido de este artículo está organizado de la siguiente forma: La sección 2 se dedica a explicar las formas de obtener e instalar LAPACK. En la sección 3 sobre álgebra lineal numérica, se explican procesos de modelamiento que conducen a problemas de álgebra lineal susceptibles de ser resueltos con LAPACK. La sección 4 se dedica a considerar dos clases de problemas que se pueden resolver con LAPACK, sistemas de ecuaciones lineales y problemas de valores y vectores propios. En la sección 5, de naturaleza más técnica, se consideran aspectos relacionados con almacenamiento de matrices y uso de la memoria. El artículo termina con algunas conclusiones presentadas en la sección INSTALACIÓN DE LAPACK 2.1 Descarga Lo primero que se necesita para instalar LAPACK, es obtener una distribución. El principal sitio del que se puede obtener LAPACK es el sitio NetLib, mantenido por la Universidad de Tennessee y ORNL (Oak Ridge National Laboratory). La página oficial es Al momento de escribir este texto, la última versión de la colección es la 3.0. Hay tres archivos principales con el código fuente: 1. lapack.tgz: Contiene la distribución completa de LAPACK 3.0 en un archivo comprimido usando tar y gzip. 2. lapack-pc.zip: archivo en formato zip con la versión 3.0 de LAPACK para Microsoft Windows. Requiere de la rutina NMAKE de Microsoft y el compilador Fortran de Digital Equipment Corporation, distribuido ahora por la propia Microsoft. 3. lapack-pc-wfc.zip: archivo en formato zip con la versión 3.0 de LAPACK para Microsoft Windows. Requiere de la rutina NMAKE de Microsoft y el compilador Fortran 77/32 de Watcom versión Adicionalmente es posible descargar las colecciones precompiladas, así como algunas de las rutinas individuales. Para estas notas se supone que se desea instalar LAPACK para LINUX, por lo que el archivo correcto que se debe descargar es lapack.tgz. El sistema operativo LINUX, que es de dominio público, puede instalarse fácilmente en microcomputadores con procesadores Intel o compatible, que son la mayoría. Lo más común, es que toda instalación de LINUX incluya compiladores de dominio público de varios lenguajes de programación, incluyendo Fortran. Así que, en principio, estamos ante una plataforma de cálculo numérico de primera calidad basada en software no comercial. Una de las distribuciones mejor mantenidas de LINUX es REDHAT LINUX, con dirección Internet
4 Compilación Lo primero que se necesita hacer es descomprimir el archivo de la distribución. Esto se realiza usando los comandos: gunzip lapack.tgz tar -xvf lapack.tar Para compilar LAPACK, se cambia al directorio./lapack/install. Desde aquí se puede configurar la instalación de LAPACK para obtener las opciones correctas para la compilación. Lo primero es usar el archivo de configuración apropiado para LINUX 3. ALGEBRA LINEAL NUMÉRICA La teoría matemática necesaria para emprender con éxito la solución numérica de problemas de álgebra lineal tiene un poco de álgebra, análisis, topología y posiblemente de otras ramas de la matemática. Toda esa teoría se denomina actualmente Álgebra Lineal Numérica (ALN). Se trata de un campo muy activo de la matemática y uno de los que más aplicaciones tiene. Su desarrollo se aceleró enormemente a partir de 1950, por la aparición y desarrollo de los computadores. Hay dos grandes familias de problemas en ALN: a. Sistemas de ecuaciones algebraicas de la forma Ax=b, con A matriz y x, b vectores. Aquí se incluyen las soluciones por mínimos cuadrados. b. Problemas de valores y vectores propios de la forma Ax=lx, con A matriz, x vector y l escalar (real o complejo). cp make.inc.linux../make.inc Antes de compilar LAPACK propiamente dicho, hay que compilar las BLAS: cd../blas/src make Si éstas funcionan sin ningún problema, se puede compilar la colección y los programas de prueba: cd../../ make Lo primero que se realiza es compilar unas cortas pruebas de desempeño que se encuentran en el directorio INSTALL, y las ejecuta, mostrando su salida en pantalla. Luego, procede a compilar la colección en sí misma, cuyo código fuente se encuentra en el directorio SRC. Luego se procederá a compilar los programas de prueba. Hay que tener en cuenta que la compilación puede tardar bastante tiempo, así que prepárese a esperar un buen rato. Si no está interesado en los programas de prueba, y sólo desea compilar la colección, ejecute el comando make directamente en el directorio SRC. Últimamente han aparecido buenos libros sobre el tema como (Demmel, 1997) y (Trefethen, 1997) que se unen a libros ya clásicos como (Golub, 1989) y (Stewart, 1973). Un libro de nivel intermedio y muy completo que se recomienda es (Noble, 1989). Para dar una idea concreta acerca de este tema de estudio, se considera un ejemplo que sirve de motivación, no solo para el estudio somero de los conceptos y métodos que se ofrecen enseguida sino también para profundizar en algunos de ellos. Este ejemplo sirve para precisar la idea de discretización. Se considera la ecuación diferencial ordinaria de segundo orden y''( t) = f( t, y, y'), a t b, ya ( ) = Ayb, ( ) = B. La función f puede depender de forma lineal o no lineal de sus variables segunda y tercera. Son muchos los problemas de la física y la ingeniería que admiten un modelo matemático de este estilo para su solución. Basta pensar, por ejemplo, en la segunda ley de Newton de la que surge un gran número. Una posible discretización de diferencias finitas para este problema es como sigue: Se genera una subdivisión en n+1 subintervalos de igual longitud del intervalo [a,b]. Los puntos extremos de los subintervalos se denotan a= t0 < t1 <... < tn+ 1 = b. Enseguida para j=1,2,...,n, aproximamos la segunda derivada y''( t j ) por el cociente 1 éy( t 2 j- 1) - 2y( tj) + y( tj+ 1) ù h ë û, y entonces resulta natural querer resolver el sistema 1 æ vj+ 1-vj-1ö v 2 j- 1 2 vj vj+ 1 f tj, vj,, j 1,2,..., n h ë é - + û ù= ç = è 2h ø REVISTA Universidad EAFIT No. 123 julio.agosto.septiembre 2001
5 76 con v(0)=a,v(n+1)=b y v(j) para j=1,...,n aproximaciones de los valores y(t(j)) utilizados para los cálculos. Se tiene entonces un sistema, en general no lineal, de n ecuaciones con n incógnitas v(1),v(2),...,v(n) que generalmente se resuelve por un método iterativo de tipo Newton. Se supone que este sistema tiene solución para preservar sencilla esta exposición. Se advierte sin embargo que estudiar condiciones suficientes o necesarias para que haya solución a sistemas no lineales es un tema de investigación actual en análisis numérico que está lleno de sutilezas y dificultades. El método de Newton para este sistema se enuncia más fácilmente con notación vectorial. Denotamos por V al vector columna cuyas componentes son v(1),v(2),...,v(n) y definimos una función de variable y valor vectorial G de la siguiente manera: n n G: R R con la componente j-ésima de GV igual al lado izquierdo de la ecuación j-ésima correspondiente. De esta forma, el sistema se puede denotar simplemente como GV=0. Su solución por el método iterativo de Newton exige en cada iteración, la solución de un sistema lineal de ecuaciones. El algoritmo numérico debe incluir un número máximo de iteraciones como criterio de parada y algún otro criterio de parada basado en lo cercanos que son iterados consecutivos. Esta noción de cercanía es una noción topológica. La teoría matemática necesaria para emprender con éxito la solución numérica de problemas de álgebra lineal tiene un poco de álgebra, análisis, topología y posiblemente de otras ramas de la matemática. Toda esa teoría se denomina actualmente Álgebra Lineal Numérica (ALN). Se trata de un campo muy activo de la matemática y uno de los que más aplicaciones tiene. Su desarrollo se aceleró enormemente a partir de 1950,por la aparición y desarrollo de los computadores. En resumen: este ejemplo presenta un problema con valores en la frontera en dos puntos, posiblemente no lineal, que generalmente debe tratarse por un método iterativo para su solución. Cuando se enuncia un método de Newton, que es uno de los más comunes y poderosos, se ve la necesidad de resolver, en cada iteración un sistema tridiagonal de ecuaciones lineales. Para resolver sistemas de ese tipo, LAPACK proporciona la rutina SGTSVX, que ofrece la máxima calidad y especialización. 4. ALGUNAS RUTINAS DE LAPACK En esta sección se presentan dos rutinas, de las muchas disponibles en LAPACK, con las que se ilustra la forma en que se pueden utilizar estos procedimientos. 4.1 Ecuaciones Lineales En primer lugar, se considera el sistema de ecuaciones lineales x=b donde A es la matriz de coeficientes, b es el lado derecho y x es la solución. El objetivo es hallar x dados A y b. Si hay varios lados derechos, se escribe AX=B, donde las columnas de B son los lados derechos individuales y las columnas de X son las soluciones correspondientes. La tarea es calcular X, dadas A y B. LAPACK cuenta con once variedades de rutinas para la solución de sistemas de ecuaciones lineales, uno para cada tipo y almacenamiento de matriz. Además, para cada tipo se cuenta con rutinas simples y expertas para cuatro clases de datos diferentes. A continuación se describe una subrutina, de las llamadas driver subroutines, que se encargan de resolver un problema completo de una matriz general. SGESV es un manejador simple para matrices generales. Soluciona un sistema real de ecuaciones lineales de la forma AX=B, donde A es una matriz N x N, y X y B son matrices de N filas y NRHS columnas. La subrutina SGESV utiliza un algoritmo de factorización LU con pivoteo parcial e intercambio de filas para resolver el sistema de ecuaciones, haciendo uso de las rutinas computacionales SGETRF y SGETRS. La primera de éstas realiza la factorización LU de la matriz A, mientras que la segunda utiliza los factores L y U para resolver el sistema de ecuaciones propiamente dicho. Con este algoritmo, una matriz de coeficientes A se descompone en 3 matrices P, L y U, tales
6 77 que A=PLU, donde P es la matriz de permutaciones debidas al pivoteo, L es la matriz triangular inferior y U es la matriz triangular superior. Vale la pena recordar que el realizar pivoteo parcial no afecta el orden de las respuestas que se obtienen del procedimiento, a diferencia de lo que ocurre cuando se utiliza pivoteo total. El prototipo de la subrutina SGESV es: SUBROUTINE SGESV ( N, NRHS, A, LDA, IPIV, B, LDB, INFO ) INTEGER INFO, LDA, LDB, N, NRHS INTEGER IPIV( * ) REAL A( LDA, * ), B( LDB, * ) donde N: Orden de la matriz A y el número de ecuaciones en el sistema. NRHS: Número de columnas de la matriz B. A: Matriz de coeficientes. Cuando la subrutina retorna, A contiene los factores L y U resultado de la factorización. Nótese que los elementos unitarios en la diagonal de L no se almacenan. LDA: Número máximo de filas que se pueden almacenar en la matriz A. Debe ser mayor o igual a N. IPIV: Arreglo de dimensión N en que se retornan las permutaciones realizadas durante el pivoteo. La fila i de la matriz A fue intercambiada con la fila IPIV(i). Pueden aparecer valores repetidos en IPIV. Por ejemplo, para n=4: é3ù é ù ê 4 ú ê ú IPIV = ê ú Þ P = ê ú ê3ú ê ú ê ú ê ú ë4û ë û B: Una matriz de LDB filas y NRHS columnas, que a la entrada contiene los conjuntos de valores del lado derecho de las ecuaciones, y a la salida la matriz X de soluciones del sistema. LDB: Número máximo de filas que se pueden almacenar en la matriz B. Debe ser mayor o igual a N. INFO: Un valor entero que indica si la subrutina tuvo éxito o no, de la siguiente forma: INFO=0: El sistema pudo ser resuelto sin ningún problema. INFO=-i: El argumento i-ésimo pasado a la subrutina era inválido. Por ejemplo, si se pasa un N menor que 0, se retorna INFO=-1, o si se pasa LDA menor que N, se retorna INFO=-4. INFO=i: La posición (i,i) del factor U calculado es 0, por lo que U es singular. Debido a esto, tratar de resolver el sistema daría como resultado una división por cero. LAPACK cuenta con once variedades de rutinas para la solución de sistemas de ecuaciones lineales, uno para cada tipo y almacenamiento de matriz. Además, para cada tipo se cuenta con rutinas simples y expertas para cuatro clases de datos diferentes. 4.2 Valores y Vectores Propios El problema estándar de valores propios o simplemente problema de valores propios para una matriz A consiste en encontrar escalares (valores propios) l y vectores propios correspondientes z no nulos tales que Az = lz. A los vectores z que cumplen esta ecuación también se les llama vectores propios derechos asociados con el valor propio l y los vectores propios izquierdos asociados con el valor propio l se definen como los vectores u no nulos H H tales que u A= lu. La matriz A puede ser real o compleja, lo mismo que los escalares l y los vectores z. Se distinguen tres problemas generalizados de valores propios para dos matrices simétricas A y B. Consisten en encontrar escalares (valores propios generalizados) l y vectores propios generalizados correspondientes z no nulos tales que: Problema 1: Az=lBz Problema 2: ABz=lz Problema 3: BAz=lz Para matrices cuadradas no simétricas A y B, se distinguen dos problemas generalizados de valores propios, a saber: 1. Encontrar escalares (valores propios generalizados) l y vectores propios generalizados correspondientes x no nulos tales que Ax = lbx. 2. Encontrar escalares (valores propios generalizados) r y vectores propios generalizados correspondientes y no nulos tales que ray = Bx. Con fines ilustrativos, se presentan algunos detalles de la definición de una subrutina para resolver problemas de valores propios. Se trata de la subrutina SSYEV que se encarga de calcular todos o solamente algunos de los valores propios y vectores propios correspondientes para una matriz real simétrica. Se basa en factorización QR. REVISTA Universidad EAFIT No. 123 julio.agosto.septiembre 2001
7 78 Esta rutina es apenas una pequeña parte de la amplia colección de rutinas que tiene LAPACK para el complicado problema de aproximar valores y vectores propios. Algunas de las rutinas disponibles se basan en algoritmos recientemente desarrollados, como RRR (Relatively Robust Representation) o en algoritmos que sólo recientemente están llamando la atención de los investigadores como el algoritmo "Divide y vencerás". La subrutina SSYEV es el manejador simple para encontrar valores y/o vectores propios de una matriz simétrica A. Se basa en los cuatro pasos siguientes: 1. Escalar la matriz, si es necesario, de acuerdo a la precisión y al epsilon de la máquina. Para los numerales 2 y 3 siguientes, se asume que, después del escalamiento, la matriz se sigue llamando A. 2. Reducir la matriz a su forma tridiagonal, mediante una llamada a SSYTRD. Es decir, se trata de encontrar una matriz T ortogonal Q tal que Q AQ = T, con T tridiagonal simétrica. 3. Calcular los valores propios de A, mediante una llamada a SSTERF, que utiliza una variante del algoritmo QR para encontrar valores propios de matrices tridiagonales simétricas. Si también se requiere de los vectores propios, entonces primero se genera la matriz ortogonal Q del paso anterior mediante una llamada a la subrutina SORGTR, y luego se computan valores y vectores propios de T, con una sola llamada a la subrutina SSTEQR. Los valores propios de A son los mismos de T. Los vectores propios correspondientes del problema con matriz A, se consiguen a partir de los que se acaban de calcular recurriendo a la transformación ortogonal del paso anterior. 4. Se reescala la matriz nuevamente si se realizó algún escalamiento en el paso 1. El prototipo de la rutina SSYEV es: SUBROUTINE SSYEV (JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, INFO ) CHARACTER INTEGER JOBZ, UPLO INFO, LDA, LWORK, N REAL A ( LDA, * ), W( * ), WORK( * ) donde: JOBZ: Un caracter de control que especifica que tareas se desea realizar. Los valores posibles son: N: Calcular sólo los valores propios V: Calcular tanto los valores como los vectores propios UPLO: Un caracter de control que indica si se almacenó la parte triangular superior U o la parte triangular inferior L de A. N: Orden de la matriz A. A: Matriz que se desea usar. LDA: Número de filas en la matriz A. Debe ser mayor o igual a N. W: Un vector de dimensión N en el que se retornan los valores propios de A, en orden ascendente. WORK: Vector real que la rutina usa como área de trabajo. LWORK: Longitud del vector WORK. INFO: Valor entero que indica si hubo éxito o no de la siguiente forma: INFO=0: Salida exitosa. INFO=-i: El argumento i-ésimo de entrada es inválido. INFO=i: El algoritmo falló. Un elemento de la diagonal i-ésima en una matriz tridiagonal intermedia no converge a cero. 5. MANEJO DE BLOQUES EN LAPACK Para mejorar el desempeño de LAPACK y comprender cómo funciona la colección, es importante entender cómo manipula las matrices sobre las que opera. Por qué bloques? Existen muchas razones por las cuales puede ser útil dividir las matrices en subbloques. La más importante es la ganancia en velocidad que se obtiene al utilizar más eficientemente la jerarquía de memoria que se encuentra en la gran mayoría de los equipos actuales. La mayor parte de las máquinas usadas actualmente tienen una jerarquía de memoria de por lo menos cinco niveles: 1. Disco (Memoria Virtual) 2. Memoria Principal (física)
8 79 3. Memoria Cache L2 4. Memoria Cache L1 5. Registros Por el momento, se hará énfasis principalmente en los 3 niveles intermedios. Cada nivel es aproximadamente un orden de magnitud más lento que el nivel inmediatamente superior, lo cual indica que hay gran potencial para incrementar el desempeño de los algoritmos si podemos mantener los datos sobre los que operan en el nivel más alto posible. Esto es conocido como el principio de localidad espacial. El dividir las matrices en subbloques y operar sobre éstos de manera independiente permite aprovechar el principio de localidad al operar sobre conjuntos de datos más pequeños cuyos elementos están relativamente muy cercanos. Como ejemplo, considérese una operación aparentemente simple como es la multiplicación de matrices. Supongamos una máquina con una configuración común hoy en día: un cache L2 de 512KB y un par de matrices relativamente grandes de 1000 x Supóngase también que cada elemento en cada una de las matrices se almacena como un valor de punto flotante de doble precisión, necesitando 8 bytes (64-bits) de almacenamiento cada uno. Según esto, en el mejor de los casos, el cache L2 no podría retener en un momento dado ni el 3% de cada una de las matrices. Si se considera el algoritmo usual (manual) para multiplicar matrices, se notará que para la matriz que aparece en el lado derecho de la operación, las referencias a memoria son a datos no contiguos, por lo que el cache se vería obligado a remplazar bloques con mucha frecuencia. El desempeño se reduciría notablemente debido al incremento en la rata de "cache misses". Es de notar además, que estos cálculos son, en realidad, optimistas, ya que usualmente, los caches L2 son unificados para datos e instrucciones, lo cual reduce aún más el espacio efectivo que se puede utilizar. Casos como éstos requieren de algoritmos capaces de operar sobre porciones relativamente pequeñas de una matriz a la vez, y así tomar plena ventaja del principio de localidad espacial. El éxito de esta estrategia depende de dos factores críticos: a. Tamaño de Bloque (NB): para aprovechar el cache de la máquina al máximo, se necesita escoger un NB tal que el cache sea capaz de acomodar fácilmente todos los bloques sobre los que se está operando. Por ejemplo, con un tamaño de bloque de 64 bytes y un cache de 512KB, se podrían acomodar tres bloques completos de 8x8 y habría espacio de sobra. Es de notar que la selección de NB también depende del proceso numérico que se desee realizar, así como de las restricciones o posibilidades que ofrezca el algoritmo desarrollado. b. Estabilidad: Es claro que, en muchos casos, de nada sirve desarrollar algoritmos super-veloces si se pierden la estabilidad y precisión numérica de las que depende el éxito de la operación que se realiza. De poco sirve obtener resultados muy rápido si éstos son completamente erróneos. Una pregunta que se presenta aquí es: Por qué no aprovechar la estructura que presentan algunas matrices para seleccionar el tamaño de bloque? Algunas matrices, notablemente las diagonales por bloques, se prestan a ser procesadas de esta manera. Sin embargo, actualmente LAPACK no posee rutinas especializadas para operar sobre este tipo de matrices, ni, como se verá en detalle más adelante, tiene en cuenta esta información en la escogencia de NB. LAPACK implementa los algoritmos por bloques en dos niveles separados: las BLAS y las rutinas computacionales. Las "Basic Linear Algebra Subprograms" o BLAS implementan las operaciones básicas de álgebra lineal de manera independiente del resto de LAPACK. Esto implica que es fácil remplazar las BLAS originales escritas en FORTRAN por versiones altamente optimizadas para la plataforma en la que se está trabajando. Las operaciones que las BLAS implementan se dividen, como se vio al principio, en operaciones de niveles 1 (vector-vector), 2 (matriz-vector) y 3 (matriz-matriz). De manera análoga, las rutinas computacionales están divididas en versiones diferentes, dependiendo de cuál categoría de funciones BLAS invoquen. Las rutinas que hacen uso del nivel 3 de las BLAS son aquellas que más relevantes resultan en esta discusión, puesto que son las que realmente implementan los algoritmos por bloques. Las rutinas computacionales y las BLAS de nivel 3 tienen una característica particular: implementan algoritmos O(n 3 ) cuando bastaría con O(n 2 ). Inicialmente, podría parecer que esto va a afectar negativamente el desempeño, pero dadas las condiciones antes mencionadas, es posible lograr todo lo contrario. Cómo? Las BLAS de nivel 3 ejecutan más operaciones, pero de forma tal que REVISTA Universidad EAFIT No. 123 julio.agosto.septiembre 2001
9 80 se realice todo el trabajo posible de una vez sobre los datos traídos de memoria. Con esto, se puede aumentar el factor q, que es la razón de operaciones realizadas sobre los datos con respecto a las referencias de memoria requeridas. En otras palabras, q es un estimador de que tan bien se puede aprovechar el principio de localidad espacial. En algunos casos, las rutinas computacionales con bloques no pueden invocar a las funciones de BLAS 3, como cuando el usuario fuerza el tamaño de bloque 1. En estos casos las rutinas invocan a sus contrapartes sin bloques para hacer el trabajo, que sólo invocan funciones de nivel 1 ó 2 de las BLAS. En LAPACK, el trabajo de seleccionar el tamaño de bloque apropiado recae sobre la rutina ILAENV. Esta rutina devuelve un valor para NB de acuerdo a dos factores: a. Subrutina que solicita el valor de NB. b. Qué tan óptimo debe ser el valor retornado para NB por esta rutina? El valor retornado por ILAENV es una constante predefinida por los autores de LAPACK, hallada de forma experimental como un valor óptimo del tamaño de bloque para las máquinas en que se escribió LAPACK originalmente. Es concebible pensar que, dados los avances en la velocidad de procesamiento y la complejidad de los procesadores en los últimos diez años, sea factible usar valores más grandes. Es una lástima, sin embargo, que ILAENV no haga parte de las BLAS sino de LAPACK propiamente dicha, ya que así al reemplazar las BLAS por una versión optimizada para nuestra máquina, se obtendría un tamaño de bloque también óptimo. 6. CONCLUSIONES LAPACK es una colección de rutinas que combina precisión y eficiencia. Usar las rutinas de LAPACK es fácil. Se considera que los cursos de métodos numéricos que se dictan en las Universidades locales pueden complementarse con el uso de algunas de las rutinas de LAPACK. De esta forma los estudiantes, y futuros profesionales, estarán conscientes de la existencia de unas rutinas de alta calidad, de uso fácil y gratuitas que podrán utilizar más adelante cuando la ocasión lo amerite. BIBLIOGRAFÍA Anderson E. et al. (1992). LAPACK Users' Guide. USA: SIAM. Demmel J.W. (1997). Applied Numerical Linear Algebra. USA: SIAM. Golub G. and Van Loan C. (1989). Matrix Computations, 2a. ed. USA: Johns Hopkins University Press. Hayes J.P. (1988). Computer Architecture and Organization. USA: McGraw Hill. Noble B. and Daniel J. W. (1989). Algebra Lineal Aplicada, 3a. ed. USA: Prentice Hall. Stewart G.W. (1973). Introduction to Matrix Computations. New York: Academic Press. Stewart G.W. (1998). Matrix Algorithms - Volume I: Basic Decompositions", USA: SIAM. Trefethen L.N. and Bau D. (1997). Numerical Linear Algebra. USA: SIAM.
Linear Algebra PACKage (LAPACK) Avances en la Generación de Bibliotecas de Álgebra Lineal Universidad Politécnica de Valencia Marzo, 2006
Linear Algebra PACKage () Avances en la Generación de Bibliotecas de Álgebra Lineal Universidad Politécnica de Valencia Marzo, 2006 Estructura Qué es la biblioteca? Organización de Funcionalidad de Sistemas
Tareas de matrices especiales
Tareas de matrices especiales Objetivos. Estudiar una clase especial de matrices. Para matrices de esta clase realizar un algoritmo rápido de multiplicación por vectores. Aplicar este algoritmo para resolver
Computación Científica Gustavo Rodríguez Gómez INAOE Agosto Dicembre 2011 1 / 46 Capítulo II 2 / 46 1 Introducción Métodos Directos Sistemas Triangulares Sustitución Hacia Atrás Invertibilidad de una Matriz
Utilización de CLAPACK Para Resolver Sistemas de Ecuaciones Lineales Mediante Paralelismo y Optimización de Memoria
Utilización de CLAPACK Para Resolver Sistemas de Ecuaciones Lineales Mediante Paralelismo y Optimización de Memoria Germán Harvey Alférez Salinas 1 y Edward Mauricio Alférez Salinas 2 1. Facultad de Ingeniería
Clase. 1. Resolución de sistemas de ecuaciones lineales: preliminares
Clase 1. Resolución de sistemas de ecuaciones lineales: preliminares 2. Método directo y exacto: Gauss 3. Método directo y exacto (II): descomposición LU 4. Métodos indirectos: Jacobi, Gauss-Seidel 2 Sistemas
Eigenvalores y eigenvectores Los dos problemas principales del álgebra lineal son: resolver sistemas lineales de la forma Ax = b y resolver el problema de eigenvalores. En general, una matriz actúa sobre
Inversión de Matrices Daniel Vaughan Es bien conocido que en diversas aplicaciones de contabilidad nacional, así como en otras áreas de la economía, es usual encontrarse con la inversión de matrices. Ejemplos
Cálculo Numérico - CO3211. Ejercicios. d ) Sabiendo que la inversa de la matriz A = es A c d
Cálculo Numérico - CO32 Ejercicios Decida cuáles de las siguientes proposiciones son verdaderas y cuáles son falsas Si una proposición es verdadera, demuéstrela, y si es falsa dé un contraejemplo: a Sea
Unidad 5 Problemas de álgebra lineal numérica.
Unidad 5 Problemas de álgebra lineal numérica. Eliminación gaussiana. Factorización LU. Sea Ax = b un sistema de n ecuaciones lineales con n incógnitas, donde A es una matriz cuadrada de orden n no singular
LABORATORIO DE COMPUTACIÓN CIENTÍFICA (Prácticas) Curso 2009-10 1 MÓDULO SE: SISTEMAS DE ECUACIONES Alumno: Lee detenidamente los enunciados. Copia las funciones y scripts que crees a lo largo de la practica,
Producto escalar. x y. x = x x y cos α =
resumen06 1 Producto escalar Vectores ortogonales y proyecciones La definición matemática de producto escalar es bastante amplia porque recoge toda expresión bilineal que sirva razonablemente para medir
Sistemas de Ecuaciones. Lineales III
Sistemas de Ecuaciones Lineales III Pivoteo: Estrategia de pivoteo parcial Adaptación a matrices con estructuras particulares: Matrices banda y tridiagonales Método de holesky 52123-1 - DIM Universidad
Solución de Ecuaciones Simultáneas por el Método de Matrices
Solución de Ecuaciones Simultáneas por el Método de Matrices Universidad de San Carlos de Guatemala Facultad de Ciencias Químicas y Farmacia Matemática IV Rony José Letona QQ 200960024 ÍNDICE ÍNDICE Índice.
I. Métodos directos para resolución de SEL. Se dice que una matriz A admite una factorización LU indirecta A = LU
I. Métodos directos para resolución de SEL 1. Factorización LU Se dice que una matriz A admite una factorización LU si dicha matriz puede escribirse como el producto de una matriz triangular inferior,
Motivación La cercha de la figura se carga con una fuerza uniforme repartida sobre el cordón superior Resolución numérica de sistemas de ecuaciones. Introducción El planteamiento del problema conduce a
4 o Grado. Multiplicación y la División. Revisión de Multiplicación. Relación entre la. Slide 2 / 105. Slide 1 / 105. Slide 4 / 105.
Slide 1 / 105 Slide 2 / 105 New Jersey Center for Teaching and Learning Iniciativa de Ciencia Progresiva Este material está disponible gratuitamente en www.njctl.org y está pensado para el uso no comercial
Análisis Numérico para Ingeniería. Clase Nro. 5
Análisis Numérico para Ingeniería Clase Nro. 5 Sistemas de Ecuaciones Lineales Temas a tratar: Tipos de Matrices Método de Triangulación de Gauss Método de Diagonalización de Gauss-Jordan Tácticas de Pivoteo
Edgar Acuña/ ESMA 6665 Lecc1-2 1 ESTADISTICA COMPUTACIONAL. Capítulo I. Matrices y solución de Ecuaciones lineales
Edgar Acuña/ ESMA 6665 Lecc1-2 1 ESTADISTICA COMPUTACIONAL Capítulo I. Matrices y solución de Ecuaciones lineales Referencias: 1. Datta, B. (1995) Numerical Linear Algebra. Brooks Cole 2. Golub, G. and
MAGMA. Matrix Algebra on GPU and Multicore Architecture. Ginés David Guerrero Hernández
PLASMA GPU MAGMA Rendimiento Trabajo Futuro MAGMA Matrix Algebra on GPU and Multicore Architecture Ginés David Guerrero Hernández gines.guerrero@ditec.um.es Grupo de Arquitecturas y Computación Paralela
1.IV Aproximación numérica de valores y vectores propios.
.IV Aproximación numérica de valores y vectores propios. Manuel Palacios Departamento de Matemática Aplicada Centro Politécnico Superior. Universidad de Zaragoza Primavera 2007 Contents Introducción 2
La importancia de las factorizaciones matriciales no negativas en la minería de datos y el procesamiento de
La importancia de las factorizaciones matriciales no negativas en la minería de datos y el procesamiento de imágenes Humberto Madrid, Irma García, Federico Garza Centro de Investigación en Matemáticas
Cálculo numérico. Sistemas de ecuaciones lineales.
José Luis Morales http://allman.rhon.itam.mx/ jmorales Departamento de Matemáticas. ITAM. 2010. Las raíces de x 2 bx + c = 0. r = b ± b 2 4c 2 b = 3.6778, c = 0.0020798 r 1 = 3.67723441190... r 2 = 0.00056558809...
INSTITUTO VALLADOLID PREPARATORIA página 93
INSTITUTO VALLADOLID PREPARATORIA página 93 página 94 CONCEPTO Una ecuación es una especie de "adivinanza numérica", o sea que se hace un planteamiento cuya respuesta debe ser un número. Por ejemplo: "
Los algoritmos numéricos
Enseñanza Los algoritmos numéricos José Guerrero Grajeda Rosa Margarita y Álvarez González UNAM, UAQ UNAM resumen Se define el concepto de algoritmo Se presentan varios casos de problemas numéricos, se
Algoritmos paralelos para la Multiplicación de Matrices
Introducción Paralelizar onclusiones Trabajo ibliografía lgoritmos paralelos para la Multiplicación de Matrices Universidad de Murcia - UM 18 de diciembre de 2008 Introducción Paralelizar onclusiones Trabajo
ARITMETICA MODULAR: Una Aritmética Divertida Luis F. Cáceres
ARITMETICA MODULAR: Una Aritmética Divertida Luis F. Cáceres La idea de número debió surgir de la necesidad que tenía el hombre de llevar registro de cosas importantes del diario vivir. Cuántas ovejas
Diseño de algoritmos paralelos Curso 2011-2012 Esquema del capítulo Visión general de algunos algoritmos serie. Algoritmo paralelo vs. Formulación paralela Elementos de un Algoritmo paralelo Métodos de
Clase No. 13: Factorización QR MAT 251. Joaquín Peña (CIMAT) Métodos Numéricos (MAT 251) / 16
Clase No 13: Factorización QR MAT 251 Joaquín Peña (CIMAT) Métodos Numéricos (MAT 251) 03102011 1 / 16 Factorización QR Sea A R m n con m n La factorización QR de A es A = QR = [Q 1 Q 2 ] R1 = Q 0 1 R
Herramienta de distribución libre para Cómputo Científico como. Científico como opción viable en la investigación científica y tecnológica
Herramienta de distribución libre para Cómputo Científico como opción viable en la investigación científica y tecnológica May 17, 2010 Contenido 1 Es indudable que el rápido avance en el desarrollo de
Ejercicios: Estimando π con métodos Monte Carlo y cálculo del mayor autovalor de una matriz
Ejercicios: Estimando π con métodos Monte Carlo y cálculo del mayor autovalor de una matriz 6 de marzo de 05. Cálculo de π con métodos Monte Carlo El objetivo de este ejercicio consiste en estimar el valor
INTRODUCCIÓN AL USO DE PAQUETES COMPUTACIONALES Profesoras Especialista Bernal Herrera Fanny Maricela Especialista Franco Valderrama Dalia Elizabeth Página 1 3.1 Objetivo particular del tema 3.0 Tema II
Capítulo 2 Resolución de sistemas lineales 2. Introducción Básico La mayor parte de este escrito trata sobre la forma en que las ecuaciones de campo son aproximadas por sistemas de ecuaciones algebraicas.
ESTIMACIÓN DE TIEMPO Y COSTO DE PRODUCTOS SOFTWARE
Análisis y costo de algoritmos Algoritmos Un algoritmo es un conjunto de pasos claramente definidos que a partir de un cierta entrada (input) produce una determinada salida (output) Algoritmos y programas
Descomposición LU de matrices José L. Vieitez IMERL, Facultad de Ingeniería, Universidad de la República 3 de agosto de 006 Abstract Descomposición LU de una matriz A = (a ij ), i, j = 1,..., n. 1 Hechos
DESARROLLO DE UN SOFTWARE PARA LA GENERACIÓN DE ESPECTROS ELASTICOS DE RESPUESTA.
DESARROLLO DE UN SOFTWARE PARA LA GENERACIÓN DE ESPECTROS ELASTICOS DE RESPUESTA. Carlos M. PISCAL. Wilson RODRIGUEZ. Andrés M. LOTERO Programa de Ingeniería Civil Universidad de la Salle Bogotá. Colombia

References: Resolución 
 Resolución 
 resolución 
 resolución 
 Resolución 
 Resolución