C   prog1a.f -- Solve Poisson Equation
C   Written Feb 25, 2003 by Eric Olson for Mathematics 484
C
C   For information about the five-point formula for solving the Poisson
C   equation see Chapter 7 in A First Course in the Numerical Analysis
C   of Differential Equations by Iserles.  For information about the
C   Gauss-Seidel iteration method used to solve the resulting system of 
C   finite differences see Chapter 10 in the same reference.

      program main

      implicit none

      include "param1.h"
C   param1.h contains the following lines uncommented
C     integer N
C     parameter(N=10)

      real u,hhf,h
      real eps,etol
      integer ipts_i,ipts_j,ipts_n,cmax
      common /soln/u(N,N),hhf(N,N),h
      common /tolr/eps,etol
      common /ipts/ipts_i(2*N),ipts_j(2*N),ipts_n(2*N),cmax

      h=8.0/N
      eps=0.001
      etol=0.001/N/N/N

      call doinit
      call gauss_seidel
      call printu
      stop
      end

C   Determine if a point (x,y) is an interior point in the domain.
C   Returns .TRUE. if an interior point .FALSE. otherwise
      logical function isint(x,y)
      implicit none
      real eps,etol
      common /tolr/eps,etol
      real x,y
      if( (x.gt.0.0+eps).and.
     *    (y.gt.0.0+eps).and.
     *    (x+y.lt.7.0-eps)) then 
        isint=.TRUE.
      else
        isint=.FALSE.
      endif
      end

C   Specify the Dirichlet condition at (x,y) on the boundary
      real function psi(x,y)
      implicit none
      real x,y
      psi=x-2.0*y
      end

C   Specify the forcing term at (x,y) in the domain
      real function f(x,y)
      implicit none
      real x,y
      f=0.0
      end

C  Print the matrix u
      subroutine printu
      implicit none
      include "param1.h"
      real u,hhf,h
      common /soln/u(N,N),hhf(N,N),h

      integer i,j;
      do i=1,N
         write(*,10) (u(i,j),j=1,N)
10       format(1X,<N>F6.1)
      enddo
      end


      subroutine doinit
      implicit none
      include "param1.h"
      real u,hhf,h
      integer ipts_i,ipts_j,ipts_n,cmax
      logical isint
      real psi,f
      common /soln/u(N,N),hhf(N,N),h
      common /ipts/ipts_i(2*N),ipts_j(2*N),ipts_n(2*N),cmax

      integer i,j,c,jmax
C  Locate the domain
      c=1
      do i=1,N
        ipts_n(c)=0;
        do j=1,N
          if(isint((i-1)*h,(j-1)*h)) then
                u(i,j)=1.0
                if(ipts_n(c).eq.0) then
                    ipts_i(c)=i
                    ipts_j(c)=j
                endif
                ipts_n(c)=ipts_n(c)+1
           else
                u(i,j)=0.0
                if(ipts_n(c).gt.0) then
                    c=c+1
                    ipts_n(c)=0
                endif
           endif
        enddo
      enddo
      cmax=c-1

C  Initialize the exterior points
      do c=1,cmax
        jmax=ipts_j(c)+ipts_n(c)-1
        i=ipts_i(c)
        do j=ipts_j(c),jmax
            if(u(i+1,j).eq.0.0) u(i+1,j)=psi(i*h,(j-1)*h)
            if(u(i-1,j).eq.0.0) u(i-1,j)=psi((i-2)*h,(j-1)*h)
            if(u(i,j+1).eq.0.0) u(i,j+1)=psi((i-1)*h,j*h)
            if(u(i,j-1).eq.0.0) u(i,j-1)=psi((i-1)*h,(j-2)*h)
        enddo
      enddo

C  Initialize the interior and forcing points
      do c=1,cmax
        jmax=ipts_j(c)+ipts_n(c)-1
        i=ipts_i(c)
        do j=ipts_j(c),jmax
            u(i,j)=0.0
            hhf(i,j)=h*h*f((i-1)*h,(j-1)*h)
        enddo
      enddo
      end


C  Solve the problem
      subroutine gauss_seidel
      implicit none
      include "param1.h"
      real u,hhf,h
      real eps,etol
      integer ipts_i,ipts_j,ipts_n,cmax
      common /soln/u(N,N),hhf(N,N),h
      common /tolr/eps,etol
      common /ipts/ipts_i(2*N),ipts_j(2*N),ipts_n(2*N),cmax
      integer k,i,j,c,jmax
      real uerr,unew,udif
      do k=1,100000
        uerr=0.0;
        do c=1,cmax
            jmax=ipts_j(c)+ipts_n(c)-1
            i=ipts_i(c);
            do j=ipts_j(c),jmax
                unew=(u(i,j-1)+u(i,j+1)+u(i+1,j)
     *              +u(i-1,j)-hhf(i,j))/4.0
                udif=unew-u(i,j);
                u(i,j)=unew
                uerr=uerr+udif*udif
            enddo
        enddo
        if(uerr.le.etol*etol) goto 20
      enddo
20    continue
      write(*,30) k 
30    format(' Convergence after',i7,' iterations')
      end
