c1t/f77
      program laplace
c     This is the Fortran 77 version of the Laplace equation solver
c     using Jacobi algorithm.
c     Y. Bonver, CIS, June 9, 1992
      parameter(max=512)
      double precision values(max,max)
      common/val/values
      common/block/nxpts,nypts,nproc
      integer job, icycle
      integer niter
      common/info/niter
      real total2
      common/time/total2
      double precision tol,tolerance
      common/itol/tol,tolerance
      integer k,maxk
      common/imaxk/k,maxk
      integer arrniter(10)
      double precision arrtolerance(10)
      common/iarr/arrniter,arrtolerance
c
10    continue
      print *
      print *
      print *,'Please select one of the following:'
      print *,'  1..to reade in bc''s and initialize problem'
      print *,'  2..to perform update cycles'
      print *,'  3..to print output'
      print *,'  4..to terminate the work'
      read *, job
      if (job.ne.4) then
        if (job.eq.1) then
          call initlz
          initst = 1
        else if (job.eq.2.and.initst.ne.0) then
          print *
          print *
          print *,'Please give the number of update cycles.'
          read *, icycle
          call update(icycle)
        else if (job.eq.3.and.initst.ne.0) then
          call output
        end if
        goto 10
      end if
      close(2)
      stop
      end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine update (icycle)
      parameter(max=512)
      double precision values(max,max), valold(max, max)
      double precision tolerance, det(max,max)
      integer icycle
      common/val/values
      common/block/nxpts,nypts,nproc
      integer niter
      common/info/niter
      real total1,total2,dtime
      common/time/total2
      double precision tol,tolerance
      common/itol/tol,tolerance
      integer k,maxk
      common/imaxk/k,maxk
      integer arrniter(10)
      double precision arrtolerance(10)
      common/iarr/arrniter,arrtolerance
c   
      niter = 0
      tolerance = 1.d0
      k = 0
      do i = 1,10
          arrniter(i) = 0.d0 
          arrtolerance(i) = 0.d0
      end do
      total1 = dtime(time)
      do while ((tolerance .gt. tol) .and. (niter. lt. icycle))
         niter = niter + 1
         do i = 1,nypts
           do j = 1,nxpts
             valold(i,j) = values(i,j)
           end do
         end do
         do i=2,nypts-1
           do j = 2, nxpts-1
             values(i,j) = (valold(i-1,j)+valold(i+1,j)
     c     +valold(i,j-1)+valold(i,j+1))*.25d0
           end do
         end do
         det(1,1)=abs(values(1,1)-valold(1,1))
         do i = 2,nypts
           do j = 2,nxpts
             det(i,j)=abs(values(i,j)-valold(i,j))
             tolerance=dmax1(det(1,1),det(i,j)) 
             det(1,1)=tolerance
           end do
          end do
c
         if (mod(niter,100) .eq. 0) then
           k = k + 1
           if (k .le. maxk) then
             print *,niter,tolerance
             arrniter(k) = niter
             arrtolerance(k) = tolerance
           end if
         end if         
      end do
      total2 = dtime(time)
      print *,'****total = ',total2
      print *,'last niter, tolerance',niter, tolerance 
      print *, 'These are the 9 points at the Northwest'
      do i=1,3
      print *, (values(i,j), j=1,3)
      end do
      print *, 'These are the 9 points at the Southeast'
      print *, 'They should be same as the Northwest after reversing 
     + the order'
      do i=nypts,nypts-2,-1
      print *, (values(i,j), j=nxpts,nxpts-2,-1)
      end do
      return 
      end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine initlz
      parameter (max=512)
      double precision values(max,max)
      common/val/values
      common/block/nxpts,nypts,nproc 
      integer niter
      common/info/niter
      double precision ledge,redge,bedge,tedge
      double precision tol,tolerance
      common/itol/tol,tolerance
      integer k,maxk
      common/imaxk/k,maxk
      integer  arrniter(10)
      double precision arrtolerance(10)
      common/iarr/arrniter,arrtolerance
c
10    continue
      print *
      print *
      print *,'Please give the number of points in x-direction'
      read *, nxpts
      if (nxpts .gt. max) then
        print *
        print *,'*** Invalid number of points ***'
        goto 10
      end if
c
      print *,'Please give the number of points in y-direction.'
      read *, nypts
      if (nypts .gt. max) then
        print *
        print *,'*** Invalid number of points ***'
        goto 10
      end if
c
      nproc = nxpts * nypts
      print *,'The whole number of points is',nproc
      do i = 1, nypts
        do j = 1, nxpts
          values(i,j) = 0.d0
        end do
      end do
      print *,'Give the left, right, bottom, and top edge values.'
      read  *, ledge, redge, bedge, tedge
      do i = 1, nypts
        values(i,1) = ledge 
        values(i,nxpts) = redge
      end do
      do j = 1, nxpts
        values(nypts,j) = bedge
        values(1,j) = tedge
      end do
c     
      print *,'Please give the value of the tolerance.'
      read *, tol
      if (tol .gt. 1.d0) then
        print *
        print *,'***Invalid value of the tolerance'
        goto 10
      end if
c
      print *,'***Please give the number of iterations to print.'
      read *, maxk
      if (maxk .gt. 10) then
        print *
        print *,'***Invalid number of iterations'
        goto 10
      end if
      return
      end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine output
      parameter(max=512)
      double precision values(max,max)
      common/val/values
      common/block/nxpts,nypts,nproc
      integer j, i, j1, j2, i1, i2
      integer niter
      common/info/niter
      double precision tol,tolerance
      common/itol/tol,tolerance
      real total2
      common/time/total2
      integer k,maxk
      common/imaxk/k,maxk
      integer arrniter(10)
      double precision arrtolerance(10)
      common/iarr/arrniter,arrtolerance
c
   10 continue    
      print *,'Please give the x-coordinate of the first',
     c         ' and the last points of the block to print'
      read(*,*) j1,j2
      print *,'Please give the y-coordinate of the first',
     c         ' and the last points of the block to print'
      read(*,*) i1, i2
      if((j2 .gt. max) .or. (i2 .gt. max)) then
        print *,'Invalid length of the block'
        goto 10
      end if
      open(unit=2,file='outf77',status='new')
      do i = 1,k
        write(2,*) arrniter(i)
        write(2,*) arrtolerance(i)
      end do
      write(2,*)
      write(2,*) niter
      write(2,*) tolerance
      write(2,*)
      do i=i1,i2
        do j=j1,j2
          print *, values(i,j)
          write (2,*) values(i,j)
        end do
      end do
      write (2,*)
      write (2,*) total2
      return
      end

