PROGRAM JACOBI
      IMPLICIT NONE
      include "mpif.h"
      INTEGER Dim,Iteration,RowPeek,ColPeek
      CHARACTER*100 buffer
      REAL*8, DIMENSION(:,:),ALLOCATABLE ::  SurfaceMatrix
      REAL*8, DIMENSION(:,:),ALLOCATABLE ::  SurfaceMatrix_t
      INTEGER i,j,ii
      REAL*8 Increment
      REAL*8 timestart,timeend
!     MPI Variables
      INTEGER ierror,NCPUS,MY_PE
      INTEGER LocalDim
      INTEGER Req1,Req2,Req3,Req4
      INTEGER STAT(MPI_STATUS_SIZE)
      
      CALL MPI_INIT(ierror)
      CALL MPI_COMM_SIZE(MPI_COMM_WORLD,NCPUS,ierror)
      CALL MPI_COMM_RANK(MPI_COMM_WORLD,MY_PE,ierror)
      
      IF(IARGC() < 4) THEN
         WRITE(6,*)"Too Few Input Parameters.. there should be four";
         CALL MPI_FINALIZE(ierror)
         CALL EXIT(2)
      ENDIF
      
      ! Get Command Line Arguments
      CALL GETARG(1,buffer)
      READ(buffer,*)Dim
      CALL GETARG(2,buffer)
      READ(buffer,*)Iteration
      CALL GETARG(3,buffer)
      READ(buffer,*)RowPeek
      CALL GETARG(4,buffer)
      READ(buffer,*)ColPeek
      
      IF((MOD(Dim,NCPUS)).ne.0) THEN
         IF(MY_PE == 0) THEN
            WRITE(6,*)"The Dimensions of the matrix need to be"
            WRITE(6,*)"divisible by the total number of cores"
         ENDIF
         CALL MPI_FINALIZE(ierror)
         CALL EXIT(4)
      ENDIF

      IF((RowPeek > Dim ).or.(ColPeek > Dim)) THEN
         IF(MY_PE == 0) THEN
            WRITE(6,*)"Cannot Peek a matrix element outside of",
     $           "the surface"
            WRITE(6,*)"Arguments 3 and 4 must be smaller than",Dim
         ENDIF
         CALL MPI_FINALIZE(ierror)
         CALL EXIT(3)
      ENDIF
      
! Initialize local matrices
      LocalDim = Dim/NCPUS
      ALLOCATE(SurfaceMatrix(Dim+2,LocalDim+2))
      ALLOCATE(SurfaceMatrix_t(Dim+2,LocalDim+2))
      DO i=1,LocalDim+2
         DO j=1,Dim+2
            SurfaceMatrix(j,i)=0.0
         ENDDO
      ENDDO

! Fill Initial Values
      DO i=2,LocalDim+1
         DO j=2,Dim+1
            SurfaceMatrix(j,i)=0.5
            ENDIF
         ENDDO
      ENDDO
      
!     Left Side
      Increment = 100.0 / (Dim + 1)
      DO i=2,LocalDim+1
        SurfaceMatrix(1,i) = (MY_PE*LocalDim+i-1)*Increment
        SurfaceMatrix(Dim+2,i)=0.0
      ENDDO
      if(MY_PE == 0) THEN
         DO i=2,Dim+1
            SurfaceMatrix(i,1)=0.0
         ENDDO
      ENDIF
      
      IF(MY_PE == (NCPUS-1)) THEN
         DO i=2,Dim+1
            SurfaceMatrix(i,LocalDim+2)=(Dim-i+2)*Increment
         ENDDO
      ENDIF
      timestart = MPI_WTIME()
      DO ii = 1,Iteration
         if(MY_PE.ne.0) THEN
            CALL MPI_IRECV(SurfaceMatrix(2,1),Dim,MPI_REAL8
     $           ,MY_PE-1,0,MPI_COMM_WORLD,req3,ierror)
         ENDIF
         if(MY_PE.ne.(NCPUS-1)) THEN
            CALL MPI_IRECV(SurfaceMatrix(2,LocalDim+2),Dim,MPI_REAL8,
     $           MY_PE+1,0,MPI_COMM_WORLD,req4,ierror)
         ENDIF
!     Send the Appropriate Rows (err columns)
         if(MY_PE.ne.(NCPUS-1)) THEN
            CALL MPI_ISEND(SurfaceMatrix(2,LocalDim+1),Dim, MPI_REAL8
     $           ,MY_PE+1,0,MPI_COMM_WORLD,req1,ierror)
         ENDIF
         IF(MY_PE.ne.0) THEN
            CALL MPI_ISEND(SurfaceMatrix(2,2),Dim,MPI_REAL8,
     $           MY_PE-1,0,MPI_COMM_WORLD,req2,ierror)
         ENDIF
         
         IF(MY_PE.ne.0) THEN
            CALL MPI_WAIT(req3,stat,ierror)
            CALL MPI_WAIT(req2,stat,ierror)
         ENDIF
         IF(MY_PE.ne.(NCPUS-1)) THEN
            CALL MPI_WAIT(req4,stat,ierror)
            CALL MPI_WAIT(req1,stat,ierror)
         ENDIF
 
!        CALL PrintSurfaceMatrix(SurfaceMatrix,LocalDim,Dim,MY_PE,NCPUS)
!         CALL MPI_FINALIZE(ierror)
!         CALL EXIT(10)
!     Actually do the Calculations
         DO i = 2,Dim+1
            DO j = 2,Dim+1
               SurfaceMatrix_t(j,i) = 0.25*(SurfaceMatrix(j-1,i) + 
     $                                     SurfaceMatrix(j,i+1) +
     $                                     SurfaceMatrix(j+1,i) +
     $                                     SurfaceMatrix(j,i-1))
            ENDDO               ! j
         ENDDO                  ! i 

!     Transfer back to the main matrix
         DO i=2,Dim
            DO j=2,Dim
               SurfaceMatrix(j,i) = SurfaceMatrix_t(j,i)
            ENDDO               ! j
         ENDDO                  ! i
      ENDDO                     ! ii
      timeend = MPI_WTIME()
      WRITE(6,'(A,F20.10)')'The Iteration Time is ',
     $     timeend-timestart
      CALL PrintSurfaceMatrix(SurfaceMatrix,Dim)
      IPRowPeek = RowPeek / (Dim/NCPUS)
      
      WRITE(6,'(A,I5,A,I5,A,F20.10)')'The Matrix Element ',ColPeek,',',
     $     RowPeek," is ",SurfaceMatrix(ColPeek+1,RowPeek+1)
      
      CALL MPI_FINALIZE(ierror)
      END PROGRAM
      
      SUBROUTINE PrintSurfaceMatrix(SM,LD,D,MP,NP)
      IMPLICIT NONE
      INCLUDE 'mpif.h'
      INTEGER LD,D,i,j,ip,ierror
      INTEGER MP,NP
      REAL*8 SM(D+2,*)
      
      WRITE(6,*)"Number of Core ",NP," My Core ",MP
      DO ip = 0,NP-1
         if(ip == MP) THEN
            WRITE(6,'(A,I10)')'Matrix from Processor',MP
            WRITE(6,*)' '
            WRITE(6,'(A)',advance='no')'       '
            DO i=1,D+2
               WRITE(6,'(I10)',advance='no'),i
            ENDDO
            WRITE(6,*)' '
            DO i=1,LD+2
               WRITE(6,'(I10)',advance='no'),i
               DO j=1,D+2
                  WRITE(6,'(F10.5)',advance='no'),SM(j,i)
               ENDDO
               WRITE(6,*)' '
            ENDDO
         ENDIF
         CALL MPI_BARRIER(MPI_COMM_WORLD,ierror)
      ENDDO
      
      RETURN
      END SUBROUTINE