maxima-src 5.32.1-1 / usr / share / maxima / 5.32.1 / src / newinv.lisp

This file is indexed.

/usr/share/maxima/5.32.1/src/newinv.lisp is in maxima-src 5.32.1-1.

This file is owned by root:root, with mode 0o644.

The actual contents of the file can be viewed below. 

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;     The data in this file contains enhancments.                    ;;;;;
;;;                                                                    ;;;;;
;;;  Copyright (c) 1984,1987 by William Schelter,University of Texas   ;;;;;
;;;     All rights reserved                                            ;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;     (c) Copyright 1980 Massachusetts Institute of Technology         ;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(in-package :maxima)

(macsyma-module newinv)

(declare-top (special *ptr* *ptc* *iar* *nonz* detl* *rr*))

(defun multbk (l ax m)
  (prog (e)
     (do ((j (1+ m) (1+ j)))
	 ((> j (* 2 m)))
       (setq e (car l) l (cdr l))
       (do ((i 1 (1+ i))) ((> i m))
	 (setf (aref ax i j) (rattimes e (aref ax i j) t))))))

(defun ctimemt (x y)
  (prog (c)
   loop (cond ((null y) (return c)))
   (setq c (nconc c (list (timesrow x (car y))))
	 y (cdr y))
   (go loop)))


(defun stora (ax m ei r)
  (declare (fixnum m r))
  (prog (det (i 0) (j 0) ro mat)
     (declare(fixnum i j))
     (setq i 0)
     loop0 (cond ((null ei) (return nil)))
     (setq mat (car ei) ei (cdr ei))
     (setq det (caar mat) mat (cdr mat))
     loop (setq j r)
     (cond ((null mat) (go loop0)))
     (setq i (1+ i) ro (car mat) mat (cdr mat))
     loop2 (cond ((null ro) (go loop)))
     (incf j)
     (setf (aref ax i (+ m j)) (ratreduce (caar ro) det))
     (setf (aref ax (aref *ptr* i) (aref *ptc* j)) nil)
     (setq ro (cdr ro))
     (go loop2)))

(defun prodhk (ax ri d r m)
  (declare (fixnum r m))
  (prog (ei e *rr* *r0 co)
     (setq *r0 r ei ri)
     loop (cond ((null ei)
		 (stora ax m (append co (list d)) r)
		 (setq detl* (cons (car d) detl*))
		 (return (cons (list d)
			       (mapcar #'(lambda (x y) (nconc x (list y)))
				       ri (nreverse *rr*))))))
     (setq e (car ei) ei (cdr ei))
     (setq co (cons (bmhk ax e d co r detl* *r0) co))
     (go loop)))

(defun obmtrx (ax r s i j)
  (declare (fixnum r s i j ))
  (prog (ans (dj 0) (ds 0) dr d)
     (declare(fixnum ds dj))
     (setq ds s dj j)
     loop (cond ((= i 0) (return ans)))
     loop1 (cond ((= j 0)
		  (setq j dj
			s ds
			ans (cons (nreverse dr) ans))
		  (setq dr nil r (1- r) i (1- i))
		  (go loop)))
     (setq s (1+ s) j (1- j))
     (setq d (aref ax (aref *ptr* r) (aref *ptc* s)))
     (cond ((or *nonz* (equal d 0)) nil)
	   (t (setq *nonz* t)))
     (setq dr (cons d dr))
     (go loop1)))

(defun bmhk (ax da b nc c0 detl *r0)
  (prog (c a sum det dy *nonz* x y)
     (setq det (car b) b (cdr b) a (car da) da (cdr da))
     (setq nc (reverse nc))
     (setq da (reverse da))
     (setq c (obmtrx ax *r0 c0 (length(cdr a)) (length b)))
     (setq *rr* (cons c *rr*))
     (cond ((null *nonz*) (return (cons '(1 . 1) c))))
     (setq sum (multmat c b))
     (setq *r0 (- *r0 (length (cdr a))))
     loop (cond ((null da) (go on)))
     (setq x (car da) y(car nc) dy (car y) y (cdr y))
     (setq x (multmat x y))
     (setq sum (addmatrix1 (ctimemt (cons (pminus (caar detl)) 1) sum) x))
     (setq det dy detl (cdr detl))
     (setq da (cdr da) nc (cdr nc))
     (go loop)
     on  (setq det (cons (ptimes (pminus (caar a)) (car det)) 1))
     (return (cons det (multmat(cdr a) sum)))))

;; tmlattice returns the block structure in the form of a list of blocks
;; each in the form of ((i1 j1) (i2 j2) etc))

(defun newinv (ax m n)
  (declare (fixnum m n ))
  (prog (j mmat bl d bm detl* dm ipdm dm2 r i ei)
     (do ((i m (1- i)))
	 ((= i 0))
       (declare (fixnum i))
       (setq mmat (cons (aref ax i (+ i m)) mmat)))
     (setq *ptr* (make-array (1+ m)))
     (setq  *ptc* (make-array (1+ m)))
     (setq bl (tmlattice ax '*ptr* '*ptc* m))
     (cond ((null bl) (merror (intl:gettext "newinv: matrix is singular."))))
     (setq bl (mapcar #'length bl))
     (setq bm (apply #'max bl))		;Chancey.  Consider mapping.
     (setq *iar* (make-array (list (1+ bm) (1+ (* 2 bm)))))
     (setq r 0)
     loop1 (cond ((null bl)
		  (tmunpivot ax '*ptr* '*ptc* m n)
		  (return (multbk mmat ax m))))
     (setq i (car bl))
     (setq dm i)
     (setq dm2 (* 2 dm))
     loop2 (cond ((= i 0) (go inv)))
     (setq j dm2 ipdm (+ i dm))
     loop3 (cond ((= j 0) (setq i (1- i)) (go loop2)))
     (setf (aref *iar* i j)
	   (cond ((> j dm)
		  (cond ((= j ipdm) '(1 . 1))
			(t '(0 . 1))))
		 (t (aref ax (aref *ptr* (+ r i)) (aref *ptc*(+ r j))))))
     (decf j)
     (go loop3)
     inv  (cond ((= r 0)
		 (setq ei (tmlin '*iar* dm dm dm))
		 (setq ei (list (cons (caar ei) (cdr ei))))
		 (stora ax m ei r)
		 (setq ei (list ei))(go next)))
     (setq d (tmlin '*iar* dm dm dm))
     (setq d (cons (caar d) (cdr d)))
     (setq ei (prodhk ax ei d r m))
     (setq d nil)
     next  (incf r (car bl))
     (setq bl (cdr bl))
     (go loop1)))

(declare-top (unspecial *nonz* detl* *rr*))

APT Browse - Built by Thomas Orozco.