src/linnew.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 1980 Massachusetts Institute of Technology         ;;;
   9 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  10
  11 (in-package :maxima)
  12 (macsyma-module linnew)
  13
  14 ;; This is a matrix package which uses minors, basically.
  15 ;; TMLINSOLVE(LIST-OF-EQUAIONS,LIST-OF-VARIABLES,LIST-OF-VARIABLES-TO-BE-OBTAINED)
  16 ;; solves the linear equation. LIST-OF-VARIABLES-TO-BE-OBTAINED can be omitted,
  17 ;; in which case all variables are obtained. TMNEWDET(MATRIX,DIMENSION)
  18 ;; computes the determinant.  DIMENSION can be omitted.  The default is
  19 ;; DIMENSION=(declared dimension of MATRIX). TMINVERSE(MATRIX) computes the
  20 ;; inverse of matrix.
  21
  22 ;; The program uses hash arrays to remember the minors if N > threshold.  If
  23 ;; $WISE is set to T, the program knocks out unnecessary elements.  But also it
  24 ;; kills necessary ones in the case of zero elements! The $WISE flag should
  25 ;; not be set to T for inverse.  The default of $WISE is NIL.
  26
  27 ;; There seem to have been a number of bugs in this code.  I changed
  28 ;; the array referencing to value cell, and some of the stuff about
  29 ;; cre form.  It now seems tminverse  and tmlinsolve, now seem to work. --wfs.
  30
  31 ;;these are arrays
  32 (declare-top (special *tmarrays* *a2* *b* *aa* *row* *col* *rowinv* *colinv* *indx*))
  33
  34 (declare-top (special n nx ix))
  35
  36 (declare-top (special $wise $fool))
  37
  38 (defvar *tmarrays* nil)
  39
  40 ;; If N < threshold declared array is used, otherwise hashed array.
  41
  42 (defparameter *threshold* 10)
  43
  44 (defun tminitialflag nil
  45   (unless (boundp '$wise) (setq $wise nil))
  46   (unless (boundp '$fool) (setq $fool nil)))
  47
  48 ;; TMDET returns the determinant of N*N matrix A2 which is in an globally
  49 ;; declared array A2.
  50
  51 (defun tmdet (a4 n)
  52   (prog (index ix)
  53      (tminitialflag)
  54      (setq ix 0 nx 0)
  55      (do ((i 1 (1+ i)))
  56          ((> i n))
  57        (push i index))
  58      (setq index (nreverse index))
  59      (tminor a4 n 1 index 0)))
  60
  61 ;; TMLIN SOLVES M SETS OF LINEAR EQUATIONS WITH N UNKNOWN VARIABLES. IT SOLVES
  62 ;; ONLY FOR THE FIRST NX UNKNOWNS OUT OF N. THE EQUATIONS ARE EXPRESSED IN
  63 ;; MATRIX FORM WHICH IS IN N*(N+M) ARRAY A2. AS USUAL , THE LEFT HAND SIDE N*N
  64 ;; OF A2 REPRESENTS THE COEFFICIENT MATRIX, AND NEXT N*M OF A2 IS THE RIGHT
  65 ;; HAND SIDE OF THE M SETS OF EQUATIONS.  SUPPOSE N=3, M=2, AND THE UNKKNOWNS
  66 ;; ARE (X1 Y1 Z1) FOR THE FIRST SET AND (X2 Y2 Z2) FOR THE SECOND. THEN THE
  67 ;; RESULT OF TMLIN IS ((DET) (U1 U2) (V1 V2) (W1 W2)) WHERE DET IS THE
  68 ;; DETERMINANT OF THE COEFFICIENT MATRIX AND X1=U1/DET, X2=U2/DET, Y1=V1/DET,
  69 ;; Y2=V2/DET ETC.
  70
  71 (defun tmlin (a4 n m nx)
  72   (prog (index r)
  73      (tmdefarray n)
  74      (tminitialflag)
  75      (do ((i 1 (1+ i)))
  76          ((> i n))
  77        (push i index))
  78      (setq index (reverse index))
  79      (setq r
  80            (do ((ix 0 (1+ ix))
  81                 (result))
  82                ((> ix nx) (reverse result))
  83              (push (do ((i 1 (1+ i)) (res))
  84                        ((> i (if (= ix 0) 1 m))
  85                         (reverse res))
  86                      (unless $wise (tmkillarray ix))
  87                      (push (tminor a4 n 1 index i) res))
  88                    result)
  89              (when (and (= ix 0) (equal (car result) '(0 . 1)))
  90                (merror (intl:gettext "tmlin: coefficient matrix is singular.")))))
  91      (tmrearray n)
  92      (return r)))
  93
  94 ;; TMINOR ACTUALLY COMPUTES THE MINOR DETERMINANT OF A SUBMATRIX OF A2, WHICH
  95 ;; IS CONSTRUCTED BY EXTRACTING ROWS (K,K+1,K+2,...,N) AND COLUMNS SPECIFIED BY
  96 ;; INDEX. N IS THE DIMENSION OF THE ORIGINAL MATRIX A2.  WHEN TMINOR IS USED
  97 ;; FOR LINEAR EQUATION PROGRAM, JRIGHT SPECIFIES A COLUMN OF THE CONSTANT
  98 ;; MATRIX WHICH IS PLUGGED INTO AN IX-TH COLUMN OF THE COEFFICIENT MATRIX FOR
  99 ;; ABTAINING IX-TH UNKNOWN. IN OTHER WORDS, JRIGHT SPECIFIES JRIGHT-TH
 100 ;; EQUATION.
 101
 102 (defun tminor (a4 n k index jright)
 103   (prog (subindx l result name aorb)
 104      (setq a4 (get-array-pointer a4))
 105      (cond ((= k n)
 106             (setq result
 107                   (if (= k ix)
 108                       (aref a4 (car index) (+ jright n))
 109                       (aref a4 (car index) k))))
 110            (t
 111             (do
 112              ((j 1 (1+ j))
 113               (sum '(0 . 1)))
 114              ((> j (1+ (- n k))) (setq result sum))
 115               (setq l (extract index j))
 116               (setq subindx (cadr l))
 117               (setq l (car l))
 118               (setq aorb (if (= k ix)
 119                              (aref a4 l (+ jright n))
 120                              (aref a4 l k)))
 121               (unless (equal aorb '(0 . 1))
 122                 (setq name (tmaccess subindx))
 123                 (setq sum
 124                       (funcall (if (oddp j) #'ratplus #'ratdifference)
 125                                sum
 126                                (rattimes
 127                                 aorb
 128                                 (if $fool
 129                                     (tminor a4 n (1+ k) subindx jright)
 130                                     (cond ((not (null (tmeval name)))
 131                                            (tmeval name))
 132                                           ((tmnomoreuse j l k)
 133                                            (tmstore name nil)
 134                                            (tminor a4 n (1+ k) subindx jright))
 135                                           (t
 136                                            (tmstore name (tminor a4 n (1+ k) subindx jright)))))
 137                                 t))))
 138               (when $wise
 139                 (when (tmnomoreuse j l k)
 140                   (tmkill subindx k))))))
 141      (return result)))
 142
 143 (defun extract (index j)
 144   (do ((ind index (cdr ind))
 145        (count 1 (1+ count))
 146        (subindx))
 147       ((null ind))
 148     (if (= count j)
 149         (return (list (car ind) (nconc subindx (cdr ind))))
 150         (setq subindx (nconc subindx (list (car ind)))))))
 151
 152 (declare-top (special vlist))
 153
 154 (defun tmratconv (bbb n m)
 155   (prog (ccc)
 156      (declare (special ccc))   ;Tell me this worked in Maclisp.  --gsb
 157                                         ;Actually, i suspect it didn't, at least ever since
 158                                         ; (sstatus punt).
 159      (setf (symbol-value 'ccc) bbb)
 160      (do ((k 1 (1+ k)))
 161          ((> k n))
 162        (do ((j 1 (1+ j)))
 163            ((> j m))
 164          (newvar1 (setf (aref *a2* k j)
 165                          (maref ccc k j)
 166                          ;; (nth j (nth k *a2*))
 167                          ;; (MEVAL (LIST (LIST 'CCC 'array) K J))  ;;just the
 168                          ))))
 169
 170      (newvar (cons '(mtimes) vlist))
 171      (do ((k 1 (1+ k)))
 172          ((> k n))
 173        (do ((j 1 (1+ j)))
 174            ((> j m))
 175          (setf (aref *a2* k j) (cdr (ratrep* (aref *a2* k j))))))))
 176
 177 (defmfun $tmnewdet (mat &optional (dim nil dim?))
 178   (prog (*aa* r vlist n)
 179      (cond (dim?
 180             (unless (integerp dim)
 181               (merror (intl:gettext "tmnewdet: second argument must be an integer; found: ~M") dim))
 182             (setq n dim))
 183            (($matrixp mat)
 184             (setq n (length (cdr mat))))
 185            (t
 186             (merror (intl:gettext "tmnewdet: first argument must be a matrix; found: ~M") mat)))
 187      (setq *aa* mat)
 188      (setq *a2* (make-array (list (1+ n) (1+ n)) :initial-element nil))
 189      (tmdefarray n)
 190      (tmratconv *aa* n n)
 191      (setq r (cons (list 'mrat 'simp varlist genvar) (tmdet '*a2* n)))
 192      (tmrearray n)
 193      (return r)))
 194
 195 (defmfun $tmlinsolve (&rest arglist)
 196   (prog (equations vars outvars result *aa*)
 197      (setq equations (cdar arglist)
 198            vars (cdadr arglist)
 199            outvars (cond ((null (cddr arglist)) vars)
 200                          (t (cdaddr arglist))))
 201      (setq vars (tmerge vars outvars))
 202      (setq nx (length outvars))
 203      (setq n (length vars))
 204      (unless (= n (length equations))
 205        (return (print 'too-few-or-much-equations)))
 206      (setq *aa*
 207            (cons '($matrix simp)
 208                  (mapcar #'(lambda (exp)
 209                              (append
 210                               '((mlist))
 211                               (mapcar #'(lambda (v)
 212                                           (prog (r)
 213                                              (setq exp ($bothcoef exp v)
 214                                                    r (cadr exp)
 215                                                    exp (meval (caddr exp)))
 216                                              (return r)))
 217                                       vars)
 218                               (list (list '(mminus) exp))))
 219                          (mapcar #'(lambda (e)
 220                                      (meval (list '(mplus)
 221                                                   ($lhs e)
 222                                                   (list '(mminus) ($rhs e)))))
 223                                  equations))))
 224      (setq result (cdr ($tmlin *aa* n 1 nx)))
 225      (return
 226        (do ((vars (cons nil outvars) (cdr vars))
 227             (labels)
 228             (dlabel)
 229             (name))
 230            ((null vars)
 231             (cons '(mlist) (cdr (reverse labels))))
 232          (setq name (makelabel $linechar))
 233          (incf $linenum)
 234          (setf (symbol-value name)
 235                (cond ((null (car vars))
 236                       (setq dlabel name)
 237                       (cadar result))
 238                      (t (list '(mequal)
 239                               (car vars)
 240                               (list '(mtimes simp)
 241                                     (cadar result)
 242                                     (list '(mexpt simp) dlabel -1))))))
 243          (push name labels)
 244          (setq result (cdr result))
 245          (when $dispflag
 246            (mtell-open "~M" (nconc (ncons '(mlabel))
 247                                    (ncons name)
 248                                    (ncons (eval name)))))))))
 249
 250 (defun tmerge (vars outvars)
 251   (append outvars
 252           (prog (l)
 253              (mapcar #'(lambda (v)
 254                          (if (member v outvars) nil (push v l)))
 255                      vars)
 256              (return (reverse l)))))
 257
 258 (defmfun $tmlin (*aa* n m nx)
 259   (prog (r vlist)
 260      (setq *a2* (make-array (list (1+ n) (+ 1 m n)) :initial-element nil))
 261      (show *a2*)
 262      (tmratconv *aa* n (+ m n))
 263      (setq r
 264            (cons '(mlist)
 265                  (mapcar
 266                   #'(lambda (res)
 267                       (cons '(mlist)
 268                             (mapcar #'(lambda (result)
 269                                         (cons (list 'mrat 'simp varlist genvar) result))
 270                                     res)))
 271                   (tmlin '*a2* n m nx))))
 272      (show *a2*)
 273      (return r)))
 274
 275 (defun tmkill (*indx* k)
 276   (prog (name subindx j l)
 277      (when (null *indx*) (return nil))
 278      (setq name (tmaccess *indx*))
 279      (cond ((not (null (tmeval name)))
 280             (tmstore name nil))
 281            (t
 282             (do ((ind *indx* (cdr ind))
 283                  (count 1 (1+ count)))
 284                 ((null ind))
 285               (setq l (extract *indx* count)
 286                     j (car l)
 287                     subindx (cadr l))
 288               (when (= j count)
 289                 (tmkill subindx (1+ k))))))))
 290
 291 (defun tmnomoreuse (j l k)
 292   (if (and (= j l) (or (> k nx) (< k (1+ ix))))
 293       t
 294       nil))
 295
 296 (defun tmdefarray (n)
 297   (prog (name)
 298      (cond ((setq *tmarrays* (get-array-pointer *tmarrays*))
 299             (tmrearray (1- (cond ((cadr (arraydims *tmarrays*)))
 300                                  (t 1))))))
 301      (setq *tmarrays* (make-array (1+ n) :initial-element nil))
 302      (do ((i 1 (1+ i)))
 303          ((> i n))
 304        (setq name (if (= i 1) (make-symbol "M") (gensym)))
 305        (cond ((< n *threshold*)
 306               (setf (symbol-value name) (make-array (1+ (tmcombi n i)) :initial-element nil))
 307               (setf (aref *tmarrays* i) (get-array-pointer name)))
 308              (t
 309               (setf (aref *tmarrays* i) (list name 'simp 'array)))))
 310      (gensym "G")))
 311
 312 ;; TMREARRAY kills the TMARRAYS which holds pointers to minors. If (TMARRAYS I)
 313 ;; is an atom, it is declared array.  Otherwise it is hashed array.
 314
 315 (defun tmrearray (n)
 316   (prog nil
 317      (do ((i 1 (1+ i)))
 318          ((> i n))
 319        (unless (atom (aref *tmarrays* i))
 320          (tm$kill (car (aref *tmarrays* i)))))))
 321
 322 (defun tmaccess (index)
 323   (prog (l)
 324      (cond ($fool (return nil)))
 325      (setq l (length index))
 326      (return
 327        (cond ((< n *threshold*)
 328               (list 'aref (aref *tmarrays* l)
 329                     (do ((i 1 (1+ i))
 330                          (x 0 (car y))
 331                          (y index (cdr y))
 332                          (sum 0))
 333                         ((> i l) (1+ sum))
 334                       (do ((j (1+ x) (1+ j)))
 335                           ((= j (car y)))
 336                         (incf sum (tmcombi (- n j) (- l i)))))))
 337              (t (cons 'aref (cons (aref *tmarrays* l) index)))))) )
 338
 339 (defun tmcombi (n i)
 340   (if (> (- n i) i)
 341       (/ (tmfactorial n (- n i)) (tmfactorial i 0))
 342       (/ (tmfactorial n i) (tmfactorial (- n i) 0))))
 343
 344 (defun tmfactorial (i j)
 345   (if (= i j)
 346       1
 347       (* i (tmfactorial (1- i) j))))
 348
 349 (defun tmstore (name x)
 350   (cond ((< n *threshold*)
 351          (eval `(setf ,name ',x)))
 352         (t
 353          (mset name (list '(mquote simp) x))
 354          x)))
 355
 356 ;; TMKILLARRAY kills all (N-IX+1)*(N-IX+1) minors which are not necessary for
 357 ;; the computation of IX-TH variable in the linear equation.  Otherwise, they
 358 ;; will do harm.
 359
 360 (defun tmkillarray (ix)
 361   (do ((i (1+ (- n ix)) (1+ i)))
 362       ((> i n))
 363     (if (< n *threshold*)
 364         (fillarray (aref *tmarrays* i) '(nil))
 365         (tm$kill (car (aref *tmarrays* i))))))
 366
 367 (defun tmeval (e)
 368   (prog (result)
 369      (return (cond ((< n *threshold*)
 370                     (eval e))
 371                    (t
 372                     (setq result (meval e))
 373                     (if (equal result e) nil (cadr result)))))))
 374
 375 (defun tm$kill (e)
 376   (kill1 e))
 377
 378 (defmfun $tminverse (*aa*)
 379   (prog (r vlist n m nx)
 380      (setq n (length (cdr *aa*)) m n nx n)
 381      (setq *a2* (make-array (list (1+ n) (+ 1 m n)) :initial-element nil))
 382      (tmratconv *aa* n n)
 383      (do ((i 1 (1+ i)))
 384          ((> i n))
 385        (do ((j 1 (1+ j)))
 386            ((> j m))
 387          (setf (aref *a2* i (+ n j))
 388                (if (= i j) '(1 . 1) '(0 . 1)))))
 389      (setq r (mapcar #'(lambda (res)
 390                          (cons '(mlist)
 391                                (mapcar #'(lambda (result)
 392                                            ($ratdisrep (cons (list 'mrat 'simp varlist genvar) result)))
 393                                        res)))
 394                      (tmlin '*a2* n m nx)))
 395      (setq r (list '(mtimes simp)
 396                    (list '(mexpt simp) (cadar r) -1)
 397                    (cons '($matrix simp) (cdr r))))
 398      (return r)))
 399
 400 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 401 ;;THIS IS A UTILITY PACKAGE FOR SPARSE
 402 ;;MATRIX INVERSION. A3 IS A N*N MATRIX.
 403 ;;IT RETURNS A LIST OF LISTS, SUCH AS
 404 ;;((I1 I2 ...) (J1 J2...) ...) WHERE (I1
 405 ;;I2 ..) SHOWS THE ROWS WHICH BELONGS TO
 406 ;;THE FIRST BLOCK, AND SO ON.  THE ROWS
 407 ;;SHOULD BE REORDERED IN THIS ORDER. THE
 408 ;;COLUMNS ARE NOT CHANGED. IT RETURNS NIL
 409 ;;IF A3 IS "OBVIOUSLY" SINGULAR.
 410
 411 ;; (DEFUN TMISOLATE (A3 N)
 412 ;;        (PROG (NODELIST)
 413 ;;           (SETQ A3 (GET A3 'ARRAY))
 414 ;;           (setq  B (*ARRAY nil 'T (1+ N) (1+ N)))
 415 ;;           (setq  ROW (*ARRAY nil 'T (1+ N)))
 416 ;;           (setq  COL (*ARRAY nil 'T (1+ N)))
 417 ;;           (DO ((I 1 (1+ I)))
 418 ;;               ((> I N))
 419 ;;               (STORE (ROW I) I)
 420 ;;               (STORE (COL I) I))
 421 ;;           (DO ((I 1 (1+ I)))
 422 ;;               ((> I N))
 423 ;;               (DO ((J 1 (1+ J)))
 424 ;;                   ((> J N))
 425 ;;                   (STORE (B I J)
 426 ;;                          (NOT (EQUAL (AREF A3 I J)
 427 ;;                                      '(0 . 1))))))
 428 ;;           (COND ((NULL (TMPIVOT-ISOLATE 1))
 429 ;;                  (SETQ NODELIST NIL)
 430 ;;                  (GO EXIT)))
 431 ;;           (DO ((I 1 (1+ I)))
 432 ;;               ((> I N))
 433 ;;               (DO ((J 1 (1+ J)))
 434 ;;                   ((> J I))
 435 ;;                   (STORE (B (ROW J) (COL I))
 436 ;;                          (OR (B (ROW I) (COL J))
 437 ;;                              (B (ROW J) (COL I))))
 438 ;;                   (STORE (B (ROW I) (COL J)) (B (ROW J) (COL I))))
 439 ;;               (STORE (B (ROW I) (COL I)) T))
 440 ;;           (DO ((I 1 (1+ I)))
 441 ;;               ((> I N))
 442 ;;               (COND ((EQ (B (ROW I) (COL I)) T)
 443 ;;                      (SETQ NODELIST
 444 ;;                            (CONS (TMPULL-OVER I N) NODELIST)))))
 445 ;;           EXIT
 446 ;;           (*TMREARRAY 'B)
 447 ;;           (*TMREARRAY 'ROW)
 448 ;;           (*TMREARRAY 'COL)
 449 ;;           (RETURN (REVERSE NODELIST)))))
 450
 451 ;; (DEFUN TMPULL-OVER (P N)
 452 ;;        (PROG (Q)
 453 ;;           (STORE (B (ROW P) (COL P)) NIL)
 454 ;;           (DO ((J 1 (1+ J)))
 455 ;;               ((> J N) (SETQ Q NIL))
 456 ;;               (COND ((EQ (B (ROW P) (COL J)) T)
 457 ;;                      (RETURN (SETQ Q J)))))
 458 ;;           (COND ((NULL Q) (RETURN (LIST (ROW P))))
 459 ;;                 (T (DO ((J 1 (1+ J)))
 460 ;;                        ((> J N))
 461 ;;                        (STORE (B (ROW Q) (COL J))
 462 ;;                               (OR (B (ROW Q) (COL J))
 463 ;;                                   (B (ROW P) (COL J))))
 464 ;;                        (STORE (B (ROW J) (COL Q))
 465 ;;                               (B (ROW Q) (COL J))))
 466 ;;                    (TMCRIP P)
 467 ;;                    (RETURN (CONS (ROW P) (TMPULL-OVER Q N)))))))
 468
 469 ;; (DEFUN TMCRIP (P)
 470 ;;        (DO ((I 1 (1+ I)))
 471 ;;         ((> I N))
 472 ;;         (STORE (B (ROW P) (COL I)) NIL)
 473 ;;         (STORE (B (ROW I) (COL P)) NIL)))
 474
 475 ;;TMPIVOT-ISOLATE CARRIES OUT PIVOTING
 476 ;;SO THAT THE ALL DIAGONAL ELEMENTS ARE
 477 ;;NONZERO. THIS GUARANTEES WE HAVE MAXIMUM
 478 ;;NUMBER OF BLOCKS ISOLATED.
 479
 480 (defun tmpivot-isolate (k)
 481   (cond ((> k n) t)
 482         (t (do ((i k (1+ i)))
 483                ((> i n) nil)
 484              (when (aref *b* (aref *row* i) (aref *col* k))
 485                (tmexchange '*row* k i)
 486                (if (tmpivot-isolate (1+ k))
 487                    (return t)
 488                    (tmexchange '*row* k i)))))))
 489
 490 (defun tmexchange (rowcol i j)
 491   (prog (dummy)
 492      (setq rowcol (get-array-pointer rowcol))
 493      (setq dummy (aref rowcol i))
 494      (setf (aref rowcol i) (aref rowcol j))
 495      (setf (aref rowcol j) dummy)))
 496
 497
 498 ;; PROGRAM TO PREDICT ZERO ELEMENTS IN
 499 ;; THE SOLUTION OF INVERSE OR LINEAR
 500 ;; EQUATION. A IS THE COEFFICIENT MATRIX.
 501 ;; B IS THE RIGHT HAND SIDE MATRIX FOR
 502 ;; LINEAR EQUATIONS. A3 IS N*N AND B IS
 503 ;; M*M. X IS AN N*M MATRIX WHERE T -NIL
 504 ;; PATTERN SHOWING THE ZERO ELEMENTS IN
 505 ;; THE RESULT IS RETURNED. T CORRESPONDS TO
 506 ;; NON-ZERO ELEMENT. IN THE CASE OF
 507 ;; INVERSE, YOU CAN PUT ANYTHING (SAY,NIL)
 508 ;; FOR B AND 0 FOR M.  NORMALLY IT RETURNS
 509 ;; T, BUT IN CASE OF SINGULAR MATRIX, IT
 510 ;; RETURNS NIL.
 511
 512 ;; (DEFUN TMPREDICT (A3 B X N M)
 513 ;;   (PROG (FLAGINV FLAG-NONSINGULAR)
 514 ;;      (SETQ A3 (GET A3 'ARRAY) B (GET B 'ARRAY) X (GET X 'ARRAY))
 515 ;;      (setq  AA (*ARRAY nil 'T (1+ N) (1+ N)))
 516 ;;      (setq  ROW (*ARRAY nil 'T (1+ N)))
 517 ;;      (SETQ FLAGINV (= M 0))
 518 ;;      (COND (FLAGINV (SETQ M N)))
 519 ;;      (DO ((I 1 (1+ I)))
 520 ;;          ((> I N))
 521 ;;          (DO ((J 1 (1+ J)))
 522 ;;              ((> J N))
 523 ;;              (STORE (AA I J)
 524 ;;                     (NOT (EQUAL (AREF A3 I J) '(0 . 1))))))
 525 ;;      (DO ((I 1 (1+ I)))
 526 ;;          ((> I N))
 527 ;;          (DO ((J 1 (1+ J)))
 528 ;;              ((> J M))
 529 ;;              (STORE (AREF X I J)
 530 ;;                     (COND (FLAGINV (EQ I J))
 531 ;;                           (T (EQUAL (AREF B I J)
 532 ;;                                     '(0 . 1)))))))
 533 ;;      (DO ((I 1 (1+ I))) ((> I N)) (STORE (ROW I) I))
 534 ;;              ;FORWARD ELIMINATION.
 535 ;;      (DO ((I 1 (1+ I)))
 536 ;;          ((> I N))
 537 ;;          (SETQ FLAG-NONSINGULAR
 538 ;;                (DO ((II I (1+ II)))
 539 ;;                    ((> II N) NIL)
 540 ;;                    (COND ((AA (ROW II) I)
 541 ;;                           (TMEXCHANGE 'ROW II I)
 542 ;;                           (RETURN T)))))
 543 ;;          (COND ((NULL FLAG-NONSINGULAR) (RETURN NIL)))
 544 ;;          (DO ((II (1+ I) (1+ II)))
 545 ;;              ((> II N))
 546 ;;              (COND ((AA (ROW II) I)
 547 ;;                     (DO ((JJ (1+ I) (1+ JJ)))
 548 ;;                         ((> JJ N))
 549 ;;                         (STORE (AA (ROW II) JJ)
 550 ;;                                (OR (AA (ROW I) JJ)
 551 ;;                                    (AA (ROW II) JJ))))
 552 ;;                     (DO ((JJ 1 (1+ JJ)))
 553 ;;                         ((> JJ M))
 554 ;;                         (STORE (AREF X (ROW II) JJ)
 555 ;;                                (OR (AREF X (ROW I) JJ)
 556 ;;                                    (AREF X (ROW II) JJ))))))))
 557 ;;      (COND ((NULL FLAG-NONSINGULAR) (GO EXIT)))       ;GET OUT  BACKWARD SUBSTITUTION
 558 ;;      (DO ((I (1- N) (1- I)))
 559 ;;          ((< I 1))
 560 ;;          (DO ((L 1 (1+ L)))
 561 ;;              ((> L M))
 562 ;;              (STORE (AREF X (ROW I) L)
 563 ;;                     (OR (AREF X (ROW I) L)
 564 ;;                         (DO ((J (1+ I) (1+ J)) (SUM))
 565 ;;                             ((> J N) SUM)
 566 ;;                             (SETQ SUM
 567 ;;                                   (OR SUM
 568 ;;                                       (AND (AA (ROW I) J)
 569 ;;                                            (AREF
 570 ;;                                                       X
 571 ;;                                                       (ROW J)
 572 ;;                                                       L)))))))))
 573 ;;             ;RECOVER THE ORDER.
 574 ;;      (TMPERMUTE 'X N M 0 0 'ROW N 'ROW)
 575 ;;    EXIT (*TMREARRAY 'ROW) (*TMREARRAY 'AA) (RETURN FLAG-NONSINGULAR)))
 576
 577 ;;TMPERMUTE PERMUTES THE ROWS OR COLUMNS
 578 ;;OF THE N*M MATRIX AX ACCORDING TO THE
 579 ;;SPECIFICATION OF INDEXLIST. THE FLAG
 580 ;;MUST BE SET 'ROW IF ROW PERMUTATION IS
 581 ;;DESIRED , OR 'COL OTHERWISE. THE RESULT
 582 ;;IS IN AX. NM IS THE DIMENSION OF
 583 ;;INDEXLIST.
 584
 585 (defun tmpermute (ax n m rbias cbias indexlist nm flag)
 586   (prog (k l)
 587      ;;      (SETQ AX (GET AX 'array)
 588      ;;            INDEXLIST (GET INDEXLIST 'array))
 589      (setq ax (get-array-pointer ax))
 590      (setq indexlist (get-array-pointer indexlist))
 591      (setf (symbol-array *indx*) (make-array (1+ nm) :initial-element nil))
 592      (do ((i 1 (1+ i)))
 593          ((> i nm))
 594        (setf (aref *indx* i) (aref indexlist i)))
 595      (do ((i 1 (1+ i)))
 596          ((> i nm))
 597        (cond ((not (= (aref *indx* i) i))
 598               (prog nil
 599                  (tmmove ax n m rbias cbias i 0 flag)
 600                  (setq l i)
 601                  loop (setq k (aref *indx* l))
 602                  (setf (aref *indx* l) l)
 603                  (cond ((= k i)
 604                         (tmmove ax n m rbias cbias 0 l flag))
 605                        (t (tmmove ax n m rbias cbias k l flag)
 606                           (setq l k)
 607                           (go loop)))))))))
 608
 609 (defun tmmove (ax n m rbias cbias i j flag)
 610   (prog (ll)
 611      (setq ax (get-array-pointer ax))
 612      (setq ll (if (eq flag '*row*)
 613                   (- m cbias)
 614                   (- n rbias)))
 615      (do ((k 1 (1+ k)))
 616          ((> k ll))
 617        (cond ((eq flag '*row*)
 618               (setf (aref ax (+ rbias j) (+ cbias k))
 619                     (aref ax (+ rbias i) (+ cbias k))))
 620              (t (setf (aref ax (+ rbias k) (+ cbias j))
 621                       (aref ax  (+ rbias k) (+ cbias i))))))))
 622
 623 ;;TMSYMETRICP CHECKS THE SYMMETRY OF THE MATRIX.
 624
 625 (defun tmsymetricp (a3 n)
 626   (setq a3 (get-array-pointer a3))
 627   (do ((i 1 (1+ i)))
 628       ((> i n) t)
 629     (cond ((null (do ((j (1+ i) (1+ j)))
 630                      ((> j n) t)
 631                    (unless (equal (aref a3 i j) (aref a3 j i))
 632                      (return nil))))
 633            (return nil)))))
 634
 635 ;;TMLATTICE CHECKS THE "LATTICE"
 636 ;;STRUCTURE OF THE MATRIX A. IT RETURNS
 637 ;;NIL IF THE MATRIX IS "OBVIOUSLY"
 638 ;;SINGULAR. OTHERWISE IT RETURNS A LIST
 639 ;;(L1 L2 ... LM) WHERE M IS THE NUMBER OF
 640 ;;BLOCKS (STRONGLY CONNECTED SUBGRAPHS),
 641 ;;AND L1 L2 ... ARE LIST OF ROW AND
 642 ;;COLUMN NUMBERS WHICH BELONG TO EACH
 643 ;;BLOCKS. THE LIST LOOKS LIKE ((R1 C1)
 644 ;;(R2 C2) ...) WHERE R R'S ARE ROWS AND
 645 ;;C'S ARE COLUMNS.
 646
 647 (defun tmlattice (a3 xrow xcol n)
 648   (setq a3 (get-array-pointer a3))
 649   (setq xrow (get-array-pointer xrow))
 650   (setq xcol (get-array-pointer xcol))
 651   (setq *b* (make-array (list (1+ n) (1+ n)) :initial-element nil))
 652   (setq *row* (make-array (1+ n) :initial-element nil))
 653   (setq *col* (make-array (1+ n) :initial-element nil))
 654   (do ((i 1 (1+ i)))
 655       ((> i n))
 656     (do ((j 1 (1+ j)))
 657         ((> j n))
 658       (setf (aref *b* i j)
 659             (not (equal (aref a3 i j) '(0 . 1))))))
 660   (do ((i 0 (1+ i)))
 661       ((> i n))
 662     (setf (aref *row* i) i)
 663     (setf (aref *col* i) i))
 664   (when (null (tmpivot-isolate 1))
 665     (return-from tmlattice nil))
 666   (do ((i 1 (1+ i)))
 667       ((> i n))
 668     (setf (aref *b* (aref *row* i) (aref *col* i)) i)
 669     (setf (aref *b* (aref *row* i) (aref *col* 0)) t))
 670   (tmlattice1 1)
 671   (tmsort-lattice xrow xcol))
 672
 673 (defun tmlattice1 (k)
 674   (cond ((= k n)
 675          nil)
 676         (t
 677          (tmlattice1 (1+ k))
 678          (do ((looppath))
 679              (nil)
 680            (if (setq looppath (tmpathp k k))
 681                (tmunify-loop k (cdr looppath))
 682                (return nil))))))
 683
 684 (defun tmpathp (j k)
 685   (cond ((equal (aref *b* (aref *row* j) (aref *col* k)) t)
 686          (list j k))
 687         (t
 688          (do ((jj k (1+ jj)) (path))
 689              ((> jj n))
 690            (when (and (equal (aref *b* (aref *row* j) (aref *col* jj)) t)
 691                       (setq path (tmpathp jj k)))
 692              (return (cons j path)))))))
 693
 694 (defun tmunify-loop (k chain)
 695   (prog (l dummyk dummyl)
 696      (setq l (car chain))
 697      (cond ((= l k) (return nil)))
 698      (setq dummyk (aref *b* (aref *row* k) (aref *col* k)))
 699      (setq dummyl (aref *b* (aref *row* l) (aref *col* l)))
 700      (setf (aref *b* (aref *row* k) (aref *col* k)) nil)
 701      (setf (aref *b* (aref *row* l) (aref *col* l)) nil)
 702      (do ((i 1 (1+ i)))
 703          ((> i n))
 704        (setf (aref *b* (aref *row* k) (aref *col* i))
 705              (or (aref *b* (aref *row* k) (aref *col* i)) (aref *b* (aref *row* l) (aref *col* i))))
 706        (setf (aref *b* (aref *row* i) (aref *col* k))
 707              (or (aref *b* (aref *row* i) (aref *col* k)) (aref *b* (aref *row* i) (aref *col* l))))
 708        (setf (aref *b* (aref *row* l) (aref *col* i)) nil)
 709        (setf (aref *b* (aref *row* i) (aref *col* l)) nil))
 710      (setf (aref *b* (aref *row* k) (aref *col* k)) dummyl)
 711      (setf (aref *b* (aref *row* l) (aref *col* l)) dummyk)
 712      (setf (aref *b* (aref *row* k) (aref *col* 0)) t)
 713      (setf (aref *b* (aref *row* l) (aref *col* 0)) nil)
 714      (tmunify-loop k (cdr chain))))
 715
 716 (defun tmsort-lattice (xrow xcol)
 717   (prog (nodelist result)
 718      (setq nodelist (tmsort1))
 719      (setq result
 720            (do ((x nodelist (cdr x)) (result))
 721                ((null x) result)
 722              (setq result
 723                    (cons (do ((next (aref *b* (aref *row* (car x)) (aref *col* (car x)))
 724                                     (aref *b* (aref *row* next) (aref *col* next)))
 725                               (res))
 726                              ((= next (car x))
 727                               (cons (list (aref *row* next) (aref *col* next)) res))
 728                            (push (list (aref *row* next) (aref *col* next)) res))
 729                          result))))
 730      (do ((list1 result (cdr list1))
 731           (i 1))
 732          ((null list1))
 733        (do ((list2 (car list1) (cdr list2)))
 734            ((null list2))
 735          (setf (aref xrow i) (caar list2))
 736          (setf (aref xcol i) (cadar list2))
 737          (incf i)))
 738      (return result)))
 739
 740 ;; (DEFUN TMLESS (I J) (B (ROW I) (COL J)))
 741
 742 (defun tmsort1 nil
 743   (do ((i 1 (1+ i))
 744        (result))
 745       ((> i n) result)
 746     (cond ((and (aref *b* (aref *row* i) (aref *col* 0)) (tmmaxp i))
 747            (do ((j 1 (1+ j)))
 748                ((> j n))
 749              (unless (= j i)
 750                (setf (aref *b* (aref *row* i) (aref *col* j)) nil)))
 751            (setf (aref *b* (aref *row* i) (aref *col* 0)) nil)
 752            (push i result)
 753            (setq i 0)))))
 754
 755 (defun tmmaxp (i)
 756   (do ((j 1 (1+ j)))
 757       ((> j n) t)
 758     (when (and (not (= i j)) (aref *b* (aref *row* j) (aref *col* i)))
 759       (return nil))))
 760
 761 ;;UNPIVOT IS USED IN PAUL WANG'S PROGRAM
 762 ;;TO RECOVER THE PIVOTING. TO GET THE
 763 ;;INVERSE OF A, PAUL'S PROGRAM COMPUTES
 764 ;;THE INVERSE OF U*A*V BECAUSE OF
 765 ;;BLOCKING. LET THE INVERSE Y. THEN
 766 ;;A^^-1=V*Y*U. WHERE U AND V ARE
 767 ;;FUNDAMENTAL TRANSFORMATION
 768 ;;(PERMUTATION). UNPIVOT DOES THIS,
 769 ;;NAMELY, GIVEN A MATRIX A3, INDEX ROW
 770 ;;AND COL ,WHICH CORRESPONDS TO THE Y , U
 771 ;; AND V, RESPECTIVELY, IT COMPUTES V*Y*U
 772 ;;AND RETURNS IT TO THE SAME ARGUMENT A.
 773
 774 (defun tmunpivot (a3 *row* *col* n m)
 775   (prog nil
 776      (setq *col* (get-array-pointer *col*))
 777      (setq *row* (get-array-pointer *row*))
 778      (setq *rowinv* (make-array (1+ n) :initial-element nil))
 779      (setq *colinv* (make-array (1+ n) :initial-element nil))
 780      (do ((i 1 (1+ i)))
 781          ((> i n))
 782        (setf (aref *rowinv* (aref *row* i)) i))
 783      (do ((i 1 (1+ i)))
 784          ((> i n))
 785        (setf (aref *colinv* (aref *col* i)) i))
 786      (tmpermute a3 n m 0 n '*colinv* n '*row*)
 787      (tmpermute a3 n m 0 n '*rowinv* n '*col*)))
 788
 789 (declare-top (unspecial n vlist nx ix))