src/invert.lisp

   1 ;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
   2
   3 (in-package :maxima)
   4
   5 (defmfun $adjoint (mat)
   6   (let* ((n ($length mat))
   7          (adj (simplify ($ident n))))
   8     (unless (like n 1)
   9       (do ((i 1 (1+ i)))
  10           ((> i n))
  11         (do ((j 1 (1+ j)))
  12             ((> j n))
  13           (maset (mul* (power -1 (+ i j))
  14                        (simplify ($determinant (simplify ($minor mat j i)))))
  15                  adj i j))))
  16     adj))
  17
  18 (defmfun $invert_by_adjoint (mat)
  19   (let*
  20     ((adj (simplify ($adjoint mat)))
  21      (det (let (($scalarmatrixp t))
  22                (ncmul2 (simplify ($row mat 1))
  23                        (simplify ($col adj 1)))))
  24      (mat1 (if (and $scalarmatrixp (= ($length mat) 1)) (maref adj 1 1) adj)))
  25     (if $detout
  26       `((mtimes) ((mexpt) ,det -1) ,mat1)
  27       (div mat1 det))))
  28
  29 (defmvar $invert_method nil)
  30 (defmvar $invert_by_adjoint_size_limit 8)
  31
  32 (defmfun $invert (&rest args)
  33   (case $invert_method
  34     (($lu) (apply #'invert-via-$invert_by_lu args))
  35     (($gausselim) (apply #'$invert_by_gausselim args))
  36     (($adjoint) (apply #'$invert_by_adjoint args))
  37     ((nil)
  38       ;; Select a method appropriate for the matrix.
  39       ;; This could be more sophisticated.
  40       (let*
  41         ((my-matrix (first args))
  42          (size (length (rest my-matrix))))
  43         (if (<= size $invert_by_adjoint_size_limit)
  44           (apply #'$invert_by_adjoint args)
  45           (apply #'$invert_by_gausselim args))))
  46     (t
  47       (mtell "invert: unrecognized invert_method=~M; assume default.~%" $invert_method)
  48       (let (($invert_method nil))
  49         (apply #'$invert args)))))
  50
  51 (defun invert-via-$invert_by_lu (m &optional (field-name (if $ratmx '$crering '$generalring)))
  52   (declare (special $ratmx $detout))
  53   ;; Call functions from package linearalgebra via MFUNCALL to autoload them if necessary.
  54   (if $detout
  55     (let*
  56       ((field (mfuncall '$require_ring field-name "$second" "$invert"))
  57        (d-i (funcall 'invert-by-lu-with-determinant m field-name))
  58        (d (first d-i))
  59        (i (second d-i))
  60        (d-times-i (multiply-matrix-elements d (funcall 'mring-mult field) i))
  61        (d^-1 (funcall (funcall 'mring-reciprocal field) d)))
  62       (list '(mtimes) d^-1 d-times-i))
  63     (mfuncall '$invert_by_lu m field-name)))
  64
  65 ;; I wonder if this function already exists somewhere. Oh well.
  66 (defun multiply-matrix-elements (a multiply m)
  67   (cons (car m) (mapcar #'(lambda (row) (cons (car row) (mapcar #'(lambda (x) (funcall multiply a x)) (cdr row)))) (cdr m))))