\documentstyle{article}

\begin{document}

\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{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}

\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}

These claims seem all-embracing---not restricted to any specific
application area.  So we can hope to adopt functional programming when
coding operating systems, telecoms applications, safety
critical embedded systems, and so on.  Our particular interest
here is in science and engineering applications.
In this article we examine the evidence for the thesis that
scientists and engineers should abandon Fortran and C or C++,
and switch to a functional programming language.

\subsection*{Claimed benefits of declarative languages}

It is claimed that compilers for imperative languages are handicapped by the
{\em assignment statement}.
\begin{verbatim}
    A = F (X)
    B = G (Y)
\end{verbatim}
Analysis of this program fragment is complicated by various issues.
{\em Aliasing} can mean that any of the variables \verb$A$, \verb$X$,
\verb$B$, \verb$Y$ may represent the same memory cell.  Common blocks or
global variables allow functions to share data---\verb$F$ or \verb$G$
could update one of the four variables as side-effect.  Potential
parallelism is therefore not obvious.  Extracting it requires detailed
interprocedural analysis.  This proves too difficult for all but the
simplest programs.

Functional languages on the other hand express a problem in terms of
well-defined, mathematically pure, functions.  The elimination of
side-effects improves the clarity of code and facilitates
software engineering.  The simplified semantics of the languages allow
formal proofs and program transformations.  The level of abstraction
provided by features like higher order polymorphic functions, together
with the liberation of the programmer from detailed concern about
memory management issues, allow rapid prototyping of applications.
Finally, since there is no notion of sequencing of operations in a
functional program, there is implicit parallelism---and the expression of
this parallelism is not tied to any architectural details.

Many of these advantages are indisputable, but one should also
ask what scientists and engineers really care about.  We suggest
\begin{itemize}
\item
Transparent relation of mathematics to programming languages.
\item
Portability and standards, especially for parallel code.
\item
Efficient programs and high performance.
%Support for data parallelism.
\end{itemize}
Functional languages are certainly more ``mathematical'' than
imperative ones; but it remains to be established that a typical
scientific applications needs much more mathematical sophistication
than ``formula translation''.  There is an emergent standard for
functional programming, in the form of Haskell;
%[MENTION ISO(?) STANDARD FOR PROLOG?]
but parallel extensions
or annotations to functional languages are still the subject of
research rather than standardisation.  The efficiency issue will be 
reviewed in some detail in the remainder of the article.

Perhaps there is a mismatch of concerns here.  {\em Formal} correctness
is rarely an issue in a scientific program.  Detractors from
declarative programming would probably claim that C and Fortran
are adequate vehicles for rapid prototyping of scientific
and engineering software.

% Ditch following section?  Just describe problem in case study section?

%\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
%

\subsection*{The remainder of this article}

Section \ref{phenomenology} collects together some evidence concerning
the relative merits of functional and imperative languages in
the domain of scientific programming.

Section \ref{conclusions} draws together our conclusions, based on this
evidence.

\section{\label{phenomenology} Phenomenology}

In this section we selectively review some results from two major
projects, which have promoted the adoption of functional languages for
scientific programming.  We also describe a small case study of our
own, in which we have applied some declarative languages to a
representative kernel algorithm.

\subsection{The LLNL Sisal project}

{\em Sisal} (Streams and Iteration in a Single-Assignment Language) is
a functional language developed at Lawrence Livermore National
Laboratory.  It differs from most other functional languages in taking
support of scientific computation as a primary design goal.
Consequently the early implementations omit many of the advanced
abstractions of purist functional languages (higher order functions,
normal order evaluation, Milner-Hindly type system, and so on).  These are
perceived as having limited usefulness in scientific programming.
Some of them incur a significant penalty in implementation efficiency,
or hinder compile-time analysis.

Sisal introduces a few extra ideas not found in conventional functional
languages, including ``parallel'' and ``sequential'' loops.  It retains
those features of functional languages which facilitate program
analysis and thus optimising transformations.  Programs are {\em
referentially transparent}.  A function passed a given set of argument
values always returns the same value, and (as a rule) within the range
of definition of a variable, its value never changes (in Sisal, this
constraint is relaxed in the eponymous iterative construct).

Dispensing with the idea of a modifiable state, the concept of ``data''
replaces that of ``memory'', and ``creation'' replaces ``update''.
Where an imperative program would update a portion of an array, a
Sisal program logically creates a new array with different values
in the modified portion.

This presents an efficiency problem if the implementation genuinely
allocates new space each time a modified array is created.  But it
is claimed that
\begin{quote}
``In place operations can be automatically compiled into most
functional programs without violating their computational semantics''
\cite{...}
\end{quote}
This optimisation is called {\em copy elimination}.  Table
\ref{copyElim} illustrates the dramatic effects of copy elimination on
two benchmarks (the Livermore Loops and a hydrodynamic code).
\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{\label{copyElim}Effect of Copy Elimination.
Cray Y-MP: one CPU.  Times in seconds}
\end{table}

The most audacious claim of Sisal is that, due to its simplified
side-effect-free semantics, compiler analysis can achieve automatic 
parallelisation more effectively than for an imperative language.
This claim is backed by data such as that in table \ref{parPerf}
\begin{table}[hbtp]
\centerline{
\begin{tabular}{|c|c|c|c|c|c|c|}
\hline
Program &
\multicolumn{2}{c|}{1 CPU} &
\multicolumn{2}{c|}{4 CPUs} &
\multicolumn{2}{c|}{8 CPUs} \\
\cline{2-7}
%       & SISAL & FORTRAN & SISAL & FORTRAN & SISAL & FORTRAN \\
       & OSC & CF77 & OSC & CF77 & OSC & CF77 \\
\hline
LOOPS  & 30    & 37      & 39    & 42      & 42    & 43      \\
\hline
SIMPLE & 20    & 21      & 61    & 21      & 85    & 21      \\
\hline
\end{tabular}
}
\caption{\label{parPerf}Parallel execution.  Cray Y-MP.
Performances in Mflop/s}
\end{table}
The comparison is between code compiled with the Sisal OSC compiler
and Fortran compiled with the Cray CF77 compiler.

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

\subsection{The FLARE Project}

The UK FLARE (Functional Languages Applied to Realistic Exemplars)
project had dual aims:
\begin{itemize}
\item
transferring of effective functional programming technology and tools
from the laboratory to serious practitioners, and
\item
demonstrating that functional languages are a suitable vehicle
for effectively programming parallel machines.
\end{itemize}
The project focussed on the Haskell language.  Haskell is a very
different language to Sisal.  It is a full-blown, lazy functional
language, with an elaborate type system and first-class, higher
order functions.  Sisal is restrictive---often more restrictive
than an imperative language.  Haskell is very free.  It retains
referential transparency, but has a powerful set of abstractions,
and for many practical purposes feels more expressive than a typical
imperative language.

Scientific Applications were investigated by a group in Department of
Computer Science at Swansea.  They have translated a sequential Fortran
finite element code, used in Computational Fluid Dynamics research.  In
two reports they have described the results of benchmarking the
sequential Haskell, and a parallel implementation for the GRIP machine
\cite{...}.

\subsubsection*{Sequential results}

Like many scientific applications, the exemplar developed
made extensive use of arrays.  In contrast to Sisal, neither of the
major Haskell compilers \cite{...} implements arrays particularly
efficiently.  The approaches currently adopted are either too slow
(non-constant access time) or too demanding of space (naive copying at
each update).  The Swansea group built their own implementation of the
Haskell array, better suited to their requirements.  (The fact that
this could be done easily is perhaps characteristic of the flexibility
of Haskell.)

Various concrete representations of the sparse arrays in the problem
were explored.  Here we only record the results with the ``fixed
bandwidth'' scheme, which is directly comparable with the original
Fortran.  Table \ref{haskellPerf} (derived from table 3 in \cite{...})
displays timings for two problem sizes---one with 200 elements and one
with 800 elements.

\begin{table}[hbtp]
\centerline{
\begin{tabular}{|c|c|c|}
\hline
Program & Fortran & Haskell \\
\hline
Problem 1 (200)   & 4.3           & 91.3                 \\
\hline
Problem 3 (800)  & 16.8           & exhausted heap space \\
\hline
\end{tabular}
}
\caption{\label{haskellPerf} Performance on Sun workstation.  Times in seconds}
\end{table}

On the smaller problem there is a factor of about 20 between Fortran
and Haskell timings.  Probably more worrying is the fact that on the
larger application, exhaustion of heap space prevented the Haskell code
from running at all.  More ``Haskell-friendly'' representations of the
sparse arrays did allow the larger problem to be run, but again showed
performance degradation by a factor 25 or more relative to the
original Fortran.

The authors conclude that
\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}

\subsubsection*{Parallel Results}

GRIP (Graph Reduction In Parallel) is an experimental
parallel architecture for executing functional programs \cite{...}.  It
contains a number of processing elements (PEs) together with
intelligent memory units and an input-output unit.  To execute a
program in parallel on GRIP it was necessary to insert explicit
primitives for ``sparking'' threads.  This was reported to be a rather
easy modification to the code.

Table \ref{GRIPPerf} is taken from section 4 of \cite{...}.  It
compares times for the sequential and parallel versions of the CFD
code.  The authors state that
\begin{quote}
``Because the GRIP system as a whole is still under development, at
present we can only run our code on an 8-element trivial test problem''
\end{quote}
\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 PE & 3 PEs & 6 PEs & 9 PEs \\
\hline
Iteration time & 0.23    & 3.4     & 86.3 & 39.0  & 30.7  & 23.0 \\
\hline
\end{tabular}
}
\caption{\label{GRIPPerf} Sequential and parallel times.  Times in seconds}
\end{table}

The authors conclude
\begin{quote}
``Parallel performance fits into an expected trend''
\end{quote}
\begin{quote}
``There is a need to investigate larger problems to establish more clearly
the advantages of parallelization''
\end{quote}
\begin{quote}
``A merit of functional programming includes relative simplicity of
parallelization''
\end{quote}


\subsection{A Case Study: Jacobi Relaxation on the 2-D Laplace equation}

As an further contribution to the emerging picture, we have taken a
couple of very simple algorithms, and coded them in Fortran 77, Sisal
and Haskell.  So far as practical, we endeavoured to use exactly the same
algorithms in the three versions.  The two problems considered
were relaxation for the Laplace equation and a matrix LU decomposition.
We report results for the former, which are representative.

The two-dimensional Laplace equation,
\begin{displaymath}
\mbox{\boldmath $\nabla$}^2 u (x, y) = 0,
\end{displaymath}
can be solved numerically through the finite difference approximation
\begin{displaymath}
u (x, y) = \frac{1}{4} (u (x + 1, y) + u (x - 1, y) +
                        u (x, y + 1) + u (x, y - 1)).
\end{displaymath}
The values of $u$ are defined on a two-dimensional grid.  In the
relaxation algorithm this system is solved iteratively through the
scheme
\begin{displaymath}
\begin{array}[t]{c}
 \\
\begin{array}[t]{r}
\nearrow \\
\mbox{``new'' value}
\end{array}
\end{array}
\begin{array}[t]{c}
u (x, y) \\
\\
\times
\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'' values}
\end{array}
\end{array}
\end{displaymath}
The detailed codes are given in the appendices.  The initial version
was Fortran.  The algorithm was transcribed to the other languages,
trying to keep as close as possible to the style of the original
program.  This was not always practical, and small changes had to be
made.

\subsubsection*{Remarks on programming style}

The Fortran program is slightly longer than the Sisal and significantly
longer than the Haskell, although the style was supposedly
Fortran-like.  To this extent the claims of declarative programs for
conciseness are vindicated.

Subjectively, it was felt Sisal was more awkward to use than
Fortran or Haskell.  Sisal 1.2 does not have very powerful facilities
for updating arrays (especially, for updating subsections of arrays).
For example, the Fortran version of the Jacobi update, leaving
boundaries unchanged, is just
\begin{verbatim}
          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
\end{verbatim}
whereas in Sisal 1.2 we have to make a new copy of the whole array,
including the unchanged boundaries.  This involves using a larger
index range and a ``masked'' loop body:
\begin{verbatim}
      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 ;
\end{verbatim}
As another example, the set of reduction operations is hard-wired into
the language, and quite restrictive.  So in the Sisal program, the
convergence tests had to return integer 0 or 1 rather than a logical,
because Sisal 1.2 has no {\em and} reduction operation.  It was
necessary to use integer product instead.

The Haskell transcription went smoothly.  The {\em while} loop had to
be converted to a tail recursive function call, but otherwise the
program looks very much like the Fortran, underneath the unfamiliar
syntax.  (We could use the built-in list-reduction function {\em and},
which has no analogue in Fortran 77).

%Haskell implementations in terms of ``array of streams'' or
%``stream of arrays''---details in paper by Carpenter and Glaser

\subsubsection*{Some timings}


The codes were run on a Sun Sparc I workstation.  The compilers were:
for Fortran, the Sun {\em f77} compiler; for Sisal, {\em OSC} version
{\em V12.9.1}, the optimising Sisal 1.2 compiler from Lawrence
Livermore; for Haskell, {\em ghc} version {\em 0.16}, the Glasgow
Haskell compiler.  We also tried the Chalmers compiler.  The Glasgow
compiler was faster on this benchmark.  Programs were timed with the
UNIX {\em time} command.

\begin{table}[hbtp]
\centerline{
\begin{tabular}{|c|c|c|c|c|}
\hline
grid & iterations to & f77 & Sisal & Haskell \\
size &  convergence  &     &       &         \\
\hline
12   & 125           & 0.2 & 0.5   & 14.1    \\
\hline
16   & 213           & 0.4 & 1.4   & 43.5    \\
\hline
24   & 443           & 1.9 & 6.2   & 207.4   \\
\hline
32   & 740           & 5.3 & 18.1  & 701.4   \\
\hline
\end{tabular}
}
\caption{Timings for Laplace Relaxation.  Times in seconds.}
\end{table}

\section{\label{conclusions}Conclusions/Questions/Perspectives}

Performance is a central concern in many scientific programs.  Sisal
does show a performance comparable to imperative languages on
scientific, array-based problems.  A suspicion lingers that the cost of
this performance is reduced expressiveness of the language.
Restricting the semantics of the language makes compiler analysis
simpler; it is questionable whether it makes for more concise
or maintainable code.  Sisal 2.0 lifts some of the restrictions of
Sisal 1.2, but implementations were not available for benchmarking.

It has been demonstrated that Sisal can be automatically
parallelised---more effectively than sequential Fortran.  A fairer
comparision might be with High Performance Fortran (HPF), which retains
all the facilities of conventional Fortran and adds parallel updates;
there is every reason to assume that HPF-like imperative languages
should speed up at least as effectively as Sisal on parallel machines.

The lazy functional languages, of which Haskell is representative, are
perhaps the most elegant and economical of all programming languages.
They certainly deserve to be more widely known amongst scientific
programmers, if only for specification and prototyping purposes.
But on the available evidence it is difficult to claim that they
are suitable vehicles for high-performance production codes.  Even our
simplest array-based programs executed a couple of orders of magnitude
more slowly than the Fortran; their profligate use of memory is another
major worry.

[TONY'S LAST WORDS...]

The case for appropriateness of functional languages for scientific and
engineering applications is ``not proven''.

Imperative languages are well matched to mass-produced microprocessor
architectures.  The compilers and tools are pretty well understood.
Functional languages should be used in applications where the are {\em
evidently} useful and superior to other technologies.

Why no mass market for ``Haskell''?

\begin{figure}[hbtp]
\centerline{\psfig{figure=standards.eps,height=2in}}
\end{figure}

\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}

%
%\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
%
