Hi everybody. Here is my routine which I use it to get an Inverse Matrix. It based on Cayley-Hamilton theorem and Boucher formula.
Routine code as below.
Thanks for your interest.
; Matrix Inverse
; based on Cayley-Hamilton theorem
; and Boucher formula
(princ "\nType mInverse to run ")
(defun c:mInverse ( / ma uma msize mp mplist spn splist fan falist
falist-1 mplist-1 faz mma rma mm i maStr okStr)
(defun *error* ( msg )
(if (and msg (not (wcmatch (strcase msg t) "*break,*cancel*,*exit*")))
(princ (strcat "\nError : " msg))
)
(princ)
) ; end *error*
(vl-load-com)
;; Matrix Transpose - Doug Wilson
;; Args: m - nxn matrix
(defun mTrp ( m )
(apply 'mapcar (cons 'list m))
) ; end mTrp
;; Matrix x Vector - Vladimir Nesterovsky
;; Args: m - nxn matrix, v - vector in R^n
(defun MxV ( m v )
(mapcar '(lambda ( r ) (apply '+ (mapcar '* r v))) m)
) ; end MxV
;; Matrix x Matrix - Vladimir Nesterovsky
;; Args: m1, m2 - nxn matrices
(defun MxM ( m1 m2 )
(mapcar (function (lambda (r) (MxV (mTrp m2) r))) m1)
;((lambda ( a ) (mapcar '(lambda ( r ) (MxV a r)) m1)) (mTrp m2))
) ; end MxM
;; matrix + matrix
;; Args: m1, m2 - nxn matrices
(defun MpM (m1 m2)
(mapcar '(lambda (a b) (mapcar '+ a b)) m1 m2)
) ; end MpM
; e.g : (mpm '((1 2) (3 4)) '((5 6) (7 8)))
; -> ((6 8) (10 12))
;; Make Vector
(defun make-vector (num n / vec)
(setq vec (list))
(repeat n
(setq vec (append vec (list num)))
)
vec
) ; end make-vector
; make vector matrix
(defun make-Vmatrix (num n / vma)
(setq vma (list))
(repeat n
(setq vma (append vma (list (make-vector num n))))
) ; end repeat
vma
) ; end make-Vmatrix
;; matrix scalar product
(defun Masp (num ma / vnum)
(setq vnum (make-vector num (length ma)))
(mapcar '(lambda (x) (mapcar '* vnum x)) ma)
) ; end masp
; e.g : (masp 2 '((1 2) (3 4)))
; -> ((2 4) (6 8))
;; make unit matrix
(defun Make-Umatrix (num n / i j nlst zlst zmatrix zrow irow Umatrix itm mastr)
(setq i 0)
(setq nlst (list))
(repeat n
(setq nlst (append nlst (list i)))
(setq i (1+ i))
)
(setq zmatrix (list))
(repeat n
(setq zlst (list))
(setq zlst (make-vector 0 n))
(setq zmatrix (append zmatrix (list zlst)))
)
(setq i 0)
(setq Umatrix (list))
(while (< i n)
(setq zrow (nth i zmatrix))
(setq j 0)
(setq irow (list))
(while (< j n)
(setq itm (nth j zrow))
(if (= j i)
(setq itm NUM)
)
(setq irow (append irow (list itm)))
(setq j (1+ j))
) ; end while
(setq Umatrix (append Umatrix (list irow)))
(setq i (1+ i))
) ; end while
Umatrix
) ; end make-Umatrix
; e.g : (make-Umatrix 1 3)
; -> ((1 0 0) (0 1 0) (0 0 1))
; sum of numbers on main cross of matrix
(defun Sum-MaCross (ma / i j n mrow numlst sum)
(setq n (length ma))
(setq numlst (list))
(setq i 0)
(while (< i n)
(setq mrow (nth i ma))
(setq j 0)
(while (< j n)
(setq itm (nth j mrow))
(if (= j i)
(setq numlst (append numlst (list itm)))
)
(setq j (1+ j))
) ; end while
(setq i (1+ i))
) ; end while
(apply '+ numlst)
) ; end Sum-MaCross
; e.g: (Sum-MaCross '((1 0 0) (0 1 0) (0 0 1)))
; -> 3
;; string to list [Lee-Mac]
(defun String-to-list ( str del / len lst pos )
(setq len (1+ (strlen del)))
(while (setq pos (vl-string-search del str))
(setq lst (cons (substr str 1 pos) lst)
str (substr str (+ pos len))
)
)
(reverse (cons str lst))
) ; end String-to-list
(defun Check-Input-String (str / lst msize i itm row ma)
(setq lst (mapcar 'atof (String-to-list str ",")))
(setq msize (fix (sqrt (length lst))))
(if (/= (expt msize 2) (length lst))
(progn
(alert (strcat "Invalid Input" "\nExit"))
(exit)
)
)
(setq ma (list) row (list) i 0)
(while lst
(setq itm (car lst))
(setq row (append row (list itm)))
(setq lst (cdr lst))
(setq i (1+ i))
(if (= (rem i msize) 0)
(progn
(setq Ma (append Ma (list Row)))
(setq row (list))
)
)
) ; end while
Ma
) ; end Check-Input-String
; e.g : (Check-Input-String "1,2,3,4")
; ->((1.0 2.0) (3.0 4.0))
;; + + + + + + + +
;; Main program
;; + + + + + + + +
(if (null maStr0)
(setq maStr0 "1,2,3,4")
)
(setq maStr (getstring (strcat "\nInput define matrix String or [eXit] <" maStr0 "> : ")))
(cond
( (= (strcase maStr) "X")
(exit)
)
( (= maStr "")
(setq maStr maStr0)
)
(T
(setq okStr maStr)
)
) ; end cond
; supplied matrix
(setq ma (Check-Input-String maStr))
; convert to floating number
(setq ma (mapcar '(lambda(x) (mapcar 'float x)) ma))
; matrix size
(setq msize (length ma))
; unit matrix 1.0
(setq uma (make-Umatrix 1.0 msize))
; a factor list
(setq falist (list))
; matrix power list
(setq mplist (list))
; specialize polynomical list
(setq splist (list))
(setq i 1)
(repeat msize
(if (null mp) ; matrix power
(setq mp ma)
(setq mp (MxM ma mp))
)
(setq mplist (append (list mp) mplist))
(setq spn (Sum-MaCross mp))
; follow Boucher formula
(if (null splist)
(setq fan (/ spn (* i -1.0)))
(setq fan (/ (+ spn (apply '+ (mapcar '* (reverse splist) falist))) (* i -1.0)))
)
(setq falist (append falist (list fan)))
(setq splist (append splist (list spn)))
(setq i (1+ i))
) ; end repeat
;; !!!
;; the end factor
;; this must be Non-Zero
(setq faz (last falist))
(if (= faz 0)
(progn
(alert (strcat "Cannot calculate"
"\nInverse Matrix"
"\n*Exit*"
)
)
(exit)
)
)
; when faz satisfied non-zero
(setq maStr0 okStr)
; add 1.0 to remain factor list
(setq falist-1 (cons 1.0 (reverse (cdr (reverse falist)))))
; add unit matrix to matrix power list
(setq mplist-1 (append (cdr mplist) (list uma)))
; multiply matrix
(setq mma (mapcar '(lambda(n m) (Masp n m))
falist-1 mplist-1
)
)
; result matrix
; unit matrix 0.0
(setq Rma (make-Vmatrix 0.0 msize))
; follow Cayley-Haminton theorem
(foreach mm mma
(setq Rma (MpM mm Rma))
) ; end foreach
(setq rma (Masp (/ -1.0 faz) rma))
(princ "[")
(mapcar 'print rma)
(princ "\n]")
(*error* nil)
(princ)
) ; end c:mInverse
; result
; inverse matrix of
; [(1.0 2.0) (3.0 4.0)] is :
; [(-2.0 1.0) (1.5 -0.5)]
