Submission of Matrix Kronecker Product for calc

["Vincent Belaïche" <vincent.b.1@hotmail.fr>]

[-- Attachment #1: Type: text/plain, Size: 4514 bytes --]


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






_________________________________________________________________
Nouveau ! Créez votre profil Messenger !
http://home.services.spaces.live.com/

[-- Attachment #2: Type: text/html, Size: 4979 bytes --]