! File:  sor_d2.f90
!
! Creation Date: 09.11.2015; Time: 19:11
! Time-stamp: <2015-11-09 20:55:26 (andrzejbaran)> 
!
! Andrzej Baran, IF UMCS
! baran@kft.umcs.lublin.pl
!
SUBROUTINE sor_d2(f,u,n,omega,iter) 
  IMPLICIT none
  INTEGER, PARAMETER :: maxit = 1000
  DOUBLE PRECISION, PARAMETER :: epsilon=1.d-5, zero=0.d0
  DOUBLE PRECISION, DIMENSION(n,n) :: u, f
  DOUBLE PRECISION :: normf, norm, omega, resid, om 
  INTEGER j, l, n, mit, nm, iter
  !
  ! Iterative solution of D^2 u(x,y) = F(x,y)
  ! SOR Chebyshev  method.
  ! Start at u(2:n-1,2:n-1) = 0;
  !
  ! f(n,n) - source term multiplied by h**2, f = h*h*F
  ! u(n,n) - solution
  ! omega  - overrelaxation parameter
  !
  nm = n - 1
  normf=zero 
  do j=2,nm 
     do l=2,nm 
        normf=normf+abs(f(j,l)) 
     end do
  end do
  om = omega
  do mit=1,maxit
     norm=zero 
     do j = 2,nm 
        do l = 2,nm 
           if(mod(j+l,2) == mod(mit,2)) then 
              resid = u(j+1,l) + u(j-1,l) + u(j,l+1)+ &
                   u(j,l-1) - 4d0*u(j,l) - f(j,l) 
              norm = norm + abs(resid) 
              u(j,l) = u(j,l) + om*resid/4d0 
           endif
        end do
     end do
     if(mit == 1) then 
        om = .25*omega
     else 
        om = .5*omega
     endif
     if( (mit > 1) .and. (norm <= epsilon*normf) ) then
        iter = mit
        return
     end if
  end do
  write(*,*) 'No convergence: maxit exceeded.'
END SUBROUTINE sor_d2
