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
Soustava lineárních rovnic
     Determinant matice
     Inverze matic
     Řešení soustavy lineárních rovnic

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


Obálka Obsah Index Hlavní stránka

změněno 1. srpna 2001
Copyright © 2000-2001 Vladimír Zábrodský, RNDr.