src/mat.lisp

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