|
From: | Vincent Belaïche |
Subject: | Submission of Matrix Kronecker Product for calc |
Date: | Sun, 2 Mar 2008 22:37:40 +0100 |
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 |
[Prev in Thread] | Current Thread | [Next in Thread] |