
\section{\label{examples}Some example parallel programs}

     We conclude this discussion of data parallel programming in
Fortran by displaying a couple of more substantial applications.
These have been chosen to illustrate the parallel features of the
Fortran 90.  The first is a parallel implementation of a standard
matrix algorithm.  The other two are implementations of more intricate
parallel algorithms.  Explanation of the working of these algorithms is
beyond the scope of this report.

     They are coded in pure Fortran 90, without HPF extensions.
They have been tested using an existing Fortran 90 compiler.

\subsection{Linear solver}

     Figures \ref{lu-1} and \ref{lu-2} give a program
for solving of linear equations by the LU-decompostion technique.
It illustrates data-parallel programming through Fortran 90
array syntax.
The implementation of \verb$EXCHANGE$, which swaps the
contents of its two arguments, is omitted.

     Section \ref{alignment} takes a simplified kernel of
this program as a basis for a distributed data programming
example.

     Note the use of the \verb$SPREAD$ intrinsic for converting
one-dimensional arrays into two-dimensional ones.  This is a
common feature of programs written using pure Fortran 90 array syntax.
The version given in section \ref{alignment} uses an HPF
\verb$FORALL$ construct, and is rather more transparent.

     The number of parallel operations needed to complete the program is
proportional to the matrix size, $N$ (whereas the total number of
arithmetic operations is proportional to $N^3$).

\begin{figure}
\begin{verbatim}
  SUBROUTINE LINEAR_SOLVER (A, B)
  REAL, DIMENSION ( : , : ) :: A, B
  ! solve A X = B.  results for X overwrite B.  A is also overwritten.
    INCLUDE 'interfaces.inc'

    CALL LU_DECOMPOSE (A, B)
    CALL SOLVE_LOWER  (A, B)
    CALL SOLVE_UPPER  (A, B)
  ENDSUBROUTINE

  SUBROUTINE LU_DECOMPOSE (A, B)
  REAL, DIMENSION ( : , : ) :: A, B
  ! Solve L U = A.  Results for L and U overwrite A.  B is
  ! pivotted as A.
    INTEGER N, R, R1, PIVOT (1)
    REAL, DIMENSION (SIZE (A, 1)) :: L_COL, U_ROW

    N = SIZE (A, 1)
    DO R = 1, N - 1
      PIVOT = R - 1 + MAXLOC (ABS (A (R : , R)))
      CALL EXCHANGE (A (R, : ), A (PIVOT (1), : ))
      CALL EXCHANGE (B (R, : ), B (PIVOT (1), : ))

      R1 = R + 1
      L_COL (R  : ) = A (R : , R)
      U_ROW (R1 : ) = A (R, R1 : ) / L_COL (R)
      A (R    , R1 : ) = U_ROW (R1 : )
      A (R1 : , R1 : ) = A (R1 : , R1 : ) - &
        SPREAD (L_COL (R1 : ), 2, N - R) * SPREAD (U_ROW (R1 : ), 1, N - R)
    ENDDO
  ENDSUBROUTINE
\end{verbatim}
\caption{\label{lu-1}Linear equation solution by LU-decomposition}
\end{figure}

\begin{figure}
\begin{verbatim}
  SUBROUTINE SOLVE_LOWER (L, B)
  REAL, DIMENSION ( : , : ) :: L, B
  ! Solve L Y = B.  Result for Y overwrites B.
    INTEGER N, M, R, R1
    REAL L_COL (SIZE (L, 1)), Y_ROW (SIZE (B, 2))
  
    N = SIZE (L, 1)
    M = SIZE (B, 2)
    DO R = 1, N
      R1 = R + 1
      L_COL (R : ) = L (R : , R)
      Y_ROW        = B (R, : ) / L_COL (R)
      B (R    , : ) = Y_ROW
      B (R1 : , : ) = B (R1 : , : ) - &
        SPREAD (L_COL (R1 : ), 2, M) * SPREAD (Y_ROW, 1, N - R)
    ENDDO
  ENDSUBROUTINE

  SUBROUTINE SOLVE_UPPER (U, Y)
  REAL, DIMENSION ( : , : ) :: U, Y
  ! Solve U X = Y.  Result for X overwrites Y.
    INTEGER N, M, R, R1
    REAL U_COL (SIZE (U, 1)), X_ROW (SIZE (Y, 2))

    N = SIZE (U, 1)
    M = SIZE (Y, 2)
    DO R = N - 1, 1, -1
      R1 = R + 1
      U_COL ( : R) = U ( : R, R1)
      X_ROW        = Y (R1, : )
      Y ( : R, : ) = Y ( : R, : ) - &
        SPREAD (U_COL (: R), 2, M) * SPREAD (X_ROW, 1, R)
    ENDDO
  ENDSUBROUTINE
\end{verbatim}
\caption{\label{lu-2}Linear equation solution (details)}
\end{figure}
\newpage  %FUDGE

\subsection{Parallel sorting}

     Parallel sorting is a more intricate problem.  The previous
example took a standard sequential algorithm and parallelised it.  Here,
we need a specialised parallel algorithm.  The one used here is odd-even merge
sort.  Parallel merge sorts work on a ``divide and conquer'' basis.
The data is split into two halves which are individually sorted
(in parallel).  Then the two sorted halves are merged by a
procedure which gives a sorted whole.  This strategy is applied recursively,
to sort the data set in $log_2 (N)$ stages, provided the array size $N$
is a power of two.  The divide and conquer strategy is handled by the
recursive function \verb$SORT$.  Merging is handled by \verb$OE_MERGE$
(which also happens to be a recursive procedure).  Detailed
discussion of the odd-even merge algorithm is omitted here.

     The program illustrates the use of recursive procedures in
Fortran 90.  An explicit interface of \verb$OE_MERGE$ is needed
(in \verb$SORT$) because it takes assumed shape array arguments---this
would be included in the file ``interface.inc''.  It also illustrates
the use of general regular array sections, with strides.

     The version given in figure \ref{sorting} is not fully
parallelised.  That requires the two invocations of \verb$SORT$
appearing as arguments to \verb$OE_MERGE$ to be evaluated in parallel,
and that the two recursive invocations of \verb$OE_MERGE$ (in the
assignments to \verb$SOD$ and \verb$SEV$) are evaluated in parallel.
Making the parallelism in these operations completely explicit, simply
using array syntax, would involve a certain amount of code
reorganisation.  Neither Fortran 90 nor full HPF provide very good
support for the divide and conquer paradigm\footnote{In full HPF one could
combine pairs of arrays into single arrays, declare the functions
\verb$PURE$, and use \verb$FORALL$ constructs with an index running
over two values to achieve the extra parallelism.  This is
slightly obscure, and it's debatable whether such convoluted
code would compile efficiently.}.

     A fully parallel implementation would involve a number of
parallel steps propertional to $(log_2 (N))^2$.
\begin{figure}[hbtp]
\begin{verbatim}
  RECURSIVE FUNCTION SORT (S) RESULT (OUT)
  INTEGER, DIMENSION (:) :: S
  INTEGER, DIMENSION (SIZE (S)) :: OUT

    INCLUDE 'interfaces.inc'
    INTEGER N
 
    N = SIZE (S)
    IF (N == 1) THEN
      OUT = S
    ELSE
      OUT = OE_MERGE (SORT (S (1 : N / 2)), SORT (S (N / 2 + 1 : N)))
    ENDIF
  ENDFUNCTION


  RECURSIVE FUNCTION OE_MERGE (S1, S2) RESULT (OUT)
  INTEGER, DIMENSION (:) :: S1
  INTEGER, DIMENSION (SIZE (S1)) :: S2
  INTEGER, DIMENSION (2 * SIZE (S1)) :: OUT

    INTEGER, DIMENSION (SIZE (S1)) :: SOD, SEV
    INTEGER N
 
    N = SIZE (S1)
    IF (N == 1) THEN
      IF (S1 (1) < S2 (1)) THEN
        OUT = (/ S1 (1), S2 (1) /)
      ELSE
        OUT = (/ S2 (1), S1 (1) /)
      ENDIF
    ELSE
      SOD = OE_MERGE (S1 (1 : N - 1 : 2), S2 (1 : N - 1 : 2))
      SEV = OE_MERGE (S1 (2 : N     : 2), S2 (2 : N     : 2))
      OUT (1) = SOD (1)
      WHERE (SOD (2 : N) < SEV (1 : N - 1))
        OUT (2 : 2 * N - 2 : 2) = SOD (2 : N)
        OUT (3 : 2 * N - 1 : 2) = SEV (1 : N - 1)
      ELSEWHERE
        OUT (2 : 2 * N - 2 : 2) = SEV (1 : N - 1)
        OUT (3 : 2 * N - 1 : 2) = SOD (2 : N)
      ENDWHERE
      OUT (2 * N) = SEV (N)
    ENDIF
  ENDFUNCTION
\end{verbatim}
\caption{\label{sorting}Odd-even Merge Sort}
\end{figure}

\subsection{Parallel list ranking}

     The list ranking problem: the entries in a linked list is of $N$
items are held in an array of size $N$.  The array \verb$NEXT$ holds the
links; the entry \verb$NEXT (I)$ is the index of the successor of element
\verb$I$ in the list.  The end of the list is flagged by the
condition \verb$NEXT (I) = I$.  The problem is to find the distance of
each element from the end of the list.

     Figure \ref{ranking} is an implementation of a standard parallel
algorithm for solving this problem.  It illustrates the use of
``irregular sections''---array sections with vector subscripts.

     The program completes in a number of parallel steps proportional to
$log_2 (N)$.
\begin{figure}[hbtp]
\begin{verbatim}
  FUNCTION RANK (NEXT)
  INTEGER NEXT (:)
  INTEGER RANK (SIZE (NEXT))

    INTEGER, DIMENSION (SIZE (NEXT)) :: P, DISTANCE, INDEX

    N = SIZE (NEXT)
    INDEX = (/ (I, I = 1, N) /)
    P = NEXT
    WHERE (P /= INDEX)
      DISTANCE = 1
    ELSEWHERE
      DISTANCE = 0
    ENDWHERE
    I = 1
    DO WHILE (I < N)
      WHERE (P /= P (P))
        DISTANCE = DISTANCE + DISTANCE (P)
        P = P (P)
      ENDWHERE
      I = 2 * I
    ENDDO
    RANK = DISTANCE
  ENDFUNCTION
\end{verbatim}
\caption{\label{ranking}List ranking by the ``Doubling'' technique}
\end{figure}

