Eliminacion Gaussiana.

UNI-CTIC

1

Eliminacion Gaussiana Mendoza Olivares, Jesus Ticse Torres, Royer Lima, 25 de Agosto del 2009

Resumen—En el siguiente trabajo, se examina uno de los algoritmos m´as conocidos para resolver un sistema lineal de ecuaciones el algoritmo de eliminaci´on gaussiana en su versi´on paralela.

I.

´ I NTRODUCCI ON

Este m´etodo trabaja directamente sobre la matriz aumentada para resolver un sistema de ecuaciones lineales, llevandola a la matriz de un sistema triangular que es eqivalente al sistema inicial. La equivalencia del sistema triangular final con el inicial se argumenta debido a que el algoritmo s´olo utilizamos operaciones elementloes entre las filas y cuya aplicaci´on individual siempre preserva la equivalencia. Implementamos el algoritmo de eliminaci´on Gaussiana en forma paralela y se realizar´a una revisi´on sobre el trabajo computacional realizado por este m´etodo. II.

O BJETIVOS

Entender el algoritmo de Eliminaci´on Gassuina para solucionar un sistema de ecuaciones. Implementar la versi´on paralela del algoritmo Eliminaci´on Gaussiana. III.

P LANTEAMIENTO DEL PROBLEMA

La principal caracteristica de la Eliminaci´on Gaussiana es transformar el sistema original AX = B en un sistema lineal equivalente UX = Y donde U es una matriz triangular superior. La raz´on de lograr e´ sta transformaci´on del sistema es que el nuevo sistema es m´as simple y resolverlo s´olo requiere un proceso hacia atr´as. Por ejemplo, tomemos A una matriz de dimenci´on n = 4.   a11 a12 a13 a14  a21 a22 a23 a24     a31 a32 a33 a34  a41 a42 a43 a44 Las opreraciones disponibles son multiplicar un vector fila por un escalar y a˜nadirlo a otra fila. 1. Nos intereza que las posiciones 21, 31, 41 de la matriz sean ceros, esto lo podemos lograr sustrayendo un multiplo adecuado de la primera fila de A. f ila2 ← f ila2 − m21 ∗ f ila1 donde:m21 = a21 /a11 . f ila3 ← f ila3 − m31 ∗ f ila1 donde:m31 = a31 /a11 .

f ila4 ← f ila4 − m41 ∗ f ila1 donde:m41 = a41 /a11 . y entonces nuestra matriz A la hemos transformado en:   a11 a12 a13 a14  0 u22 u23 u24     0 u32 u33 u34  0 u42 u43 u44 2. Ahora, enfoquemos nuestra atenci´on al subproblema de una matriz de 3x3   u22 u23 u24  u32 u33 u34  u42 u43 u44 En forma similar, hacemos ceros los elementos por debajo de u22 , para lograrlo operemos f ila3 ← f ila3 − m32 ∗ f ila2 donde:m32 = u32 /u22 . f ila4 ← f ila4 − m42 ∗ f ila2 donde:m42 = u42 /u22 . 3. Luego, nuestra matriz  a11  0   0 0

se ha trasformado en:  a12 a13 a14 u22 u23 u24   0 v33 v34  0 v43 v44

ahora solo nos resta operar con una fila para reducir v43 a cero. f ila4 ← f ila4 − m43 ∗ f ila3 donde:m43 = v43 /v33 . Lo que nos da lugar a la matriz triangular superior:   a11 a12 a13 a14  0 u22 u23 u24     0 0 v33 v34  0 0 0 w44 Contando con el  a11  0   0 0

sistema: a12 u22 0 0

a13 u23 v33 0

    a14 x1 y1 x2  y2  u24    =   v34  x3  y3  w44 x4 y4

Introduzcamos la idea de matriz aumentada: [AkB] y sobre ella, realicemos las operaciones indicadas en los pasos 1, 2, ..., n − 1. Al hacerlo, estaremos trabajando sobre la parte activa de A pero de igual forma sobre las filas correspondientes a B. El resultado ser´a [U kY ] donde Y es el resultado de los cambios ocurridos a B al aplicarle las operaciones fila que hemos considerado.

UNI-CTIC

2

IV.

A LGORITMO

El metodo de eliminaci´on de gauss es un m´etodo estandar para resolver el sistema de ecuaciones Ax = b. Se comienza por transformar el sistema a la forma equivalente U x = b, donde U es una matriz triangular superior (es decir todos los elemntos debajo de la diagonal son iguales a cero). El sistema U x = b es ahora resuelto por sustituci´on hacia atraz. Dada una matriz A = (aij ), un algor´ıtmo para llevar a la practica el metodo de eliminaci´on de gaussiana que acabamos de escribir es el siguiente. Algorithm 1 input n, (aij ) for k = 1, 2, ..., n − 1 do for i = k + 1, k + 2, ..., n do z ← aik /akk aik ← 0 for j = k + 1, k + 2, ..., n do aij ← aij − zakj end for end for end for output(aij ) El algoritmo que se aplico para el caso del problema esta paralelisado se muestra acontinuaci´on. 1. Primero reservamos memoria en base a la variable definida DIM EN SION , que nos indica la dimension de la matriz A que vamos triangularizar. 2. Si llamamos p el numero de procesos que ejecutaran el programa y myRank la variable que representa el id de cada proceso, si myRank = 0 entonces genera y llena los valore de la matriz de forma aleatoria, cumpliendo las siguientes restricciones: Genera valores de A que sean diferentes de cero en la diagonal. Genera valores de A que sean diferentes de cero en las 6 sub-diagonales arriba y las 6 abajo, el resto es cero. Los valores en la diagonal deben ser mayores en valor absoluto que los de las otras 12 subdiagonales. 3. El modelo de programacion paralela que se aplica esta basado en la particion de data en este caso de las columnas de la matriz: La columna 0 los va ejecutar el proceso con myRank = 0. La columna 1 los va ejecutar el proceso con myRank = 1 Y asi sucesivamente, siempre que la dimension de la matriz sea divisible por el numero de procesos (esto para que todos los procesos trabajen igual). 4. Acontinuaci´on cada proceso trabaja sobre su correspondiente columna: Estos son todos los elementos en la k−esima columna debajo de la diagonal

for i = k + 1, k + 2, ..., DIM EN SION − 1 do Aik Aik = −1,0. A kk buf f er(i − k − 1) = Aik end for Donde buf f er es una variable para la transmision de data entre procesos. 5. Difundimos todos los datos desde proceso root a todos los procesos. 6. El proceso 0 ejecuta lo siguiente : for i = 0, 1, ..., DIM EN SION do if i mod p <> 0 then recibe la data buf f er de los procesos con id (i mod p). end if for j = 0, 1, ..., DIM EN SION do Aji = buf f er end for end for for i = 0, 1, ..., DIM EN SION do for i = 0, 1, ..., DIM EN SION do Aij = 0,0 end for end for 7. Los demas procesos con Id distinto de cero: for i = 0, 1, ..., DIM EN SION do if myRank = (i mod p) then for i = 0, 1, ..., DIM EN SION do Enviamos buf f er al proceso 0. end for end if end for Este algortimo presenta los mismos bucles que se muestran tambien en el algoritmo desarrollado en serie, con la unica diferencia con este ultimo algoritmo cada columna es computada por un proceso y estos se comunican todos los pasos de actualizacion que puedan realizar en las entradas de la matriz.

UNI-CTIC

3

V.

P RUEBA

Utilizamos una matriz de 4x4 para probar el programa. matriz aleatoria de entrada:   8 1 2 0  −2 10 1 −3     4 −4 7 2  −5 4 −2 11 Llevandolo  8 1 −2 10   4 −4 −5 4 → 8 1 0 41/4  0 0 0 0 → 

8 0  0 0

1 10,25 0 0

a su forma diagonal   2 0 8 0 1 −3 → 0 7 2 −2 11 0

 1 2 0 41/4 3/2 −3 → −9/2 6 2 37/8 −6/8 11 

→

Obtenido por el programa eliminacion gaussiana.c obtenemos:   8 1 2 0   0 10,25 1,5 −3    0 0 6,658537 0,682927  0 0 0 12,5 Para este caso nos resulta la misma matriz calculada, con error cero. R ESULTADOS

Luego, hacemos las pruebas con dos matrices aleatorias A de tama˜no 120x120 y 240x240 cuyos u´ nicos elementos distintos de cero se encuentran en la diagonal principal y en las subdiagonales arriba y las 6 debajo de la diagonal principal. Los valores en la diagonal son mayores en valor absoluto que los de las otras 12 sub-diagonales Matriz aleatoria de 120x120 N´umero de procesadores = 1 Prueba 1 2 3

Tiempo(seg) 0.294493 0.374795 0.304104

N´umero de procesadores = 2 Prueba Tiempo(seg) 1 0.419840 2 0.359263 3 0.308416

Prueba 1 2 3

Tiempo(seg) 7.035414 6.640708 6.566449

N´umero de procesadores = 2

triangular :

2 0 3/2 −3   273/41 28/41  117/82 −1013/82   8 1 2 0 0 41/4 3/2 −3    0 0 273/41 28/41 0 0 25/2 0 2 0 1,5 −3   6,6585 0,6829 0 12,5

VI.

Matriz aleatoria de 240x240 N´umero de procesadores = 1

Prueba 1 2 3

Tiempo(seg) 6.554764 6.51030 6.277527 VII.

C ONCLUSIONES

En este trabajo hemos analizado el algoritmo Eliminaci´on Gaussianaen su version en paralela, aprovechando los procesadores m´as eficientemente.. Con el desarrollo del programa hecho con la herramienta MPI se pudo hallar valores aproximados del tiempo de procesado para hallar la matriz triangular superior y podemos notar que el tiempo de c´omputo no mejora considerablemente al usar 2 procesadores. Probablemente este tiempo se reduce al usar m´as procesadores. Es evidente la importancia del m´etodo de Eliminaci´on Gaussiana para resolver sistema de ecuaciones lineales, ya que transforma en un sistema m´as simple de resolver. Se puede mejorar el algoritmo aprovechando las caracteristicas de una matriz particular para economizar operaciones y espacio. R EFERENCIAS [1] [MPI: A Message Passing Interface Standar ] http://wwwunix.mcs.anl.gov/mpi/mpi-standard/mpi-report-1.1/mpi-report.html [2] [Blas Barney, Lawrence Livermore National Laboratory :MPI ] https://computing.llnl.gov/tutorials/mpi/ [3] [Wikipedia] http://es.wikipedia.org/wiki/Eliminaci´on de Gauss-Jordan

UNI-CTIC

4

´ A NEXO :C ODIGO UTILIZADO eliminaci´on gaussiana.c

# include # include # include # include # include

< s t d i o . h> < s t d l i b . h> <math . h> <mpi . h>

# d e f i n e DIMENSION 10 # d e f i n e NUMDIAG 4 i n t main ( i n t a r g c , c h a r ∗ a r g v [ ] ) { MPI Status s t a t u s ; i n t p , my rank , i , j , k , d e s t ; d o u b l e ∗∗A, ∗∗ A s e r ; double ∗ b u f f e r ; double l ; double s t a r t p a r , end par , s t a r t s e r , end ser , s e r i a l t i m e , p a r a l l e l t i m e , speedup ; double diff , error ; M P I I n i t (& a r g c , &a r g v ) ; MPI Comm size (MPI COMM WORLD, &p ) ; MPI Comm rank (MPI COMM WORLD, &my rank ) ; i f ( DIMENSION % p != 0 ) { i f ( my rank == 0 ) { p r i n t f ( ” ( %d ) no s e p u e d e d i v i d i r en un numero de p r o c e s o s ( %d ) \ n ” ,DIMENSION , p ) ; } fflush ( stdout ); fflush ( stderr ); M P I B a r r i e r (MPI COMM WORLD ) ; MPI Finalize ( ) ; exit (1); } / / r e s e r v a n d o memoria A = ( d o u b l e ∗ ∗ ) m a l l o c ( DIMENSION∗ s i z e o f ( d o u b l e ∗ ) ) ; f o r ( i = 0 ; i = i − NUMDIAG && j <= i + NUMDIAG) { i f ( i == j ) { A[ i ] [ j ] = ( d o u b l e ) ( rand ( ) %1 0 0 0 0 ) + 5000 ; }else{ A[ i ] [ j ] = ( d o u b l e ) ( rand ( ) %1 0 0 0 0 ) − 5000 ; } }else{ A[ i ] [ j ] = } } }

0.0;

f o r ( i = 0 ; i
UNI-CTIC

i f ( my rank == 0 ) { s t a r t s e r = MPI Wtime ( ) ; f o r ( k = 0 ; k
5