I was trying to implement the simple function: Given two 1-dimensional
single-float arrays, a and b, and three fixnums i, j and n, form the
following sum:
sum_{k=0}^{n-1} a[(i*n)+k] * k[(j*n)+k]
where the notation is mathematical (non-Lisp).
Here is a first attempt:
(defun obvious (a i b j n)
(declare (type (simple-array single-float (*)) a b)
(fixnum i j n)
(optimize (speed 3) (safety 0) (space 0) (debug 0)))
(let ((acc 0.0)
(n*i (the fixnum (* n i)))
(n*j (the fixnum (* n j))))
(declare (single-float acc)
(fixnum n*i n*j))
(do ((ak n*i (1+ ak))
(bk n*j (1+ bk)))
((>= ak (the fixnum (+ n n*i))) acc)
(declare (fixnum ak bk))
(incf acc (* (aref a ak) (aref b bk))))))
Unfortunately ACL5.0.1(Solaris) does not allocate acc in a register
even though it *seems* (via the :explain declaration) to be doing so.
Now here is the weird apart. If I define an intermediate function to
do the simpler sum
sum_{k=0}^{n-1} a[i+k] * b[j+k]
(defun intermediate (a i b j n)
(declare (type (simple-array single-float (*)) a b)
(fixnum i j n)
(optimize (speed 3) (safety 0) (debug 0) (space 0)))
(let ((acc 0.0))
(declare (single-float acc))
(do ((ak i (1+ ak))
(bk j (1+ bk)))
((>= ak (+ i n)) acc)
(declare (fixnum ak bk))
(incf acc (* (aref a ak) (aref b bk))))))
and then make the definition:
(defun nonobvious (a i b j n)
(intermediate a (* i n) b (* j n) n))
everything works out fine. When ACL5.0.1 compiles intermediate it
happily assigns acc to a register! (I tried a bizzillion other
combinations before I hit on this one.) This is completely bizarre to
me; both do loops are almost identical!!!!!
Unfortunately, I have a lot of similar code and ACL's treatment of
loops is hard to predict. Does anyone (especially from Franz) know of
a method in this (highly frustrating) madness? For example suppose I
want to implement this sum:
sum_{k=0}^{n-1} a[i + (k*da)] * b[j + (k*da)]
Then one would think that
(defun second-sum (a i da b j db n)
(declare (type (simple-array single-float (*)) a b)
(fixnum i j n da db))
(let ((acc 0.0)
(last-ak (the fixnum (+ i (* n da)))))
(declare (single-float acc)
(fixnum last-ak))
(do ((ak i (+ ak da))
(bk j (+ bk db)))
((>= ak last-ak) acc)
(declare (fixnum ak bk))
(incf acc (* (aref a ak) (aref b bk))))))
would do the trick, since it is so similar to "intermediate".
However, it doesn't. ACL does not assign acc to a register :( Any
help will be appreciated. Thanks,
--shiv--