Another Example - Laplace's Equation

The problem addressed here is solution of the two-dimensional Laplace equation with Dirichlet boundary conditions by the Jacobi relaxation method. An HPJava implementation:

  Procs2 p = new Procs2() ;
  on(p) {
    Range x = new ExtBlockRange(N, p.dim(0)) ;
    Range y = new ExtBlockRange(N, p.dim(1)) ;

    double [[-,-]] a = new double [[x, y]] ;

    // Initialize `a' - set boundary values.

    overall(i = x for :)
      overall(j = y for :)
        if(i` == 0 || i` == N - 1 || j` == 0 || j` == N - 1)
          a [i, j] = i` * i` - j` * j` ;
        else
          a [i, j] = 0.0 ;
    
    
    // Main loop.

    double [[-,-]] b = new double [[x, y]], r = new double [[x, y]] ;

    do {
      Adlib.writeHalo(a) ;

      overall(i = x for 1 : N - 2)
        overall(j = y for 1 : N - 2) {

          b [i, j] = 0.25 * (a [i - 1, j] + a [i + 1, j] +
                             a [i, j - 1] + a [i, j + 1]) ;

          r [i, j] = Math.abs(b [i, j] - a [i, j]) ;
        }

      HPutils.copy(a, b) ;

    } while(Adlib.maxval(r) > EPS) ;
  }

If i is a distributed index scoped by an overall construct, the expression i` (written with a backquote, and read "i-primed") represents the integer global index associated with the loop. In this example - because no other explicit upper and lower bounds were specified in the overall header - the value of i` ranges between 0 and N-1 - all global index values associated with the range x. The initialization is done with a pair of nested overall constructs. Inside, a conditional tests if we are on an edge of the array. If so, the element values are set to some chosen expression - the boundary function. Interior elements are zeroed. Actually the assignment to zero could be have been omitted, because zero is the default value for double array elements in Java.

Notice that one can freely use ordinary Java constructs like if inside an overall construct. One can also include local variable declarations, and so on. Declarations of collective objects like distributed arrays and process grids are allowed as well, so long as one understands that the APG is restricted within the overall body (see the section on the Active Process Group). In this sense the HPJava distributed control construct are true, compositional control constructs. This is different from the HPF forall construct, which has restrictive rules about what kind of statement can appear in its body.

A single iteration of Jacobi relaxation involves replacing each interior element by the average of its neighbouring values. We make use of the technique of ghost regions mentioned in the previous section. Because the arrays - and in particular the array a - are created using ranges of class ExtBlockRange, extra space is allocated around the edges of the locally held segment of the array. In this case the "halo area" has the default width, 1. We illustrate in the case where the process grid is 3 by 2 and N is 8:

Here the uninitialized ghost cells are represented by the yellow zeroes. The communication method Adlib.writeHalo() copies element values from neighboring processes into the ghost regions. We can visulize the state of the array after the call to this method as:

The ghost-cells on the extreme edges of the array are left unchanged by this call to Adlib.writeHalo(). They are not used in the present algorithm (but if the algorithm needs it one can give options to the writeHalo() method that cause these cells to be updated with wrap-around edge conditions). As mentioned earlier, the exact interpreation of the ghost-cells is determined by the program, not the HPJava language. In fact one could initialize the unused outer ghost regions to hold the fixed boundary condition values for our equation; we didn't do that here because it adds some complexity.

The computation of the elements of b is done in overall constructs that specify non-default lower and upper bounds for the iteration - in this case 1 and 6. So the area of b that is actually computed is the red area below:

Ghost cell values from a can only accessed from within an HPJava program by the special syntax of shifted indexes. We said before that with one exception a subscript in a distributed dimension of an array must be a distributed index symbol. The exception occurs in the case where the array has ghost regions. In this case the subscript can have the form i + expression or i - expression, where i is a distributed index symbol - bound to a location in the range of the array as usual - and +/-expression is an int expression whose value lies between minus the lower ghost region width and plus the upper ghost region width. A shifted index always subscripts into the extended local segment of the distributed array. It is always the programmer's responsibility to ensure that, if a shifted access implies reading the value of a ghost cell, the accessed cells have been initialized properly (typically through some communication operation).

The properties of shifted indexes are relatively complex for a feature appearing as a language primitive. But the stencil pattern of update is so important that this "impure" feature was admitted into the language early on, and it has proven very useful.

After computing the whole of b, all its elements are copied to a using the HPutils.copy() library method (this is the correct order of operations to implement Jacobi relaxation.) The main loop terminates when the largest change in any element is smaller than some predefined value EPS. The collective library function Adlib.maxval() finds the largest element of a distributed array, and broadcasts its value to all processes in the APG. Because this involves a relatively expensive global synchronization, in practise one may not wish to test the convergence criterion in every iteration.

Next: Red-Black Relaxation