3.3 Coding the Algorithm

In order to represent Jacobi's method in Fortran 77 syntax (i.e. with DO loops), the program must explicitly define a temporary array to hold the results of the intermediate computations of each iteration.[1] At the end of each iteration, this temporary array must be copied back onto the main array, as in the following code: 

DO k = 1, number_of_iterations
 DO i = 2, n-1                ! Update non-edge
  DO j = 2, n-1               ! elements only
   temp(i, j) = (slab(i, j-1)+slab(i-1, j)+slab(i+1, j)+slab(i, j+1))/4
  END DO
 END DO
 DO i = 2, n-1
  DO j = 2, n-1
   slab(i, j) = temp(i, j)
  END DO
 END DO
END DO

In order to make use of the compiler's nearest-neighbor optimization, Digital recommends coding nearest-neighbor algorithms using Fortran 90/95 syntax, rather than using INDEPENDENT DO loops. To achieve good speed-up in our example problem, you should code the algorithm using FORALL, as follows: 

DO k = 1, number_of_iterations
  FORALL (i=2:n-1, j=2:n-1)      ! Non-edge elements only
    slab(i, j) = (slab(i, j-1)+slab(i-1, j)+slab(i+1, j)+slab(i, j+1))/4
  END FORALL
END DO

There is no need to explicitly define a temporary array to hold intermediate results, because a FORALL structure computes all values on the right side of the assignment statement before making any changes to the left side. (For a full comparison between FORALL structures and DO loops, refer to Section 2.2.3).

For More Information:

[1] Algorithms that modify the array "in place", without the use of temporaries, actually accelerate the convergence. We have chosen Jacoby's method only because of the simplicity of the algorithm.