LU Rozklad

Pole jsou neuvěřitelně silný programovací prostředek. Umožnují jako jediné objekty provádět rozsáhle operace nad jejími prvky najednou a tak umožňují využití počítačů naplno. Je to samozřejmě společně s jiným velmi silným prostředkem cykly. Fortran ovšem umožňuje zacházet s poli mnohem elegantněj.

LU rozklad je určený na řešení soustav lineárních rovnic. Na rozdíl od Gaussovy eliminace je ovšem přímo vyvinuta na použití maticových operací. Uvedeme jeho modul Croutuv rozklad.

Zajimavé je použití proměnné ix. Ta není viditelná pro vnejší svět mimo modulu. Slouží na uložení mezivýsledků indexů při přehazování řádků během pivotování. Pro většinu použití není potřeba a proto je viditelná pouze uvnitř modulu. Pokud je z nějakých důvodů potřeba jej přece jen použít pak je na její získání třeba uvedena subroutina getix. Důležité je, správně alokovat a dealokovat pamět tak aby v průběhu běhu programu nezůstávaly v paměti její nepoužitelné zbytečné kusy.

Module Crout

  integer, allocatable, private :: ix(:)

contains 

subroutine lu (n,a,d)

  implicit none
  integer :: n, i, k, k1
  real :: d,a(:,:),s(n),piv

  deallocate(ix)
  allocate(ix(n))
  d = 0.0
  do k = 1, n 
    piv = 0.0

! Vypocet L        
    do i = k, n
      a(i,k) = a(i,k) - sum(a(i,1:k-1)*a(1:k-1,k))
      if( abs( a(i,k) ) > piv )then
         k1 = i 
         piv = abs( a(i,k) )
      endif
    enddo
! Singularni matice
    if( piv == 0.0 ) return
! Vymena radku
    if( k /= k1 )then 
       s = a(k,:)
       a(k,:) = a(k1,:)
       a(k1,:) = s
    endif
    ix(k) = k1
    
! Vypocet U
    do j = k + 1, n 
        a(k,j) = (a(k,j) - sum(a(k,1:k-1)*a(1:k-1,j)))/a(k,k)
    enddo
enddo

d = 1.0
end subroutine

subroutine lusol(n,a,b)

Integer :: n, i, k
real :: a(:,:),b(:),s

if( .not. allocated(ix) ) stop 'Matrix not factorized.'

! prehozeni indexu
do i = 1, n - 1
   if( ix(i) /= i )then
      s = b(i) 
      b(i) = b(ix(i)) 
      b(ix(i)) = s 
   endif
enddo
! zpetna substitice do spodni matice
do i = 1, n 
  b(i) = (b(i) - sum(a(i,1:i-1)*b(1:i-1)) / a(i,i)
enddo
do i = n - 1, 1, -1 
  b(i) = b(i) - sum(a(i,i+1:n)*b(i+1:n))
enddo
end subroutine

subroutine getix(n,ixx)

real :: ixx(:)

ixx(1:n) = ix(1:n)
end subroutine

! jaka bude inverze?

end module