CÁLCULO NUMÉRICO DE AUTOVALORES USANDO ALGORITMOS PARALELOS EN PROBLEMAS DE QUÍMICA CUÁNTICA.
Descripción
UNIVERSIDAD DE CARABOBO ÁREA DE ESTUDIOS DE POSTGRADO FACULTAD DE INGENIERÍA MAESTRÍA EN MATEMÁTICAS Y COMPUTACIÓN
CÁLCULO NUMÉRICO DE AUTOVALORES USANDO ALGORITMOS PARALELOS EN PROBLEMAS DE QUÍMICA CUÁNTICA.
Autor: Li . JOSÉ LUIS RAMÍREZ B. CI: 12.035.391 Tutor: Dr. GERMÁN LARRAZÁBAL CI: 6.131.761
Valen ia, O tubre de 2005.
UNIVERSIDAD DE CARABOBO ÁREA DE ESTUDIOS DE POSTGRADO FACULTAD DE INGENIERÍA MAESTRÍA EN MATEMÁTICAS Y COMPUTACIÓN
CÁLCULO NUMÉRICO DE AUTOVALORES USANDO ALGORITMOS PARALELOS EN PROBLEMAS DE QUÍMICA CUÁNTICA.
Autor: Li . JOSÉ LUIS RAMÍREZ B. CI: 12.035.391 Trabajo de Grado presentado ante el Área de Estudios de Postgrado de la Universidad de Carabobo para optar al título de Magister en Matemáti as y Computa ión.
Valen ia, O tubre de 2005. ii
UNIVERSIDAD DE CARABOBO ÁREA DE ESTUDIOS DE POSTGRADO FACULTAD DE INGENIERÍA MAESTRÍA EN MATEMÁTICAS Y COMPUTACIÓN
CÁLCULO NUMÉRICO DE AUTOVALORES USANDO ALGORITMOS PARALELOS EN PROBLEMAS DE QUÍMICA CUÁNTICA.
Autor: Li . JOSÉ LUIS RAMÍREZ B. CI: 12.035.391 Aprobado en el Área de Estudios de Postgrado de la Universidad de Carabobo por Miembros de la Comisión Coordinadora del Programa: Ing. Ri ardo Olivero Prof. Katty Ordaz Prof. Melba Rodríguez
Valen ia, O tubre de 2005. iii
UNIVERSIDAD DE CARABOBO ÁREA DE ESTUDIOS DE POSTGRADO FACULTAD DE INGENIERÍA MAESTRÍA EN MATEMÁTICAS Y COMPUTACIÓN
VEREDICTO Nosotros, Miembros del Jurado designado para la evalua ión del Trabajo de Grado titulado:
CÁLCULO NUMÉRICO DE AUTOVALORES USANDO ALGORITMOS
PARALELOS EN PROBLEMAS DE QUÍMICA CUÁNTICA, presentado por el Li . José Luis Ramírez B., para optar al título de Magíster en Matemáti as y Computa ión, estimamos que el mismo reúne los requisitos para ser onsiderado omo:
Miembros del Jurado: Apellido, Nombre C.I.
Apellido, Nombre C.I.
Apellido, Nombre C.I.
Valen ia, O tubre de 2005. iv
DEDICATORIA
Este
trabajo
está
dedi ado
prin ipal-
mente a mi hija Fran ys y mi esposa Fran ia por estar a mi lado apoyandome.
Igualmente
está
dedi ado
a
mi
Padre
José Luis y a mi Madre Leisy, por apoyarme en todas mis etapas de estudio.
Y
a
mi
hermana
Móni a
y
mi
ahijado
Daniel por ser grandes amigos.
Finalmente
a
toda
mi
familia
por
ayu-
darme uando más lo ne esito, por ser una familia in ondi ional.
v
ÍNDICE GENERAL
ÍNDICE DE CUADROS ÍNDICE DE FIGURAS RESUMEN
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix ...........................................................x
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi
INTRODUCCIÓN
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xii
CAPÍTULO I: EL PROBLEMA
................................................. 1
1 Planteamiento del Problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 2 Formula ión del Problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 3 Objetivos de la Investiga ión . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 3.1 Objetivo General . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 3.2 Objetivos Espe í os . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 4 Justi a ión de la Investiga ión . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
CAPÍTULO II: MARCO TEÓRICO
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1 Ante edentes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2 Bases Teóri as . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.1 Forma Diagonal de una Matriz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.2 Método de Ja obi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.3 Itera ión
QR
y Autovalores . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.3.1 Tridiagonaliza ión de Matri es . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 2.3.2 Algoritmo
QR
y
QL
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
2.3.3 GIVENS y la Fa toriza ión
QR
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
2.3.4 A elera ión de la Convergen ia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 2.3.4.1 Desplazamiento Simple . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 2.3.4.2 Desplazamiento Implí ito . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 2.3.5 Criterio de Parada . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 2.4 Computa ión Paralela . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 2.4.1 Computador Paralelo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
vi
2.4.2 Programa ión Multithread . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 2.5 Estru tura de Alma enamiento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 2.5.1 Compressed Row Storage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 2.5.2 Compressed Column Storage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 3 De ni ión de Términos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
CAPÍTULO III: MARCO METODOLÓGICO
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
1 Tipo y Diseño de la Investiga ión . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 2 Modalidad de la Investiga ión . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 3 Té ni as e Instrumentos para Re ole tar la Informa ión . . . . . . . . . . . . . . . . . . . . . . . . 57 4 Fases Metodológi as . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 4.1 Fase I: Fo alizar . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 4.2 Fase II: Analizar . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 4.3 Fase III: Diseñar . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 4.4 Fase IV: Comparar . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 5 Plan de Trabajo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 5.1 Re ursos Humanos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 5.2 Re ursos Finan ieros . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 5.3 Re ursos Institu ionales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 6 Resultados de las Fases Metodológi as . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 6.1 Fase I . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 6.2 Fase II . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 6.3 Fase III . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 6.3.1 Plataforma Computa ional . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 6.3.1.1 Bibliote a UCSparseLib . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 6.3.1.2 El Estándar OPENMP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 6.3.2 Crea ión del Módulo EIGENVALUES para UCSparseLib . . . . . . . . . . . . 70 6.3.2.1 Considera iones de Implementa ión . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
vii
6.3.3 Tridiagonaliza ión de Householder y uso del paralelismo . . . . . . . . . . . . . . 75
CAPÍTULO IV: RESULTADOS
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
1 FASE IV . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 1.1 Resultados Numéri os . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 1.1.1 Cál ulo Se uen ial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 1.1.2 Cál ulo Paralelo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
CAPÍTULO V: CONCLUSIONES Y RECOMENDACIONES
. . . . . . . . . . . . . . 89
1 Con lusiones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 2 Re omenda iones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
ANEXOS
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
A. MÓDULO EIGENVALUES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 B. PROGRAMA PRINCIPAL . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
BIBLIOGRAFÍA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
viii
ÍNDICE DE CUADROS 1.- Fases Metodológi as . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 2.- Re ursos Humanos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 3.- Re ursos Finan ieros . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 4.- Re ursos Institu ionales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 5.- Modelos de Programa ión Paralela . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 6.- Matri es de Prueba. Matrix Market . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 7.- Matri es de Prueba. Matrix Market . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 8.- Matri es de Prueba. Matrix Market . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 9.- Matri es de Prueba. Matrix Market . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 10.- Matri es de Prueba. CATIVIC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 11.- Tiempo Se uen ial. Caso: Matrix Market . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 12.- Tiempo Se uen ial. Caso: CATIVIC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 13.- Se uen ial vs Paralelo (Autovalores). Caso: Matrix Market . . . . . . . . . . . . . . . . . . . . . . 87 14.- Se uen ial vs Paralelo (Autovalores y Autove tores). Caso: Matrix Market . . . . . . . 87 15.- Se uen ial vs Paralelo (Autovalores). Caso: CATIVIC . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 16.- Se uen ial vs Paralelo (Autovalores y Autove tores). Caso: CATIVIC . . . . . . . . . . . . 88
ix
ÍNDICE DE FIGURAS 1.- Clasi a ión General de las Matri es . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.- Bibliote a UCSparseLib y módulos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 3.- Algoritmo
QR
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
4.- Tridiagonaliza ión y Diagonaliza ión de una Matriz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 5.- Algoritmo Se uen ial de Tridiagonaliza ión . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 6.- Algoritmo
QR
on Desplazamiento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
7.- Modelo de paralelismo SISD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 8.- Modelo de paralelismo SIMD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 9.- Modelo de paralelismo MISD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 10.- Modelo de paralelismo MIMD on Memoria Compartida . . . . . . . . . . . . . . . . . . . . . . . . 45 11.- Modelo de paralelismo MIMD on Memoria Distribuida . . . . . . . . . . . . . . . . . . . . . . . . . 45 12.- Modelo de paralelismo MIMD on Memoria Compartida Distribuida . . . . . . . . . . . . 46 13.- Aproxima ión de la E ua ión de S hrödinger . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 14.- Modelo Fork-Join . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 15.- Modelo de Paraleliza ión In remental . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
x
UNIVERSIDAD DE CARABOBO ÁREA DE ESTUDIOS DE POSTGRADO FACULTAD DE INGENIERÍA MAESTRÍA EN MATEMÁTICAS Y COMPUTACIÓN CÁLCULO NUMÉRICO DE AUTOVALORES USANDO ALGORITMOS PARALELOS EN PROBLEMAS DE QUÍMICA CUÁNTICA. Autor: Li . José Luis Ramírez B. Tutor: Dr. Germán Larrazábal. Fe ha: Julio, 2005. RESUMEN En este trabajo son implementados los algoritmos más utilizados para el ál ulo de autovalores de una matriz
A simétri a de n×n on la ara terísti a prin ipal de ser una matriz
dispersa, además de que el orden de la matriz supera los miles, lo ual ha e ne esario utilizar una estru tura e iente para el manejo de matri es dispersas y adaptar los métodos de
ál ulo de autovalores a tales estru turas, se emplearán algoritmos que permitan obtener de la matriz todos sus autovalores. En bus a de redu ir los tiempos de eje u ión, estos algoritmos están implementados apli ando paralelismo. El método
QR
o
QL
de Fran is
tiene la ara terísti a de diagonalizar la matriz a través de matri es similares, por lo que los autovalores de la matriz original son preservados en ada itera ión, además de ser el método por ex elen ia para el ál ulo de los mismo y sus autove tores. En este ontexto, el estudio es realizado en modernas super omputadoras que permiten eje utar más de una instru
ión y que a la vez permiten pro esar múltiples datos, y en onjunto on la bibliote a UCSparseLib de la Universidad de Carabobo, la ual ya tiene las estru turas ne esarias y e ientes para el manejo de estru turas dispersas, se bus a mejorar el tiempo de eje u ión del algoritmo serial, apli ando
multithreading, usando la bibliote a OpenMP,
de tal manera de obtener resultados ables a problemas de gran magnitud.
Palabras laves: Autovalores, Algoritmo QR, OpenMP, UCSparseLib.
xi
INTRODUCCIÓN En este trabajo son desarrollados los algoritmos
QR
y
QL
para el ál ulo de autovalores
y autove tores, de forma tal que puedan sean apli ados en problemas de la quími a uánti a, en este aso, para resolver de forma aproximada la e ua ión de S hrödinger a través de las e ua iones de Hartree-Fo k por medio del método paramétri o
Self-Consistent-Field.
Esta apli a ión quími a tiene la ara terísti a de que al ser aproximada mediante un modelo matri ial, la matriz de oe ientes sobre la ual se va a realizar el ál ulo de autovalores, es una matriz simétri a y dispersa. Por lo tanto se ne esita que el desarrollo de tales algoritmos se reali e sobre una estru tura que permita manejar de forma e iente las matri es dispersas.
Además se realiza un estudio del algoritmo que ha sido el método lási o para el ál ulo de autovalores y sus respe tivos autove tores uando la matriz es simétri a, este es el método
QR
estudiado por Fran is (1961) y desde enton es ha sido desarrollado por diversos
autores, tales omo Wilkinson, Reins h y Bowdler (1968), entre otros. Dado que el método
QR
es un algoritmo iterativo, se han presentado diversas versiones del mismo, en bus a
de a elerar la onvergen ia, tal omo se plantea en Trefethen y Baun (1997).
Otra ara terísti a de la apli a ión quími a, es que las matri es pueden ser de gran tamaño, de este modo, junto on la ara terísti a de ser dispersa, pare e no ser onveniente alma enar toda la estru tura en un formato denso, más bien se bus a utilizar un formato que permita el alma enamiento de aquellos elementos que realmente son signi ativos para la apli a ión, en este aso, los elementos no nulos. La bibliote a UCSparseLib posee los formatos ne esarios para el alma enamiento de estas matri es, así se logra un ahorro signi ativo en espa io omo lo plantea Larrazabal (2004).
xii
Sin embargo, al utilizar estos formatos de alma enamiento, se presenta un
overhead
que
debe ser disminuido, ya que los algoritmos se uen iales están pensados sobre estru turas densas y no sobre estru turas dispersas, por lo que este
overhead
es ausado por a
esos
no uniformes sobre la data. Para ello, los algoritmos serán implementados sobre arquite turas paralelas, de modo tal de disminuir el
overhead
presentado, en este trabajo se
utiliza la bibliote a OpenMP sobre una arquite tura de multipro esadores on memoria
ompartida. Se realiza el estudio de los distintos algoritmos involu rados para obtener de forma e iente los autovalores de una matriz a través del método
QR,
evaluando uales
son los i los ríti os de los algoritmos.
Este trabajo realizado onsta de in o apítulos; en el primero se presenta el planteamiento del problema, los objetivos propuestos y la justi a ión de la investiga ión. El segundo
apítulo presenta el mar o teóri o sobre el ual se basa la des omposi ión de autovalores por medio de la diagonaliza ión de matri es ha iendo mayor énfasis en el método
QR
y
sus variantes, en este apítulo también se introdu en los on eptos bási os de la omputa ión paralela y de forma similar se plantea el alma enamiento de estru turas dispersas. En el ter er apítulo se exponen las ara terísti as del enfoque metodológi o, expli ando
ada una de las fases involu radas en la investiga ión. El uarto apítulo presenta los resultados obtenidos, tanto en el desarrollo de la investiga ión omo en la apli a ión de los algoritmos paralelos implementados, realizando una ompara ión on el mismo algoritmo en su forma se uen ial. Por último, en el quinto apítulo se presentan las on lusiones y re omenda iones obtenidas al ulminar la investiga ión.
xiii
CAPÍTULO I: EL PROBLEMA
2
CAPÍTULO I EL PROBLEMA 1 Planteamiento del Problema Problemas que involu ren el ál ulo de autovalores de una matriz pueden derivarse de un gran número de dis iplinas de las ien ias y la ingeniería. Ellos onstituyen la herramienta bási a empleada para el diseño de onstru
iones, puentes, y turbinas resistentes a vibra iones, permiten modelar el en olamiento de una red y analizar la estabilidad de redes elé tri as o de ujo de uidos, también permiten entender fenómenos físi os o estudiar patrones de bifur a ión en sistemas dinámi os.
Por lo tanto la solu ión de mu hos problemas físi os, generalmente es redu ida enton es a resolver un problema de autovalores. Por lo tanto, estos problemas tienen onsiderable aten ión y representan una de las más importantes áreas del álgebra lineal numéri a. La Línea de Investiga ión del Programa de Maestría en Matemáti a y Computa ión en la ual se enmar a este trabajo es la que se presenta en el Departamento de Mátemáti as de la Maestría, entre las uales en ontramos la Línea de Investiga ión de Métodos Numéri os, ya que on los mismos se plantea la solu ión a diversos tipos de problemas físi os on ayuda del omputador.
En uanto al problema de autovalores, este puede ser des rito mediante el siguiente sistema de e ua iones:
Ax = λx Donde el es alar triz
(1)
λ es denominado autovalor, y x es el orrespondiente autove tor de la ma-
A. Si A es una matriz simétri a de orden n, enton es A posee exa tamente n autovalo-
3
res reales, y los autove tores orrespondientes forman una base del espa io
n−dimensional.
Por ello el uso del omputador en la resolu ión de problemas de autovalores y valores singulares ha ondu ido a dos aspe tos importantes de investigar: velo idad y exa titud. Debido a la ne esidad de resolver grandes problemas, el primer aspe to de investiga ión es obtener algoritmos más rápidos y el análisis de la velo idad de onvergen ia del método en estudio. Mientras que el problema de la exa titud puede ser estable ido omo la antidad de dígitos signi ativos que posee el valor al ulado, omo es señalado por Kin aid y Cheney (1994) en su texto.
El estudio del autovalores o espe tro (para matri es y operadores lineales, respe tivamente), revelan informa ión importante a er a del omportamiento del sistema, ya sea lineal o no lineal. En el trabajo de Trefethen (1999) se muestran resultados a er a de la evolu ión del ál ulo de autovalores, sobre todo, en este trabajo se reseña omo en 1990 la obten ión del espe tro de una matriz
30 × 30
tomó una importante antidad de minu-
tos, mientras que a tualmente se espera obtener similares resultados pero on matri es de mayor orden.
Con la apari ión de omputadoras paralelas, mu hos algoritmos han sido rediseñados para las nuevas arquite turas. Mu hos trabajos desarrollados en el álgebra lineal numéri a están basados en la implementa ión de
software
para omputadoras de arquite tura avanzada,
omo es planteado por Dongarra y Eijkout (2000), el ual se entra en uatro he hos bási os: 1.- motiva ión del trabajo; 2.- desarrollos estándar para su uso en el álgebra lineal; 3.- aspe tos referentes al diseño de algoritmos e implementa iones paralelas; y 4.- futuras dire
iones de investiga ión.
El problema de obtener los autovalores de una matriz uadrada
A,
es la determina ión de
4
solu iones no triviales de la e ua ión (1). El ál ulo de autovalores
λ
de manera explí ita
onstruyendo la e ua ión ara terísti a.
det(A − λI) = 0
(2)
La resolu ión de la e ua ión (2) no es la mejor op ión, ex epto en asos espe iales, dado que los oe ientes de la e ua ión ara terísti a no pueden ser al ulados de manera estable. Y aún, si la e ua ión ara terísti a (2) pudiese ser determinada de manera exa ta, el
ál ulo de sus raí es, en una pre isión nita, puede ser altamente inestable a pequeñas perturba iones en los oe ientes lo ual puede ondu ir a grandes errores en los resultados nales. Así pues, la des omposi ión de autovalores se en uentra entre las más difí iles de al ular, y es en general obtenida omo un resultado aproximado proveniente de un algoritmo iterativo. Por lo tanto diferentes onsidera iones omputa ionales deben ser tomadas en uenta al momento de resolver un problema aso iado a autovalores. Mu hos de los métodos han sido desarrollados aprove hando las ara terísti as de la matriz de
oe ientes. En la gura 1 se muestran algunas de las ara terísti as, que podrían afe tar el algoritmo a desarrollar.
Figura 1: Clasi a ión General de las Matri es
Los elementos de interés también forman parte de las ara terísti as de los algoritmos a desarrollar, entre estos elementos tenemos los siguientes:
Todos los autovalores y posibles autove tores.
5
Solo el autovalor más grande en magnitud y posible autove tor.
Solo el autovalor más pequeño en magnitud y posible autove tor.
Aquellos autovalores tal que
Real(λj ) < 0.
Aquellos autovalores tal que
|λj − α| < β .
Métodos para omputadores de arquite turas diferentes.
La resolu ión de problemas rela ionados on autovalores, fre uentemente involu ran matri es de gran tamaño. Saad (1996) ha e notar que la frase "gran tamaño", es relativa y varía rápidamente en onjunto on el avan e de la omputa ión. De esta manera una matriz formada por ientos de las y ientos de olumnas puede ser onsiderada grande si se usa un omputador personal, mientras que, de manera similar, una matriz formada por miles de las y miles de olumnas puede ser onsiderada grande si se trabaja en un super omputador.
Mu has de las matri es que apare en en problemas prá ti os poseen una estru tura dispersa, es de ir, en ellas se en uentran una gran antidad de elementos iguales a ero, o lo que es igual, po a antidad de elementos signi ativos, omo es referido en Saad (1996). De esta de ni ión se puede de ir que una matriz tridiagonal posee una estru tura dispersa y a la vez triangular, lo ual no pare e ser muy onvin ente.
University de Carabobo Sparse Li-
Por otra parte la bibliote a numéri a UCSparseLib (
brary ), es una bibliote a que permite resolver sistemas lineales dispersos, la ual fue desarrollada por Larrazábal (2004). En este trabajo se de nen los formatos de alma enamiento
Compressed Row Stora-
para una matriz dispersa, formatos basados en estru turas CRS (
ge ) y CCS (Compressed Column Storage ).
6
La bibliote a está ompuesta por un onjunto de módulos. La siguiente gura muestra la organiza ión general de la bibliote a:
Figura 2: Bibliote a UCSparseLib y módulos
Fuente: Larrazábal (2004).
Según espe i a Larrazábal (2004), esta bibliote a está siendo utilizada omo
template
y
omo herramienta de desarrollo en varías apli a iones industriales, tal omo son los asos del
Computa ional S ien e Resear h Center
de San Diego University y del Departamento
de Modelaje y Simula ión de Ya imientos de INTEVEP S.A. Se puede resaltar que la bibliote a UCSparseLib está siendo usada de manera e iente tanto a nivel industrial omo
ientí o. Pero es de notar que en esta bibliote a ha e falta la existen ia de un módulo que permita resolver problemas de autovalores por diferentes té ni as, de modo tal de tener una bibliote a más robusta, que permita la apli a ión de la misma en un mayor número de apli a iones.
El problema prin ipal de quími a uánti a es en ontrar solu ión aproximada a la e ua ión de S hrödinger no relativista y esta ionaria, la ual viene dada por la siguiente rela ión:
b >= E|φ > H|φ En la e ua ión (3),
E
es la energía total del sistema en el estado
Hamiltoniano del sistema de nú leos y ele trones.
(3)
|φ >
y
b H
es el operador
En el Instituto Venezolano de Investiga iones Cientí as (IVIC), a tualmente se desa-
7
rrolla un software que permite resolver de manera aproximada la e ua ión (3) apli ando las e ua iones de Hartree-Fo k las uales a su vez son aproximadas por las e ua iones no lineales de Roothan en una forma auto onsistente por medio del pro eso iterativo deno-
Self-Consistent-Field ).
minado SCF (
Este pro edimiento SCF, onsiste en la siguiente se uen ia de pasos:
1. Espe i ar una molé ula a partir de un onjunto de parámetros base.
2. Cal ular todas las integrales mole ulares requeridas.
3. Diagonalizar la matriz de solapamiento para obtener matriz de transforma ión
4. Obtener una matriz
5. Cal ular la matriz
6. Sumar
G
P
G
X.
mediante suposi ión de un onjunto ini ial de oe ientes.
para la matriz
P
y las integrales de dos ele trones.
al ore-Hamiltoniano para obtener la matriz de Fo k.
7. Cal ular la matriz de Fo k transformada.
8. Diagonalizar
9. Cal ular
F′
para obtener
C′
y
E.
C = XC ′ .
10. Formar una nueva matriz
P
usando
C.
11. Determinar si el pro eso onverge.
En uanto a los Métodos Paramétri os, estos son muy e ientes en tiempo de ál ulo para evaluar propiedades ele tróni as de sistemas omplejos y amorfos, donde los métodos ab-initio son imprá ti os. Sin embargo, los fundamentos de los Métodos Paramétri os, hasta ha e po o, eran muy os uros. El futuro del desarrollo de los Métodos Paramétri os
8
es onvin ente, debido a la apa idad de reprodu ir, en forma ualitativa, la lo aliza ión de puntos ríti os del Lapla iano de la Densidad. Previas apli a iones de métodos semiempíri os, usando esquemas no estándar de parametriza ión, han resultado exitosos en el análisis de rea
iones atalíti as.
Por lo tanto, la apli a ión de los Métodos Paramétri os en la Catálisis Computa ional puede ser un fa tor de isivo en el desarrollo de la industria quími a en países del ter er mundo, debido al mayor poder de ompetitividad, es de ir:
Po o osto de inversión en investiga ión.
Ahorro de tiempo y dinero en desarrollos te nológi os.
Self-Consistent-Field ) se puede
En el pro eso SCF (
notar que en el paso 8 se ha formado
la forma matri ial de las e ua iones de Roothan, del ual se desea obtener la diagonaliza ión de la matriz de Fo k obtenida en el paso 7, de modo tal que el pro eso SCF
Self-Consistent-Field ) pueda
(
ser llevado a abo de forma efe tiva. Un pro eso de diago-
naliza ión de una matriz no es más que la obten ión de los autovalores de la matriz.
La importan ia de obtener on una buena pre isión los autovalores y respe tivos autove tores de la matriz transformada de Fo k, es que estos des riben las intera
iones entre las distintas partí ulas.
Lo anteriormente planteado, permite apli ar métodos numéri os para el ál ulo del espe tro de una matriz, tomando en uenta el desarrollo de arquite turas de avanzada que permiten la apli a ión de algoritmos paralelos para su eje u ión.
9
2 Formula ión del Problema ¾Se podrán obtener óptimos resultados en tiempo de eje u ión al al ular mediante Métodos Numéri os usando algoritmos paralelos los autovalores de la matriz de Fo k empleada en el pro eso SCF que aproxima la e ua ión de S hrödinger en problemas de quími a
uánti a?
3 Objetivos de la Investiga ión 3.1 Objetivo General Cal ular mediante métodos numéri os usando algoritmos paralelos los autovalores de la matriz de Fo k transformada que aproxima la e ua ión de S hrödinger en problemas de quími a uánti a.
3.2 Objetivos Espe í os 1. De nir las bases prin ipales de la té ni a de des omposi ión de autovalores.
2. Determinar las bases prin ipales del algoritmo de des omposi ión de autovalores.
3. Estable er el algoritmo omputarizado paralelo en la bibliote a UCSparseLib sobre la base del algoritmo estudiado.
4. Comparar los resultados numéri os que se obtienen al apli ar algoritmos paralelos que permitan obtener los autovalores de la matriz de Fo k la ual aproxima a la e ua ión de S hrödinger en problemas de quími a uánti a.
10
4 Justi a ión de la Investiga ión La importan ia de esta investiga ión es realizar un aporte al Departamento de Quími a Computa ional, donde se ha estado desarrollando un software que permita resolver problemas de quími a uánti a, en espe ial el problema proveniente de apli ar Métodos Paramétri os en la atálisis omputa ional.
Por ello el valor agregado de la misma será bajo las siguientes razones:
Teóri o: El desarrollo de algoritmos que vayan a la par on el desarrollo de las arquite turas de avanzada de las omputadoras, permitirá explotar métodos numéri os altamente es alables, omo se ha estado estudiando en el método de Ja obi por Demmel, Heath y van der Vorst (1993).
Prá ti o: El desarrollo de algoritmos e ientes permitirá que la bibliote a UCSparseLib se vea fortale ida mediante la in orpora ión de un módulo que permita el ál ulo de autovalores por diversos métodos.
Metodológi o: Obtener la diagonaliza ión de la matriz de Fo k, o lo que es igual, los autovalores de di ha matriz, permitirá que los estudios que se realizan en el IVIC tomen mayor importan ia, ya que al obtener resultados ables en una apli a ión numéri a y on un buen tiempo de eje u ión les permitirá realizar pruebas de mayor envergadura.
So ial: Esta investiga ión permitirá ampliar el espe tro de ono imientos aso iados
on tópi os avanzados del Álgebra Lineal, lo ual aportará a los estudiantes de la materia nuevas té ni as teóri as dentro de su pro eso de aprendizaje.
De forma que la presente investiga ión estará entrada en la re ole
ión de informa ión, estudio y análisis de té ni as que permitan obtener todos los autovalores de una matriz
11
simétri a que posea estru tura dispersa o densa. Tomando en uenta esta onsidera ión resulta interesante explotar la paraleliza ión de estas té ni as.
CAPÍTULO II: MARCO TEÓRICO
13
CAPÍTULO II MARCO TEÓRICO 1 Ante edentes El ál ulo de autovalores es un problema de onsiderable interés teóri o y tiene un an ho rango de apli a iones. Por ejemplo, este problema es ru ial en la resolu ión de sistemas de e ua iones diferen iales, análisis de modelos de re imiento de pobla ión, ál ulo de poten ias de matri es, entre otras. Áreas tales omo físi a, so iología, quími a, biología, e onomía y estadísti a, han enfo ado onsiderablemente su aten ión sobre el uso y ál ulo de los autovalores y autove tores.
La mayoría de los métodos numéri os empleados para obtener todos los autovalores y posibles autove tores de una matriz, están basados bajo la siguiente de ni ión, omo lo estable e Strang (1982) en su texto.
De ni ión 1: Se di e que dos matri es una matriz no singular es similar a similitud.
A
B
S
de tamaño
y por lo tanto
B
y
n×n
es similar a
B
es unitariamente similar a
Proposi ión 1: Supóngase que
A
A
y
B
si
B
de tamaño
tal que
A; S
n×n
A = S −1 BS .
mientras que
S
son similares, si existe Se di e enton es que
A
es una transforma ión de
es unitaria.
son matri es similares, es de ir,
A = S −1 BS ,
en-
ton es las siguientes propiedades se satisfa en: 1.
A
2.
Ax = λx ⇒ B(P −1x) = λ(P −1 x).
3.
Bw = λw ⇒ B(P w) = λ(P w).
y
B
poseen los mismos autovalores. De he ho
PA (λ) ≡ PB (λ).
En fun ión a la Proposi ión 1 y la De ni ión 1, el Teorema 1 estable e que los autovalores
14
de dos matri es son iguales si las matri es son similares.
Teorema 1: Supongamos que
A
on autove tor
enton es
Sx
x
A
y
B
son matri es similares y que
aso iado. Enton es
λ
es un autove tor aso iado a
Demostra ión: Supongamos que
x 6= 0
λ
es también un autovalor de
λ
para la matriz
es un autovalor de
B,
y si
A = S −1 BS
B.
es tal que
S −1 BSx = Ax = λx Al multipli ar en la izquierda por la matriz
S,
obtenemos:
BSx = λSx Puesto que
x 6= 0 y S
al autovalor
es no singular
Sx 6= 0. Por tanto, Sx es un autove tor de B
aso iado
λ. 2
Como señalan Golub y van der Vorst (2000) el ál ulo numéri o de autovalores y autove tores es deli ado, en parti ular uando los autove tores ha en pequeños ángulos on otros autove tores. En asos extremos, uando la matriz es defe tuosa,
A
puede ser redu ida
a la forma anóni a de Jordan, pero pequeñas perturba iones arbitrarias en
A
pueden
ondu ir a una matriz no defe tuosa. Esto nos lleva a la pregunta prin ipal: ¾Cómo se puede al ular los autovalores y autove tores de manera e iente y uan exa tos son?.
15
2 Bases Teóri as 2.1 Forma Diagonal de una Matriz La aproxima ión estándar para el ál ulo numéri o del problema de autovalores es redu ir los operadores involu rados en una forma más simple, que onlleve a los autovalores y a los autove tores de manera dire ta, por ejemplo, una forma diagonal. El siguiente teorema, el ual en ontramos en Strang (1982) muestra omo la forma diagonal de una matriz nos da los autovalores de la misma.
Teorema 2: Supóngase que una matriz
A
de
n×n
tiene
n
autove tores linealmente inde-
S,
pendientes. Enton es si se eligen estos ve tores omo las olumnas de una matriz sigue que
S −1 AS
es una matriz diagonal
Λ,
on los autovalores de
λ1
S −1 AS = Λ =
Demostra ión: De niendo al ve tor olumna autove tores formarán las olumnas de
S
A
en la diagonal.
λ2 ..
.
λn
xi
omo autove tor de la matriz
y olo ando el produ to
AS
.. .
.. .
.. .
. ..
. ..
. ..
.. .
.. .
.. .
. ..
. ..
. ..
A.
Tales
por olumna queda
que:
se
x1 x2 · · · xn λ1 x1 λ2 x2 · · · λn xn AS = A . = . . . . . . . . . . . . . . . ..
16
La matriz del lado izquierdo se puede des omponer omo el siguiente produ to:
. . .
. . .
. . .
λ1 x1 λ2 x2 · · · λn xn . . . . .. .. . . ..
. ..
. ..
. . .
. . .
. . .
λ1 x1 x2 · · · xn λ2 = . . . . .. .. . . ..
. ..
. ..
..
.
λn
2
De esta manera se tiene el siguiente resultado, el ual es válido dado que se supuso que los autove tores que forman las olumnas de
S
S
son linealmente independientes, por lo tanto
es una matriz no singular.
AS = SΛ ⇔ S −1 AS = Λ ⇔ A = SΛS −1 Las matri es simétri as
(A = At )
surgen de manera natural en la investiga ión, y omo
ara terísti a espe ial se tiene que toda matriz simétri a es diagonalizable. Los siguientes teoremas estable en tal vera idad.
Teorema 3: Todas la raí es del polinomio ara terísti o de una matriz simétri a son números reales.
Demostra ión: En prin ipio se tiene que para toda matriz simétri a y ualquier es alar
λ
se umple la siguiente rela ión:
(λI − A)t = (λI)t − At = λI − A Sea ahora
A
una matriz simétri a y supóngase que existe un número omplejo
que es raíz de la e ua ión ara terísti a de
A.
Si
b = 0,
enton es
λ
λ = a + bi,
es un número real.
17
Se tiene que
det ((a + bi) I − A) = 0
((a + bi) I − A) S
y enton es
es singular para toda matriz
S.
(a + bi) I − A
es singular. Por lo tanto
En parti ular:
((a + bi) I − A) ((a − bi) I − A) = a2 + b2 I − 2aA + A2 + b2 I = (aI − A)2 + b2 I es singular. Puesto que
(aI − A)2 + b2 I
es singular, existe un ve tor
v
no nulo tal que:
(aI − A)2 + b2 I v = 0 = v t (aI − A)2 + b2 I v enton es se tiene que
0 = v t (aI − A) (aI − A) + b2 I v = v t (aI − A) (aI − A) v + b2 v t Iv = v t (aI − A)t (aI − A) v + b2 |v|2 = ((aI − A) v)t (aI − A) v + b2 |v|2 = |(aI − A) v|2 + b2 |v|2 En esta rela ión se puede notar que
b = 0.
λ
Por lo tanto,
|(aI − A) v|2 ≥ 0
y dado que
v 6= 0,
es una raíz real de la e ua ión ara terísti a de
nos queda que
A. 2
Teorema 4: Sean
λ1
y
λ2
autovalores distintos de una matriz simétri a
autove tores que orresponden a
λ1
y
λ2 ,
respe tivamente, enton es
A. Si v1 v2
y
v2
son los
v1
y
A
es una matriz de
son ve tores
ortogonales.
Demostra ión: Si tamaño
n × n,
u
y
v
enton es
son ve tores olumnas de tamaño
u · v = ut v .
n×1
y
Por lo tanto:
Au · v = (Au)t v = ut At v = ut At v = u · At v
18
y enton es
(Au) · v = u · At v Sean ahora
v1
y
autovalores, on
v2
autove tores de la matriz simétri a
λ1 6= λ2 .
(4)
A,
y
λ1
Por la e ua ión (4) y la simetría de
y
λ2
sus orrespondientes
A,
λ1 (v1 · v2 ) = (λ1 v1 ) · v2 = (Av1 ) · v2 = v1 · At v2
= v1 · (Av2 ) = v1 · (λ2 v2 ) = λ2 (v1 · v2 ) Por lo tanto, que
(λ1 − λ2 ) (v1 · v2 ) = 0. Puesto que λ1 6= λ2 , enton es λ1 − λ2 6= 0. Se on luye
v1 · v2 = 0,
y de esta manera
v1
y
v2
son ortogonales.
2 En resumidas uentas toda matriz simétri a es siempre diagonalizable por medio de matri es ortogonales, se di e enton es que
A
es ortogonalmente diagonalizable. El siguiente
pro edimiento des rito por Kolman (1999) permite diagonalizar una matriz simétri a mediante una matriz ortogonal
P.
1. Formar el polinomio ara terísti o
f (λ) = det (λIn − A).
2. Determinar las raí es del polinomio ara terísti o de
3. Para ada valor
λj
de
A.
Las uales son reales.
A de multipli idad kj , se determina una base de kj
para el espa io solu ión de
(λj In − A) x = 0
(el espa io propio de
autove tores
λj ).
4. Para ada espa io propio, se transforma la base obtenida en el paso 3 en una base ortonormal mediante el pro eso de Gram-S hmidt. La totalidad de estas bases ortonormales determina un onjunto ortonormal de independientes.
n
autove tores de
A
linealmente
19
5. Sea
P
la matriz uyas olumnas son los
hallados en el paso 4. Enton es
D,
P
n
autove tores linealmente independientes
es una matriz ortogonal y
P −1 AP = P t AP =
una matriz diagonal, uyos elementos de la diagonal son los autovalores de
orrespondientes a las olumnas de
A
P.
2.2 Método de Ja obi El método de Ja obi fue ini ialmente propuesto en 1846, este mismo redu e una matriz simétri a real a una matriz diagonal mediante una se uen ia de rota iones planares. Para matri es hasta de orden 10, el método es ompetitivo on aquellos más so sti ados. Si la velo idad no es prioridad, es también a eptable para matri es de orden 20, omo lo señalan Mathews y Fink (1999).
Dado que el método de Ja obi obtiene una matriz diagonal la ual es similar a la matriz
A,
enton es, por las de ni iones dadas anteriormente, los autovalores de la matriz
diagonal oin iden on los autovalores de la matriz
A.
Una rota ión en el plano es una transforma ión lineal diseñada para an elar los elementos fuera de la diagonal prin ipal, sea transforma ión lineal
y = Rx,
donde
R
x
un ve tor de dimensión
es una matriz
n × n,
dada por:
n
y onsidérese la
20
1 ··· 0 . .. 0 · · · cos(φ) .. . R= 0 · · · − sin(φ) . .. 0 ··· 0
···
0
· · · sin(φ) · · · 0 . .. · · · cos(φ) · · · 0 .. . ··· 0 ··· 1 . . .
El efe to de la transforma ión
p
y = Rx
yj = xj
← la p ← la q
↑
↑
ol
··· 0
ol
q
es la siguiente:
uando
j 6= p
y
j 6= q
yp = xp cos(φ) + xq sin(φ) yq = −xp sin(φ) + xq cos(φ) la transforma ión es vista omo una rota ión del ángulo
φ.
n−dimensional
en el plano
por medio del ángulo
−φ.
x = R−1 y
en la
realiza la misma rota ión en el plano
xp xq
Se puede observar que
R−1 = Rt
o
R
yp = 0
por medio
yq = 0
Sele
ionando el ángulo apropiado, se puede ha er que
imagen. La transforma ión inversa
xp xq y
es una matriz ortogonal, es de ir:
Rt R = I
Las siguientes rela iones estable en el he ho de que dada una matriz simétri a, la apli a ión de matri es ortogonales produ e una matriz simétri a similar a la matriz original. Consideremos el problema (1) y supóngase que de nida por
K
es una matriz no singular y que
B
está
21
B = K −1 AK
(5)
multipli ando ambos miembros por la dere ha por el término
K −1 x,
produ e el siguiente
resultado
BK −1 x = K −1 AKK −1 x = K −1 Ax = K −1 λx = λK −1 x
(6)
Si se de ne el siguiente ambio de variable
y = K −1 x
o
Ky = x
(7)
empleando la e ua ión (7) en la e ua ión (6), el problema ini ial pasa a tomar la siguiente forma.
By = λy
(8)
Comparando la e ua ión (1) y la e ua ión (8) se observa que la transforma ión de similitud preserva los autovalores
Supóngase enton es que
R
λ.
es una matriz ortogonal y que
D
está de nida por
D = Rt AR multipli ando ambos miembros por el lado dere ho por el término
(9)
Rt x, produ e el siguiente
resultado
DRt x = Rt ARRt x = Rt Ax = Rt λx = λRt x Si se de ne el siguiente ambio de variable
(10)
22
y = Rt x
o
Ry = x
(11)
empleando la e ua ión (11) en la e ua ión (10), el problema ini ial pasa a ser:
Dy = λy
(12)
Como se estable ió anteriormente, los autovalores de la e ua ión (1) y los de la e ua ión (12) son iguales. Sin embargo, en la e ua ión (11) es mu ho más sen illo realizar la onversión de
x
a
y
y regresar
y
a
x,
debido a que
R−1 = Rt .
Si además, suponemos de entrada que la matriz
D t = Rt AR es de ir,
D
t
A
= Rt A Rt
es una matriz simétri a. Enton es si
matriz ortogonal, enton es la transforma ión de
es simétri a, enton es tenemos que:
t
= Rt AR = D
A
A
a
es una matriz simétri a y
D
R
es una
dada por la e ua ión (9) preserva
la simetría de la matriz así omo sus autovalores.
El he ho de que el método de Ja obi fuera exitoso en el pro eso de diagonaliza ión de una matriz simétri a mediante transforma iones de similitud ortogonales inspiró a varios investigadores a hallar métodos similares para matri es no simétri as. Lo ual fue rápidamente realizado mediante la forma de S hur, según es señalado en Golub y van der Vorst (2000).
Como estable e Giménez, Hernández y Vidal (1995) uando se apli a el algoritmo se uen ial del método de Ja obi a matri es simétri as se puede explotar la simetría a tualizando solamente la parte triangular superior de la matriz. Sin embargo el mismo autor di e que no es tan fá il explotar la simetría uando se implementa Ja obi en forma paralela.
23
2.3 Itera ión QR y Autovalores El método
QR
fue ini iado por Kublanovskaya y Fran is (1961), alrededor de 1961. Este
es un método iterativo on ebido para al ular los autovalores de
A mediante transforma-
iones de similitud. Como su nombre lo indi a, este algoritmo ha e uso de la fa toriza ión
QR, de modo tal que la matriz A se puede des omponer omo el produ to de dos matri es, una matriz
Q
ortogonal y una matriz
R
triangular superior, es de ir:
A = QR
(13)
Tal fa toriza ión se lleva a abo de diversas formas omo es planteado por Golub y Van Loan (1996), las formas más utilizadas es a través de las rota iones de Givens y las transforma iones de Householder.
Las rota iones de Givens tienen omo objetivo, al igual que las rota iones de Ja obi, eliminar un elemento de la matriz. De este modo si denotamos
Gj
la
j−ésima
rota ión de
Givens, enton es
Qt A = R es una matriz triangular superior donde
Q = G1 · · · Gk
y
k
es el total de rota iones que
hay que realizar para hallar tal fa toriza ión. Si la matriz es simétri a, enton es al apli ar las transforma iones de Givens por olumnas pueden ser apli adas también por las y de este modo llegar a una matriz tridiagonal, la ual es similar a la matriz
A,
en aso de
no ser simétri a, al apli ar las matri es de Givens por olumna se obtiene una matriz de Hessenberg.
24
Las transforma iones de Householder umplen on el mismo objetivo tanto para las matri es simétri as omo para las matri es no simétri as, on la ventaja on respe to a Givens que elimina una olumna en una sola itera ión. El siguiente teorema dado en Burden y Faires (1998) de ne laramente omo se obtienen las transforma iones de Householder.
De ni ión 2: Sea
w ∈ Rn
on
w t w = 1.
La matriz de
n×n
P = I − 2ww t
(14)
re ibe el nombre de transforma ión de Householder.
Teorema 5: Si
P = I − 2ww t
es una transforma ión de Householder, enton es
métri a y ortogonal; por lo tanto,
P −1 = P .
Demostra ión: Dado que
ww t
t
t = w t w t = ww t
se dedu e que
P t = I − 2ww t Más aún, dado que
t
= I − 2ww t = P
w t w = 1,
P P t = I − 2ww t
I − 2ww t = I − 2ww t − 2ww t + 4ww tww t = I − 4ww t + 4ww t = I
Así que:
P
es si-
25
P −1 = P t = P 2 Al apli ar el algoritmo
QR
de forma iterativa, la matriz
A
se va transformado en una
matriz diagonal similar, esto viene justi ado por las siguientes rela iones.
Sea la e ua ión dada por la fa toriza ión (13), multipli ando di ha e ua ión por
Qt ,
por
la izquierda, queda que:
Qt A = Qt QR
omo
Q
es una matriz ortogonal, se tiene que
Q = Q−1 = Qt ,
(15)
por lo tanto, la e ua ión
(15) se transforma en
Qt A = R
(16)
B = RQ
(17)
si se de ne
y empleando la e ua ión (16), se tiene que
B = Qt AQ por lo que
B
es semejante a
A,
y por lo tanto se preservan los autovalores.
Si este pro eso se realiza de manera iterativa, la su esión de matri es onverge a una matriz diagonal la ual es semejante a es llevado a abo.
A.
En la gura 3 se puede ver omo este pro eso
26
Figura 3: Algoritmo Entrada:
An×n
Salida: Matriz
real y
An×n
M
QR.
número máximo de itera iones.
tal que
diag(A) = λi , i = 1, . . . , n.
Ini io
A1 = A Paso 2: Para k = 1 : M ha er Paso 3: Ak = Qk Rk Paso 4: Ak+1 = Rk QK Paso 1:
pasos 3 y 4
n-para Fin
2.3.1 Tridiagonaliza ión de Matri es El método de Ja obi elimina solo un par de elementos a la vez en ada paso, las transforma iones de Householder elimina los elementos de una la on ex ep ión de los primeros elementos. Este método fué estudiado por Wilkinson, Reins h y Martin (1968), y en el mismo do umento se presentan distintas versiones del algoritmo de tridiagonaliza ión, las
uales se diferen ian bási amente, en ual es el método es usado para diagonalizar la matriz tridiagonal resultante.
Tridiagonalizar la matriz y luego diagonalizarla por medio del método
QR
o
QL
estu-
diado por Wilkinson, Reins h y Bowdler (1968), resulta ser un método más e iente para resolver el problema de autovalores, según es señalado por Dubrulle (1972), al igual omo es planteado en un trabajo previo realizado por Wilkinson (1960).
El método de Householder tiene omo objetivo redu ir una matriz simétri a
A
n × n,
transforma-
a una matriz del mismo orden on forma tridiagonal, apli ando
n−2
de orden
iones ortogonales. Cada transforma ión elimina la parte requerida de una olumna y su
orrespondiente la. El ingrediente prin ipal del método son las ono idas matri es
P
de
Householder, las uales poseen la forma dada en la e ua ión (14). Rees ribiendo la matriz
27
P
de la siguiente manera:
P =I− donde el es alar
H
viene dado por la siguiente rela ión
H≡ Sea
x
u · ut H
1 2 |u| 2
el ve tor onformado por la primera olumna de la matriz
ve tor
u
A,
tomando ento es el
mediante la siguiente rela ión
u = x ∓ |x| e1 donde
e1
es el ve tor anóni o
forma ión
P
al ve tor
x,
[1, 0, . . . , 0]t .
De esta manera, apli ando la matriz de trans-
queda que:
u · (x ∓ |x|e1 )t · x H 2u · (|x|2 ∓ |x|x1 ) = x− 2|x|2 ∓ 2|x|x1
P ·x = x−
= x−u
= ∓|x|e1 Lo ual omprueba que al apli ar la matriz
P
de Householder sobre un ve tor
x,
ésta
ha e ero los elementos del ve tor a ex ep ión del primero. Bajo esta idea, para redu ir la matriz
A a una
forma tridiagonal, se sele
iona omo ve tor
de la primera olumna, de modo tal que los últimos
n−2
x los últimos n − 1
elementos
elementos serán he hos ero.
28
... 0 a1,1 a1,2 a1,3 . . . a1,n 1 0 0 0 a2,1 PA = 0 a3,1 . . (n−1) .. .. P 1 an,1 0 a1,1 a1,2 a1,3 . . . a1,n k = 0 . .. 0 Donde
(n−1)
P1
denota la matriz de Householder de dimensión
k
viene dado más o menos la magnitud del ve tor
el valor
(n−1)×(n−1), mientras que [a2,1 , . . . , an,1]t .
Para nalizar
la transforma ión se multipli a la matriz de transforma ión por el lado dere ho, sabiendo que
Pt = P,
y queda que:
a1,1 k 0 . . . 0 k ′ A = P AP = 0 . .. 0 Sele
ionando omo ve tor
onstruye la matriz
P2 ,
x
los
n−2
últimos elementos de la segunda olumna y se
de la siguiente forma:
29
··· 0 1 0 0 0 1 0 ··· 0 P2 ≡ 0 0 . . (n−2) .. .. P 2 0 0 El bloque identidad en la parte superior izquierda asegura que la tridiagonaliza ión obtenida en los pasos previos no sea destruida en los pasos posteriores. Claramente, una se uen ia de
n−2
transforma iones, redu irán la matriz
Para no llevar a abo el produ to matri ial
p≡
P · A · P,
A
a su forma tridiagonal.
se al ula el ve tor
A·u H
enton es
u · ut = A − p · ut A·P =A· 1− H A′ = P · A · P = A − p · ut − u · pt + 2Ku · ut donde el es alar
K
está de nido por la siguiente rela ión
K=
ut · p 2H
Si se es ribe que
q ≡ p − Ku se tiene enton es que
30
A′ = A − q · ut − u · q t Estos resultados son mostrados en Wilkinson, Reins h y Martin (1968), en el mismo la rutina de redu
ión de Householder omienza en la
n-esima
olumna de
A
y no en la
primera omo se expli ó anteriormente. De forma más detallada, las e ua iones son las siguientes: En la etapa
m (m = 1, 2, . . . , n − 2)
el ve tor
ut = (ai,1 , ai,2 , . . . , ai,i−2 , ai,i−1 ±
√
u
posee la siguiente estru tura:
σ, 0, . . . , 0)
(18)
donde
i ≡ n − m + 1 = n, n − 1, . . . , 3 De esta forma el valor de
σ
viene dado por la siguiente rela ión
σ = |x|2 = (ai,1 )2 + · · · + (ai,i−1 )2 mientras que el signo de
ai,i−1 ,
σ
en la e ua ión (18) viene determinado por el signo del elemento
de esta manera se bus a minimizar el error de redondeo, lo ual onlleva a obtener
una buena aproxima ión, omo es planteado por Wilkinson (1960).
Para que el algoritmo
QR sea e az, la matriz simétri a debe poseer una forma tridiagonal,
omo lo plantea Burden (1998), es por ello que este algoritmo es un pro eso de dos fases. A la matriz
A
se le apli a una fase 1 que onsiste en la tridiagonaliza ión de la matriz. Y
luego se apli a una fase 2 que es diagonalizar la matriz obtenida en la fase 1, por medio de la fa toriza ión
QR, tal omo se presenta
De esta manera, dado que la matriz
T
en la gura 4, dada en Trefethen y Baun (1997).
es también simétri a, lo ual se asegura debido a
31
Figura 4: Tridiagonaliza ión y Diagonaliza ión de una Matriz
x x x x x
x x x x x
x x x x x
x x x x x
x x x x x
A = At
Fase 1
−→
x x x x x x x x x x x x x
T
Fase 2
−→
x x x x x
D
Fuente: Trefethen y Baun (1997).
las transforma iones de similitud realizadas, enton es apli ar el método
QR
a esta matriz
es más e onómi o, ya que solo se alma ena la diagonal prin ipal y la subdiagonal. Y diagonalizar la matriz, onsiste enton es en anular los de ir, la fa toriza ión
QR de
n−1
elementos de la subdiagonal, es
una matriz tridiagonal onsiste en
n−1
transforma iones de
Givens. En la gura 5 se muestra el algoritmo que permite obtener una matriz tridiagonal a travéz del pro eso anterior.
Figura 5: Algoritmo Se uen ial de Tridiagonaliza ión. Entrada:
An×n
real y simétri a.
Salida: Matri es
Q
Q
y
T
tales que
T = Qt AQ,
donde
T
es tridiagonal y
es ortogonal.
Ini io
Q=I k = 1 : n − 2 ha er pasos 3 al 9 t Paso 3: u = A(k, k + 1 : n) Paso 4: v = u + ||u||2 signo(u1 )e1 t Paso 5: p = 2A(k + 1 : n, k + 1 : n)v/v v t t Paso 6: w = p − (p v)v/v v t t Paso 7: A(k + 1 : n, k + 1 : n) = A(k + 1 : n, k + 1 : n) − vw − wv Paso 8: y = −2Qv t Paso 9: Q = Q + yv
Paso 1:
Paso 2: Para
n-para Fin
La matriz
Z = Pn−2 · · · P1
de transforma ión es ne esaria para el ál ulo de autove tores
32
empleando el algoritmo
QR
al ser apli ado sobre una matriz simétri a que es tridiagona-
lizada por las transforma iones de Householder, tal omo es planteado por Reeve y Heath (1999).
2.3.2 Algoritmos QR y QL Como se mostró anteriormente, la idea bási a detras del algoritmo
QR
es que ualquier
matriz real puede ser des ompuesta mediante el produ to dado en la e ua ión (13), donde
Q
es una matriz ortogonal y
R
es triangular superior. En general esta des omposi ión es
llevada a abo apli ando transforma iones de Householder sobre las olumnas de la matriz
A.
Considerando la matriz generada por la e ua ión (17) se on luye que ambas matri es
son similares y por lo tanto se preservan las siguientes propiedades de la matriz original: simetria, forma tridiagonal y forma de Hessenberg.
Por la rela iones dadas anteriormente, no existe una rela ión espe ial a er a de la ele
ión de que uno de los fa tores de
A sea
triangular superior, se pudiese tomar este fa tor omo
triangular inferior. Realizando este ambio se da paso al algoritmo
QL
dado que:
A = QL donde
L es
triangular inferior y
Q sigue
(19)
siendo ortogonal. En Press, Vetterling, Teukolsky
y Flannery (2002) se estable e que el error de redondeo es menor que en el mismo algoritmo
QR,
y por lo tanto es una té ni a re omendable para el ál ulo de autovalores.
Press et al (2002) estable e que el algoritmo
QL
posee un orden de eje u ión
O(n3 )
por
itera ión para una matriz genéri a, lo ual es pra ti amente prohibitivo. Sin embargo, el orden de eje u ión es nal y
O(n2)
O(n)
por itera ión uando la matriz posee una estru tura tridiago-
uando la matriz posee una forma de Hessenberg, lo que lo onvierte en un
33
algoritmo altamente e iente sobre estas formas.
2.3.3 GIVENS y la Fa toriza ión QR Suponiendo a partir de ahora, de que se posee una matriz tridiagonal
ión
QR
A,
la fa toriza-
es apli ada a una matriz on la siguiente forma:
a b 0 ··· 0 1 2 . .. . b2 a2 b3 . . .. A = 0 b3 a3 . 0 .. . . . . . . . . . . bn 0 · · · 0 bn an Si
b2 = 0
A.
Cuando
o
bn = 0, bj = 0,
enton es
donde
a1
an
o bien
2 < j < n,
produ en automáti amente un autovalor de
el problema puede ser minimizado onsiderando
matri es más pequeñas, omo por ejemplo:
a b 0 ··· 0 1 2 . . b2 a2 b3 . . . . .. 0 b . a3 0 3 .. . . . . .. . . . . bj−1 0 · · · 0 bj−1 aj−1 Si ninguna bj es ero, el método
y
a bj+1 0 ··· 0 j . .. . bj+1 aj+1 bj+2 . . .. 0 b . 0 j+2 aj+2 .. .. .. .. . . . . bn 0 ··· 0 bn an
QR forma una su esión de matri es A = A(1) , A(2) , A(3) , . . .,
así:
1.
2.
A(1) = A
se fa toriza omo un produ to
R(1)
es triangular superior.
A(2)
se de ne omo
A(2) = R(1) Q(1) .
A(1) = Q(1) R(1) ,
donde
Q(1)
es ortogonal y
34
De esta forma, siguiendo los pasos mostrados en la gura 3, se obtiene sobre la diagonal de
A(M )
los autovalores de la matriz
A,
eliminando los elementos de la subdiagonal.
La idea es obtener esta matriz diagonal, en un número de pasos menor a los ne esarios si se empleara el algoritmo
QR,
omo es dado en la gura 3.
2.3.4 A elera ión de la onvergen ia del Método QR Para que la diagonaliza ión de la matriz tridiagonal obtenida se reali e de una forma más a elerada, se emplean té ni as ono idas omo desplazamiento de la matriz de oe ientes. El algoritmo presentado en la gura 3 no presenta tales desplazamientos.
La idea de la té ni a de desplazamiento onsiste en apli ar la fa toriza ión
T − σI ,
en vez de apli arla solo a la matriz
matriz identidad y
σ
T,
donde
T
es la matriz tridiagonal e
es el valor que ambia el espe tro entero de la matriz a
La idea prin ipal de esta té ni a de desplazamiento es la siguiente, si original y
Ai
QR a la matriz
es la matriz a tual en la itera ión
hallar la fa toriza ión
QR
i,
de la matriz desplazada
miento realizado para de nir la matriz que es muy er ano a un autovalor de
Ai+1 .
A,
A
se sele
iona un es alar
Ai − σI
Si el es alar
σ
I
es la
λi − σ . es la matriz
σ,
para luego
y luego devolver el desplazaes sele
ionado de modo tal
esta té ni a a elerará la rapidez de onvergen ia
del método.
(Ai − σI) = Qi Ri
(20)
y luego se reversa el desplazamiento realizado, y se tiene que:
Ai+1 = Ri Qi + σI
(21)
35
Para omprobar que el desplazamiento realizado en (20) y (21) preservan la similitud, multipli ando la e ua ión (21) por
Q
por el lado izquierdo, tenemos que:
Ai+1 = Ri Qi + σI ⇒ Qi Ai+1 = Qi Ri Qi + σQi I tenemos enton es que
Qi Ri = Ai − σI , in orporando esta
e ua ión en (22) nos queda que:
Qi Ai+1 = (Ai − σI)Qi + σQi I ⇒ Qi Ai+1 = Ai Qi − σQi + σQi multipli ando ahora por
Qt
(22)
(23)
por el lado dere ho llegamos a que:
Qi Ai+1 Qti = Ai Qi Qti − σQi Qti + σQi Qti
an elando los dos términos similares de la dere ha, y dado que
(24)
Qt Q = QQt = I
se llega
a que:
Qi Ai+1 Qti = Ai
(25)
es más, si la e ua ión (25) se quiere mostrar omo una transforma ión de similitud, se multipli a por
Q por el lado dere ho y por Qt
por el lado izquierdo, on lo ual se on luye
que:
Qti Qi Ai+1 Qti Qi = Qti Ai Qi ⇒ Ai+1 = Qti Ai Qi
De modo tal que
Ai+1
es similar a la matriz
(26)
Ai , por lo tanto en ada itera ión se preservan
los autovalores de la matriz original.
Existen diversas estrategias de desplazamientos, entre las uales tenemos la de desplazamiento simple y desplazamiento implí ito, que pueden ser apli ados al método
QR.
36
2.3.4.1 Desplazamiento Simple Como lo plantea Meyer (2000), una forma es tomar en ada paso del método
σk ,
QR
al valor
omo una aproxima ión real de un autovalor. Por lo tanto, esta té ni a está basada
en el he ho de que un desplazamiento ausado por un autovalor provo a una de a ión,
σk
omo un autovalor aproximado.
A(k)
omo valor de desplazamiento.
por lo tanto es sugerido sele
ionar el desplazamiento Se puede sele
ionar el elemento
(k, k)
de la matriz
En Tisseur (1996) y Golub y Van Loan (1996), se ha e referen ia a que esta estrategia produ e una onvergen ia uadráti a en el aso no simétri o y onvergen ia úbi a en el
aso simétri o.
Sin embargo, en el aso simétri o, Wilkinson propor ionó razones heurísti as por las uales se puede sele
ionar un desplazamiento alternativo, ono ido omo
Wilkinson-Shift,
que
puede ser una mejor op ión. La idea de Wilkinson es sele
ionar omo valor de desplazamiento al autovalor de la matriz de orden
2×2
ak−1,k−1 ak−1,k A(k) = ak,k−1 ak,k que sea más er ano al valor
(k)
Ak,k ,
si ambos autovalores están igual de er a, enton es uno
de los autovalores es sele
ionado de forma arbitraria. Una fórmula numéri amente estable para obtener el desplazamiento de Wilkinson está de nida por
a2k−1,k q σ = ak,k − sign(δ) |δ| + δ 2 + a2k,k−1
37
donde
δ = (ak−1,k−1 − ak,k )/2.
Si
δ = 0,
enton es
sign(δ)
es olo ado arbitrariamente
a 1 o -1.
En Tisseur (1996) se plantea que si los valores
σk
son jados, es de ir,
σ ≡ σk
por ada
itera ión y si los autovalores son ordenados tales que
|λ1 − σ| > |λ2 − σ| > · · · > |λn − σ| enton es la
i−esima
entrada subdiagonal en la se uen ia de matri es
A(k)
onverge a ero
linealmente on una rata de
|λi+1 − σ| |λi − σ|
De esta manera, el algoritmo presentado en la gura 3, es modi ado, introdu iendo la idea del desplazamiento
σk
por alguna de las té ni as men ionadas anteriormente.
QR on Desplazamiento A0 = A Para i = 1, . . . , hasta onvergen ia determinar un es alar σi fa torizar Ai−1 − σi I = Qi Ri
al ular Ai = Ri Qi + σi I
Figura 6: Algoritmo Ini io on
Fin Para Fin
2.3.4.2 Desplazamiento Implí ito El algoritmo
QR
on desplazamiento implí ito, es una té ni a similar al desplazamien-
to de Wilkinson, se bus a el autovalor de la matriz
2×2
más er ano al elemento de la
diagonal, esta té ni a fue presentada en Wilkinson, Reins h y Bowdler (1968), y la diferen ia radi a en llevar a abo la transforma ión
QR
sin substraer el desplazamiento de los
38
elementos de la diagonal, omo es llevado a abo en los desplazamientos anteriores.
La lave prin ipal del desplazamiento implí ito es produ ir
A,
nes unitarias dire tamente a la matriz
Ak+1
apli ando transforma io-
en vez de apli arlas a las matri es desplazadas
omo en la e ua ión (21). Miminis y Paige (1991) estable iron teorías más sólidas para el uso del desplazamiento implí ito en el algoritmo
El algoritmo
QR
o en sus variantes.
QR on desplazamiento implí ito para el ál ulo de autovalores es usualmen-
te utilizado asumiendo que ninguno de los desplazamientos son autovalores. Sin embargo, re ientemente, una lase de algoritmos basados en algunas de las teorís del algopritmo
QR
han propor ionado algoritmos dire tos (más que iterativos) para la obten ión de autovalores, tal omo lo plantea Miminis y Paige (1991), en este mismo do umento se presenta el teorema prin ipal del desplazamiento implí ito.
Teorema 6: Sea triz simétri a
Q ∈ Rn×n
y
U ∈ Rn×n
matri es ortogonales tales que, dada una ma-
M ∈ Rn×n T = Qt MQ S = U t MU
son ambas matri es tridiagonales. Si redu ida,
Qe1 = ±Ue1 ,
y
Qe1 = Ue1
|eti T ei−1 | = |eti Sei−1 |
y
para
T
es no es redu ida enton es
S
es no
2 ≤ i ≤ n.
2.3.5 Criterio de Parada El algoritmo triz
A
QR
de orden
presenta un pro eso iterativo para hallar los autovalores de una ma-
n × n,
pero en el pro eso se debe indi ar ual es el riterio a seguir, para
indi ar que se ha hallado un autovalor.
39
En Ahues y Tisseur (1997) se estable en riterios para indi ar uando se ha obtenido un autovalor, por ejemplo, sea pequeño, er a de ero. Si
i
el entero más grande tal que
ai,i−1
es lo su ientemente
i = n se ha hallado un autovalor. Si i = n − 1 se han hallado dos
autovalores. De esta manera el riterio estable ería que un autovalor es obtenido uando
|Aki,i−1 | ≤ u k A kF donde
(27)
u representa la unidad de redondeo de la máquina y k · kF
representa la norma de
Frobenius.
Con el riterio dado en la e ua ión (27), en Ahues y Tisseur (1997) se men iona que el algoritmo
QR
es estable, en sentido de lo estable ido por Tisseur (1996). Sin embargo,
este riterio no es lo su ientemente satisfa torio para ierto tipos de matri es, omo es mostrado por Ahues y Tisseur (1997), mediante la siguiente matriz.
1 10−2 A= −4 10 10−6
10
−2
10
−4
10−6
10
−4
10
−6
10
−6
10
−8
10−8
10−10
10−8 10−10 10−12
En donde algunos autovalores pueden ser menores que
onvergiendo a un autovalor que es menor que
olo ar a
(k)
an,n−1
u k A kF , por lo que, si a(k) n,n
está
u k A kF , el riterio dado en (27) podría
en ero uando este es más grande que
(k)
an,n ,
en uyo aso
(k)
an,n
no es una
buena aproxima ión al autovalor en forma relativa.
Un riterio que se emplea en mu has de las implementa iones del algoritmo
QR,
es el
propuesto por Wilkinson, en el que se ompara el elemento de la subdiagonal on los
40
elementos ve inos de la diagonal. De esta manera, en la itera ión
k,
el elemento
(k)
ai,i−1
es
olo ado en ero si se umple que
(k) (k) (k) |ai,i−1 | ≤ u |ai−1,i−1 | + |ai,i |
(28)
El riterio dado en la e ua ión (28) es es en ialmente una prueba heurísti a, ya que a pesar de ser utilizada y uyo omportamiento ha sido el orre to, no hay una teoría matemáti a que la sostenga. Sin embargo si el riterio (28) se umple, enton es automáti amente se satisfa e el riterio (27).
2.4 Computa ión Paralela Conven ionalmente, el diseño de omputadores involu ra un simple ujo de instru
iones. Estas instru
iones son pro esadas de forma se uen ial lo ual onlleva a un inter ambio de la data que está en memoria a unidades fun ionales, para luego regresarla a la memoria. De forma más espe í a:
Una instru
ión es alar es onseguida y se de odi a.
Se al ulan las dire
iones de la data que van a ser usadas omo operandos.
Los operandos son tomados de la memoria.
El ál ulo es llevado a abo en una unidad fun ional.
El operando resultante es es rito nuevamente en memoria.
Se puede de ir enton es, que un omputador se uen ial onsiste de una memoria one tada a un pro esador por medio de un bus de datos. Sin embargo, estos omponentes (pro esador, memoria y bus de datos) presentan uellos de botella sobre el pro esamiento general del sistema omputarizado. Como es señalado por Grama et al (2003) a través de
41
los años se han realizado distintas innova iones sobre la arquite tura de los omputadores para manejar de forma e iente estos uellos de botella. Una de las más importantes innova iones es la multipli idad en unidades de pro esamiento, buses de datos y unidades de memoria.
El progreso de la omputa ión ha estado motivado a simula iones numéri as de sistemas
omplejos tales omo lima, dispositivos me áni os, ir uitos ele tróni os, pro esos de fabri a ión, rea
iones quími as, entre otros. Sin embargo, realizar tales simula iones, a un buen nivel, han llevado a que el omputador sea apaz de pro esar grandes antidades de informa ión de forma so sti ada. Estas apli a iones pueden in luir video onferen ia, ambientes de trabajo olaborativo, diagnósti os en medi ina on la asisten ia del omputador, bases de datos paralelas empleadas para la toma de de isiones y omputa ión ientí a y realidad virtual.
De esta manera surgió una demanda por un rendimiento más e iente, por lo que se fueron realizando modi a iones para mejorar el diseño de los omputadores. Sin embargo, se hizo evidente que un número de fa tores fueron limitando tales mejoras, omo por ejemplo, la velo idad de inter ambio entre dispositivos, retardos de inter one
ión, et . Como señala Dongarra et al (1998), si se realizara una mejora dramáti a en estos fa tores, todavía existe un fa tor que limita el rendimiento, y este es la velo idad de la luz, sin embargo los omputadores de hoy en día poseen tiempos de i lo que están sobre el orden de los nanosegundos.
Enfrentando esta limita ión fundamental, los diseñadores de omputadores se han movilizado en dire
ión del paralelismo. De esta forma, la omputa ión paralela no es solo una estrategia para al anzar un alto rendimiento, es una visión pre isa de omo la omputa ión puede es alar de un sen illo pro esador a una poten ia omputa ional que no posee
42
límites. Sin embargo omo plantea Dongarra et al (2003), el amino al paralelismo es alable ha sido duro, pero no se puede ver a la omputa ión paralela omo un éxito in ompetente. La verdad es que la omputa ión paralela ha permitido re er las velo idades máximas de super omputadores, en una tarifa que ex ede la ley de Moore, la ual estable e que el fun ionamiento del pro esador es doblado ada 18 meses.
2.4.1 Computador Paralelo Se puede visualizar a un omputador paralelo, omo un onjunto de pro esadores que están disponibles para trabajar ooperativamente, y de esta manera resolver un problema
omputa ional. Esta de ni ión es lo su ientemente amplia para in luir super omputadores paralelos que poseen ientos o miles de pro esadores, redes de esta iones de trabajo y redes de esta iones de trabajo on múltiples pro esadores.
La se uen ia de instru
iones de odi adas y eje utadas por el CPU es ono ido omo ujo de instru
iones. De manera similar, un ujo de data, se re ere a la se uen ia de operandos espe i ados para tales instru
iones. En Dongarra et al (1998) en ontramos la
ono ida taxonomía de Flynn, la ual es la manera lási a de organizar las omputadoras, y aunque no ubre todas las posibles arquite turas, propor iona una importante penetra ión en varias arquite turas de omputadoras. Esto da lugar a uatro tipos de omputadoras, empleando se uen ias de instru
iones sen illas o múltiples y se uen ias de datos sen illas o múltiples, las uales son:
Single Instru tion stream, Single Data stream ): Este es el modelo tradi ional
SISD (
de omputa ión se uen ial, donde una unidad de pro esamiento re ibe una sola se uen ia de instru
iones que operan en una se uen ia de datos. Un esquema de su fun ionamiento es mostrado en la gura 7.
Single Instru tion stream, Multiple Data stream ):
SIMD (
A diferen ia de SISD, en
43
Figura 7: Modelo de paralelismo SISD
este aso se tienen múltiples pro esadores que sin ronizadamente eje utan la misma se uen ia de instru
iones, pero en diferentes datos. El tipo de memoria que estos sistemas utilizan es distribuida. Un esquema de este modelo es visualizado en la gura 8.
Figura 8: Modelo de paralelismo SIMD
Multiple Instru tion stream, Single Data stream ): En este modelo, se uen ias
MISD (
de instru
iones pasan a través de múltiples pro esadores. Diferentes opera iones son realizadas en diversos pro esadores. Cada pro esador posee su propia unidad de ontrol y omparten una memoria en omún. La gura 9 representa una esquematiza ión de este modelo.
Multiple Instru tion stream, Multiple Data stream ): Este tipo de omputador
MIMD (
es paralela al igual que las SIMD, la diferen ia es que los MIMD son asín ronos, es de ir, no hay un reloj entral. Cada pro esador en un sistema MIMD puede eje utar su propia se uen ia de instru
iones y tener sus propios datos. Esta ara terísti a es la más general y poderosa de esta lasi a ión.
Desde algún tiempo, el aprove hamiento del multipro esamiento ha sido a tivamente in-
44
Figura 9: Modelo de paralelismo MISD
vestigado sobre las arquite turas SIMD y MIMD, aunque, omo estable e Dongarra et al (1998), la úni a máquina en el mer ado on arquite tura SIMD es desarrollada por la ompañia MasPar. El mismo autor estable e que aparentemente las máquinas SIMD tienen mu has limita iones sobre un número importante de problemas ientí os y por lo tanto no han sido onsideradas para propósito general en la omputa ión.
Así, el prin ipal desarrollo de la omputa ión ientí a paralela re ae en los omputadores pertene ientes a la lase MIMD. Estos sistemas MIMD se lasi an a su vez de la siguiente manera:
Sistemas de Memoria Compartida: En este tipo de sistemas ada pro esador tiene a
eso a toda la memoria, es de ir, hay un espa io de dire
ionamiento ompartido. Se tienen tiempos de a
eso a memoria uniformes ya que todos los pro esadores se en uentran igualmente omuni ados on la memoria prin ipal y las le turas y es rituras de todos los pro esadores tienen exa tamente las mismas laten ias; y además el a
eso a memoria es por medio de un du to omún. En esta on gura ión, debe asegurarse que los pro esadores no tengan a
eso simultáneamente a regiones de memoria de una manera en la que pueda o
urrir algún error. La gura 10 esquematiza este tipo de arquite tura.
45
Figura 10: Modelo MIMD on Memoria Compartida
Sistemas de Memoria Distribuida: Estos sistemas tienen su propia memoria lo al. Los pro esadores pueden ompartir informa ión solamente enviando mensajes, es de ir, si un pro esador requiere los datos ontenidos en la memoria de otro pro esador, deberá enviar un mensaje soli itándolos. Esta omuni a ión se le ono e omo Paso de Mensajes. Estos sistemas pueden ser representados omo en la gura 11
Figura 11: Modelo MIMD on Memoria Distribuida
Sistemas de Memoria Compartida Distribuida: Es un luster o una parti ión de pro esadores que tienen a
eso a una memoria ompartida omún pero sin un anal
ompartido. Esto es, físi amente ada pro esador posee su memoria lo al y se inter one ta on otros pro esadores por medio de un dispositivo de alta velo idad, y todos ven las memorias de ada uno omo un espa io de dire
iones globales. El a
eso a
46
la memoria de diferentes lusters se realiza bajo el esquema de A
eso a Memoria No Uniforme (NUMA), la ual toma menos tiempo en a
esar a la memoria lo al de un pro esador que a
esar a memoria remota de otro pro esador. Un esquema de esta arquite tura está representada en la gura 12.
Figura 12: Modelo MIMD on Memoria Compartida Distribuida
2.4.2 Programa ión Multithread Multithreading
es un paradigma moderno de programa ión para implementar apli a iones
paralelas sobre multipro esadores de memoria ompartida, bus ando mejorar el rendimiento del omputador. La programa ión
multithreading
permite que un programa pueda ser
eje utado en múltiples piezas de ódigo simultáneamente. Evidentemente, al utilizar este tipo de programa ión hay que tener uidado en omo la data es a
esada, esto debido a que la memoria es ompartida, existen datos ompartidos que pueden ser utilizados por más de un pro esador al mismo tiempo.
La programa ión
multithread es llavada a abo sobre arquite turas multithreaded, las uales
apare en debido a las limita iones que imponen las dependen ias de datos y la alta laten ia de memoria en los pro esadores onven ionales. Estas arquite turas se ara terizan por mantener varios
threads
eje utándose a la vez dentro del pro esador. Di hos
threads
47
pueden ser de una misma apli a ión, on lo que se redu iría su tiempo de eje u ión, o de apli a iones diferentes, on lo que se aumentaría la produ tividad del sistema.
El problema prin ipal de las arquite turas ra la genera ión de los
threads,
multithread
es el soporte que requieren pa-
espe ialmente si lo que se pretende es redu ir el tiempo
de eje u ión de una apli a ión extrayendo lo que se ono e omo paralelismo a nivel de
threads. A este respe to existen varias alternativas: el ompilador puede realizar un análisis del ódigo y determinar que
threads
se pueden generar, otra posibilidad es que sea el
programador quien, a través de dire tivas (por ejemplo OpenMP, Pthreads), espe i que el número de
threads de que onsta la apli a ión y en que segmentos de ódigo se generarán.
El desarrollo de la programa ión
multithread
en sistemas UNIX omenzó a mediados de
1980, pero a pesar de existir un a uerdo de lo que es son ne esarias para soportar
multithreading,
multithreading
y que ara terísti as
las interfa es usadas para implementar este
tipo de programa ión han variado grandemente.
Algoritmos numéri os y apli a iones on un alto grado de paralelismo, omo por ejemplo la multipli a ión de matri es, son eje utados mu ho más rápido uando son implementados mediante el uso de
threads
sobre multipro esadores.
Al momento de sele
ionar un paradigma de programa ión paralela se debe tomar en
uenta los siguientes objetivos prin ipales:
Explotar el máximo rendimiento de la máquina sobre una plataforma parti ular.
Portabilidad del ódigo y rendimiento a través de varias plataformas de alto rendimiento.
Una programa ión que permita la rea ión sen illa de orre tos, on ables y e ien-
48
tes programas.
2.5 Estru tura de Alma enamiento Aunque existen diversas formas de alma enar una matriz, la idea prin ipal de los formatos de alma enamiento es alma enar y manejar de forma e iente los elementos no nulos de una matriz, alma enándolos de forma ontinua empleando un arreglo lineal. Entre estos formatos, los más ono idos son
Compressed Row Storage (CRS), Compressed
Column Storage (CCS), Jagged Diagonal Storage (JDS), Transpose Jagged Diagonal Storage (TJDS).
La bibliote a UCSparseLib está basada en formatos de alma enamiento tales
omo CRS ó CCS.
2.5.1 Compressed Row Storage En el formato CRS, los elementos no nulos son alma enados por la, para ellos se emplean tres arreglos lineales, uno que permita alma enar los elementos reales no nulos y dos arreglos de números enteros, uno para alma enar los índi es de las olumnas de los
omponentes del arreglo de no nulos, mientras que el otro arreglo alma ena la posi ión ini ial de ada la.
Para ilustrar omo los elementos no nulos son alma enados bajo este esquema, supongamos una matriz
A
de tamaño
5 × 5,
on elementos no nulos
ai,j .
49
a 0 0 0 a1,4 1,1 0 a2,2 a2,3 0 0 A= 0 a3,2 a3,3 a3,4 0 0 0 a4,3 a4,4 a4,5 0 0 0 a5,4 a5,5
(29)
Empleando el formato CRS des rito anteriormente, los elementos no nulos quedan alma enados en los siguientes arreglos.
lista_de_valores
a1,1
a1,4
a2,2
índi e_de_ olumnas
a2,3 1
4
posi ión_ini io
a3,2 2
1
3
3
a3,3 2
5
a3,4 3
8
4
11
a4,3 3
4
a4,4 5
a4,5 4
a5,4
a5,5
5
13
2.5.2 Compressed Column Storage El formato de alma enamiento CCS, es similar al formato CRS. Este posee tres ve tores, un ve tor real para alma enar los elementos no nulos y dos ve tores de enteros, uno para alma enar los índi es de las las de los elementos no nulos, mientras que el otro alma ena la posi ión ini ial de ada olumna.
El alma enamiento de los elementos no nulos de la matriz dada en (29) en el formato CCS, es llevado a abo mediante el uso de los siguientes arreglos.
lista_de_valores
a1,1
a2,2
a3,2
a2,3
a3,3
a4,3
a3,4
a4,4
a5,4
a1,4
a4,5
a5,5
50
índi e_de_ las
1
2
3
posi ión_ini io
2
1
3
2
4
4
3
7
4
5
10
1
4
5
13
El formato CCS es el equivalente al formato CRS de la matriz
AT , tal omo es estable ido
en Bai et al (2000).
3 De ni ión de Términos A elera ión: El valor
S
es un multipli ador que indi a uantas ve es es más rápido el
programa paralelo on respe to al programa se uen ial.
S
viene dado por la siguiente
rela ión:
S(P ) = T (1)/T (P ) donde
T (n)
es el tiempo total de eje u ión sobre un sistema que posee
n
elemen-
tos de pro esamiento. Cuando la a elera ión es igual al número de elementos de pro esamiento, se di e que la a elera ión es perfe tamente lineal.
Ben hmark :
Pro eso sistemáti o y ontinuo para omparar resultados en términos
de produ tividad y rendimiento sobre iertos patrones.
Cluster : Un onjunto de omputadoras independientes inter one tadas, usadas omo un re urso uni ado de ómputo.
Convergen ia: Dirigirse dos o más líneas a unirse en un punto. Con urrir varias opiniones al mismo n.
Densidad de una matriz: Número real determinado omo la antidad de elementos no nulos de una matriz dispersa dividida entre el produ to de el número de las y el número de olumnas de la matriz.
51
Es alabilidad: Un sistema se di e es alable si es apaz de in rementar sus re ursos y rendimiento a las ne esidades soli itadas de manera efe tiva.
Espe tro de una matriz: Se llama espe tro de una matriz todos los autovalores de
A ∈ Cn×n
al onjunto de
A:
σ(A) = {λ ∈ C : λ
es autovalor de
A} = {λ ∈ C : ker(A − λI) 6= {0}}.
Estru turas dispersas: Una estru tura dispersa es una matriz poblada prin ipalmente
on elementos nulos.
Flags :
Conjunto de op iones de ompila ión que permiten optener un ódigo más
robusto.
Fork-join : Pro eso de genera ión de threads. En donde un pro eso maestro se bifur a en un onjunto de pro esos es lavos en una región paralela. Al nal de la región paralela estos pro esos es lavos nalizan quedando solamente el pro eso maestro.
Hamiltoniano: El Hamiltoniano uánti o
H
es un observable que orresponde a la
energía total del sistema. El Hamiltoniano genera la evolu ión temporal de los estados
uánti os a través de la e ua ión de S hrödinger.
Hardware : Se denomina hardware o soporte físi o al onjunto de elementos materiales que omponen un omputador. En di ho onjunto se in luyen los dispositivos ele tróni os y ele trome áni os, ir uitos, ables, tarjetas, armarios o ajas, periféri os de todo tipo y otros elementos físi os.
Matriz de Hessenberg:
H
es una matriz de Hessenberg si es triangular on una
diagonal adya ente distinta de ero.
Matriz diagonal: Matriz
D
que posee elementos
di,j = 0
si
i 6= j .
52
Matriz ortogonal:
Q
males, de modo que ángulos,
es una matriz ortogonal si es uadrada on olumnas ortonor-
Qt Q = I
k Qx k=k x k
y
impli a que
Qt = Q−1 .
(Qx)t (Qy) = xt y .
Conserva las longitudes y los
Ejemplos: Rota ión, Re exión, Permu-
ta ión.
Matriz tridiagonal: Matriz
T
que posee elementos
ti,j = 0
si
|i − j| > 1. T −1
es de
rango 1 por en ima y por debajo de la diagonal.
Métodos Numéri os: Son té ni as mediante las uales es posible formular problemas que se pueden resolver usando opera iones aritméti as y que nalmente son resueltos usando una omputadora.
Multipro esamiento: Té ni a en la que se emplean varios pro esadores en un solo sistema omputa ional, on la nalidad de in rementar el poder de pro esamiento de la máquina.
Multithreading :
Conjunto de pro esos eje utandose en una se uen ia rápida dentro
de un mismo programa, sin importar ual es el método lógi o de multitarea que está siendo usado por el sistema operativo. Capa idad de poder eje utar múltiples
threads
ompletamente transparente al usuario.
Overhead : Pro esamiento adi ional resultante de una arga inne esaria de datos. Paralelismo: Forma e az de pro eso de la informa ión, que favore e la explota ión de su esos on urrentes en el pro eso de ómputo.
Rendimiento: Es la efe tividad del desempeño de una omputadora, sobre una apli a ión en parti ular. En las medi iones de rendimiento están involu rados velo idad,
osto y e ien ia.
Software :
Es la parte lógi a del omputador, esto es, el onjunto de programas que
puede eje utar el
hardware
para la realiza ión de las tareas de omputa ión a las
que se destina. Es el onjunto de instru
iones que permite la utiliza ión del equipo.
53
Super omputador: Es un omputador de propósito general apaz de resolver problemas individuales on desempeños muy altos respe to a otros omputadores onstruidos en la misma épo a. Según esta de ni ión, siempre han existido super omputadores. A tualmente todo super omputador es paralelo; algunos tienen po os pero rápidos y aros pro esadores y otros mu hos pro esadores baratos.
Template :
Es un modelo, forma o patrón utilizado omo guía para elaborar nuevas
unidades de ódigo.
Thread : Es la entidad bási a por la ual el sistema operativo asigna tiempo de CPU. Un
thread
puede eje utar una parte del ódigo de la apli a ión, in luyendo alguna
parte que está siendo eje utada por otro thread. Todos los
threads
de un pro eso
omparten el espa io virtual de dire
ionamiento, variables globales y re ursos del sistema operativo de los pro esos.
Unidad de redondeo: Se de ne a la unidad de redondeo o presi ión de la máquina al número más pequeño de la forma
ε = 2−k
tal que
f l(1 + ε) 6= 1.
ara terísti a de la máquina, su sitema operativo, et .
Este valor es una
CAPÍTULO III: MARCO METODOLÓGICO
55
CAPÍTULO III MARCO METODOLÓGICO 1 Tipo y Diseño de la Investiga ión El nivel de profundidad de la investiga ión es de tipo exploratoria y des riptiva:
1. Exploratoria: el estudio está orientado a de nir un mar o on eptual rela ionado
on las teorías de diagonaliza ión de una matriz simétri a, mediante el método de itera iones
QR, así omo el estudio de su omportamiento. En tal sentido, se realizará
una mez la de on eptos rela ionados on la nalidad de analizar el impa to que tiene di has teorías sobre los autovalores, lo uál permitirá obtener informa ión relevante para onstruir los algoritmos omputarizados ne esarios para desarrollar el estudio planteado.
2. Des riptiva: el estudio está orientado a realizar el análisis de las teorías de la forma diagonal de una matriz, on el método de Ja obi, y teoría itera ión
QR,
así omo
de los algoritmos involu rados. Posteriormente, el estudio ontempla el análisis omparativo de los resultados obtenidos. Como produ to de estos estudios se obtendrá informa ión que soportará la onstru
ión de los algoritmos omputarizados que permitan la resolu ión de problemas de autovalores.
En tal sentido, Hernández, Fernández y Baptista (1998) indi a que "los estudios des riptivos miden de manera más bien independiente los on eptos o variables a los que se re eren. Aunque pueden integrar las medi iones de ada una de di has variables para de ir ómo es y ómo se mani esta en fenómeno de interés...".
56
2 Modalidad de la Investiga ión La investiga ión que indi ará el nivel de profundidad que se desea obtener en este estudio y que está re ejada en los objetivos de este trabajo estará on ebida dentro de la modalidad de Investiga ión de Campo, omo lo muestra el manual de Trabajos de Grado y Maestría de la Universidad Pedagógi a Experimental Libertador (1990), ya que ontempla la investiga ión bibliográ a de on eptos rela ionados on las teorías ya men ionadas, las
uales serán sometidas a un pro eso de estudio y evalua ión, lo que sentará las bases para la onstru
ión de los algoritmos omputarizados que apoyen la resolu ión de problemas de autovalores.
Algunas ve es una investiga ión puede in luir elementos de distintos tipos de estudios,
omo es señalado por Hernández et al (1998). Ésta, estará basada en un tipo de estudio exploratorio, dado que se ne esita avanzar en el ono imiento de las té ni as existentes para el ál ulo de los autovalores de una matriz, esto on el n de reunir informa ión para posteriores desarrollos.
También se utilizarán elementos de una investiga ión omparativa, lo ual permitirá estable er diferen ias y semejanzas entre las distintas té ni as para diagonalizar una matriz, sin estable er orrela iones o ausalidades entre ellos. Finalmente se hará uso de elementos de una investiga ión evaluativo, on la ual se estudiarán los resultados obtenidos de las des rip iones, análisis y ompara iones he has a las distintas té ni as para el ál ulo de autovalores, entre ellas las que son motivo de este estudio, el Método de las Itera iones
QR
y sus variantes.
Además la realiza ión de esta investiga ión impli a el desarrollo e implementa ión de rutinas de
software
espe í as, en el ontexto de una bibliote a numéri a, desarrollada
57
para manipular omputa ionalmente problemas rela ionados on el álgebra lineal, entre ellos, el problema del ál ulo de autovalores.
3 Té ni as e Instrumentos para Re ole tar la Informa ión La re ole
ión de la presente investiga ión se realizará por medio de onsultas a fuentes de informa ión primarias y se undarias:
1. Informa ión primaria: Para este n se utilizarán las té ni as de observa ión y la de re opila ión y análisis bibliográ o.
a)
Observa ión: Esta té ni a sirve de apoyo en todas las fases de la investiga ión para re opilar, analizar, onstruir y omparar los resultados obtenidos.
b)
Re opila ión y análisis bibliográ o: onsidera la búsqueda de informa ión relativa a las diferentes teorías planteadas omo mar o referen ial para el uso de los autovalores, así omo su análisis, on la nalidad de sele
ionar la informa ión lave referen ial al momento de profundizar el estudio y desarrollo de los algoritmos omputarizados.
)
Entrevistas a expertos en el área de quími a uánti a, los uales ofre erán informa ión relativa sobre la apli a ión de los autovalores en esta área de interés.
4 Fases Metodológi as La metodología a utilizar estará enfo ada a la resolu ión de problemas de índole matemáti o y omputa ional, ya que se plantea la resolu ión del problema de autovalores apli ando métodos numéri os, omo lo es, el método de las itera iones para el ál ulo de autovalores.
QR
y sus variantes
58
En tal sentido, se han estru turado uatro fases, ada una de las uales se ompleta una vez que se tienen iertos produ tos, los uales son utilizados omo elementos de la siguiente fase. En el uadro 1 se muestra una des rip ión general de las uatro fases y los produ tos que se obtendrán en ada una de ellas.
FASES
Cuadro 1: Fases Metodológi as
Fo alizar : Identi ar el problema, y • des ribirlo.
PRODUCTOS
Considera el desarrollo de una de lara ión es rita del problema, la obten ión de ante edentes teóri os y análisis de experien ias similares en el manejo de las teorías aso ia-
Analizar : Estudiar las teorías de forma diagonal de una matriz, on el método de Ja obi, y teoría itera ión
QR
Diseñar : Desarrollar los algoritmos
omputarizados que permitan la
das a la investiga ión.
•
Teoría de forma diagonal de una
• • •
Método de Ja obi.
resolu ión de problemas de autovalores
Comparar :
Desarrollar el análisis
respe tivo sobre el desarrollo de los
matriz. Teoría de Itera ión
QR.
Diseñar los algoritmos omputarizados
basados
en
los
elementos
men ionados.
•
Probar los algoritmos omputariza-
•
Determinar el grado de exa titud de
dos en la bibliote a UCSpaseLib. las respuestas dadas por los algorit-
algoritmos basados en ál ulo de
mos paralelos diseñados previamen-
autovalores, apli ados en los proble-
te.
mas de quími a uánti a
Seguidamente, se presentan las uatro fases de la metodología a seguir para el desarrollo de la presente investiga ión:
4.1 Fase I: Fo alizar Esta fase tiene omo nalidad enfo ar el problema y de nirlo on laridad. Comprenderá las siguientes a tividades:
59
Revisión do umental: Esta etapa onsiste en la re opila ión de informa ión sobre las teorías men ionadas anteriormente. La re ole
ión de la informa ión se llevará a abo por medio de los siguientes medios:
Revisión de Fuentes Externas:
Visitas a bibliote as y diversos entros de informa ión.
Revisión de textos e investiga iones a nes a los temas en análisis.
Consultas a páginas web, vía Internet.
Estudios previos realizados sobre los autovalores.
Revisión de Fuentes expertas:
Se realizarán diversas onsultas al personal experto en el estudio de los autovalores, su apli a ión y respuestas obtenidas. Se onsultará a personal experto en otras áreas de investiga ión, a n de validar los on eptos presentados.
4.2 Fase II: Analizar La nalidad de esta fase es entender el problema, profundizar en los aspe tos teóri os men ionados. Dentro de esta fase se onsiderarán las siguientes a tividades:
1. Estudio de on eptos laves: onsidera el estudio en profundidad de todo el mar o
on eptual de las teorías aso iadas on los autovalores en términos de: sus bases
on eptuales, su algoritmo de resolu ión, asos de ex ep ión, entre otros aspe tos.
2. Determinar los elementos a onsiderar en el diseño de los algoritmos omputarizados.
60
4.3 Fase III: Diseñar Esta fase onsidera el diseño de un
software
orientado a la resolu ión de autovalores.
Para ello se eje utarán las siguientes a tividades:
1. Diseño del
software
para la resolu ión de autovalores
a)
Evaluar la herramienta a utilizar para desarrollar el
b)
Diseño general de los algoritmos a onsiderar en el
)
Programa ión del
software
en uestión.
software.
software.
2. Eje u ión de pruebas al
software
diseñado.
4.4 Fase IV: Comparar En esta fase se ompararán los resultados numéri os en fun ión al tiempo de eje u ión ne esitado para obtener los autovalores de una matriz, a través de los algoritmos, tanto seriales omo paralelos, los uales fueron implementados en la fase anterior.
5 Plan de Trabajo A ontinua ión se des riben los re ursos involu rados dentro del presente trabajo de investiga ión, los uales umplen un papel primordial dentro de ada una de las fases metodológi as a ser ubiertas:
5.1 Re ursos Humanos El personal humano que intervendrá en el desarrollo de la presente investiga ión, es presentado en el uadro 2 en onjunto on una breve des rip ión de la parti ipa ión en el
61
trabajo.
Cuadro 2: Re ursos Humanos
RECURSO HUMANO
PARTICIPACIÓN
Personal de institu iones de
Servirán de apoyo para iden-
edu a ión superior e institu-
ti ar bibliografía a tualiza-
iones rela ionadas on el te-
da y ompleta, rela ionada
ma: UCV, UC, IVIC.
on el tema
Peronal experto en el estu-
Servirán de usuarios y eva-
dio del ál ulo de autovalo-
luadores del desempeño del
res on la apli a ión del mé-
software
todo de Ja obi y las itera io-
on lusiones obtenidas en la
nes
QR.
a diseñar y de las
investiga ión.
5.2 Re ursos Finan ieros Para el desarrollo del trabajo de investiga ión no se requerirá una inversión importante en re ursos, ya que se uenta, on anti ipa ión, de los mismos, o son fa ilitados por alguna institu ión superior omo lo es, en este aso, la Universidad de Carabobo, los mismos son mostrados en el uadro 3.
Cuadro 3: Re ursos Finan ieros
RECURSOS
Computadora. Impresora. A
esorios para omputadoras.
5.3 Re ursos Institu ionales Los prin ipales re ursos institu ionales on los que se ontará, para el desarrollo de esta investiga ión son mostrados en el uadro 4, en el mismo se presenta una breve des rip ión de la parti ipa ión de la institu ión en la investiga ión.
62
Cuadro 4: Re ursos Institu ionales
RECURSO INSTITUCIONAL
PARTICIPACIÓN
Tutores de la Universidad de Carabobo.
Servirán de guía y apoyo para desarrollar
Re ursos de omputa ión.
Servirán de punto de partida para estable-
el trabajo de investiga ión
er onexiones on Internet y poder obtener informa ión lave y a tualizada sobre el trabajo a realizar y el
software
a desa-
rrollar.
6 Resultados de las Fases Metodológi as A ontinua ión se muestran los resultados obtenidos en las primeras tres etapas, resultados que permitirán realizar las pruebas de nitivas y de esta forma omparar numéri amente los valores obtenidos en la siguiente unidad.
6.1 FASE I En esta primera fase se estable e, on ayuda del Profesor Fernando Ruette investigador del IVIC, la ne esidad del ál ulo de autovalores y autove tores en una apli a ión de la quími a uánti a. Des ribiendo la modela ión del problema quími o en un problema matemáti o al ual se le dará solu ión omputa ionalmente. Al nal de esta se
ión se indi an los resultados obtenidos en esta fase, los uales servirán de apoyo para de nir
laramente el método numéri o a desarrollar en la fase II.
El Profesor en su exposi ión planteó las diferentes áreas de investiga ión que se llevan a abo en el Departamento de Quími a Computa ional, entre ellas se en uentra el área de quími a uánti a en espe ial atálisis. En esta área se está desarrollando un
software
de-
nominado CATIVIC, el ual es un programa quími o- uánti o desarrollado espe ialmente para modelar rea
iones en el área de atálisis heterogénea. Sin embargo, puede utilizarse
63
también para sistemas orgáni os y organometáli os. Este programa permite el estudio de los pro esos que o urren en una rea
ión atalíti a, tales omo adsor ión, rompimiento y forma ión de enla es, forma ión de omplejos pre ursores, et . CATIVIC está fundamentado en la teoría de Simula ión de Fun ionales Paramétri os suponiendo bases mínimas óptimas transformadas, omo se estable e en Ruette et al (2004). Esto permite obtener resultados teóri os razonables en rela ión on los datos experimentales, en un tiempo de
ál ulo relativamente orto.
En tal exposi ión se des ribieron los on eptos fundamentales, tanto del Análisis Fun ional omo de la Quími a Cuánti a, de modo tal de estable er la idea entral del problema.
De ni ión 3: Espa io de Bana h: Es un espa io ve torial real o omplejo que posee norma y es ompleto en la métri a de nida por la norma.
De ni ión 4: Espa io de Hilbert (H): Un espa io de Bana h que tiene de nido un produ to interno se llama espa io de Hilbert, en donde el produ to interno es la extensión del produ to es alar en un espa io de fun iones.
< ΨI , ΨJ >=
De esta forma se tiene que el operador donde
H
Z
Ψ∗I ΨJ dτ
HΨI = EI ΨI
es el operador hamiltoniano y
Ψ ∈ H,
es un mapeo de fun iones en otro,
así que
H
es un operador hermíti o
(autovalor real). Mientras que por otro lado, el fun ional de nido por
f (ΨI ) = EI
< ΨI , HΨI >=
es un mapeo entre un espa io de fun iones en el espa io de los reales (R).
Para sistemas mole ulares,
H
tiempo, de la siguiente forma:
está de nido para el aso no relativista e independiente del
64
De esta forma, y de niendo un onjunto de fun iones bases para un ele trón
{ϕ}n , enton es
una fun ión para mu hos ele trones estaría de nida por:
ΨI = det φI1 , φI2 , . . . , φIa , . . . , φIn
donde
φa =
n X
dak ϕk
k=1
De este modo, la apli a ión del Método Varia ional usando
Ψ = det {φ1 , φ2 , . . . , φa , . . . , φn }
resulta en las e ua iones de Hartree-Fo k.
Por lo expli ado por el profesor Ruette en la entrevista, se puede on luir que la e ua ión de S hrödinger derivan en las e ua iones de Roothan, que son aproxima iones de las e ua iones de Hartree-Fo k, lo ual se puede mostrar en la siguiente gura
De esta forma se llega a las e ua iones matri iales, en donde se supone las e ua iones de Roothan se pueden es ribir omo autovalores a resolver.
F ′ C ′ = C ′ ε,
E
diagonal, y
la ual es la e ua ión de
65
Figura 13: Aproxima ión de la E ua ión de S hrödinger
El pro eso indi ado es llevado a abo por medio del
software
CATIVIC, apli ando el
método SCF des rito en el apítulo I. En este pro eso se bus a diagonalizar la matriz
F ′,
matriz ono ida omo matriz de Fo k.
Como resultados de esta etapa se tiene que:
El Instituto Venezolano para la Investiga ión Cientí a (IVIC) ha desarrollado un
software
que permite resolver de forma aproximada la e ua ión mole ular de S hrö-
dinger, por medio del método SCF.
En di ho
software, existe la ne esidad de diagonalizar la matriz de Fo k, pro eso que
puede ha er lento el método SCF.
La matriz de Fo k, es una matriz que posee las ara terísti as prin ipales de ser simétri a, on estru tura dispersa y puede poseer un gran tamaño.
66
6.2 FASE II En esta segunda fase se realizó una revisión teóri a de los distintos métodos para el ál ulo de autovalores, bus ando sus prin ipales ventajas y desventajas en fun ión al objetivo general planteado en el apítulo I y tomando en uenta los resultados obtenidos en la fase I, el ual onsiste en hallar todos los autovalores de la matriz de Fo k y sus posibles autove tores.
De esta forma, realizando un estudio omparativo de los distintos métodos que permiten
al ular los autovalores de una matriz, se pretende estable er ual es el método numéri o que permite obtener de una matriz simétri a sus autovalores. Entre los mismos, se en uentran:
Método de Poten ias.
Arnoldi reini ializado.
Ja obi-Davidson
Método
QR
En Li (2004) se realizó un estudio en donde son apli ados estos métodos y se obtuvieron las siguientes on lusiones. El método de Arnoldi reinizializado es sensible al ve tor ini ial, en tanto que Ja obi-Davidson onverge lentamente para iertos asos. Por otro lado el método de Poten ias solo permite obtener el autovalor más grande en magnitud, por lo que para obtener todos los autovalores de la matriz, habría que apli ar un pro eso de de a ión de la misma, una vez obtenido el autovalor, mientras que el método
QR,
a pesar
de que demanda mayor antidad de alma enamiento en memoria y tiempo de CPU uando la matriz es grande, este método es el que propor iona los resultados más robustos además de que permite obtener tanto los autovalores omo autove tores de la matriz.
67
Otro método para el ál ulo de Autovalores, en el que se pueden obtener tanto autovalores omo autove tores, es el ini iado por Ja obi en 1846, en el ual se bus a diagonalizar la matriz a través de matri es ortogonales, de modo tal que se obtiene una su esión de matri es que en teoría debe onverger a una matriz diagonal, pero este método tiende a ser lento y por lo tanto es preferible emplear otros métodos omo
QR,
omo es itado
por Gentle, Härdle y Mori (2004). La idea prin ipal de Ja obi fue tomada por Givens y Householder, es de ir, una serie de transforma iones de similitud son apli adas a la matriz
A,
para transformarla en una matriz tridiagonal y luego apli ando el método
QR
o
QL,
se obtiene la diagonaliza ión de la misma.
En Bai et al (2000) se reseña que el algoritmo más exitoso para resolver el problema de autovalores de una matriz
A
uadrada es el algoritmo
ito. Como se puede notar, el método
QR
QR
on desplazamiento implí-
es el método numéri o por ex elen ia para el
ál ulo de autovalores y autove tores, por lo que en este trabajo, se ha sele
ionado implementar el algoritmo
QR,
uyos detalles teóri os son presentados en el apítulo II.
De esta manera, los resultados obtenidos en la fase II, des riben el método numéri o a implementar en la fase III, de modo tal de que los resultados numéri os que se obtengan sean ables. Así, en la fase III, la tarea prin ipal será paralelizar el algoritmo se uen ial en fun ión de obtener resultados equivalentes pero en una medida menor de tiempo.
6.3 FASE III En fun ión a los resultados obtenidos en la fase I y la fase II, se pro edió a la implementa ión del algoritmo
QR
para matri es simétri as según las de ni iones y algoritmos
dados en el apítulo II. Los resultados más notables de las dos fases anteriores son las
68
siguientes:
La matriz generada por CATIVIC, es una matriz simétri a que puede superar el orden de los miles, y además posee una estru tura dispersa.
El algoritmo a paralelizar, es el método numéri o
QR
o
QL
para el ál ulo de auto-
valores y posibles autove tores.
Dado enton es que la matriz es simétri a, los algoritmos implementados solo alma enan el triángulo inferior de la matriz, pero además omo tal matriz posee una estru tura dispersa, se utilizó la bibliote a UCSpaseLib para el alma enamiento de la misma bajo los formatos CRS/CCS, según lo señalado en el apítulo II. La herramienta omputa ional que permitirá apli ar paralelismo a los algoritmos implementados será OpenMP.
A ontinua ión se realiza una breve des rip ión de la plataforma omputa ional empleada para el desarrollo y pruebas de los algoritmos en la bibliote a UCSparseLib, así mismo se des ribe la bibliote a UCSpaseLib y la bibliot a OpenMP. Luego se indi a omo los algoritmos fueron a oplados a la bibliote a UCSpaseLib.
6.3.1 Plataforma Computa ional Todas las pruebas se desarrollaron en máquinas on arquite tura pe í amente en el servidor
Sun re V440.
SUN mi rosystems,
es-
Este servidor posee uatro pro esadores Ul-
traSPARC IIIi de 1.28 GHz, on memoria entral de 8GB, ontrolado mediante el Sistema Operativo Solaris versión 9.
Toda la bibliote a UCSparseLib, en donde fue implementado el algoritmo de ál ulo de autovalores, fué desarrollada en ANSI C. Los
ags
de optimiza ión utilizados, los uales
69
permiten un ódigo óptimo, fueron los siguientes: -fast, -xar h=v9 y -xopenmp, este último
ag
permite que el ompilador re onoz a que existen dire tivas de la bibliote a OpenMP
en el ódigo.
6.3.1.1 Bibliote a UCSparseLib La bibliote a UCSparseLib, es una bibliote a reada para manejar matri es dispersas, en prin ipio, esta bibliote a fué reada para resolver sistemas lineales de e ua iones que involu raran matri es de gran tamaño, pero on la ara terísti a prin ipal de ser estru turas dispersas.
En esta bibliote a, para no tener la ne esidad de alma enar los elemento no signi ativos, se emplearon estru turas de alma enamiento, de modo tal de poder manejar matri es de gran tamaño que posean una baja densidad. Las estru turas bási as de alma enamiento de la bibliote a son CRS y CCS, aunque puede manejar otros formatos de alma enamiento distintos a los men ionados.
En nuestro aso, dado que las matri es on las uales trabajaremos poseen una estru tura dispersa, esta bibliote a fué empleada tanto para el alma enamiento de di has matri es, así omo para la elabora ión de los algoritmos prin ipales de ál ulo de autovalores mediante el método
QR
y
QL
y el método de tridiagonaliza ión.
6.3.1.2 El Estándar OPENMP OpenMP es un modelo de programa ión paralela para máquinas multipro esadores de memoria ompartida, el ual está basado en un onjunto de dire tivas de ompilador. OpenMP soporta un modelo de programa ión de variables ompartidas, en onjunto on
70
una bibliote a de rutinas y variables de ambiente.
OpenMP está disponible en mu has máquinas de memoria ompartida, in luyendo las
-NUMA, y es el estándar industrial soportado por la mayoría de los fabri antes, omo estable en Fowler y Greenough (2002). Estas dire tivas permiten al programador de nir
uales i los pueden ser paralelizables y omo las variables lo ales deben ser tratadas, es de ir, el ódigo serial puede ser onservado en mu hos asos, de esta manera se puede ha er un mantenimiento entre la versión serial y la versión paralela de una manera más sen illa, aunque sin embargo, a ve es se debe alterar el ódigo para explotar el paralelismo.
El programador puede on entrarse solo en aquellos i los que poseen mayor tiempo de
onsumo, y de esta forma el programador puede hallar aquellos i los que sean más seguros de paralelizar para ha er más e iente al ompilador. OpenMP también puede apli ar paralelismo de alto nivel, omo por ejemplo eje utando mu has subrutinas de forma paralela. De esta forma se pueden paralelizar ódigos existentes sin ne esidad de rees ribirlo ompletamente, lo ual no es posible on mu hos lenguajes de programa ión paralela, omo se puede eviden iar en el uadro 5, el ual podemos en ontrar en OpenMP (1997). En el mismo podemos notar que solo OpenMP y X3H5 permiten paraleliza ión in remental (pro eso de onvertir un algoritmo se uen ial en paralelo po o a po o), pero de estos dos modelos solo OpenMP es es alable.
6.3.2 Crea ión del Módulo EIGENVALUES para UCSparseLib A ontinua ión se des riben las rutinas prin ipales que onforman el módulo EIGENVALUES en la bibliote a UCSparseLib, de modo tal de llevar a abo las pruebas ne esarias tanto de forma se uen ial omo de forma paralela.
71
Cuadro 5: Modelos de Programa ión Paralela
X3H5 MPI Pthreads
HPF
OpenMP
Es alabilidad
no
si
a ve es
si
si
Paraleliza ión In remen-
si
no
no
no
si
si
si
si
si
si
si
si
no
si
si
si
no
no
si
si
de
si
no
no
si
si
Orientado al rendimiento
no
si
no
intenta
si
tal Portabilidad Enlazamiento
on
For-
tran Alto Nivel Soporta
paralelismo
datos Fuente: OpenMP (1997).
Las rutinas en argadas de llevar a abo el método
QR
o
QL
para el ál ulo de auto-
valores y posibles autove tores sobre una matriz simétri a tridiagonal poseen la siguiente des rip ión.
QR_Eig :
Cal ula solo los autovalores de una matriz simétri a tridiagonal. En esta
rutina no se utiliza ningún tipo de desplazamiento para a elerar la onvergen ia.
QR_Eig_Shift : Cal ula solo los autovalores de una matriz simétri a tridiagonal. La rutina plantea el uso del desplazamiento simple tomando
QR_Eig_Wilk_Shift :
(k)
σk = Ak,k .
Cal ula solo los autovalores de una matriz simétri a tridia-
gonal. La rutina plantea el uso del desplazamiento simple de nido por Wilkinson, estrategia ono ida omo
Wilkinson-shift.
tred2 : Transforma una matriz simétri a a una forma tridiagonal apli ando las transforma iones de Householder, esta rutina debe ser utilizada para tridiagonalizar la matriz si solo se desean sus autovalores.
tred3 : Transforma una matriz simétri a a una forma tridiagonal apli ando las transforma iones de Householder, esta rutina debe ser utilizada para tridiagonalizar la matriz si se desean de ella sus autovalores y autove tores mediante la rutina
tql2.
72
tql2 :
Obtiene los autovalores y autove tores de una matriz simétri a tridiagonal
empleando el algoritmo
tql1 : Obtiene todo
QL
QL
on desplazamiento implí ito.
los autovalores de una matriz simétri a tridiagonal empleando el mé-
on desplazamiento implí ito.
Triangular_Lower :
Alma ena el triángulo inferior de una matriz simétri a bajo el
mismo formato disperso de la matriz de entrada.
Los parámetros de entrada prin ipales y de salida en las rutinas men ionadas y des ritas anteriormente son los siguientes:
QR_Eig : La matriz tridiagonal viene suministrada en dos arreglos, uno en el que se alma ena la diagonal prin ipal y otro para alma enar la subdiagonal prin ipal. Los autovalores son retornados en un ter er arreglo.
QR_Eig_Shift : La matriz tridiagonal viene suministrada en dos arreglos, uno en el que se alma ena la diagonal prin ipal y otro para alma enar la subdiagonal prin ipal. Los autovalores son retornados en un ter er arreglo.
QR_Eig_Wilk_Shift : La matriz tridiagonal viene suministrada en dos arreglos, uno en el que se alma ena la diagonal prin ipal y otro para alma enar la subdiagonal prin ipal. Los autovalores son retornados en un ter er arreglo.
tred2 : Esta rutina tiene un parámetro de entrada el ual es el triángulo inferior de la matriz simétri a y retorna en dos ve tores la matriz tridiagonal, un ve tor en el que se alma ena la diagonal prin ipal y en otro se alma ena una de las dos subdiagonales prin ipales.
tred3 : Esta rutina tiene un parámetro de entrada el ual es el triángulo inferior de la matriz simétri a y retorna en dos ve tores la matriz tridiagonal, un ve tor en el que
73
se alma ena la diagonal prin ipal y en otro se alma ena una de las dos subdiagonales prin ipales y retorna la matriz de transforma ión para el ál ulo de autove tores.
tql2 : Esta rutina tiene omo parámetros de entrada la matriz tridiagonal y la matriz de transforma ión, obtenidas por la rutina
tred3, la matriz tridiagonal es alma enada
en dos arreglos, uno para la diagonal prin ipal y otro para la primera subdiagonal. Los autove tores son retornados en la matriz de transforma ión, mientras que los autovalores son retornados en la diagonal prin ipal de la matriz tridiagonal, la subdiagonal se pierde en el pro eso.
tql1 : Esta rutina tiene omo parámetros de entrada la matriz tridiagonal, obtenidas por la rutina
tred2, la matriz tridiagonal es alma enada en dos arreglos, uno para la
diagonal prin ipal y otro para la primera subdiagonal. Los autovalores son retornados en la diagonal de la matriz tridiagonal mientras que la subdiagonal se pierde en el pro eso.
Triangular_Lower : Esta rutina posee omo parámetro de entrada la matriz dispersa bajo algún formato (CRS/CCS) y omo parámetro de salida se obtiene bajo el mismo formato el triángulo inferior de la matriz de entrada, esta matriz es utilizada por las rutinas
tred2
En las rutinas
y
tred3
para su tridiagonaliza ión.
QR_Eig, QR_Eig_Shift
y
QR_Eig_Wilk_Shift,
se tienen otros paráme-
tros omo lo son el nivel de toleran ia a medir y el número máximo de itera iones que empleará el algoritmo
Las rutinas
QR
para el ál ulo de los autovalores.
QR_Eig, QR_Eig_Shift
y
QR_Eig_Wilk_Shift
fueron derivadas del pseudo-
ódigo presentado por Burden y Faires (1998) en su respe tivo texto, mientras que las rutinas
tred2, tred3, tql2
y
tql1
fueron derivadas de los pro edimientos es ritos en el lenguaje
ALGOL 60 presentados por Wilkinson, Reins h y Martin (1968) y por Wilkinson, Reins h
74
y Bowdler (1968).
En el anexo A se en uentran las abe eras de las rutinas que onforman el módulo EIGENVALUES de la bibliote a UCSparseLib, entre las que en ontramos las men ionadas anteriormente y de igual manera aquellas rutinas que sirven de soporte para el módulo en general, y que por no poseer una importan ia primordial, no fueron men ionadas en este tópi o. Mientras que en el anexo B se en uentra el programa prin ipal on el que se realizaron las pruebas, en el mismo se puede observar omo es la invo a ión de las rutinas planteadas para la resolu ión del problema de autovalores.
6.3.2.1 Considera iones de la Implementa ión En este trabajo se ha empleado el riterio dado en la e ua ión (28) para determinar uando
onsiderar que se ha obtenido un autovalor y proseguir on la búsqueda de los siguientes autovalores. En di ha e ua ión se ha tomado omo valor
u a la unidad de redondeo de la
máquina sobre la ual se está trabajando, de este modo aseguramos que los autovalores poseen una buena antidad de dígitos signi ativos.
Otra onsidera ión a tomar en uenta en la implementa ión de los algoritmos planteados tiene que ver on la simetría de la matriz a la ual se le desean obtener sus autovalores, debido a esto, no es ne esario alma enar toda la matriz, solo basta alma enar el triángulo inferior o superior de la matriz. En este aso, la implementa ión fué realizada tomando en
uenta que la matriz es triangular inferior. Es de notar que debido a esta onsidera ión la matriz dispersa resulta ser más dispersa aún, por lo que se redu e a así la mitad la densidad de la matriz.
75
6.3.3 Tridiagonaliza ión de Householder y uso del paralelismo El método para el ál ulo de autovalores por ex elen ia es el método ono ido omo Método
QR
QL
o
on desplazamiento implí ito, siempre y uando la matriz tenga una
forma tridiagonal para asegurar la onvergen ia del método, de este modo, antes de apli ar el algoritmo
QL,
se pro ede a obtener la matriz similar tridiagonal mediante el algoritmo
dado en la gura 5, en el mismo se puede notar que el ter er paso onsiste en generar el ve tor
u
i−esima
tomando los la de
A,
n−k
últimos elementos de la la
en el uarto paso se genera el ve tor
v
k
de
A
y son alma enados en la
el ual posee la siguiente forma.
√ v = ai,1 , ai,2 , . . . , ai,i−2 , ai,i−1 ∓ σ, 0, . . . , 0 Ahora en el quinto paso se genera el ve tor
p,
en donde ada elemento del ve tor es un
produ to matriz-ve tor, por lo que en este paso el algoritmo onsume mayor antidad de tiempo. En este paso del algoritmo se pueden usar
threads
de modo tal de minimizar el
tiempo de rea ión del ve tor. La onstru
ión de ada elemento de
p
involu ra un doble
i lo que permita re orrer la matriz, para a elerar el re orrido de la matriz, el i lo externo que genera a ada elemento de es dividido entre el número de
p
es la región a paralelizar y ada i lo interno al mismo
threads
Más expli itamente, se tiene un
lanzados al ini io de la eje u ión.
thread
maestro que eje uta la apli a ión, y uando se
en uentra on la región paralela, se generan un onjunto de
threads
el problema de forma más e iente, es lo que se llama el modelo región paralela el de los
threads
thread
que permiten resolver
fork-join,
al llegar a la
maestro se bifur a y al nalizar la región se realiza una re ole
ión
lanzados, lo ual se puede mostrar en la gura 14.
Estas regiones paralelas son levantadas por el omputador a través de las dire tivas de
ompilador de la bibliote a OpenMP, realizando una paraleliza ión in remental, omo se
76
Figura 14: Modelo Fork-Join.
muestra en la gura 15, de modo tal de determinar uales son los i los ríti os en esta parte de la apli a ión.
Figura 15: Modelo Paraleliza ión In remental.
A ontinua ión se muestra un segmento del algoritmo que permite tridiagonalizar una matriz simétri a, en el mismo se puede notar el uso de la bibliote a UCSparseLib mediante el llamado de rutinas bási as y para la le tura y es ritura en la matriz en el formato de alma enamiento. Así mismo, en este trozo de ódigo se puede observar omo iertos i los son paralelizados on la ayuda de dire tivas de prepro esador de la bibliote a OpenMP.
for( ii = nbrows-1; ii >0; ii- - ) { ll = ii - 1; hh = 0.0; s ale = 0.0; dSet( nbrows, wai, 0.0 ); For_TDMatrix_Row( AA, ii, row, ACCESS_READ ) { nzi = row.nz;
77
for( nn = 0; nn 0 ) { for( kk = 0; kk
Lihat lebih banyak...
Comentarios
