
\section{\label{implementations}Lazy implementation of Red-black
relaxation}

     With the geometry sorted out, we can tackle the main program.
Define a data type which captures the local
state in a relaxation calculation:
\begin{verbatim}
  type Obs    = Float

  data State  = Interior Float Obs |
                Boundary Obs

  siteObs :: State -> Obs
  siteObs (Interior _  u) = u
  siteObs (Boundary    u) = u
   
  siteTrans :: State -> [Obs] -> State
  siteTrans (Interior c _) nbs =
             Interior c ((sum nbs + c) /
                     fromIntegral (length nbs))
  siteTrans  b             _   = b
\end{verbatim}
There are two types of site: normal sites on the ``interior''
of the lattice, and those on some designated boundary (typically
the sites on the edge of the array).
Sites on the interior
have a local charge density and a potential value.  Those on the
boundary just have a potential value, which will remain fixed
throughout the iteration.  The function \verb$siteObs$
extracts the ``observable'', i.e. the potential, from a
\verb$State$ variable.  \verb$siteTrans$ performs the update of
equation (\ref{relaxation}) on the interior sites and leaves boundary
sites unchanged.  It is passed the old value of the site and a list
of potential values at nearest neighbours.

     The history of a single site can be regarded as a stream of
\verb$State$ values\footnote{By definition, a ``stream'' is
an infinite list}.  The site evolves through application of
a transition function which depends on the (changing) environment of
the site.  This pattern of evolution can be captured through
the function\footnote{{\bf process} can also be implemented
by the standard function {\bf scanl}.  Note the lazy pattern-match on
the last parameter of {\bf process}.
This laziness is required for the implementations below.}
\begin{verbatim}
  process :: (s -> e -> s) -> s -> [e] -> [s]
  process trans st ~(env : envs) =
      st : process trans (trans st env) envs
\end{verbatim}
When the transition function is \verb$siteTrans$, the ``current environment''
is a list of neighbour
potential values, and the last argument to \verb$process$ will
be a stream of such lists---the history of the environment.

     The result of the computation will be an array of lists
of \verb$State$ values, defining the history of the whole lattice.
Suppose we are given this array and need to calculate
the array of environments from it.  This can be done by
\begin{verbatim}
  env :: Array Coord [State] ->
           Array Coord [[Obs]]
  env = amap transpose .
          nns (\_ _ -> map siteObs)
\end{verbatim}
The application of \verb$nns$, in which the mapping function
extracts a stream of potential values from the stream of
\verb$State$ values, gives an array containing a list of
streams of neighbour histories.  Mapping \verb$transpose$
over it turns these into streams of neighbour lists.
These are just what \verb$siteTrans$ needs as its argument.
So the implementation of relaxation is
\begin{verbatim}
  evolve :: Array Coord State ->
              Array Coord [State]
  evolve init = x
    where x = azipWith (process siteTrans)
                init e
          e = env x
\end{verbatim}
The argument \verb$init$ is an array of initial states (\verb$azipWith$ is
defined to pick up the lattice from its first array parameter).

     As well as the lattice size, the array of initial states defines:
which array elements are on the boundary of the physical system being
described; the charge distribution; the initial guess for the potential
on \verb$Interior$ sites; and the fixed values of the potential at
\verb$Boundary$ sites.  It's required that no \verb$Interior$ site
should be on the edges of the array (only \verb$Boundary$ sites are
allowed there).  Because \verb$siteTrans$ doesn't evaluate its neighbours
argument for \verb$Boundary$ sites, laziness ensures automatically
that no neighbour access violates array bounds.

     The code above implements Jacobi relaxation.
More general schemes, like Gauss-Siedel,
can easily be implemented in the same framework.
Suppose a site $A$ has site $B$ as a neighbour.
The scheme above calculates the new value at $A$ in terms
of the previous generation value at $B$.  To calculate
the value of $A$ in terms of the {\em latest} value of $B$,
we should look one place forward in the history of $B$---in
other words the $B$ stream should be tailed before
feeding it into the environment of $A$.
To organise this, it is convenient to attribute a colour
to each site.  If the colour of $B$, say, is less than
the colour $A$, $A$ is updated in terms of the latest
value of $B$.  Otherwise it is updated in terms of the
previous generation value of $B$.  The function \verb$env$
is changed to
\begin{verbatim}
  env = amap transpose .
          nns (\i j -> map siteObs .
            if colour j < colour i
              then tail
              else id)
\end{verbatim}
Red-black Relaxation in 2 dimensions can be implemented by setting
\begin{verbatim}
  colour (x, y) = (x + y) `mod` 2
\end{verbatim}

