
%\documentstyle{csart}

%\begin{document}

\section{\label{examples}Examples}

     All but one of the examples given here assume the ``hidden domain''
implementation outlined in section \ref{static}, where the index grid
is identified with the problem grid.  This allows algorithms to be
specified more straightforwardly, but doesn't always permit the
most efficient implementation.

\subsection{Simple Grid Problems}

     As a simple prototype of a grid-oriented problem 
take solution of the Laplace equation by red-black
relaxation.  A straightforward and fairly detailed implementation is
given in figure \ref{red-black}.
\begin{figure}
\begin{verbatim}

    [n + 1] INDEX x, y :
    RANK XY IS {x, y} :
    XY REAL32 U :
    BOOL converged :
    SEQ
      WHERE XY
        IF
          (x = 0) OR (x = n) OR (y = 0) OR (y = n)
            U := boundary (x, y)
          TRUE
            U := 0.0
      parity    := 0
      converged := FALSE
      WHILE NOT converged
        XY BOOL ok :
        XY REAL32 neighb.xn, neighb.xp, neighb.yn, neighb.xp :
        SEQ
          rotate (x, neighb.xn, U,  1)
          rotate (x, neighb.xp, U, -1)
          rotate (y, neighb.yn, U,  1)
          rotate (y, neighb.yp, U, -1)
          WHERE XY
            IF
              (((x + y) REM 2) = parity) AND
                  (0 < x) AND (x < n) AND (0 < y) AND (y < n)
                REAL32 new :
                SEQ
                  new := 0.25 * ((neighb.xn + neighb.xp)
                               + (neighb.yn + neighb.yp))
                  ok  := ABS (U - new) < epsilon
                  U   := new
              TRUE
                ok := TRUE
          parity    := (parity + 1) REM 2
          converged := and.2 (ok)

\end{verbatim}
\caption{\label{red-black}Red-black Relaxation for Laplace Equation}
\end{figure}
The ordinary function $boundary$ defines the fixed edge values of $U$.
The ``predefine'' $and.2$ performs logical ``and'' across the elements of a
two-dimensional boolean
distribution (a possible definition is given in appendix \ref{primitives}).

    One practical issue that arises is the requirement for whole
distributions of neighbours.  This is unlikely to be a problem for
Laplace's equation, where memory requirements are modest.  In more
memory-intensive applications it could be awkward.  There are
techniques for significantly reducing the amount of memory used by
temporary arrays.  For example in figure \ref{reduced-mem} the $y$
neighbours are declared as just \verb${y}$ distributions---normally
such arrays require less memory than \verb${x, y}$ distributions.
\begin{figure}
\begin{verbatim}

    ...
    SEQ
      ...
      WHILE NOT converged
        ...
        {x,y} REAL32 neighb.xn, neighb.xp :
        SEQ
          rotate (x, neighb.xn, U,  1)
          rotate (x, neighb.xp, U, -1)
          WHERE {x}
            {y} REAL32 neighb.yn, neighb.xp :
            SEQ
              rotate (y, neighb.yn, U,  1)
              rotate (y, neighb.yp, U, -1)
              WHERE {y}
                ... body of `WHERE XY' loop
          ...

\end{verbatim}
\caption{\label{reduced-mem}Red-black Relaxation with Reduced Memory}
\end{figure}

     There is another style of geometric programming
which can also be expressed in the current notation.  In this
case the grid defined by the indicies is regarded as the
process grid, and domain decomposition is managed explicitly
by the programmer.  For a slightly
simplified version of the Laplace equation an example is
given in figure \ref{domain}.
\begin{figure}
\begin{verbatim}

    [n / s] INDEX x, y :
    RANK XY IS {x, y} :
    XY [s] [s] REAL32 U :
    SEQ
      ... initialisation of 'U'
      WHILE NOT converged
        XY BOOL ok :
        SEQ
          WHERE XY
            ok := TRUE
          SEQ i = 0 FOR s
            SEQ j = 0 FOR s
              XY REAL32 neighb.xn, neighb.xp, neighb.yn, neighb.xp :
              SEQ
                WHERE XY
                  SEQ
                    neighb.xn := U [(i - 1) MOD s] [j]
                    neighb.xp := U [(i + 1) MOD s] [j]
                    neighb.yn := U [i] [(j - 1) MOD s]
                    neighb.yp := U [i] [(j + 1) MOD s]
                IF
                  i = 0
                    rotate (x, neighb.xn, U,  1)
                  TRUE
                    SKIP
                IF
                  i = (s - 1)
                    rotate (x, neighb.xp, U, -1)
                  TRUE
                    SKIP
                IF
                  j = 0
                    rotate (y, neighb.yn, U,  1)
                  TRUE
                    SKIP
                IF
                  j = (s - 1)
                    rotate (y, neighb.yp, U, -1)
                  TRUE
                    SKIP
                WHERE XY
                  REAL32 new :
                  SEQ
                    new := 0.25 * ((neighb.xn + neighb.xp)
                                 + (neighb.yn + neighb.yp))
                    ok  := ok AND (ABS (U [i] [j] - new) < epsilon)
                    U [i] [j] := new
          converged := and.2 (ok)
                
\end{verbatim}
\caption{\label{domain}Relaxation with Explicit Domain Decomposition}
\end{figure}
We ignore here the issue alternating parity and
treatment of boundaries.  This implementation is clearly more
complicated than the first one (almost as complicated as an ordinary
{\em occam} implementation).  We have explicitly kept an $s$ by $s$
sequential
segment of the array on each of $n/s$ by $n/s$ ``processes''.
Note that in this way the space occupied by the neighbours is much smaller
than the data array (if $s$ is moderately large) ---we only
need four words per processor instead of some arrays the same
size as $U$.

     Incidentally the \verb$MOD$ operator is a piece of shorthand---it
means roughly the same as the {\em occam} \verb$REM$ operator, redefined
to return the ``correct'' value for negative arguments.

\subsection{Simple Matrix Operations}

     Figure \ref{transpose} implements matrix transpose for a two-dimensional
distribution (of \verb$INT$s) in terms of the simpler operation $rotate$.
In practise both $in.situ.rotate$ and (a polymorphic form of)
$transpose.matrix$ would probably be provided in predefined libraries,
implemented more efficiently.
\begin{figure}
\begin{verbatim}

    PROC in.situ.rotate ([] INDEX x, {x} \/ RANK r TYPE t A,
                          VAL INT shift)

      -- Rotate input data in 'A', returning result also in 'A'.

      WHERE r
        {x} t temp :
        SEQ
          rotate (x, temp, A, shift)
          WHERE {x}
            temp := shift
    :

    PROC transpose.matrix ([] INDEX x, y, {x, y} INT result, A)

      -- Set 'result' to transpose of 'A'.
      -- Must have SIZE x = SIZE y.

      SEQ i = 0 FOR SIZE x
        SEQ
          WHERE {x, y}
            IF
              x = y
                result := A
              TRUE
                SKIP
          in.situ.rotate (x, result, 1)
          in.situ.rotate (y, A     , 1)
    :

\end{verbatim}
\caption{\label{transpose}Matrix Transpose on a Grid}
\end{figure}

     In a similar vein figure \ref{multiply} embodies a standard
distributed algorithm for multiplying two square matrices.  Notice
the call to $broadcast$.  Remember that at this point in the program
$x$ is an index, $y$ is simply an integer value (the ``loop counter''
for the $y$ loop).
\begin{figure}
\begin{verbatim}

    [n] INDEX x, y :
    {x, y} INT A, B, result :
    SEQ
      ... initialisation of 'A', 'B'
      WHERE {x, y}
        result := 0
      SEQ i = 0 FOR n
        SEQ
          WHERE {y}
            INT diag :
            SEQ
              diag := broadcast (x, A, y)
              WHERE {x}
                result := result + (diag * B)
          in.situ.rotate (x, A, -1)
          in.situ.rotate (y, B, -1)

\end{verbatim}
\caption{\label{multiply}Matrix Multiplication on a Grid}
\end{figure}

\subsection{LU Decomposition}

     The example given here for LU or Cholesky decomposition
is based on the recurrence
\\ \\
\begin{math}
L = \left(
\begin{array}{cc}
l_{11} & 0   \\
l_C    & L'
\end{array} \right),\; U = \left(
\begin{array}{cc}
1      & u_R \\
0      & U'
\end{array} \right),\; A = \left(
\begin{array}{cc}
a_{11} & a_R \\
a_C    & A'
\end{array} \right)
\end{math}
\\ \\
with $LU = A$, i.e.
\\ \\
\begin{math}
l_{11} = a_{11}, \; l_C = a_C, \; u_R = l_{11}^{-1}a_R
\end{math}
\\ \\
and
\\ \\
\begin{math}
L'U' = A' - l_C u_R
\end{math}
\\ \\
If $A$, $L$ and $U$ are $n \times n$ matrices, 
$A'$, $L'$ and $U'$ are $(n-1) \times (n-1)$ matrices, so $n$
applications of this operation reduce $A$ to LU form.
This isn't quite the way LU decomposition is presented in
most textbooks, but it suffices for present purposes.

     Figure \ref{LU} gives a ``column-oriented'' distribution.  The
function $locate.largest$ is an ordinary sequential routine for
locating the largest element of an array.  Note that the broadcast
of $n$ numbers in $reduction$ takes time dominated by a term
proportional to $n$.  The arithmetic in one reduction stage also
takes parallel time proportional to $n$---in fact considerably longer
if there are several columns per process.  Thus the communications
overhead is at most a constant factor.

     The algorithm is not particularly well load-balanced if the
distribution of data is the most direct one, because at any point
one is only operating on a slice of the whole array.  Again this
leads to a constant inefficiency, one that is often
alleviated by ``wrap-around'' mapping of data to processors.
The specification here is at too high a level to express this
mapping, although it is possible that the implementation might
adopt it.
\begin{figure}
\begin{verbatim}

    INT FUNCTION find.pivot ({[n] INDEX x} [n] REAL32 mat, VAL INT col)

      -- Calc pivot for column 'col', and broadcast it.

      INT pivot :
      VALOF
        {x} INT temp :
        SEQ
          WHERE {x}
            IF
              x = col
                temp := col + locate.largest ([mat FROM col FOR n - col])
              TRUE
                SKIP
          broadcast (x, pivot, temp, col)
        RESULT pivot 
    :

    PROC reduction ({[n] INDEX x} [n] REAL32 mat, VAL INT stage)

      INT pivot :
      [n] REAL32 first.col :
      SEQ
        pivot := find.pivot (mat, stage)
        WHERE {x}
          mat [stage], mat [pivot] := mat [pivot], mat [stage]
        broadcast (x, first.col, mat, stage)
        WHERE {x}
          IF
            x > stage
              SEQ
                mat [stage] := mat [stage] / first.col [stage]
                SEQ i = stage + 1 FOR n - (stage + 1)
                  mat [i] := mat [i] - (first.col [i] * mat [stage])
            TRUE
              SKIP
    :

    [n] INDEX x :
    {x} [n] REAL32 A :

    SEQ r = 0 FOR n
      reduction (A, r)

\end{verbatim}
\caption{\label{LU}Column-oriented LU Decomposition}
\end{figure}

%\end{document}
