QR-Zerlegung mit dem Gram-Schmidt-Verfahren

module mysubs

contains

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

    integer, parameter  :: fileid = 1
    character(len=256)  :: fileName
    integer             :: status, i

    write (*,*) 'Gram-Schmidtsches Orthogonalisierungsverfahren'
    write (*,*)
    write (*,*) 'Welche Datei enthaelt die Matrix?&
                 & (In grams.dat liegt ein Beispiel.)'
    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) m, n
    if (status /= 0) then
       write (*,*) 'Konnte Dimension nicht lesen.'
       stop
    end if

    allocate(a(n, m))
    allocate(r(n, m))

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

    close(unit=fileid)
    return
  end subroutine init


  subroutine result(a, r)
    use kinds
    implicit none
    real(kind=REAL8), dimension(:,:), pointer :: a, r

    integer   :: n, m, i
    character(len=80) :: format

    n = size(a, 1)
    m = size(a, 2)
    write (format, '(A,I4,A)')   '(', n, '(F12.3, " "))'

    write (*,*) 'Man erhaelt die Matrizen'
    write (*,*)
    write (*, format) (a(i,1:m), i=1,n)
    write (*,*)
    write (*, format) (r(i,1:m), i=1,n)

    return
  end subroutine result

end module mysubs


program GramS
  ! Gram-Schmidtsches Orthogonalisierungsverfahren}

  use kinds
  use mysubs
  implicit none

  real(kind=REAL8), dimension(:,:), pointer :: a, r
  real(kind=REAL8)  :: sum
  integer           :: n, m
  integer           :: i, j, k

  call init(a, r, n, m)

  do k = 1, m
     sum = 0
     do i = 1, n
        sum = sum + a(i, k)*a(i, k)
     end do
     r(k, k) = sqrt(sum)
     do i = 1, n
        a(i, k) = a(i, k) / r(k, k)
     end do
     do j = k + 1, m
        sum = 0
        do i = 1, n
           sum = sum + a(i, k) * a(i, j)
        end do
        r(k, j) = sum
        do i = 1, n
           a(i, j) = a(i, j) - a(i, k) * r(k, j)
        end do
     end do
  end do
  
  call result(a, r)
  stop
end program GramS