program sor
implicit none
integer,parameter :: nx = 65, ny = 65
real*8,dimension(nx,ny) :: u,rhs,res,uxct
real*8 :: l2,h,error
real*8,parameter :: eps = 0.00000001d0

integer i,j,l

interface
real*8 function relaxc(nx, ny,u,rhs,res,h)
implicit none
integer :: nx,ny
real*8,dimension(nx,ny) :: u, rhs,res
real*8 :: h
integer i,j
real*8 h2,jacel,sum
end function

end interface


h = 1.0d0 / (nx - 1 )

	call set_up_u(nx,ny,u,uxct,h)
   call set_up_rhs(nx,ny,rhs,uxct)

	do i = 1 , nx
		u(i,1) = uxct(i,1)
		u(i,ny) = uxct(i,ny)
	end do
	do j = 1, ny
		u(1,j) = uxct(1,j)
		u(nx,j) = uxct(nx,j)
	end do

   l2 = relaxc(nx,ny,u,rhs,res,h)
   do while (l2 > eps) 
         l  = l + 1 
			if (l.eq.4944) exit
         l2 = relaxc(nx,ny,u,rhs,res,h)
         if (mod(l,100).eq.0) then
	         write(6,*) 'Iteration ' , l , 'norm = ',l2
				error = 0.0d0
				do i = 1, nx
				do j = 1, ny
					error = error + abs(u(i,j)-uxct(i,j))
				end do
				end do
				error = error/ nx / ny
				write(6,*) 'Error ',error
			end if
	end do
	write(6,*) 'Iteration ' , l , 'norm = ',l2

stop
end

real*8 function relaxc(nx, ny,u,rhs,res,h)
implicit none
integer :: nx,ny
real*8,dimension(nx,ny) :: u, rhs,res
real*8 :: h
integer i,j
real*8 h2,jacel,sum,sum2,sum3

      h2    = 1.0/( h * h)
      jacel = (h * h)/ 4.0
      sum = 0.0d0

		do i = 2 , nx-1
		do j = 2 , ny-1
            res(i,j) = h2 * ( u(i+1,j) + u(i-1,j) + u(i,j+1) + u(i,j-1) - &
                               4.0 * u(i,j) )  - rhs(i,j)

            u(i,j) = u(i,j) + jacel * res(i,j)
            sum = sum + res(i,j) * res(i,j) 
		end do
		end do
		sum2 = sum / nx/nx/ny/ny
		sum3 = sqrt(sum2)
		relaxc = sum3

return
end

subroutine set_up_u(nx, ny, u, uxct,h)
implicit none
integer nx,ny
real*8,dimension(nx,ny) :: u,uxct
real*8 :: h
real*8,parameter :: pi = 3.1415926535897932d0

integer i,j

do j = 1 , ny
do i = 1 , nx
	uxct(i,j) = sin(pi*(i-1)*h)
	u(i,j) = 0.0d0
end do
end do
 
return
end

subroutine set_up_rhs(nx,ny,rhs,uxct)
implicit none
integer :: nx,ny
real*8,dimension(nx,ny) :: rhs,uxct
real*8,parameter :: pi = 3.1415926535897932d0

integer i,j


do j = 1 ,ny
do i = 1 ,nx

	rhs(i,j) = -pi*pi * uxct(i,j)

end do
end do

