(defun mmul (A B)
(let* ((m (car (array-dimensions A)))
(n (cadr (array-dimensions A)))
(l (cadr (array-dimensions B)))
(C (make-array `(,m ,l) :initial-element 0)))
(loop for i from 0 to (- m 1) do
(loop for k from 0 to (- l 1) do
(setf (aref C i k)
(loop for j from 0 to (- n 1)
sum (* (aref A i j)
(aref B j k))))))
C))
;; Cofactor i j
(defun matrix-cofactor (A i j)
(let* (
(tmp 0)
(n (car (array-dimensions A)))
(ans (make-array `(,(- n 1) ,(- n 1)) :initial-element 0)))
(loop for ii from 0 to (- i 1) do
(loop for jj from 0 to (- j 1) do (setf (aref ans ii jj) (aref A ii jj)))
(loop for jj from (+ j 1) to (- n 1) do (setf (aref ans ii (- jj 1)) (aref A ii jj))))
(loop for ii from (+ i 1) to (- n 1) do
(loop for jj from 0 to (- j 1) do (setf (aref ans (- ii 1) jj) (aref A ii jj)))
(loop for jj from (+ j 1) to (- n 1) do (setf (aref ans (- ii 1) (- jj 1)) (aref A ii jj))))
ans))
;; Calculate det value of matrix A
(defun matrix-det (A)
(let* ((n (car (array-dimensions A))) (tmp 0))
(if (= n 1)
(aref A 0 0)
(let ((tmp 0) (MA 0))
(loop for j from 0 to (- n 1) do
(setf MA (matrix-cofactor A 0 j))
(if (= (mod j 2) 0)
(setf tmp (+ tmp (* (aref A 0 j) (matrix-det MA))))
(setf tmp (- tmp (* (aref A 0 j) (matrix-det MA))))))
tmp))))
;; Calculate inverse matrix of matrix A
(defun matrix-inverse (A)
(let* ((n (car (array-dimensions A))) (tmp 0) (inv-A (make-array `(,n ,n) :initial-element 0)) (det-val 0))
(setf det-val (matrix-det A))
(loop for i from 0 to (- n 1) do
(loop for j from 0 to (- n 1) do
(setf (aref inv-A j i) (/ (matrix-det (matrix-cofactor A i j)) det-val))))
inv-A))
(let ((g 0) (g0 0))
(setf g (make-array '(3 3) :initial-contents '((1 0 0) (2 2 0) (3 4 5))))
(setf g0 (make-array '(1 1) :initial-contents '((2))))
; (matrix-cofactor g 1 1)
(format t "~A~%" (matrix-inverse g)))