[GSLL-devel] [Gsll-devel] Efficient access to externally generated double-float arrays?

[Sebastian Sturm]

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