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)))))))))