\documentstyle{article}

\begin{document}

\def\gap{\vspace{0.3in}}

\title{Procedural versus Declarative Programming for Scientific
       and Engineering Applications: The Half-Time Score}

\author{Tony Hey, Bryan Carpenter and Hugh Glaser}
\date{}
\maketitle

\section*{Outline}

\begin{enumerate}
  \item Introduction
  \item Phenomenology
  \begin{itemize}
     \item LNLL Sisal project
     \item FLARE project
     \item Case study
  \end{itemize}
  \item Conclusions/Questions/Perspectives
\end{enumerate}

\newpage

\section{Introduction}

Declarative language enthusiasts claim that there are real benefits to
be gained from programming in a declarative style.

\begin{quote}
``We believe, as did McGraw in 1984, that increased productivity,
generality, utility, portability, and performance are only possible
if programmers avoid the constraints of imperative languages and
adopt a higher level of abstraction.''  {\em David Cann, CACM 1992}
\end{quote}
{\em or}

%\newpage

\begin{quote}
``Declarative functional programming has been shown to be a promising
approach in overcoming the difficulties of writing large correct
software systems and of exploiting the power of parallel computing.''
{\em P.\ Grant et al, FLARE Report 1992}
\end{quote}
{\em or}

\begin{quote}
``Hudak points out that functional languages are not oriented to
any particular machine model and they are easier to compile for
parallel machine execution than conventional imperative languages
because imperative languages are biased towards a sequential machine
model.''  {\em P.\ Hudak, 1991}
\end{quote}

\newpage

\begin{itemize}
\item
These claims seem all embracing and {\em not} restricted to any
specific application area
\begin{list}{\bf ?}{}
\item Operating systems
\item Telecoms applications
\item Safety critical embedded systems
\item Science and Engineering applications
\end{list}
\end{itemize}

\begin{list}{$\Rightarrow$}{}
\item
In this talk we examine the evidence for the thesis that
scientists and engineers should abandon Fortran and C(++)
and switch to a functional programming language.
\end{list}

\newpage

\section*{Claimed benefits of declarative languages}

\begin{quote}
\em Compilers for imperative languages {\em handicapped} by the
{\em assignment statement}
\end{quote}

\begin{verbatim}
    A = F (X)
    B = G (Y)
\end{verbatim}

\begin{list}{$\Rightarrow$}{}
\item
Common blocks allow functions to share data
\item
Aliasing can meant that \verb$A$, \verb$X$, \verb$B$, \verb$Y$, may
represent same memory cell
\item
Parallelism not obvious and requires detailed interprocedural analysis
\item
{\em too difficult for all but simple programs}
\end{list}

%\newpage

\gap

\begin{quote}
\em Functional languages express problem in terms of well-defined
functions
\end{quote}

\begin{list}{$\Rightarrow$}{}
\item
Referential transparency
\end{list}

\begin{itemize}
\item
no side effects, better software engineering
\item
simpler semantics allows formal proofs and program transformation
\item
higher level of abstraction allows rapid prototyping
\item
implicit parallelism, no need to worry about architectural details
\end{itemize}

\newpage

\subsection*{What do scientists and engineers care about?}

\begin{list}{---}{}
\item
Transparent relation of mathematics to programming languages
\item
Portability and standards, especially for parallel code
\item
Efficient programs and high performance
\item
Support for data parallelism.
(Only really useful type of parallelism for MPP's)
\end{list}

\subsection*{? Mismatch of concerns ?}

\begin{list}{---}{}
\item
Formal correctness rarely an issue
\item
C, Fortran adequate for rapid prototyping of Sci/Eng problems
\end{list}

\newpage

\section*{Data parallelism}

\begin{list}{---}{}
\item
Common problem---solution of PDE
\begin{displaymath}
\mbox{\boldmath $\nabla$}^2 u (x, y) = - c (x, y).
\end{displaymath}
{\em 2-D Poisson's equation}
\item
Numerical solution may be obtained by discretization
\begin{displaymath}
u (x, y) = \frac{1}{4} (
\begin{array}[t]{l}
u (x + 1, y) + u (x - 1, y) + \\
u (x, y + 1) + u (x, y - 1) + c (x, y) \,\,)
\end{array}
\end{displaymath}
\end{list}
Finite difference approximation
\begin{list}{$\Rightarrow$}{}
\item
2-D grid of points
\item
Relaxation algorithm
\end{list}

\begin{displaymath}
\begin{array}[t]{c}
u (x, y) \\
\uparrow \\
\mbox{``new'' value}
\end{array}
= \frac{1}{4} (\sum
\begin{array}[t]{c}
\mbox{neighbours} \\
\times \\
\times \, \, \circ \, \, \times \\
\times
\end{array}
\begin{array}[t]{c}
+ \\
\begin{array}[t]{l}
\nwarrow \\
\mbox{``old'' value}
\end{array}
\end{array}
c (x, y))
\end{displaymath}

\newpage

\centerline{\psfig{figure=dataPar.eps,height=3in}}

\gap

Sequential version $\Rightarrow$ Nested Do Loops

\begin{list}{$\Rightarrow$}{}
\item
Parallel implementations ``easy'' on Vector, SIMD and MIMD
\end{list}

\gap

\centerline{\Large \em ``SPMD style''}
\centerline{\Large \em Single Program Multiple Data}

\newpage

\section{Phenomenology}

\subsection*{LNLL SISAL PROJECT}

\begin{list}{}{}
\item
Concept of ``data'' replaces ``memory''
\item
Concept of ``creation'' replaces ``update''
\end{list}

\begin{description}
\item{Problem:} Without side-effects an array update must define an entirely
new array
\end{description}

\begin{quotation}
``In place operations can be automatically compiled into most
functional programs without violating their computational semantics''

{\em COPY ELIMINATION}
\end{quotation}

\begin{table}[hbtp]
\centerline{
\begin{tabular}{|c|c|c|}
\hline
PROGRAM & NO OPTIMIZATION & OPTIMIZATION \\
\hline
LOOPS   & 681.6           & 16.0         \\
\hline
SIMPLE  & 177.7           &  9.4         \\
\hline
\end{tabular}
}
\caption{Cray Y-MP: One CPU.  Seconds}
\end{table}

\newpage

\begin{table}[hbtp]
\centerline{
\begin{tabular}{|c|c|c|c|c|c|c|}
\hline
PROGRAM &
\multicolumn{2}{c|}{ONE CPU} &
\multicolumn{2}{c|}{4 CPUs} &
\multicolumn{2}{c|}{8 CPUs} \\
\cline{2-7}
       & SISAL & FORTRAN & SISAL & FORTRAN & SISAL & FORTRAN \\
\hline
LOOPS  & 30    & 37      & 39    & 42      & 42    & 43      \\
\hline
SIMPLE & 20    & 21      & 61    & 21      & 85    & 21      \\
\hline
\end{tabular}
}
\caption{PARALLEL EXECUTION.  Mflop/s}
\end{table}

Cann concludes:  Functional \ldots
\begin{enumerate}
\item
Programs easier to write and maintain and are more concise
\item
Programs amenable to analysis
\item
Programs exhibit implicit parallelism and favour automatic parallelization
\item
Programs can run as fast as conventional languages.
\end{enumerate}

\begin{list}{$\Rightarrow$}{}
\item
But most of Sci/Eng community remain unconvinced
\end{list}

\newpage

\subsection*{UK FLARE Project}

\begin{quote}
Functional Languages Applied to Realistic Exemplars
\end{quote}
Dual aims:
\begin{list}{---}{}
\item
transferring effective functional programming technology and tools
from laboratory to serious practitioners
\item
demonstrating that functional languages are a suitable vehicle
for effectively programming parallel machines
\end{list}

\gap

Scientific Applications being investigated by group in Swansea
CS Department
\begin{list}{$\Rightarrow$}{}
\item
Haskell for CFD application
\begin{enumerate}
\item
Sequential: Fortran $\rightarrow$ Haskell
\item
Parallel: version for GRIP
\end{enumerate}
\end{list}

\newpage

\subsubsection*{Sequential results}

\begin{itemize}
\item
Both Haskell compilers produce inadequate optimization of array operations
\begin{list}{---}{}
\item either too slow or too space demanding
\end{list}
\end{itemize}

\begin{table}[hbtp]
\centerline{
\begin{tabular}{|c|c|c|}
\hline
Program & Fortran & Haskell \\
\hline
FB Problem 1 (200)   & 4.3           & 91.3                 \\
\hline
FB Problem 3 (800)  & 16.8           & exhausted heap space \\
\hline
\end{tabular}
}
\caption{Timings on Sun (seconds)}
\end{table}
\begin{list}{$\Rightarrow$}{}
\item ``Heap space problems prevent us attempting larger applications''
\item Haskell/Fortran $\approx$ 25/1
\end{list}

\begin{quote}
``The way forward is to introduce some parallelism \ldots
which is commonly held to be easier in a functional setting than
in an imperative one''
\end{quote}

\newpage

\subsubsection{Parallel Results}

\begin{list}{$\Rightarrow$}{}
\item
Need to insert explicit parallelism primitives for GRIP implementation
\end{list}

\begin{table}[hbtp]
\centerline{
\begin{tabular}{|c|c|c|cccc|}
\hline
           & Fortran & Haskell & \multicolumn{4}{c|}{GRIP} \\
\cline{2-7}
           & Sparc 1 & Sparc 1 & 1    & 3    & 6    & 9    \\
\hline
Speed(sec) & 0.23    & 3.4     & 86.3 & 39.0 & 30.7 & 23.0 \\
\hline
\end{tabular}
}
\end{table}

\begin{list}{$\Rightarrow$}{}
\item
``Can only run our code on an 8-element trivial test problem''
\end{list}

\gap

Conclude?

\begin{list}{---}{}
\item
``Parallel performance fits expected trend''
\item
``Need to investigate larger problems to establish more clearly
the advantages of parallelization''
\item
``Merit of functional programming includes relative simplicity of
parallelization''
\end{list}

\centerline{\Large ?\ ?\ ?}

\newpage

\subsection*{Case Study: Jacobi Relaxation on 2-D Laplace's equation}

\begin{itemize}
\item
Fortran version $\rightarrow$ SISAL, HASKELL versions
\item
Comments
\begin{list}{---}{}
\item
Length of code:
\begin{displaymath}
\mbox{Fortran} > \mbox{Sisal} > \mbox{Haskell}
\end{displaymath}
\item
Sisal (v1.2):
\begin{list}{}{}
\item
lacks facilities for updating sections of arrays
\item
has a rather restrictive set of reduction operations
\item
is less easy to customize than Haskell
\end{list}
\item
Haskell implementations in terms of ``array of streams'' or
``stream of arrays''---details in paper by Carpenter and Glaser
\end{list}
\end{itemize}

\newpage

\begin{verbatim}
      program main
      implicit none

c     Simple Jacobi relaxation for Laplace equation.
        
        integer n
        real eps
        parameter (n = 64, eps = 1.0e-6)

        real u (0 : n, 0 : n), uold (0 : n, 0 : n)
        real bdry
        integer i, j, niter
        logical convgd

c       Zero interior and initialise boundaries.

        do i = 0, n
          do j = 0, n
            if (i .eq. 0 .or. i .eq. n .or. j .eq. 0 .or. j .eq. n) then
	      u (i, j) = bdry (float (i) / n, float (j) / n)
            else
	      u (i, j) = 0.0
            end if
          end do
        end do

        convgd = .false.
        niter  = 0
        do while (.not. convgd)

c         Save old value.

          do i = 0, n
	    do j = 0, n
              uold (i, j) = u (i, j)
            end do
          end do

c         Update.

          do i = 1, n - 1
	    do j = 1, n - 1
              u (i, j) = (uold (i, j - 1) + uold (i, j + 1) +
     &                    uold (i - 1, j) + uold (i + 1, j)) / 4
            end do
          end do

c         Check convergence.

          convgd = .true.
          do i = 1, n - 1
	    do j = 1, n - 1
              convgd = convgd .and.
     &                 (abs (u (i, j) - uold (i, j)) .lt. eps)
            end do
          end do

          niter = niter + 1
        end do

        print *, niter
c        print *, u
      end

      real function bdry (x, y)
      implicit none
      real x, y

        bdry = x**2 - y**2
      end
\end{verbatim}

\newpage

\begin{verbatim}
define main

% Simple Jacobi relaxation for Laplace equation.
 
function bdry (x, y : real returns real)
  x * x - y * y
end function

%function main (returns integer, array [array [real]])
function main (returns integer)
  let
    n     := 64 ;
    eps   := 1.0E-6
  in
    for initial

      % Zero interior and initialise boundaries.

      u      := for i in 0, n cross j in 0, n
		returns array of
		  if i = 0 | i = n | j = 0 | j = n then
		    bdry (real (i) / real (n), real (j) / real (n))
		  else
		    0.0
		  end if
	        end for ;

      convgd := 0 ;
      niter  := 0

    while convgd = 0 repeat

      % Update.

      u := for i in 0, n cross j in 0, n
	   returns array of
             if i ~= 0 & i ~= n & j ~= 0 & j ~= n then
               (old u [i] [j - 1] + old u [i] [j + 1] +
		old u [i - 1] [j] + old u [i + 1] [j]) / 4.0
             else
               old u [i] [j]
             end if
	   end for ;

      % Check convergence.

      convgd := for i in 1, n - 1 cross j in 1, n - 1
                returns value of product
                  if abs (u [i] [j] - old u [i] [j]) < eps then 1 else 0 end if
                end for ;

      niter  := old niter + 1
%    returns value of niter value of u
    returns value of niter
    end for
  end let
end function
\end{verbatim}

\newpage

\begin{verbatim}
module Main (main) where

-- Simple Jacobi relaxation for Laplace equation.

n   = 6
eps = 1.0e-6

main = print (let

                -- Zero interior and initialise boundaries.

                u = array ((0, 0), (n, n))
                          [(i, j) :=
                             if i == 0 || i == n || j == 0 || j == n then
                               bdry (fromIntegral i / fromIntegral n)
                                    (fromIntegral j / fromIntegral n)
                             else
                               0.0
                          | i <- [0..n], j <- [0..n]]

              in
                loop False 0 u)

loop convgd niter u =
  if convgd then
    --(niter, u)
    niter
  else
    let

      -- Update.

      u'      = u // [(i, j) :=
                          (u ! (i, j - 1) + u ! (i, j + 1) +
                           u ! (i - 1, j) + u ! (i + 1, j)) / 4.0
                     | i <- [1..n - 1], j <- [1..n - 1]]

      -- Check convergence.

      convgd' = and [abs (u' ! (i, j) - u ! (i, j)) < eps
                    | i <- [1..n - 1], j <- [1..n - 1]]

      niter'  = niter + 1
    in 
      loop convgd' niter' u'
    
bdry x y = x * x - y * y
\end{verbatim}

\newpage

\begin{verbatim}
?- op(700, xfx, ->).

isN(16).
isEps(1.0E-6).

main(Res) :- gridDomain(UInit), initialVals(UInit), loop(0, UInit, Res).

% `gridDomain' asserts that a mapping (``array'') has an N x N grid as
% domain (``index space'').  NB: it imposes a particular ordering---
% reverse lexicographic---on sites in the mapping list.

gridDomain(U) :- isN(N), enumDomain(N, N, U).

enumDomain(0, 0, [site(0, 0) -> _])     :- !.
enumDomain(I, 0, [site(I, 0) -> _ | U]) :-
			!, IM is I - 1, isN(N), enumDomain(IM, N, U).
enumDomain(I, J, [site(I, J) -> _ | U]) :-
			JM is J - 1, enumDomain(I, JM, U).

% `initialVals' asserts that a mapping defines a suitable initial state.

initialVals([]).
initialVals([site(I, J) -> B | UTail]) :-
    isN(N),
    (I =:= 0 ; I =:= N ; J =:= 0 ; J =:= N), !,
    X is float(I) / float(N),
    Y is float(J) / float(N),
    bdry(X, Y, B),
    initialVals(UTail).
initialVals([_ -> 0.0 | UTail]) :- initialVals(UTail).

% The main loop...

loop(Niter, U, Res) :- gridDomain(UNew),
                       update(U, UNew),
                       Miter is Niter + 1,
                       loopTest(Miter, U, UNew, Res).

% `update' asserts that all mappings in 2nd arg equal the average of
% values of neighbour mappings in 1st arg.

update(_, []).
update(Uold, [site(I, J) -> X | UTail]) :-
    isN(N),
    (I =:= 0 ; I =:= N ; J =:= 0 ; J =:= N), !,
    element(site(I, J) -> X, Uold),
    update(Uold, UTail).
update(Uold, [site(I, J) -> X | UTail]) :-
    JM is J - 1, element(site(I, JM) -> A, Uold),
    JP is J + 1, element(site(I, JP) -> B, Uold),
    IM is I - 1, element(site(IM, J) -> C, Uold),
    IP is I + 1, element(site(IP, J) -> D, Uold),
    X  is (A + B + C + D) / 4.0,
    update(Uold, UTail).

% Check for convergence.

%loopTest(Niter, OldU, U, res(Niter, U)) :- converged(OldU, U), !.
loopTest(Niter, OldU, U, Niter) :- converged(OldU, U), !.
loopTest(Niter, _,    U, Res)           :- loop(Niter, U, Res).

converged([], []).
converged([_ -> X | U], [_ -> Y | V]) :- Dif is abs(X - Y),
					 isEps(Eps), Dif < Eps,
					 converged(U, V).

bdry(X, Y, B) :- B is X * X - Y * Y.

element(X, [X | _]) :- !.
element(X, [_ | L]) :- element(X, L).
\end{verbatim}

\newpage

\subsubsection*{Some timings}

\begin{table}[hbtp]
\centerline{
\begin{tabular}{|c|c|c|c|c|c|}
\hline
n  & iter & f77 & Sisal & Haskell & Prolog \\
\hline
12 & 125  & 0.2 & 0.5   & 14.1    & 280    \\
\hline
16 & 213  & 0.4 & 1.4   & 43.5    & 1080   \\
\hline
24 & 443  & 1.9 & 6.2   & 207.4   &        \\
\hline
32 & 740  & 5.3 & 18.1  & 701.4   &        \\
\hline
\end{tabular}
}
\end{table}

\subsubsection*{Comments}

\begin{itemize}
\item
Difficult to optimize functional versions because of complex relation
between execution speed and memory utilisation
$\Rightarrow$ cf. Profiling tools of Runciman and Wakeling
\item
Functional versions capture ``essence'' of algorithm in a way that
Fortran probably never can $\Rightarrow$ could have a role in
specification or prototyping
\end{itemize}

\section{Conclusions/Questions/Perspectives}

Conclude: Case for appropriateness of functional languages for
scientific and engineering applications ``not proven''.

Imperative languages well matched to mass-produced microprocessor
architectures $\Rightarrow$ compilers, tools pretty well understood

\begin{list}{$\Rightarrow$}{}
\item
Functional languages should be used in applications where
the are {\em evidently} useful and superior to other technologies
\end{list}

\newpage

Why no mass market for ``Haskell''?

\gap

\centerline{\psfig{figure=standards.eps,height=2in}}

\gap

\begin{itemize}
\item
Standards vital
\begin{quotation}
SQL, POSIX, HPF, MPI
\end{quotation}
\item
Best parallel programming tool
\begin{quotation}
FORGE90 from APR

$\rightarrow$ sold $< 100$ licenses!
\end{quotation}
\end{itemize}

\end{document}
