Hi,
I'm having a problem with a macro and CLISP. The code works in
CMU, but CLISP refuses it, and I'm not confident enough to know who/what is
right.
The problem relates to iterating over banded matrices. I've defined
some routines for manipulating band diagonal matrices - matrices for
which all elements are zero except for those that are within a certain
distance of the main diagonal e.g
1 2 3 0 0 0
4 5 6 7 0 0
0 8 9 1 2 0
0 0 3 4 5 6
0 0 0 7 8 9
is a 5x6 band diagonal matrix, with upper bandwidth (number of diags
above main diagonal) of 2 and a lower bandwidth of 1. To store the
matrix you need only store the non zero diags, so the above is stored
as
_ 1 2 3
4 5 6 7
8 9 1 2
3 4 5 6
7 8 9 _
where underscores denote wasted elements. There is then a simple
translation scheme from the indexing of one storage scheme to the
other. I've implemented a simple macro bmref that does precisely
that.
So given a banded matrix A, i want to be able to loop over only the
non zero elements. So for the example, rather than looping over the m
by n matrix A as
(loop for i fixnum from 1 upto m do
(loop for j fixnum from 1 upto n do
(bmref A i j)))
I need
(loop for i fixnum from 1 upto (min m (+ n lb)) do
(loop for j fixnum from (max 1 (- i lb) upto (min n (+ i ub)) do
(bmref A i j)))
where lb and ub represent the upper and lower bandwidths of A. This
way you only loop over those elements that have actually been stored.
So I wanted a macro for which you could say
(iterate-banded-matrix (A m n lb ub)
... code ...
(loop-rows-A (i 1 upto m)
... more code ...
(loop-cols-A (j 1 upto n)
... more code ...)
... more code ...)
... more code...)
and the "upto"s could be "downto"s, and you could also swap the order
of the loop-rows-A and loop-cols-A in a sensible fashion. The A is
appended to the loop names so that you could have nested calls to
iterate-banded-matrices for different matrices.
My solution was to define a macro iterate-banded-matrix that defines
loop-cols-A and loop-rows-A in a macrolet. When loop-rows-A is
expanded, it shadows loop-cols-A (again with a macrolet) so that
loop-cols-A expands correctly, since the expansion of loop-cols-A
depends on whether we are inside or outside of loop-rows-A.
loop-cols-A is defined similarly, shadowing loop-rows-a with a
macrolet.
Thats the idea. The following implementation works in CMU, but CLISP
refuses the ,@(mk-loop-range ...) in the innermost macrolets. So is
this OK and CLISP has a bug or have I made a subtle error? Secondly,
is there a better way?
Thanks in advance,
Sp
(defun mk-loop-range (dirn min max)
(case dirn
(upto `(,min upto ,max))
(downto `(,max downto ,min))))
(defmacro iterate-banded-matrix ((A &optional m n lb ub) &body body)
"Iterate over the elements of a banded matrix."
(let ((row-macro (intern (concatenate 'string "LOOP-ROWS-" (symbol-name A))))
(col-macro (intern (concatenate 'string "LOOP-COLS-" (symbol-name A)))))
`(let ((,m (banded-matrix-m ,A))
(,n (banded-matrix-n ,A))
(,lb (banded-matrix-lower-bandwidth ,A))
(,ub (banded-matrix-upper-bandwidth ,A)))
(macrolet
((,row-macro ((i i-min i-dirn i-max) &body row-body)
`(macrolet
((,',col-macro ((j j-min j-dirn j-max) &body col-body)
`(loop for ,j fixnum from
,@(build-loop-range
j-dirn
`(max ,j-min (- ,',i ,',',lb))
`(min ,j-max (+ ,',i ,',',ub)))
do
,@col-body)))
(loop for ,i fixnum from
,@(build-loop-range i-dirn i-min `(min ,i-max (+ ,',n ,',lb)))
do
,@row-body)))
(,col-macro ((j j-min j-dirn j-max) &body col-body)
`(macrolet
((,',row-macro ((i i-min i-dirn i-max) &body row-body)
`(loop for ,i fixnum from
,@(build-loop-range
i-dirn
`(max ,i-min (- ,',j ,',',ub))
`(min ,i-max (+ ,',j ,',',lb)))
do
,@row-body)))
(loop for ,j fixnum from
,@(build-loop-range j-dirn j-min `(min ,j-max (+ ,',m ,',ub)))
do
,@col-body))))
,@body))))
Other stuff you might want to look at.
;;;------------------------------------------------------------------------------
;;;------------------------------------------------------------------------------
;;;Implementation - Banded matrices.
;;;
;;;Banded matrices are implemented as vectors of diagonals of length
;;;min(m, n+lb), where m,n are the dimesnions of the matrix, and lb
;;;and ub are the upper and lower bandwidths. Thus the matrix
;;;
;;;
;;; a b 0 0
;;; c d e 0
;;; 0 f g h
;;; 0 0 i j
;;;
;;;is stored as
;;;
;;; _ a b
;;; c d e
;;; f g h
;;; i j _
;;;
;;;where underscores denote unused locations. The unused locations
;;;could be eleminated, but at the expense of complicating the
;;;indexing. Under the current system,
;;;
;;; a(i,j) => i-1 th element of the lb+j-i th vector
;;;
(defvar *band-matrix-max-print-dimension* 25
"Banded matrices with dimension bigger than this are printed in compact form.")
(deftype double-vector ()
"A vector of double floats."
`(simple-array double-float (*)))
(defstruct (banded-matrix
(:print-function print-banded-matrix))
(m 0 :type fixnum)
(n 0 :type fixnum)
(lower-bandwidth 0 :type fixnum)
(upper-bandwidth 0 :type fixnum)
(diagonals (make-array 0 :element-type 'double-vector)
:type (simple-array double-vector (*))))
;;;------------------------------------------------------------------------------
;;;------------------------------------------------------------------------------
;;;Macros.
;;;
;;;We provide a macro BMREF for accessing the elements of a banded
;;;matrix. This macro just translates from the standard matrix
;;;indexing syste, to the sparse format, it does no error checking.
;;;
(defmacro bmref (A i j &optional lb)
`(aref
(aref (banded-matrix-diagonals ,A)
(+ (- ,j ,i) ,(or lb `(banded-matrix-lower-bandwidth ,A))))
(1- ,i)))
;;;------------------------------------------------------------------------------
;;;------------------------------------------------------------------------------
;;;Printing
;;;
;;;We print in one of two forms depending on the size of the matix.
;;;Compact form just show dims and upper and lower bandwidth, while
;;;the full form display the elements.
(defun print-banded-matrix (A &optional (strm t) k)
(declare (type banded-matrix A) (stream strm) (ignore k))
(let ((m (banded-matrix-m A))
(n (banded-matrix-n A))
(lb (banded-matrix-lower-bandwidth A))
(ub (banded-matrix-upper-bandwidth A)))
(declare (fixnum m n lb ub))
(if (< (max m n) *band-matrix-max-print-dimension*)
(loop for i fixnum from 1 upto m do
(format strm "~%")
(loop for j fixnum from 1 upto n do
(format strm "~10,4G"
(if (and (<= (- j i) ub) (<= (- i j) lb))
(bmref A i j)
0.0D0))))
(format strm "<Banded Matrix (~S ~S) (~S ~S)>" m n lb ub))))
;;;------------------------------------------------------------------------------
;;;------------------------------------------------------------------------------
;;;Constructors.
;;;
;;;Functions for constructing banded matrices of various types
(defun banded-matrix (m n lb ub &optional diagonal-elements)
"Construct an m by n zero matrix with lower bandwidth lb and upper
bandwidth ub."
(declare (fixnum m n lb ub))
(let ((diags (make-array (+ 1 lb ub) :element-type 'double-vector))
(len (min m (+ n lb))))
(declare (fixnum len)
(type (simple-array double-vector (*)) diags))
(dotimes (i (+ 1 lb ub))
(setf (aref diags i)
(make-array len :element-type 'double-float
:initial-element (etypecase diagonal-elements
(null 0.0D0)
(real (coerce diagonal-elements 'double-float))
(vector (coerce (aref diagonal-elements i) 'double-float))))))
(make-banded-matrix :m m :n n
:lower-bandwidth lb :upper-bandwidth ub
:diagonals diags)))
--
"You need as many clues as you can get as to how these things work
when you're a buffoon." - D. Row