Gaußalgorithmus

Die Vorwärtselimination im Gaußalgorithmus erfordert ca. , die Rücksubstitution ca. Operationen.

module mysubs

contains

  subroutine init(a, c, x, n)
    use kinds
    implicit none
    real(kind=REAL8), dimension(:,:), pointer :: a
    real(kind=REAL8), dimension(:), pointer :: c, x
    integer, intent(out)   :: n

    integer, parameter  :: fileid = 1
    character(len=256)  :: fileName
    integer             :: status
    integer             :: i
    
    write (*,*) 'Gauss-Elimination ohne Pivotisierung'
    write (*,*)
    write (*,*) 'In welcher Datei stehen Matrix und Loesungsvektor?'
    write (*,*) '(Beispiel in gauss.dat)'
    read (*,*) fileName

    open(unit=fileid, file=fileName, action='read', iostat=status)
    if (status /= 0) then
       write (*,*) 'Konnte Datei ', fileName, 'nicht öffnen.'
       stop
    end if

    read (fileid, *, iostat=status) n
    if (status /= 0) then
       write (*,*) 'Konnte Dimension nicht lesen.'
       stop
    end if

    allocate(a(n, n))
    allocate(c(n), x(n))

    read(fileid, *, iostat=status) (a(i,1:n), i=1,n)
    if (status /= 0) then
       write (*,*) 'Konnte Array a nicht lesen.'
       stop
    end if

    read(fileid, *, iostat=status) c(1:n)
    if (status /= 0) then
       write (*,*) 'Konnte Lösungsvektor c nicht lesen.'
       stop
    end if

    close(unit=fileid)
    return
  end subroutine init


  subroutine result(x)
    use kinds
    implicit none
    real(kind=REAL8), dimension(:), pointer :: x
    
    integer           :: i, n
    
    n = size(x, 1)
    write (*,*) 'Das Ergebnis lautet'
    write (*, '(F10.3)') (x(i), i=1,n)

    return
  end subroutine result

end module mysubs


program Gauss
  ! Gauss-Elimintation ohne Pivotisierung

  use kinds
  use mysubs
  implicit none

  real(kind=REAL8), dimension(:,:), pointer :: a
  real(kind=REAL8), dimension(:), pointer :: c, x
  integer             :: n, i, j, k
  real(kind=REAL8)    :: sum, factor
  
  call init(a, c, x, n)

  ! Vorwaertselimination
  do k = 1, n - 1
     do i = k + 1, n
        ! Hier muesste die Pivotisierung kommen
        if (a(k,k) == 0) then
           write (*,*) 'Matrix muss pivotisiert werden.'
           stop
        end if
        factor = a(i, k) / a(k, k)
        do j = k + 1, n
           a(i, j) = a(i, j) - factor * a(k, j)
        end do
        c(i) = c(i) - factor * c(k)
     end do
  end do

  ! Rueckwaertsubstitution
  x(n) = c(n) / a(n, n)
  do i = n - 1, 1, -1
     sum = 0
     do j = i + 1, n
        sum = sum + a(i, j) * x(j)
     end do
     x(i) = (c(i) - sum) / a(i, i)
  end do

  call result(x)
  stop
end program Gauss