LU - ROZKLAD
PROBLÉM
Zapsat čtvercovou matici A. řádu N jako součin matice L. (lower triangular, tj. dolní trojúhelníková matice) a matice U. (upper triangular, tj. horní trojúhelníková matice).
ALGORITMUS
Následující algoritmus (Press a kol. jej přisuzují Croutovi) je rozkladem v místě (in situ). Předpokládá se, že diagonální prvky matice L. jsou rovny 1: L.J.J=1 for J=1,...,N. Ostatní prvky matice L. budou uloženy pod hlavní diagonálou matice A.; Prvky matice U. budou uloženy v A. na hlavní diagonále a nad ní. Algoritmus neprovádí permutaci řádků. Zaznamenává ji však ve vektoru Indx., který je jedním z výsledků výpočtu. LU rozklad provedený funkcí LUDCMP vyžaduje okolo (1/3)* N** 3 operací násobení a stejně tolik operací sečítání.
IMPLEMENTACE
Jednotka: vnitřní funkce
Globální proměnné: vstupní čtvercová matice A. řádu N, výstupní vektor Indx.
Parametr: přirozené číslo N
Vrací: +1 nebo -1 v závislosti na tom, zda počet výměn řádků byl sudý nebo lichý (využívá se při výpočtu determinantu)
Výsledek: Indx. zaznamená permutaci řádků; LU rozklad matice A. v místě
Poznámka:
W. H. Press a kol. ve své knize uvádí: "Je-li v daném kroku výpočtu hlavní prvek nulový, matice je singulární. Pro práci se singulárními maticemi může být vhodné nahradit nulovou hodnotu prvku hodnotou Tiny" (doporučují volit Tiny=1E-20).
LUDCMP: procedure expose A. Indx.
parse arg N
D = 1; Tiny = 1E-20
do I = 1 to N
Max = 0
do J = 1 to N
W = ABS(A.I.J)
if W > Max then Max = W
end
if Max = 0 then
call ERROR "LUDCMP: Chyba - singulární matice"
Vv.I = 1 / Max
end
do J = 1 to N
do I = 1 to J - 1
Sum = A.I.J
do K = 1 to I - 1
Sum = Sum - A.I.K * A.K.J
end
A.I.J = Sum
end
Max = 0
do I = J to N
Sum = A.I.J
do K = 1 to J - 1
Sum = Sum - A.I.K * A.K.J
end
A.I.J = Sum
W = Vv.I * ABS(Sum)
if W >= Max then do; Max = W; Imax = I; end
end
if J <> Imax then do
do K = 1 to N
W = A.Imax.K; A.Imax.K = A.J.K; A.J.K = W
end
D = -D; Vv.Imax = Vv.J
end
Indx.J = Imax
if A.J.J = 0 then A.J.J = Tiny
if J <> N then do
W = 1 / A.J.J
do I = J + 1 to N
A.I.J = A.I.J * W
end
end
end
return D
ERROR: say ARG(1); exit
|
SOUVISLOSTI
Literatura
Press W.H., Teukolsky S.A., Vetterling W.T., Flannery B.P. Numerical Recipes in C : the art of scientific computing
- 2nd ed. University Press, Cambridge, 1992
