module mysubs

contains

  subroutine getsrc(q)
    use kinds
    implicit none
    
    real(kind=REAL8), dimension(:,:)  :: q
    integer  :: n
    integer  :: d1, d2
    integer  :: x1, x2, y1, y2

    !
    ! elektronische Linse:
    !
    !                             |         |
    !                   U = 1  -> |   D2    | <-  U = -1
    !          Y2  ---------------|         |----------------
    !                                                      D1
    !
    !          Y1  ---------------|         |----------------
    !                   U = -1 -> |         | <-  U =  1
    !                             |         |  
    !             0               X1        X2
    !

    n  = size(q, 1)
    d1 = n/2
    d2 = n/4
    
    q = 0.0
    
    x1 = (n - d2)/2     
    x2 = (n + d2)/2
    y1 = (n - d1)/2     
    y2 = (n + d1)/2

    q(1:x1, y1) = -1.0
    q(1:x1, y2) =  1.0
    
    q(x2:n, y1) =  1.0
    q(x2:n, y2) = -1.0
      
    q(x1, 1:y1) = -1.0
    q(x2, 1:y1) =  1.0
    
    q(x1, y2:n) =  1.0
    q(x2, y2:n) = -1.0

    return
  end subroutine getsrc

  subroutine store(u)
    use kinds
    implicit none
    
    real(kind=REAL8), dimension(:,:)  :: u
    character(len=80), parameter      :: ofile = 'poisson.out'
    integer                           :: n, status, i, j

    n  = size(u, 1)
    open(1, file=ofile,iostat=status)
    if (status /= 0) then
       write (*,*) 'Datei ', ofile, ' konnte nicht geoeffnet werden !'
       stop
    end if
    write (1, *) n-2
    write (1,'(F10.5)') ((u(i,j), j=2, n-1), i=2, n-1)

    close(1)
    return
  end subroutine store

  function sweep(a, f, b)
    use kinds
    implicit none

    real(kind=REAL8)                  :: sweep
    real(kind=REAL8), dimension(:,:)  :: a, b, f
    real(kind=REAL8)                  :: resid
    integer                           :: n, i, j

    n  = size(a, 1)
    sweep = 0.0

    !$OMP DO
    do i = 2, n-1
       do j = 2, n-1
          resid = a(i-1,j) + a(i+1,j) + a(i,j+1) + a(i,j-1) + &
                 (-4.0)*a(i,j) - f(i,j)
          sweep = sweep + abs(resid)
          b(i,j) = a(i,j) + 0.25*resid
       end do
    end do
    !$OMP END DO

    return
  end function sweep

end module mysubs


program poisson3
  ! Loesen der Poisson-Gleichung mit Jacobi-Relaxation

  use kinds
  use mysubs
  implicit none

  integer, parameter          :: n      = 200
  integer, parameter          :: maxits = 1000
  real(kind=REAL8), parameter :: eps    = 1.0e-4

  real(kind=REAL8), dimension(n,n) :: q, u1, u2
  real(kind=REAL8)                 :: error, total_error
  integer                          :: it

  call getsrc(q)        ! Quellen-Matrix einlesen
  u1 = 0.0              ! Start-Loesung setzen
  u2 = 0.0

  !$OMP PARALLEL PRIVATE(error)
  do it = 1, maxits
     error = sweep(u1, q, u2)
     error = sweep(u2, q, u1)

     ! sum local errors
     !$OMP SINGLE
     total_error = 0.0
     !$OMP END SINGLE

     !$OMP CRITICAL
     total_error = total_error + error
     !$OMP END CRITICAL
     !$OMP BARRIER

     if (total_error/(n*n) < eps) exit
  end do

  !$OMP SINGLE
  print *, 'Fehler nach ', it, ' Iterationen: ', total_error/(n*n)
  !$OMP END SINGLE

  !$OMP END PARALLEL

  call store(u1)          ! Loesung ausgeben

  stop
end program poisson3