src/newinv.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
  13 (macsyma-module newinv)
  14
  15 (declare-top (special *ptr* *ptc* *iar* *nonz* detl* *rr*))
  16
  17 (defun multbk (l ax m)
  18   (prog (e)
  19      (do ((j (1+ m) (1+ j)))
  20          ((> j (* 2 m)))
  21        (setq e (car l) l (cdr l))
  22        (do ((i 1 (1+ i))) ((> i m))
  23          (setf (aref ax i j) (rattimes e (aref ax i j) t))))))
  24
  25 (defun ctimemt (x y)
  26   (prog (c)
  27    loop (cond ((null y) (return c)))
  28    (setq c (nconc c (list (timesrow x (car y))))
  29          y (cdr y))
  30    (go loop)))
  31
  32
  33 (defun stora (ax m ei r)
  34   (declare (fixnum m r))
  35   (prog (det (i 0) (j 0) ro mat)
  36      (declare(fixnum i j))
  37      (setq i 0)
  38      loop0 (cond ((null ei) (return nil)))
  39      (setq mat (car ei) ei (cdr ei))
  40      (setq det (caar mat) mat (cdr mat))
  41      loop (setq j r)
  42      (cond ((null mat) (go loop0)))
  43      (setq i (1+ i) ro (car mat) mat (cdr mat))
  44      loop2 (cond ((null ro) (go loop)))
  45      (incf j)
  46      (setf (aref ax i (+ m j)) (ratreduce (caar ro) det))
  47      (setf (aref ax (aref *ptr* i) (aref *ptc* j)) nil)
  48      (setq ro (cdr ro))
  49      (go loop2)))
  50
  51 (defun prodhk (ax ri d r m)
  52   (declare (fixnum r m))
  53   (prog (ei e *rr* *r0 co)
  54      (setq *r0 r ei ri)
  55      loop (cond ((null ei)
  56                  (stora ax m (append co (list d)) r)
  57                  (setq detl* (cons (car d) detl*))
  58                  (return (cons (list d)
  59                                (mapcar #'(lambda (x y) (nconc x (list y)))
  60                                        ri (nreverse *rr*))))))
  61      (setq e (car ei) ei (cdr ei))
  62      (setq co (cons (bmhk ax e d co r detl* *r0) co))
  63      (go loop)))
  64
  65 (defun obmtrx (ax r s i j)
  66   (declare (fixnum r s i j ))
  67   (prog (ans (dj 0) (ds 0) dr d)
  68      (declare(fixnum ds dj))
  69      (setq ds s dj j)
  70      loop (cond ((= i 0) (return ans)))
  71      loop1 (cond ((= j 0)
  72                   (setq j dj
  73                         s ds
  74                         ans (cons (nreverse dr) ans))
  75                   (setq dr nil r (1- r) i (1- i))
  76                   (go loop)))
  77      (setq s (1+ s) j (1- j))
  78      (setq d (aref ax (aref *ptr* r) (aref *ptc* s)))
  79      (cond ((or *nonz* (equal d 0)) nil)
  80            (t (setq *nonz* t)))
  81      (setq dr (cons d dr))
  82      (go loop1)))
  83
  84 (defun bmhk (ax da b nc c0 detl *r0)
  85   (prog (c a sum det dy *nonz* x y)
  86      (setq det (car b) b (cdr b) a (car da) da (cdr da))
  87      (setq nc (reverse nc))
  88      (setq da (reverse da))
  89      (setq c (obmtrx ax *r0 c0 (length(cdr a)) (length b)))
  90      (setq *rr* (cons c *rr*))
  91      (cond ((null *nonz*) (return (cons '(1 . 1) c))))
  92      (setq sum (multmat c b))
  93      (setq *r0 (- *r0 (length (cdr a))))
  94      loop (cond ((null da) (go on)))
  95      (setq x (car da) y(car nc) dy (car y) y (cdr y))
  96      (setq x (multmat x y))
  97      (setq sum (addmatrix1 (ctimemt (cons (pminus (caar detl)) 1) sum) x))
  98      (setq det dy detl (cdr detl))
  99      (setq da (cdr da) nc (cdr nc))
 100      (go loop)
 101      on  (setq det (cons (ptimes (pminus (caar a)) (car det)) 1))
 102      (return (cons det (multmat(cdr a) sum)))))
 103
 104 ;; tmlattice returns the block structure in the form of a list of blocks
 105 ;; each in the form of ((i1 j1) (i2 j2) etc))
 106
 107 (defun newinv (ax m n)
 108   (declare (fixnum m n ))
 109   (prog (j mmat bl d bm detl* dm ipdm dm2 r i ei)
 110      (do ((i m (1- i)))
 111          ((= i 0))
 112        (declare (fixnum i))
 113        (setq mmat (cons (aref ax i (+ i m)) mmat)))
 114      (setq *ptr* (make-array (1+ m)))
 115      (setq  *ptc* (make-array (1+ m)))
 116      (setq bl (tmlattice ax '*ptr* '*ptc* m))
 117      (cond ((null bl) (merror (intl:gettext "newinv: matrix is singular."))))
 118      (setq bl (mapcar #'length bl))
 119      (setq bm (apply #'max bl))         ;Chancey.  Consider mapping.
 120      (setq *iar* (make-array (list (1+ bm) (1+ (* 2 bm)))))
 121      (setq r 0)
 122      loop1 (cond ((null bl)
 123                   (tmunpivot ax '*ptr* '*ptc* m n)
 124                   (return (multbk mmat ax m))))
 125      (setq i (car bl))
 126      (setq dm i)
 127      (setq dm2 (* 2 dm))
 128      loop2 (cond ((= i 0) (go inv)))
 129      (setq j dm2 ipdm (+ i dm))
 130      loop3 (cond ((= j 0) (setq i (1- i)) (go loop2)))
 131      (setf (aref *iar* i j)
 132            (cond ((> j dm)
 133                   (cond ((= j ipdm) '(1 . 1))
 134                         (t '(0 . 1))))
 135                  (t (aref ax (aref *ptr* (+ r i)) (aref *ptc*(+ r j))))))
 136      (decf j)
 137      (go loop3)
 138      inv  (cond ((= r 0)
 139                  (setq ei (tmlin '*iar* dm dm dm))
 140                  (setq ei (list (cons (caar ei) (cdr ei))))
 141                  (stora ax m ei r)
 142                  (setq ei (list ei))(go next)))
 143      (setq d (tmlin '*iar* dm dm dm))
 144      (setq d (cons (caar d) (cdr d)))
 145      (setq ei (prodhk ax ei d r m))
 146      (setq d nil)
 147      next  (incf r (car bl))
 148      (setq bl (cdr bl))
 149      (go loop1)))
 150
 151 (declare-top (unspecial *nonz* detl* *rr*))