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