;;; -*-  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* *R0 MUL* $SPARSE *DET* *RR* AX))
(declare-top (FIXNUM *R0))

(DEFUN MULTBK (L AX M)
       (PROG (E)
		 (DO ((J (f1+ M) (f1+ J))) ((> J (f* 2 M)))
		 (SETQ E (CAR L) L (CDR L))
	     (DO ((I 1 (f1+ I))) ((> I M))
		     (STORE (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 (f1+ I) RO (CAR MAT) MAT (CDR MAT))
       LOOP2 (COND ((NULL RO) (GO LOOP)))
             (SETQ J (f1+ J))
	     (STORE (AREF AX I (f+ M J)) (RATREDUCE (CAAR RO) DET))
	     (STORE (AREF AX (AREF *PTR* I) (AREF *PTC* J)) NIL)
	     (SETQ RO (CDR RO))
	     (GO LOOP2)))

(DEFUN PRODHK (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 (FUNCTION
					   (LAMBDA (X Y)
					      (NCONC X (LIST Y))))
				      RI (NREVERSE *RR*))))))
	     (SETQ E (CAR EI) EI (CDR EI))
	     (SETQ CO (CONS (BMHK E D CO R DETL*) 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 (f1- R) I (f1- I))
		    (GO LOOP)))
             (SETQ S (f1+ S) J (f1- 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 (DA B NC C0 DETL)
       (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 (f- *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)))))

(DECLARE-TOP(SPECIAL BL))
(COMMENT 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)
  ;(DECLARE (FIXNUM  J R BM DM DM2 ))
	(declare (special bl))			;Why?  I don't know why.  --gsb
	(DO ((I M (f1- I))) ((= I 0))
	    (DECLARE (FIXNUM I))
	    (SETQ MMAT(CONS (AREF AX I (f+ I M)) MMAT)))
	(setq *ptr* (*ARRAY nil T  (f1+ M)))
	(setq  *PTC* (*ARRAY nil T  (f1+ M)))
	(SETQ BL (TMLATTICE AX '*PTR* '*PTC* M))
	(COND ((NULL BL) (merror "SINGULAR")))
	(SETQ BL (MAPCAR 'LENGTH BL))
	(SETQ BM (APPLY 'MAX BL))		;Chancey.  Consider mapping.
	(setq *iar* (*ARRAY nil T (f1+ BM) (f1+ (f* 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 (f* 2 DM))
     LOOP2 (COND ((= I 0)(GO INV)))
	(SETQ J DM2 IPDM(f+ I DM))
     LOOP3 (COND ((= J 0) (SETQ I (f1- I)) (GO LOOP2)))
	(STORE (AREF *IAR* I J)
		(COND ((> J DM)
		       (COND ((= J IPDM) '(1 . 1))
			     (T '(0 . 1))))
		      (T (AREF AX (AREF *PTR* (f+ R I)) (AREF *PTC*(f+ R J))))))
	(SETQ J (f1- 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 EI D R M))
	(SETQ D NIL)
     NEXT  (SETQ R(f+ R (CAR BL)))
	(SETQ BL (CDR BL))
	(GO LOOP1)))

#-NIL
(DECLARE-TOP(UNSPECIAL BL *NONZ* DETL* *R0 MUL* *DET* *RR* AX))