src/mat.lisp

   1 ;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
   2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   3 ;;;     The data in this file contains enhancements.                   ;;;;;
   4 ;;;                                                                    ;;;;;
   5 ;;;  Copyright (c) 1984,1987 by William Schelter,University of Texas   ;;;;;
   6 ;;;     All rights reserved                                            ;;;;;
   7 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
   8 ;;;     (c) Copyright 1982 Massachusetts Institute of Technology         ;;;
   9 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  10
  11 (in-package :maxima)
  12
  13 (macsyma-module mat)
  14
  15 ;; this is the mat package
  16
  17 (declare-top (special *ech* *tri*
  18                       mul* *det*
  19                       xm* xn* ax))
  20
  21 ;;these are arrays.
  22 (defvar *row*)
  23 (defvar *col*)
  24 (defvar *colinv*)
  25
  26
  27 (defun solcoef (m *c varl flag)
  28   (prog (cc answer leftover)
  29      (setq cc (cdr (ratrep* *c)))
  30      (if (or (atom (car cc))
  31              (not (equal (cdar cc) '(1 1)))
  32              (not (equal 1 (cdr cc))))
  33      ;; NOTE TO TRANSLATORS: NOT CLEAR WHAT IS "UNACCEPTABLE" HERE
  34          (merror (intl:gettext "solve: unacceptable variable: ~M") *c))
  35      (setq answer (ratreduce (prodcoef (car cc) (car m)) (cdr m)))
  36      (if (not flag) (return answer))
  37      (setq leftover
  38            (rdis (ratplus m (rattimes (ratminus answer) cc t))))
  39      (if (or (not (freeof *c leftover))
  40              (dependsall (rdis answer) varl))
  41          (rat-error "`non-linear'"))
  42      (return answer)))
  43
  44 (defun formx (flag nam eql varl)
  45   (prog (b ax x ix j)
  46      (setq varlist varl)
  47      (mapc #'newvar eql)
  48      (and (not $algebraic)
  49           (some #'algp varlist)
  50           (setq $algebraic t))
  51      (setf (symbol-value nam) (make-array (list (1+ (setq xn* (length eql)))
  52                                                 (1+ (setq xm* (1+ (length varl)))))))
  53      (setq nam (get-array-pointer nam))
  54      (setq ix 0)
  55      loop1
  56      (cond ((null eql) (return  varlist)))
  57      (setq ax (car eql))
  58      (setq eql (cdr eql))
  59      (incf ix)
  60      (setf (aref nam ix xm*) (const ax varl))
  61      (setq j 0)
  62      (setq b varl) (setq ax (cdr (ratrep* ax)))
  63      loop2
  64      (setq x (car b))
  65      (setq b (cdr b))
  66      (incf j)
  67      (setf (aref nam ix j) (solcoef ax x varl flag))
  68      (cond (b (go loop2)))
  69      (go loop1)))
  70
  71 (defun dependsall (exp l)
  72   (cond ((null l) nil)
  73         ((or (not (freeof (car l) exp)) (dependsall exp (cdr l))) t)
  74         (t nil)))
  75
  76 (setq *det* nil *ech* nil *tri* nil)
  77
  78 (defun ptorat (ax m n)
  79   (prog (i j)
  80      (setq ax (get-array-pointer ax))
  81      (setq i (1+ m))
  82      (incf n)
  83      loop1
  84      (when (equal i 1) (return nil))
  85      (decf i)
  86      (setq j n)
  87      loop2
  88      (when (equal j 1) (go loop1))
  89      (decf j)
  90      (setf (aref ax i j) (cons (aref ax i j) 1))
  91      (go loop2)))
  92
  93 (defun meqhk (z)
  94   "If Z is of the form lhs = rhs, return the expression lhs - rhs.
  95   Otherwise, return Z unchanged."
  96   (if (and (not (atom z)) (eq (caar z) 'mequal))
  97       (simplus (list '(mplus) (cadr z) (list '(mtimes) -1 (caddr z))) 1 nil)
  98       z))
  99
 100 (defun const (e varl)
 101   (setq varl (mapcar #'(lambda(x) (caadr (ratrep* x))) varl))
 102   (setq e (cdr (ratrep* e)))
 103   (let ((zl (make-list (length varl) :initial-element 0)))
 104     (ratreduce (pctimes -1 (pcsubsty zl varl (car e)))
 105                (pcsubsty zl varl (cdr e)))))
 106
 107 (defvar *mosesflag nil)
 108
 109 (defun make-param ()
 110   (let ((param (intern (format nil "~A~D" '$%r (incf $%rnum)))))
 111     (tuchus $%rnum_list param)
 112     param))
 113
 114 (defun ith (x n)
 115   (if (atom x) nil (nth (1- n) x)))
 116
 117 (defun polyize (ax r m mul)
 118   (declare (fixnum m))
 119   (do ((c 1 (1+ c)) (d))
 120       ((> c m) nil)
 121     (declare (fixnum c))
 122     (setq d (aref ax r c))
 123     (setq d (cond ((equal mul 1) (car d))
 124                   (t (ptimes (car d)
 125                              (pquotientchk mul (cdr d))))))
 126     (setf (aref ax r c) (if $sparse (cons d 1) d))))
 127
 128 ;; TWO-STEP FRACTION-FREE GAUSSIAN ELIMINATION ROUTINE
 129
 130 (defun tfgeli (ax n m &aux ($sparse (and $sparse (or *det* *inv*))))
 131   ;;$sparse is also controlling whether polyize stores polys or ratforms
 132   (setq ax (get-array-pointer ax))
 133   (setq mul* 1)
 134   (do ((r 1 (1+ r)))
 135       ((> r n) (cond ((and $sparse *det*)(sprdet ax n))
 136                      ((and *inv* $sparse)(newinv ax n m))
 137                      (t (tfgeli1 ax n m))))
 138     (do ((c 1 (1+ c))
 139          (d)
 140          (mul 1))
 141         ((> c m)
 142          (and *det* (setq mul* (ptimes mul* mul)))
 143          (polyize ax r m mul))
 144       (cond ((equal 1 (setq d (cdr (aref ax r c)))) nil)
 145             (t (setq mul (ptimes mul (pquotient d (pgcd mul d)))))))))
 146
 147 ;; The author of the following programs is Tadatoshi Minamikawa (TM).
 148 ;; This program is one-step fraction-free Gaussian elimination with
 149 ;; optimal pivoting.  DRB claims the hair in this program is not
 150 ;; necessary and that straightforward Gaussian elimination is sufficient,
 151 ;; for sake of future implementers.
 152
 153 ;; To debug, delete the comments around PRINT and BREAK statements.
 154
 155 (declare-top (special permsign a rank delta nrow nvar n m variableorder
 156                       dependentrows inconsistentrows l k))
 157
 158 (defun tfgeli1 (ax n m)
 159   (prog (k l delta variableorder inconsistentrows
 160          dependentrows nrow nvar rank permsign result)
 161      (setq ax (get-array-pointer ax))
 162      (setq *col* (make-array (1+ m) :initial-element 0))
 163      (setq *row* (make-array (1+ n) :initial-element 0))
 164      (setq *colinv* (make-array (1+ m) :initial-element 0))
 165      ;; (PRINT 'ONESTEP-LIPSON-WITH-PIVOTING)
 166      (setq nrow n)
 167      (setq nvar (cond (*rank* m) (*det* m) (*inv* n) (*ech* m) (*tri* m) (t (1- m))))
 168      (do ((i 1 (1+ i)))
 169          ((> i n))
 170        (setf (aref *row* i) i))
 171      (do ((i 1 (1+ i)))
 172          ((> i m))
 173        (setf (aref *col* i) i) (setf (aref *colinv* i) i))
 174      (setq result
 175            (cond (*rank* (forward t) rank)
 176                  (*det* (forward t)
 177                         (cond ((= nrow n) (cond (permsign  (pminus delta))
 178                                                 (t delta)))
 179                               (t 0)))
 180                  (*inv* (forward t) (backward) (recoverorder1))
 181                  (*ech* (forward nil) (recoverorder2))
 182                  (*tri* (forward nil) (recoverorder2))
 183                  (t (forward t) (cond ($backsubst (backward)))
 184                     (recoverorder2)
 185                     (list dependentrows inconsistentrows variableorder))))
 186      (return result)))
 187
 188 ;;FORWARD ELIMINATION
 189 ;;IF THE SWITCH *CPIVOT IS NIL, IT AVOIDS THE COLUMN PIVOTING.
 190 (defun forward (*cpivot)
 191   (setq delta 1)                  ;DELTA HOLDS THE CURRENT DETERMINANT
 192   (do ((k 1 (1+ k))
 193        (nvar nvar)   ;PROTECTS AGAINST TEMPORARAY RESETS DONE IN PIVOT
 194        (m m))
 195       ((or (> k nrow) (> k nvar)))
 196     (cond ((pivot ax k *cpivot) (return nil)))
 197     ;; PIVOT IS T IF THERE IS NO MORE NON-ZERO ROW LEFT. THEN GET OUT OF THE LOOP
 198     (do ((i (1+ k) (1+ i)))
 199         ((> i nrow))
 200       (do ((j (1+ k) (1+ j)))
 201           ((> j m))
 202         (setf (aref ax (aref *row* i) (aref *col* j))
 203                (pquotient (pdifference (ptimes (aref ax (aref *row* k) (aref *col* k))
 204                                                (aref ax (aref *row* i) (aref *col* j)))
 205                                        (ptimes (aref ax (aref *row* i) (aref *col* k))
 206                                                (aref ax (aref *row* k) (aref *col* j))))
 207                           delta))))
 208     (do ((i (1+ k) (1+ i)))
 209         ((> i nrow))
 210       (setf (aref ax (aref *row* i) (aref *col* k)) 0))
 211     (setq delta (aref ax (aref *row* k) (aref *col* k))))
 212   ;; UNDOES COLUMN HACK IN PIVOT.
 213   (or *cpivot (do ((i 1 (1+ i))) ((> i m)) (setf (aref *col* i) i)))
 214   (setq rank (min nrow nvar)))
 215
 216 ;; BACKWARD SUBSTITUTION
 217 (defun backward ()
 218   (declare(special ax delta m rank))
 219   (do ((i (1- rank) (1- i)))
 220       ((< i 1))
 221     (do ((l (1+ rank) (1+ l)))
 222     ((> l m))
 223       (setf (aref ax (aref *row* i) (aref *col* l))
 224     (let ((mess1  (pdifference
 225              (ptimes (aref ax (aref *row* i) (aref *col* l))
 226                  (aref ax (aref *row* rank) (aref *col* rank)))
 227              (do ((j (1+ i) (1+ j)) (sum 0))
 228                  ((> j rank) sum)
 229                (setq sum (pplus sum (ptimes (aref ax (aref *row* i) (aref *col* j))
 230                             (aref ax (aref *row* j) (aref *col* l))))))) )
 231           (mess2 (aref ax (aref *row* i) (aref *col* i))  ))
 232       (cond ((equal 0 mess1) 0)
 233         ((equal 0 mess2) 0)
 234         (t   ;;   (pquotient mess1 mess2) ; fixed by line below. RJF 1/12/2017
 235
 236          (car (ratreduce mess1 mess2))
 237          )
 238         ))))
 239     (do ((l (1+ i) (1+ l)))
 240     ((> l rank))
 241       (setf (aref ax (aref *row* i) (aref *col* l)) 0)))
 242   ;; PUT DELTA INTO THE DIAGONAL MATRIX
 243   (setq delta (aref ax (aref *row* rank) (aref *col* rank)))
 244   (do ((i 1 (1+ i)))
 245       ((> i rank))
 246     (setf (aref ax (aref *row* i) (aref *col* i)) delta))  )
 247
 248 ;;RECOVER THE ORDER OF ROWS AND COLUMNS.
 249
 250 (defun recoverorder1 ()
 251   ;;(PRINT 'REARRANGE)
 252   (do ((i nvar (1- i)))
 253       ((= i 0))
 254     (setq variableorder (cons i variableorder)))
 255   (do ((i (1+ rank) (1+ i)))
 256       ((> i n))
 257     (cond ((equal (aref ax (aref *row* i) (aref *col* m)) 0)
 258            (setq dependentrows (cons (aref *row* i) dependentrows)))
 259           (t (setq inconsistentrows (cons (aref *row* i) inconsistentrows)))))
 260   (do ((i 1 (1+ i)))
 261       ((> i n))
 262     (cond ((not (= (aref *row* (aref *colinv* i)) i))
 263            (prog ()
 264               (moverow ax n m i 0)
 265               (setq l i)
 266               loop
 267               (setq k (aref *row* (aref *colinv* l)))
 268               (setf (aref *row* (aref *colinv* l)) l)
 269               (cond ((= k i) (moverow ax n m 0 l))
 270                     (t (moverow ax n m k l)
 271                        (setq l k)
 272                        (go loop))))))))
 273
 274 (defun recoverorder2 ()
 275   (do ((i nvar (1- i)))
 276       ((= i 0))
 277     (setq variableorder (cons (aref *col* i) variableorder)))
 278   (do ((i (1+ rank) (1+ i)))
 279       ((> i n))
 280     (cond ((equal (aref ax (aref *row* i) (aref *col* m)) 0)
 281            (setq dependentrows (cons (aref *row* i) dependentrows)))
 282           (t (setq inconsistentrows (cons (aref *row* i) inconsistentrows)))))
 283   (do ((i 1 (1+ i)))
 284       ((> i n))
 285     (cond ((not (= (aref *row* i) i))
 286            (prog ()
 287               (moverow ax n m i 0)
 288               (setq l i)
 289               loop
 290               (setq k (aref *row* l))
 291               (setf (aref *row* l) l)
 292               (cond ((= k i) (moverow ax n m 0 l))
 293                     (t (moverow ax n m k l)
 294                        (setq l k)
 295                        (go loop)))))))
 296   (do ((i 1 (1+ i)))
 297       ((> i nvar))
 298     (cond ((not (= (aref *col* i) i))
 299            (prog ()
 300               (movecol ax n m i 0)
 301               (setq l i)
 302               loop2
 303               (setq k (aref *col* l))
 304               (setf (aref *col* l) l)
 305               (cond ((= k i) (movecol ax n m 0 l))
 306                     (t (movecol ax n m k l)
 307                        (setq l k)
 308                        (go loop2))))))))
 309
 310 ;;THIS PROGRAM IS USED IN REARRANGEMENT
 311 (defun moverow (ax n m i j)
 312   (do ((k 1 (1+ k))) ((> k m))
 313     (setf (aref ax j k) (aref ax i k))))
 314
 315 (defun movecol (ax n m i j)
 316   (do ((k 1 (1+ k))) ((> k n))
 317     (setf (aref ax k j) (aref ax k i))))
 318
 319 ;;COMPLEXITY IS DEFINED AS FOLLOWS
 320 ;; COMPLEXITY(0)=0
 321 ;; COMPLEXITY(CONSTANT)=1
 322 ;; COMPLEXITY(POLYNOMIAL)=1+SUM(COMPLEXITY(C(N))+COMPLEXITY(E(N)), FOR N=0,1 ...M)
 323 ;; WHERE POLYNOMIAL IS OF THE FORM
 324 ;;    SUM(C(N)*X^E(N), FOR N=0,1 ... M)     X IS THE VARIABLE
 325
 326 (defun complexity (exp)
 327   (cond ((null exp) 0)
 328         ((equal exp 0) 0)
 329         ((atom  exp) 1)
 330         (t (+ (complexity (car exp)) (complexity (cdr exp))))))
 331
 332 (defun complexity/row (ax i j1 j2)
 333   (do ((j j1 (1+ j)) (sum 0))
 334       ((> j j2) sum)
 335     (incf sum (complexity (aref ax (aref *row* i) (aref *col* j))))))
 336
 337 (defun complexity/col (ax j i1 i2)
 338   (do ((i i1 (1+ i)) (sum 0))
 339       ((> i i2) sum)
 340     (incf sum (complexity (aref ax (aref *row* i) (aref *col* j))))))
 341
 342 (defun zerop/row (ax i j1 j2)
 343   (do ((j j1 (1+ j)))
 344       ((> j j2) t)
 345     (cond ((not (equal (aref ax (aref *row* i) (aref *col* j)) 0)) (return nil)))))
 346
 347 ;;PIVOTING ALGORITHM
 348 (defun pivot (ax k *cpivot)
 349   (prog (row/optimal col/optimal complexity/i/min complexity/j/min
 350          complexity/i complexity/j complexity/det complexity/det/min dummy)
 351      (setq row/optimal k complexity/i/min 1000000. complexity/j/min 1000000.)
 352      ;;TEST THE SINGULARITY
 353      (cond ((do ((i k (1+ i)) (isallzero t))
 354                 ((> i nrow) isallzero)
 355              loop (cond ((zerop/row ax i k nvar)
 356                          (cond (*inv* (merror (intl:gettext "solve: singular matrix.")))
 357                                (t (exchangerow i nrow)
 358                                   (decf nrow)))
 359                          (unless (> i nrow) (go loop)))
 360                         (t (setq isallzero nil))))
 361             (return t)))
 362
 363      ;; FIND AN OPTIMAL ROW
 364      ;; IF *CPIVOT IS NIL, (AX I K) WHICH IS TO BE THE PIVOT MUST BE NONZERO.
 365      ;; BUT IF *CPIVOT IS T, IT IS UNNECESSARY BECAUSE WE CAN DO THE COLUMN PIVOT.
 366      findrow
 367      (do ((i k (1+ i)))
 368          ((> i nrow))
 369        (cond ((or *cpivot (not (equal (aref ax (aref *row* i) (aref *col* k)) 0)))
 370               (cond ((> complexity/i/min
 371                         (setq complexity/i (complexity/row ax i k m)))
 372                      (setq row/optimal i complexity/i/min complexity/i))))))
 373      ;; EXCHANGE THE ROWS K AND ROW/OPTIMAL
 374      (exchangerow k row/optimal)
 375
 376      ;; IF THE FLAG *CPIVOT IS NIL, THE FOLLOWING STEPS ARE NOT EXECUTED.
 377      ;; THIS TREATMENT WAS DONE FOR THE LSA AND ECHELONTHINGS WHICH ARE NOT
 378      ;; HAPPY WITH THE COLUMN OPERATIONS.
 379      (cond ((null *cpivot)
 380             (cond ((not (equal (aref ax (aref *row* k) (aref *col* k)) 0))
 381                    (return nil))
 382                   (t (do ((i k (1+ i))) ((= i nvar))
 383                        (setf (aref *col* i) (aref *col* (1+ i))))
 384                      (setq nvar (1- nvar) m (1- m))
 385                      (go findrow)))))
 386
 387      ;;STEP3 ... FIND THE OPTIMAL COLUMN
 388      (setq col/optimal 0
 389            complexity/det/min 1000000.
 390            complexity/j/min 1000000.)
 391
 392      (do ((j k (1+ j)))
 393          ((> j nvar))
 394        (cond ((not (equal (aref ax (aref *row* k) (aref *col* j)) 0))
 395               (cond ((> complexity/det/min
 396                         (setq complexity/det
 397                               (complexity (aref ax (aref *row* k) (aref *col* j)))))
 398                      (setq col/optimal j
 399                            complexity/det/min complexity/det
 400                            complexity/j/min (complexity/col ax j (1+ k) n)))
 401                     ((equal complexity/det/min complexity/det)
 402                      (cond ((> complexity/j/min
 403                                (setq complexity/j
 404                                      (complexity/col ax j (1+ k) n)))
 405                             (setq col/optimal j
 406                                   complexity/det/min complexity/det
 407                                   complexity/j/min complexity/j))))))))
 408
 409      ;; EXCHANGE THE COLUMNS K AND COL/OPTIMAL
 410      (exchangecol  k col/optimal)
 411      (setq dummy (aref *colinv* (aref *col* k)))
 412      (setf (aref *colinv* (aref *col* k)) (aref *colinv* (aref *col* col/optimal)))
 413      (setf (aref *colinv* (aref *col* col/optimal)) dummy)
 414      (return nil)))
 415
 416 (defun exchangerow (i j)
 417   (prog (dummy)
 418      (setq dummy (aref *row* i))
 419      (setf (aref *row* i) (aref *row* j))
 420      (setf (aref *row* j) dummy)
 421      (cond ((= i j) (return nil))
 422            (t (setq permsign (not permsign))))))
 423
 424 (defun exchangecol (i j)
 425   (prog (dummy)
 426      (setq dummy (aref *col* i))
 427      (setf (aref *col* i) (aref *col* j))
 428      (setf (aref *col* j) dummy)
 429      (cond ((= i j) (return nil))
 430            (t (setq permsign (not permsign))))))
 431
 432 ;; Displays list of solutions.
 433
 434 (defun solve2 (llist equations)
 435   (setq $multiplicities nil)
 436   (map2c #'(lambda (equatn multipl)
 437              (setq equations
 438                    (nconc equations (list (displine equatn))))
 439              (push multipl $multiplicities)
 440              (if (and (> multipl 1) $dispflag)
 441                  (mtell (intl:gettext "solve: multiplicity ~A~%") multipl)))
 442          llist)
 443   (values (setq $multiplicities (cons '(mlist simp) (nreverse $multiplicities)))
 444           equations))
 445
 446 ;; Displays an expression and returns its linelabel.
 447
 448 (defun displine (exp)
 449   (let ($nolabels (tim 0))
 450     (elabel exp)
 451     (cond ($dispflag (remprop *linelabel* 'nodisp)
 452                      (setq tim (get-internal-run-time))
 453                      (mterpri)
 454                      (displa (list '(mlabel) *linelabel* exp))
 455                      (timeorg tim))
 456           (t (putprop *linelabel* t 'nodisp)))
 457     *linelabel*))
 458
 459 (declare-top (unspecial permsign a rank delta nrow nvar n m variableorder
 460                         dependentrows inconsistentrows l k))