
!2345678901234567890123456789012345678901234567890123456789012345678901234567890

      program

      integer n, niter
      parameter(n = 100, niter = 10000)

      integer i, j

      real h 

      real u (n, n), f (n, n), sol (n, n)

      real fac

      real t1, t2

      integer iter, parity, i1, j1

      h = 1.0 / (n - 1)
  
! Initialization:
!    sol = exact solution
!    f   = forcing term
!    u   = approximation solution
  
      fac = 4.0 - 2.0 * (cos(h) + cos(3 * h))
  
      do i = 1, N
        do j = 1, N
          sol (i, j) = sin(3 * (h * (i - 1)) + h * (j - 1))
          if(i .eq. 1 .or. i .eq. N .or.
     &       j .eq. 1 .or. j .eq. N) then
            u (i, j) = sol (i, j)
          else
            u (i, j) = 0.0
            f (i, j) = fac * sol (i, j)
          endif
        enddo
      enddo

!      print *, 'initial u:'
!      do i = 1, N
!        print '(11f7.3)', (u (i, j), j = 1, N)
!      enddo

! Main loop

      call etime(t1)

      do iter = 1, niter
        do parity = 0, 1
          do i1 = 2, N - 1
            do j1 = 2 + mod(i1 - 1 + parity, 2), N - 1, 2

              u (i1, j1) = 0.25 * (f (i1, j1) +
     &                             u (i1 + 1, j1) + u (i1 - 1, j1) +
     &                             u (i1, j1 + 1) + u (i1, j1 - 1))
            enddo
          enddo

!          print *, 'iteration =', iter
!          do i = 1, N
!            print '(11f7.3)', (u (i, j), j = 1, N)
!          enddo
        enddo
      enddo

      call etime(t2)

!      print *, 'exact solution:'
!      do i = 1, N
!        print '(11f7.3)', (sol (i, j), j = 1, N)
!      enddo

      print *, 'elapsed time =', t2 - t1

      end
