share/linearalgebra/lu.lisp

   1 ;;  Copyright 2005, 2006, 2021 by Barton Willis
   2
   3 ;;  This is free software; you can redistribute it and/or
   4 ;;  modify it under the terms of the GNU General Public License,
   5 ;;  http://www.gnu.org/copyleft/gpl.html.
   6
   7 ;; This software has NO WARRANTY, not even the implied warranty of
   8 ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
   9
  10 (in-package :maxima)
  11 ($put '$lu 3 '$version)
  12
  13 ;; Return the i,j entry of the Maxima matrix m. The rows of m have been permuted according
  14 ;; to the Maxima list p.
  15
  16 (defun m-elem (m p i j)
  17   (if (arrayp m) (aref m (aref p i) j) (nth j (nth (nth i p) m))))
  18
  19 ;; Set the i j entry of the Maxima matrix or an array to x.
  20
  21 (defun setmatelem (m x i j)
  22   (if (arrayp m) (setf (aref m i j) x) (setf (nth j (nth i m)) x)))
  23
  24 ;; Return m[perm[i], k] - sum(m[perm[i],s] * m[perm[s],k],s,0,n)
  25
  26 (defun partial-matrix-prod (m p i k n fadd fsub fmult add-id fname)
  27   (cond ((eq fname '$floatfield)
  28          (partial-matrix-prod-float m p i k n))
  29         ((eq fname '$complexfield)
  30          (partial-matrix-prod-complex-float m p i k n))
  31         (t
  32          (let ((l (aref p i)))
  33            (loop for s from 0 to n do
  34              (setq add-id (funcall fadd add-id (funcall fmult (aref m l s) (aref m (aref p s) k)))))
  35            (setf (aref m l k) (funcall fsub (aref m l k) add-id))))))
  36
  37 (defun partial-matrix-prod-float (m p i k n)
  38   (let ((add-id 0.0))
  39     (declare (type flonum add-id))
  40     (let ((l (aref p i)))
  41       (loop for s from 0 to n do
  42         (setq add-id (+ add-id (* (aref m l s) (aref m (aref p s) k)))))
  43         (setf (aref m l k) (- (aref m l k) add-id)))))
  44
  45 (defun partial-matrix-prod-complex-float (m p i k n)
  46   (let ((add-id 0.0))
  47     (let ((l (aref p i)))
  48       (loop for s from 0 to n do
  49         (setq add-id (+ add-id (* (aref m l s) (aref m (aref p s) k)))))
  50         (setf (aref m l k) (- (aref m l k) add-id)))))
  51
  52 ;; Return the infinity norm (the largest row sum) of the r by c array mat. The function
  53 ;; fn coerces matrix elements into flonum floats. The argument 'mat' is a Maxima
  54 ;; style matrix; thus mat = (($matrix) ((mlist) a b c) etc).
  55
  56 (defun array-infinity-norm (mat fn)
  57   (reduce #'max (mapcar #'(lambda (s) (reduce #'+ s))
  58                         (array-to-row-list mat #'(lambda (s) (abs (funcall fn s)))))))
  59
  60 (defun $lu_factor (mat &optional (fld '$generalring))
  61   ($require_nonempty_matrix mat '$first '$lu_factor)
  62   ($require_square_matrix mat '$first '$lu_factor)
  63   (setq fld ($require_ring fld '$second '$lu_factor))
  64
  65   (let* ((c ($length mat)) (perm (make-array c)) (cnd) (fn))
  66     ;(setq mat (full-matrix-map (mring-maxima-to-mring fld) mat))
  67     (setq mat (maxima-to-array mat (mring-maxima-to-mring fld)))
  68
  69     (decf c)
  70     (loop for i from 0 to c do (setf (aref perm i) i))
  71
  72     ;; When the matrix elements can be converted to CL floats, find
  73     ;; the infinity norm of m. The norm is used to bound the condition
  74     ;; number.
  75
  76     (cond ((not (eq nil (mring-coerce-to-lisp-float fld)))
  77            (setq fn (mring-coerce-to-lisp-float fld))
  78            (setq cnd (array-infinity-norm mat fn)))
  79           (t (setq cnd 0.0)))
  80     (lu-factor mat perm c fld cnd)))
  81
  82 (defun $get_lu_factors (x)
  83   (let ((mat ($first x)) (mp) (p ($second x)) (perm) (r) (c) (id) (lower) (upper) (zero))
  84
  85     (setq r ($matrix_size mat))
  86     (setq c ($second r))
  87     (setq r ($first r))
  88     (setq id ($args ($identfor mat)))
  89     (loop for i from 1 to c do
  90       (push (nth (nth i p) id) perm)
  91       (push (nth (nth i p) mat) mp))
  92     (setq perm (reverse perm))
  93     (setq mp (reverse mp))
  94     (push '($matrix) perm)
  95     (push '($matrix) mp)
  96
  97     (setq lower (copy-tree mp))
  98     (setq upper (copy-tree mp))
  99     (setq id ($identfor ($first ($first mat))))
 100     (setq zero ($zerofor ($first ($first mat))))
 101     (loop for i from 1 to r do
 102       (loop for j from 1 to c do
 103         (if (= i j) (setmatelem lower id i j))
 104         (if (> i j) (setmatelem upper zero i j))
 105         (if (< i j) (setmatelem lower zero i j))))
 106     `((mlist) ,($transpose perm) ,lower ,upper)))
 107
 108 (defun lu-factor (m perm c fld &optional (cnd 1.0))
 109   (let ((pos) (kp1) (mx) (lb) (ub) (save) (fname (mring-name fld))
 110         (add-id (funcall (mring-add-id fld))))
 111     (flet
 112         ((fzerop (a) (funcall (mring-fzerop fld) a))
 113          (fadd (a b) (funcall (mring-add fld) a b))
 114          (fsub (a b) (funcall (mring-sub fld) a b))
 115          (fabs (a) (funcall (mring-abs fld) a))
 116          (fmult (a b) (funcall (mring-mult fld) a b))
 117          (fdiv (a b) (funcall (mring-div fld) a b))
 118          (fgreat (a b) (funcall (mring-great fld) a b)))
 119
 120   (loop for k from 0 to c do
 121           (loop for i from k to c  do (partial-matrix-prod m perm i k (- k 1)
 122                                                          #'fadd #'fsub #'fmult add-id fname))
 123           (setq mx (fabs (m-elem m perm k k)))
 124           (setq pos k)
 125           (loop for s from k to c do
 126             (if (fgreat (fabs (m-elem m perm s k)) mx) (setq pos s mx (fabs (m-elem m perm s k)))))
 127
 128         (setq save (aref perm k))
 129           (setf (aref perm k) (aref perm pos))
 130           (setf (aref perm pos) save)
 131
 132         (setq kp1 (+ 1 k))
 133           (loop for i from kp1 to c do
 134               (when (not (fzerop (m-elem m perm k k)))
 135                   (setmatelem m (fdiv (m-elem m perm i k) (m-elem m perm k k)) (aref perm i) k)
 136                   (partial-matrix-prod m perm k i (- k 1) #'fadd #'fsub #'fmult add-id fname))))
 137
 138       (cond ((not (eq nil (mring-coerce-to-lisp-float fld)))
 139              (multiple-value-setq (lb ub) (mat-cond-by-lu m perm c (mring-coerce-to-lisp-float fld)))
 140              (setq m (array-to-maxima-matrix m (mring-mring-to-maxima fld)))
 141              (setq lb ($limit (mul lb cnd)))
 142              (setq ub ($limit (mul ub cnd)))
 143              (setq perm (array-to-maxima-list perm #'(lambda (s) (+ s 1))))
 144              `((mlist) ,m ,perm ,(mring-name fld) ,lb ,ub))
 145             (t
 146              (setq perm (array-to-maxima-list perm #'(lambda (s) (+ s 1))))
 147              (setq m (array-to-maxima-matrix m (mring-mring-to-maxima fld)))
 148              `((mlist) ,m ,perm ,(mring-name fld)))))))
 149
 150 ;; The first argument should be a matrix in packed LU form. The Maxima list perm
 151 ;; specifies the true row ordering. When r is false, reflect the matrix horizontally
 152 ;; and vertically.
 153
 154 (defun m-elem-reflect (mat perm n r i j)
 155   (cond ((and r (= i j)) 1)
 156         ;;(r (m-elem mat perm (+ n (- i) 1) (+ n (- j) 1)))
 157         (r (m-elem mat perm (+ n (- i)) (+ n (- j))))
 158         (t (m-elem mat perm i j))))
 159
 160 ;; The first argument mat should be a matrix in the packed LU format.
 161 ;; When l-or-u is lower, return the i,j element of the lower triangular
 162 ;; factor; when l-or-u is upper, return the j, i element of upper triangular
 163 ;; factor. The first argument mat should be a matrix in the packed LU format.
 164
 165 (defun mat-cond-by-lu (mat perm n fn)
 166   (let ((lb0) (ub0) (lb1) (ub1))
 167     (multiple-value-setq (lb0 ub0) (triangular-mat-cond mat perm n fn t))
 168     (multiple-value-setq (lb1 ub1) (triangular-mat-cond mat perm n fn nil))
 169     (values ($limit (mul lb0 lb1)) ($limit (mul ub0 ub1)))))
 170
 171 ;; Return lower and upper bounds for the infinity norm condition number of the lower or
 172 ;; upper triangular part of the matrix mat. The function fn coerces the matrix
 173 ;; elements to flonum floats.  When the matrix is singular, return infinity.
 174 ;; This code is based on pseudo-code (algorithm 2.1) in ``Survey of condition
 175 ;; number estimation,'' by Nicholas J. Higham, SIAM Review, Vol. 29, No. 4, December,
 176 ;; 1987.  The lower and upper bounds can differ from the true value by arbitrarily
 177 ;; large factors.
 178
 179 (defun triangular-mat-cond (mat perm n fn r)
 180   (let ((z) (d-max 0.0) (z-max 0.0) (s) (d))
 181     (setq z (make-array (+ 1 n)))
 182
 183     (catch 'singular
 184       (loop for i from n downto 0 do
 185         (setq d (abs (funcall fn (m-elem-reflect mat perm n r i i))))
 186         (if (= 0.0 d) (throw 'singular (values '$inf '$inf)) (setq d (/ 1 d)))
 187         (setq d-max (max d-max (abs (funcall fn d))))
 188         (setq s 1.0)
 189         (loop for j from (+ 1 i) to n do
 190           (incf s (* (abs (funcall fn (m-elem-reflect mat perm n r i j))) (aref z j))))
 191         (setf (aref z i) (* d s))
 192         (setq z-max (max z-max (aref z i))))
 193       (values d-max z-max))))
 194
 195 (defun $lu_backsub(m b1)
 196   ($require_list m '$first '$lu_backsub)
 197   (when (< ($length m) 3)
 198       (merror (intl:gettext "The first argument to 'lu_backsub' must be a list with at least 3 members")))
 199
 200   (let* ((mat) (n) (r) (c) (bb) (acc) (perm) (id-perm) (b)
 201          (fld (get (mfuncall '$ev (fourth m)) 'ring)) (cc)
 202          (fadd (mring-add fld))
 203          (fsub (mring-sub fld))
 204          (fmult (mring-mult fld))
 205          (fdiv (mring-rdiv fld))
 206          (add-id (funcall (mring-add-id fld))))
 207
 208     (setq mat (copy-tree ($first m)))
 209     (setq perm ($second m))
 210     (setq n ($matrix_size mat))
 211     (setq r ($first n))
 212     (setq c ($second n))
 213
 214     (setq mat (matrix-map (mring-maxima-to-mring fld) mat))
 215     (setq b (copy-tree b1))
 216     (setq c ($second ($matrix_size mat)))
 217
 218     (setq cc ($second ($matrix_size b)))
 219     (setq b (matrix-map (mring-maxima-to-mring fld) b))
 220
 221     (setq bb (copy-tree b))
 222     (loop for i from 1 to r do
 223       (loop for j from 1 to cc do
 224               (setmatelem bb (m-elem b perm i j) i j)))
 225
 226     (setq id-perm (mfuncall '$makelist 'i 'i 1 r))
 227
 228     (loop for q from 1 to cc do
 229       (loop for i from 2 to c do
 230         (setq acc add-id)
 231         (loop for j from 1 to (- i 1) do
 232           (setq acc (funcall fadd acc (funcall fmult (m-elem mat perm i j) (m-elem bb id-perm j q)))))
 233         (setmatelem bb (funcall fsub (m-elem bb id-perm i q) acc) i q)))
 234
 235     (loop for q from 1 to cc do
 236       (loop for i from r downto 1 do
 237         (setq acc (m-elem bb id-perm i q))
 238         (loop for j from (+ 1 i) to c do
 239           (setq acc (funcall fsub acc (funcall fmult (m-elem mat perm i j) (m-elem bb id-perm j q)))))
 240         (setmatelem bb (funcall fdiv acc (m-elem mat perm i i)) i q)))
 241
 242     (setq bb (matrix-map (mring-mring-to-maxima fld) bb))
 243     bb))
 244
 245 (defun $invert_by_lu (m  &optional (fld '$generalring))
 246   ($require_square_matrix m '$first '$invert_by_lu)
 247   ($lu_backsub ($lu_factor m fld) ($identfor m)))
 248
 249 ;; Return a Lisp list of two elements, the determinant, and the inverse of M.
 250 (defun invert-by-lu-with-determinant (m fld-name)
 251   (let*
 252    ((i ($identfor m))
 253     (m1 ($lu_factor m fld-name))
 254     (fld (get fld-name 'ring))
 255     (d (determinant-by-lu-factors m1 fld)))
 256    (list d ($lu_backsub m1 i))))
 257
 258 (defun $determinant_by_lu (m &optional (fld-name '$generalring))
 259   ($require_square_matrix m '$first '$determinant_by_lu)
 260
 261   (let ((fld ($require_ring fld-name '$second '$determinant_by_lu)))
 262     (setq m ($lu_factor m fld-name))
 263     (determinant-by-lu-factors m fld)))
 264
 265 ;; Assume that M has already been factored by $LU_FACTOR
 266 ;; and FLD is some field (not a field name).
 267 (defun determinant-by-lu-factors (m fld)
 268
 269   (let* ((acc (funcall (mring-mult-id fld)))
 270          (fmult (mring-mult fld))
 271          (fconvert (mring-maxima-to-mring fld))
 272          (n ($first ($matrix_size ($first m))))
 273          (perm) (d))
 274
 275     (setq perm ($second m))
 276
 277     (setq m ($first m))
 278     (loop for i from 1 to n do
 279       (setq d (funcall fconvert (m-elem m perm i i)))
 280       ;;(if ($matrixp d) (setq d ($determinant_by_lu d fld)))
 281       (setq acc (funcall fmult acc d)))
 282     (multiple-value-bind (perm-sorted perm-sign) (bbsort1 (cdr perm))
 283       (declare (ignore perm-sorted))
 284       (funcall (mring-mring-to-maxima fld)
 285                (if perm-sign (funcall (mring-negate fld) acc) acc)))))
 286
 287 (defun $mat_cond (m p)
 288   ($require_square_matrix m '$first '$mat_cond)
 289   (mul (mfuncall '$mat_norm m p) (mfuncall '$mat_norm ($invert_by_lu m) p)))
 290
 291 (defun $linsolve_by_lu (m b &optional (fld '$generalring))
 292   ($require_square_matrix m '$first '$linsolve_by_lu)
 293   (setq b (if ($listp b) ($transpose b) b))
 294   ($require_matrix b '$second '$linsolve_by_lu)
 295   ($require_ring fld '$third '$linsolve_by_lu)
 296   (if (= ($second ($matrix_size m)) ($first ($matrix_size b))) t
 297     (merror (intl:gettext "Incompatible matrix sizes")))
 298
 299   (setq m ($lu_factor m fld))
 300   `((mlist) ,($lu_backsub m b) ,(if (floatp ($last m)) ($last m) nil)))
 301
 302 (defun invertible-matrixp (mat fld)
 303   (let ((OK t) (n))
 304       (cond (($matrixp mat)
 305               (setq n ($matrix_size mat))
 306               (cond ((not (eql ($first n) ($second n))) ;nonsquare matrices aren't invertible
 307                  nil)
 308               (t
 309                 (setq mat (fourth ($get_lu_factors ($lu_factor mat fld))))
 310                 (setq n ($first n))
 311                 (while (and OK (> n  0))
 312                   (setq OK (and OK (not ($zeromatrixp ($inpart mat n n)))))
 313                  (decf n 1))
 314                 OK)))
 315             (t (not ($zeromatrixp mat))))))