PROGRAM jacobi
	IMPLICIT NONE

	integer, parameter                  ::  dim = 64
	real*8, parameter                   ::	epsilon = 1.0e-6
	integer                             ::  iter, time1, time2, rate

	real*8, dimension(dim, dim)         ::  old, new

! ------------------------ initializations ----------------------------

	CALL RANDOM_NUMBER(old)		! random initial data ...

	old(1,:)   = 1.0		! ... but unit boundary conditions
	old(dim,:) = 1.0
	old(:,1)   = 1.0
	old(:,dim) = 1.0

	new = 1.0			! anything with boundary = 1 is fine

	iter = 0			! we count the iterations ...

	CALL system_clock(time1)	! ... and measure the running time

! ------------------------ Jacobi iteration ---------------------------

	DO
		new(2:dim-1, 2:dim-1) = 0.25 * (old(1:dim-2, 2:dim-1) + &
						old(3:dim, 2:dim-1) + &
						old(2:dim-1, 1:dim-2) + &
						old(2:dim-1, 3:dim))
		iter = iter + 1
		IF (MAXVAL(ABS(new-old)) < epsilon) EXIT
		old = new
	END DO

! ----------------------- output of results ---------------------------

	CALL system_clock(time2, rate)

! output the diagonal of the array : all values are within epsilon of 1
	DO iter=1,64
		print *, new(iter,iter)
	END DO

	print *, "running time", (time2 - time1) / rate, " seconds"
	print *, "iterations", iter
	print *, "epsilon", MAXVAL(new-old)

END