Pivotisierung

Dieses Programm ist sowohl für das Gauß-Verfahren als auch für das Gauß-Jordan-Verfahren einsetzbar.

module mysubs

contains

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

    integer, parameter  :: fileid = 1
    character(len=256)  :: fileName
    integer             :: status
    integer             :: i
    
    write (*,*) 'Darstellung der Pivotisierung der Matrix A und des Vektors c'
    write (*,*) 'des Gleichungssystems Ax = c'
    write (*,*)
    write (*,*) 'In welcher Datei liegen die Daten (Beispiel in pivot.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))

    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(a, c)
    use kinds
    implicit none
    real(kind=REAL8), dimension(:,:), pointer :: a
    real(kind=REAL8), dimension(:), pointer :: c
    
    integer           :: i, n
    character(len=80) :: format
    
    n = size(c, 1)
    write (format, '(A,I4,A)')   '(', n, 'F12.3, A, F12.3)'
    write (*, format) (a(i,1:n), ':', c(i),  i=1,n)

    return
  end subroutine result

end module mysubs


program Pivot
  ! Pivotisierung der Matrix A und des Vektors c des Gleichungssystems
  !     Ax = c

  use kinds
  use mysubs
  implicit none

  real(kind=REAL8), dimension(:,:), pointer :: a
  real(kind=REAL8), dimension(:), pointer   :: c
  integer          :: piv
  real(kind=REAL8) :: maxa, dummy
  integer          :: i, j, k, n

  call init(a, c, n)

  ! piv = maxloc(abs(a(k:n, k)))
  do k=1, n
     piv = k
     maxa  = abs(a(k,k))
     do i=k+1, n
        dummy = abs(a(i, k))
        if (dummy > maxa) then
           maxa = dummy
           piv = i
        end if
     end do

     if (piv /= k) then
        do j=k, n
           dummy       = a(piv, j)
           a(piv, j) = a(k, j)
           a(k, j)     = dummy
        end do
        dummy    = c(piv)
        c(piv) = c(k)
        c(k)     = dummy
     end if
  end do

  call result(a, c)
  stop
end program Pivot