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