Here is something similar to what I suggested, although still very unfinished. It's probably not the proper way to write a macro, but I guess other people on the list will know how to do it correctly. Efficient linearization can be added in the same manner; also, the macro should be modified s.t. it accepts several array specifications at once, i.e. with-foreign-array ((array-1 :double) (array-2 :int) (array-3 :double) ...), etc.
Of course, if you can figure out a way to incorporate the optimizations into gref without such a clumsy workaround, I'd be all for it.
best regards,
Sebastian
(defun mapcons (fn x)
(if (atom x)
x
(funcall fn (let ((a (mapcons fn (car x)))
(d (mapcons fn (cdr x))))
(if (and (eql a (car x)) (eql d (cdr x)))
x
(cons a d))))))
(defmacro with-fast-access-to-single-foreign-array ((array element-type) &body body)
(alexandria:with-unique-names (array-fptr)
`(let ((,array-fptr (grid::foreign-pointer ,array)))
,@(mapcons
(lambda (expr)
(if (and (consp expr)
(eq (first expr) 'grid:gref*)
(eq (second expr) array))
(list 'cffi:mem-aref array-fptr element-type (elt expr 2))
expr)) body))))
(defun macro-force-function (dim)
"Given an integer dim, this constructs a function that, when supplied with a
N-dimensional vector Z and some output vector (-> pointer?), yields the
corresponding forces"
(declare (fixnum dim))
(let ((temp-values (make-array 2 :element-type 'double-float :initial-element 0.0d0)))
(lambda (zvector output)
(with-fast-access-to-single-foreign-array (output :double)
(with-fast-access-to-single-foreign-array (zvector :double)
(do ((i 0 (1+ i))) ((= i dim)) (declare (fixnum i))
(setf (aref temp-values 0) 0.0d0)
(do ((m 0 (1+ m))) ((> m i)) (declare (fixnum m))
(do ((n i (1+ n))) ((= n dim)) (declare (fixnum n))
(setf (aref temp-values 1) 0.0d0)
(do ((k m (1+ k))) ((> k n)) (declare (fixnum k))
(incf (aref temp-values 1) (grid:gref* zvector k)))
(incf (aref temp-values 0) (expt (aref temp-values 1) -2))))
(setf (grid:gref* output i)
(- (grid:gref* zvector i)
(aref temp-values 0)))))))))