;;; This codes generates and compiles stencil application functions for
;;; arbitrary stencils.  For simplicity, we assume Dirichlet boundary
;;; conditions.

(defun compress-stencil (width stencil)
  "Transforms a 2D-stencil given as array in a sparse form where entries
with equal value are put together."
  (let ((sparse-stencil ()))
    (dotimes (i 3)
      (dotimes (j 3)
	(let ((value (coerce (aref stencil i j) 'double-float)))
	  (unless (zerop value)
	    (let ((pos (+ (1- i) (* (1- j) width)))
		  (entry-group (assoc value sparse-stencil)))
	      (if entry-group
		  (push pos (cdr entry-group))
		  (push (list value pos) sparse-stencil)))))))
    sparse-stencil))

(defun apply-stencil-code (width height stencil)
  "Code for applying a stencil to an array of given width and height."
  (let ((cstencil (compress-stencil width stencil)))
    `(lambda (result values)
      (declare (type (simple-array double-float (*)) result values))
      (declare (optimize (speed 3) (safety 0) (debug 0) (compilation-speed 0)))
      (loop
       for pos1 of-type (integer ,width ,(* width height))
       from ,width below ,(* width (1- height)) by ,width do
       (loop
	for pos2 of-type (integer ,width ,(* width height))
	from (1+ pos1) below (+ pos1 ,(1- width)) do
	(incf (aref result pos2)
	      (+
	       ,@(loop
		  for (value . indices) in cstencil collecting
		  `(* ,value
		    (+ ,@(loop for index in indices collecting
			       `(aref values (+ pos2 ,index)))))))))))))

(defun stencil-application-function (width height stencil)
  "Returns a function for a stencil application."
  (compile nil (apply-stencil-code width height stencil)))

(defun test (stencil n &optional (count (/ 100000000 (* n n))))
  "Applies stencil count times to an n x n grid function."
  (let ((entries (make-array (expt n 2) :element-type 'double-float
			     :initial-element 1.0d0))
	(stencil-function (stencil-application-function n n stencil)))
    (time 
     (dotimes (i count (aref entries (1+ n)))
       (funcall stencil-function entries entries)))))

(defparameter *nine-point-stencil*
  #2A((-1/3  -1/3  -1/3)
      (-1/3   8/3  -1/3)
      (-1/3  -1/3  -1/3)))

(test *nine-point-stencil* 25)