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 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 $linenum $dispflag $linechar $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 PLUGED 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 varlist genvar)) 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;;SHOUD 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 PIVOTTING 476;;SO THAT THE ALL DIAGONAL ELEMENTS ARE 477;;NONZERO. THIS GARANTIES 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 NUBERS 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 COLUMMS. 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 PIVOTTING. 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)) 790