PROGRAM JACOBI_MAIN
*
* Jacobi iterations: Sequential code
*
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER ( NSIZE = 1024*1024*15/8 )
      DIMENSION F(NSIZE), U(NSIZE), UNEW(NSIZE), UA(NSIZE)
      REAL SECOND
*
      PRINT *, 'Enter A, B, N, M, NIT: '
      READ(*,*) A, B, N, M, NIT
      IF ( (N+1)*(M+1) .GT. NSIZE ) THEN
         PRINT *, 'Insufficient buffer size (increase NSIZE).'
         STOP
      ENDIF
*
      CALL INIT(A, B, N, M, U, UA, F)
      T0=SECOND()
      CALL JACOBI(NIT, A, B, N, M, U, UNEW, F)
      T1=SECOND()
      PRINT *, 'Wall time = ', T1-T0
      CALL CHKSOL(N, M, U, UA)
*
      STOP
      END
*
*-----------------------------------------------------------------------
*
      SUBROUTINE INIT(A, B, N, M, U, UA, F)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION F(0:N,0:M), U(0:N,0:M), UA(0:N,0:M)
*
      HX=A/DBLE(N)
      HY=B/DBLE(M)
*
      DX=1.D0/(HX*HX)
      DY=1.D0/(HY*HY)
      D=1.D0/(DX+DX+DY+DY)
*
* Analytical solution
      DO J=0,M
      DO I=0,N
         X = I*HX
         Y = J*HY
         UA(I,J) = -0.25D0*( X**2 + Y**2 )
      ENDDO
      ENDDO
*
* Initial approximation and RHS
      DO J=1,M-1
      DO I=1,N-1
         F(I,J)=1.D0 * D
         U(I,J)=0.D0
      ENDDO
      ENDDO
*
* Boundary condition
      DO I=0,N
         U(I,0)=UA(I,0)
         U(I,M)=UA(I,M)
      ENDDO
      DO J=0,M
         U(0,J)=UA(0,J)
         U(N,J)=UA(N,J)
      ENDDO
*
      RETURN
      END
*
*-----------------------------------------------------------------------
*
      SUBROUTINE CHKSOL(N, M, U, UA)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION U(0:N,0:M), UA(0:N,0:M)
*
* Compute Linfty(error)
      ERR=0.D0
      DO J=1,M-1
      DO I=1,N-1
         ERR = MAX( ERR, ABS( U(I,J)-UA(I,J) ) )
      ENDDO
      ENDDO
      PRINT 1000, ERR
1000  FORMAT(' Error = ', 1PG23.16)
*
      RETURN
      END
*
*-----------------------------------------------------------------------
*
      SUBROUTINE JACOBI(NIT, A, B, N, M, U, UNEW, F)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION F(0:N,0:M), U(0:N,0:M), UNEW(0:N,0:M)
*
      HX=A/DBLE(N)
      HY=B/DBLE(M)
*
      DX=1.D0/(HX*HX)
      DY=1.D0/(HY*HY)
      D=1.D0/(DX+DX+DY+DY)
      DX=DX*D
      DY=DY*D
*
      DO IT=1, NIT
         DO J=1,M-1
         DO I=1,N-1
            UNEW(I,J) = F(I,J) + DX*( U(I-1,J) + U(I+1,J) ) +
     +                           DY*( U(I,J-1) + U(I,J+1) )
         ENDDO
         ENDDO
*
         DO J=1,M-1
         DO I=1,N-1
            U(I,J)=UNEW(I,J)
         ENDDO
         ENDDO
*
* Compute residual
		if (.false.) then
         RESIDU=0.D0
         DO J=1,M-1
         DO I=1,N-1
            RES = F(I,J) - U(I,J) + DX*( U(I+1,J) + U(I-1,J) ) +
     +                              DY*( U(I,J+1) + U(I,J-1) )
            RESIDU = MAX( RESIDU, ABS(RES) )
         ENDDO
         ENDDO
         PRINT 1000, IT, RESIDU/D
1000     FORMAT(' It = ', I5, '   Residual = ', 1PG23.16)
		endif
*
      ENDDO
*
      RETURN
      END