Paraleln´ı LU rozklad Luk´aˇs Michalec ´ ı n.L. Katedra fyziky, Pˇr´ırodovˇedeck´a fakulta Univerzity J.E. Purkynˇe v Ust´ roˇcn´ık, specializace
Abstract Semin´ arn´ı pr´ ace se zab´ yv´a ˇreˇsen´ı soustavy line´arn´ıch rovnic pomoc´ı LU dekompozice a zrychlen´ı v´ ypoˇctu pomoc´ı paralelizace t´eto metody.
1
´ Uvod
K jedn´ım z nejstarˇs´ım matematick´ ym probl´em˚ um patˇr´ı ˇreˇsen´ı soustavy line´arn´ıch rovnic. S t´ımto probl´emem se m˚ uˇzeme napˇr´ıklad setkat v line´arn´ım programov´an´ı, pˇredpovˇed´ıch ˇci odhadov´an´ım. ˇ sen´ı tohoto probl´emu je nˇekolik, napˇr´ıklad Gaussova eliminace, Cramerov Reˇ pravidlo, pomoc´ı inverzn´ı matice, metodou nejmenˇs´ıch ˇcterc˚ u a LU rozklad.
1
2
Teorie
Takovou soustavou line´arn´ıch rovnic m˚ uˇze b´ yt napˇr´ıklad soustava: 3x + 5y − z = 2 4x + 1y + 2z = 1 1x − 2y + 1z = 5
(1)
V t´eto soustavˇe m´ame tˇri promnˇen´e (x,y,z). Graficky si tuto soustavu m˚ uˇzeme pˇredstavi jako soustavu tˇr´ı ploch a my hled´ame bod, ve kt´er´em se dan´e plochy prot´ınaj´ı. M´ısto rovnicov´eho z´apisu, m˚ uˇzeme pouˇz´ıt i maticov´ y z´apis: A~x = ~b
(2)
kde A je matice koeficinet˚ u, x je vektor promnˇen´ ych a b je vektor poˇzadavk˚ u. Pokud pˇrep´ıˇseme soustavu 1, dostaneme: 3 5 −1 x 2 4 1 2 × y = 1 (3) 1 −2 1 z 5
2.1
LU rozklad
Metoda LU dekompozice je zaloˇzn´a takov´em rozkladu matice A, aby se zjednoduˇsil v´ ypoˇcet v´ ysledn´eho ˇreˇsen´ı. LU rozklad jde pouˇz´ıt jen na regul´arn´ı ˇctvercov´e ˇreˇsen´ı, to znamen´a, ˇze rovnice maj´ı pr´avˇe jedno ˇreˇsen´ı a to ˇreˇsen´ı existuje. Naˇs´ım u ´kol je rozloˇzit matici A na matice L a U, tak aby: A=L∗U
(4)
kde matice L je doln´ı trojuheln´ıkov´a a matice U je horn´ı trojuheln´ıkov´a. T´ımto rozkladem dostaneme rovnici: L ∗ U ∗ ~x = ~b
(5)
souˇcin matice U a vektoru x si onaˇc´ıme U ∗ ~x = ~z, takˇze dostaneme: L ∗ ~z = ~b 2
(6)
Prvn´ım krokem je potom uˇz dopˇredn´a substituce v rovnici: L ∗ ~z = ~b a spoˇcteme ~z. Dalˇs´ım krokem potom bude zpˇetn´a substituce, dosah´ıme za ~z a spoˇcteme ~x: U ∗ ~x = ~z
3
V´ ypoˇ cet LU rozkladu
Algoritmus rozkladu si uk´aˇzeme na konkr´etn´ım pˇr´ıkladˇe. Mˇejme matici: 1 2 3 A = 2 3 1 (7) 4 2 0 z t´eto matice, potˇrebujeme 1 L = l21 l31
z´ıskat matice L a U: 0 0 u11 u12 u13 1 0 , U = 0 u22 u23 0 0 u33 l32 1
(8)
Pokud si rozep´ıˇseme n´asoben´ı matic L a U, tak aby v´ ysledkem byla matice A, tak dostaneme soustavu jednoduch´ yh rovnic: u11 ∗ 1 = 1 → u11 = 1 u12 ∗ 1 = 2 → u12 = 2 u13 ∗ 1 = 3 → u13 = 3 l21 ∗ u11 = 2 → l21 = 2 .. . coˇz m˚ uˇzeme zapsat algoritmem jako: DO k = 1 , n -1 DO i = k +1 , n a (i , k ) = a (i , k ) / a (k , k ) END DO DO i = k +1 , n a (i , j ) = a (i , j ) - a (i , k ) * a (k , j ) END DO END DO Pro uˇsetˇren´ı m´ısta, zap´ıˇseme matice L a U do matice A. 3
(9)
3.1
Dopˇ redn´ a substituce
LU rozkladem jsme dostali 1 L = 4 2
matice L a U: 0 0 2 −1 7 1 0 , U = 0 2 −23 2.5 1 0 0 41.5
(10)
a d´ale budeme ˇreˇsit: L ∗ ~y = ~b
(11)
Zde opˇet dostaneme soustavu rovnic: y1 = 11 4y1 + y2 = 25 2y1 + 2.5y2 + y3 = 16
(12)
opˇet, ˇreˇsen´ı t´eto soustavy je trivi´aln´ı a pokud to pˇrep´ıˇseme do fortranovsk´eho k´odu: DO i = 2 , n DO j = 1 , i - 1 y ( i ) = y ( i ) - a (i , j ) * y ( j ) ENDDO ENDDO
3.2
Zpˇ etn´ a substituce
Posledn´ım krokem je zpˇetn´a substituce, kde ˇreˇs´ıme rovnici: U ∗ ~x = ~y
(13)
2 −1 7 x1 11 0 2 −23 x2 = −19 0 0 41.5 x3 41.5
(14)
opˇet, kdyˇz jsi rozep´ıˇseme n´asoben´ı matice a vektoru, tak dostaneme trivi´aln´ı soustavu rovnic: x3 = 1 2x2 − 23x3 = −19 2x1 − x2 + 7x3 = 11 4
(15)
Fortranovsk´ y k´od vypad´a n´asledovnˇe: DO i = n , 1 , -1 DO j = i + 1 , n x ( i ) = y ( i ) - a (i , j ) * x ( j ) ENDDO x ( i ) = x ( i ) / a (i , i ) ENDDO A t´ımto jsme z´ıskali ˇreˇsen´ı soustavy line´arn´ıch rovnic.
3.3
Paralelizace metody LU
Paralelizace t´eto metody spoˇc´ıv´a v distribuci cykl˚ u mezi vˇsemy procesory. Takovou uk´azkou distribuce je n´asleduj´ıc´ı pˇr´ıklad: DO i = myrank , n -1 , nproc print * , " procesor " , myrank , " dela " ,i , " . iteraci " ENDDO v´ ystupem potom je: procesor procesor procesor procesor procesor ...
0 1 2 0 1
dela dela dela dela dela
0. 1. 2. 3. 4.
iteraci iteraci iteraci iteraci iteraci
A jelikoˇz vˇsechny tˇri kroky LU dekompozice jsou iteraˇcn´ı, tak lze distribuci cykl˚ u prov´est u vˇsech krok˚ u.
5
3.3.1
Paralelizace rozkladu LU
DO k = 1 , n IF ( map ( k ) == myrank ) THEN DO i = k +1 , n a (i , k ) = a (i , k ) / a (k , k ) ENDDO ENDIF
CALL MPI_BCAST ( a (k , k ) , n - k +1 , MPI_REAL , map ( k ) , MPI_COMM_WORLD , ier DO j = k +1 , n IF ( map ( j ) == myrank ) THEN DO i = k +1 , n a (i , j ) = a (i , j ) - a (i , k ) * a (k , j ) ENDDO ENDIF ENDDO ENDDO Zde byla distribuce cykl˚ u zajiˇstˇena pomoc´ı podm´ınek v cyklu, kv˚ uli synchronizaci. Kaˇzd´ y procesor si nejdˇr´ıve spoˇc´ıt´a vlastn´ı ˇc´ast L matice, pak si ty procesory mezi s sebou ty hodnoty vymˇen´ı a pak pokraˇcuj´ı u ˇc´asti U matice.
3.3.2
Paralelizace dopˇ redn´ e distribuce
DO i = 2 , n s = 0.0 DO j = 1 + myrank , i - 1 , nprocs s = s + a (i , j ) * y ( j ) ENDDO
CALL MPI_ALLREDUCE (s , ss , 1 , MPI_REAL , MPI_SUM , MPI_COMM_WORLD , ier y ( i ) = y ( i ) - ss ENDDO V t´eto ˇca´sti, kdy poˇc´ıt´ame vektor ~y , tak jsou iterace na sobˇe z´avisl´e, protoˇze postupnˇe dosazujeme vypoˇcten´e promnˇen´e do dalˇs´ıch rovnic. Proto jsme zde 6
jen zparalelizovali sumu n´asobk˚ u, kter´a je potˇreba k v´ ypoˇctu. Toto urychlen´ı se projev´ y hlavnˇe u hodnˇe velk´ ych matic.
3.3.3
Paralelizace zpˇ etn´ e distribuce
DO i = n , 1 , -1 s = 0.0 IF ( map ( i +1) <= myrank ) THEN ii = i + 1 + myrank - map ( i +1) ELSE ii = i + 1 + myrank - map ( i +1) + nprocs ENDIF DO j = ii , n , nprocs s = s + a (i , j ) * b ( j ) ENDDO
CALL MPI_ALLREDUCE (s , ss , 1 , MPI_REAL , MPI_SUM , MPI_COMM_WORLD , ier b ( i ) = ( b ( i ) - ss ) / a (i , i ) ENDDO v t´eto ˇca´sti jsme opˇet rozparalelizovali sumu n´asobk˚ u,kter´a je potˇrebn´a k dosazov´an´ı do rovnic. Jen bylo sloˇzitˇejˇs´ı zpoˇcten´ı, co kter´ y proces m´a n´asobit a od jak´eho indexu zaˇc´ıt.
4
V´ ysledky a diskuze
V´ ysledkem m´ame program, kter´ y je schopn´ y pomoc´ı LU dekompozice spoˇc´ıtat ˇreˇsen´ı soustavy line´arn´ıch rovnic.
7
Figure 1: Uk´azka v´ ystupu programu Pokud by jsme zmˇeˇrili rychlost v´ ypoˇctu u velk´ ych matic, zjistili bychom, ˇze jsme uˇsetˇrili v´ıce ˇcasu.
5
Z´ avˇ er
Tato u ´loha n´am uk´azala, jak´ ym lehk´ ym zp˚ usobem se d´a urychlit iteraˇcn´ı v´ ypoˇcet a jak v˚ ubec funguje metoda LU dekompozice.
8
