Submission of Matrix Kronecker Product for calc

[Vincent Belaïche]

Hello,

I am new to this list (I have just subscribed). I have just written a
function for CALC that computes the Kronecker product of two matrices.

Maybe this can be useful to somebody else. I am not sure that this is
the correct forum to submit that. Sorry if I bothered anybody.

BR,
Vincent.

PS-1 : I think that this is perfectible to the following extent
a) I use reverse instead of nreverse, which is more memory consuming.
Not knowing the internals of Emacs I was unsure how the garbage
collector free the memory if you nreverse a list not passing to nreverse
the pointer of the first element (ie the point to the list itself). This
is why I used reverse.

b) I started from a template created by the "ZP" key for creating a
permanent user function. The starting point was a simple multiplication.
Then I modified the function to make a Kronecker product. However I
don't know how to set the 'calc-user-defn property (what is there is
still what there was originally after "ZP" keystroke in calc mode.


PS-2 : The function works on two dimensional matrices. I had no need to
make it generic for n-dimensional (and maybe I was not proficient enough
in Lisp to do it... ;-/ ).

PS-3 : Code below the --- line :
--------------------------------------------------
;;; Definition of Kronecker matrix product
(put 'calc-define 'calc-kron '(progn
(defun calc-kron nil (interactive) (calc-wrapper (calc-enter-result 2
"kron" (cons (quote calcFunc-kron) (calc-top-list-n 2)))))
(put 'calc-kron 'calc-user-defn 't)
(defun calcFunc-kron (x y)
"Kronecker product of matrices x and y.
If x and y are not matrices then they are first embedded into matrices
before computation.
After computation the result may be desembedded from matrix so that the
Kronecker product
of two scalars is a scalar,
of one scalar and a vector or vice verse is a vector, or
of two vectors is a vector.
"
;; (matrix-backet-level x) returns 0 if x is not a matrix,
;; 1 if x is a vector that is not a
matrix
;; 2 if x is a matrix
(defun matrix-backet-level (x)
(let ((returned-value 0) (x-sub x))
(while
(and (Math-vectorp x-sub) (< returned-value 2))
(setq x-sub (elt x-sub 1))
(setq returned-value (1+ returned-value))
)
returned-value
))
;; (matrix-embed x n) return x embeded in n level of brackets
;; (matrix-embed x 0) returns x
;; (matrix-embed x 1) returns (vec x)
;; (matrix-embed x 2) returns (vec (vec x))
(defun matrix-embed (x n)
(let ((returned-value x) (n-ctr 0))
(while
(< n-ctr n)
(setq returned-value (list 'vec returned-value))
(setq n-ctr (1+ n-ctr))
)
returned-value
))
;; evaluate level of bracketting of arguments x and y
(let
(
(bl-x (matrix-backet-level x))
(bl-y (matrix-backet-level y))
prod-val
)
(let
(
(emb-x (matrix-embed x (- 2 bl-x)))
(emb-y (matrix-embed y (- 2 bl-y)))
prod-next-row
)
(dolist (x-row (reverse (cdr emb-x)) prod-val)
(dolist (y-row (reverse (cdr emb-y)))
(setq prod-next-row nil)
(dolist (x-col (reverse (cdr x-row)) prod-next-row)
(dolist (y-col (reverse (cdr y-row)))
(setq prod-next-row (cons (math-mul x-col y-col)
prod-next-row))
)
)
(setq prod-next-row (cons 'vec prod-next-row))
(setq prod-val (cons prod-next-row prod-val))
)
)
(setq prod-val (cons 'vec prod-val))
)
(cond
((and (eq bl-x 0) (eq bl-y 0)) (elt (elt prod-val 1) 1))
((and (< bl-x 2) (< bl-y 2)) (elt prod-val 1))
(t prod-val)
)
)
)
;; I don't know how to set this property ... this statement does not
implement the Kronecker product
(put 'calcFunc-kron 'calc-user-defn '(* (var x var-x) (var y var-y)))
(define-key calc-mode-map "zk" 'calc-kron)
))

---

Windows Live Messenger 2008 vient de sortir, encore plus de fun ! Téléchargez gratuitement Messenger 2008