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