
The sequential code for this example:

A is the ideal matrix.  The boundary values on the 4 sides of A
do not change.


Example 2.11 Jacobi iteration - sequential code
REAL A(0:n+1,0:n+1), B(1:n,1:n)
...
!Main Loop
DO WHILE(.NOT.converged)
   ! perform 4 point stencil
   DO j=1, n
      DO i=1, n
         B(i,j)=0.25*(A(i-1,j)+A(i+1,j)+A(i,j-1)+A(i,j+1))
      END DO
   END DO

   ! copy result back into array A
   DO j=1, n
      DO i=1, n
         A(i,j) = B(i,j)
      END DO
   END DO
...
! convergence test omitted
END DO
The code fragment describes the main loop of an iterative solver where
,at each iteration, the value at a point is replaced by the average of
the North , South, East and West neighbours(a fout point pencil is used
to keep example simple).Boundary values do not change.We focus on the 
inner loop, where most of the computation is done, and use fortran 90
syntax, for clarity.



Parallel code:

The matrix will be distributed by columns.

p = numProcs
m = number of columns of ideal matrix on this proc
    ex: 5, 5, 5, 4, 4  (for n=23, p=5).
The local matrix A has dimensions (n+2) x (m+2), indexed by
     0..n+1 x 0..m+1. The top and bottom rows appear to be constant.

The initial assumption seems to be that A is initialized correctly
(including ghost cells) on all procs.  No assumption seems to be made 
about B.


E.g.

  x

Example2.16 Using nonblocking communications in Jacobi computation.
...
REAL, ALLOCATABLE A(:,:), B(:,:)
INTEGER req(4)
INTEGER status(MPI_STATUS_SIZE, 4)
...
! Compute number of processes and myrank
CALL MPI_COMM_SIZE(comm, p, ierr)
CALL MPI_COMM_RANK(comm, myrank, ierr)

! compute size of local block
m = n/p
IF (myrank.LT.(n-p*m)) THEN
   m = m+1
END IF

! Compute neighbors
IF (myrank.EQ.0) THEN
   left = MPI_PROC_NULL
ELSE
   left = myrank - 1
END IF
IF (myrank.EQ.p-1)THEN
  right = MPI_PROC_NULL
ELSE
  right = myrank+1
END IF

! Allocate local arrays
ALLOCATE (A(0:n+1,0:m+1), B(n,m))
...
!Main Loop
DO WHILE(.NOT.converged)
   ! compute
   DO i=1, n
         B(i,1)=0.25*(A(i-1,j [sic])+A(i+1,j [sic])+A(i,0)+A(i,2))
         B(i,m)=0.25*(A(i-1,m)+A(i+1,m)+A(i,m-1)+A(i,m+1))
   END DO

  ! Communicate
      CALL MPI_ISEND(B(1,1),n, MPI_REAL, left, tag, comm, req(1), ierr)
      CALL MPI_ISEND(B(1,m),n, MPI_REAL, right, tag, comm, req(2), ierr)
      CALL MPI_IRECV(A(1,0),n, MPI_REAL, left, tag, comm, req(3), ierr)
      CALL MPI_IRECV(A(1,m+1),n, MPI_REAL, right, tag, comm, req(4), ierr)

  ! Compute interior
   DO j=2, m-1
      DO i=1, n
         B(i,j)=0.25*(A(i-1,j)+A(i+1,j)+A(i,j-1)+A(i,j+1))
      END DO
   END DO
   DO j=1, m
      DO i=1, n
         A(i,j) = B(i,j)
      END DO
   END DO

  ! Complete communication
   DO i=1, 4
      CALL MPI_WAIT(req(i), status(1.i), ierr)
   END DO
...
END DO
