.TITLE	MATPAK V1.01 FOR BASIC-11 V02-01
	.IDENT	/V1.01/
	.PSECT	MATPAK,RO
 
; MATRIX PACKAGE FOR DEC RSX-11M BASIC
;
; AUTHOR: C.J. DORAN
;
; COPYRIGHT (C) 1979 SIRA INSTITUTE LTD.,
;		SOUTH HILL, CHISLEHURST, KENT, ENGLAND.
;
; IMPLEMENTS THE NORMAL BASIC MATRIX OPERATIONS THROUGH THE CALL INTERFACE:
;
; MATZER(A(),M%,N%)		ZEROES A(M%,N%)
; MATIDN(A(),M%)		SETS A(M%,M%) TO THE UNIT MATRIX
; MATCON(A(),M%,N%)		SET A(M%,N%) TO A MATRIX OF ALL ONES
; MATTRN(A(),B(),M%,N%)		TRANSPOSES MATRIX B(N%,M%) INTO A(M%,N%)
; MATADD(A(),B(),C(),M%,N%)	A(M%,N%) = B(M%,N%) + C(M%,N%)
; MATSUB(A(),B(),C(),M%,N%)	A(M%,N%) = B(M%,N%) - C(M%,N%)
; MATSCA(A(),B(),M%,N%,S)	MULTIPLIES B(M%,N%) BY SCALAR S GIVING A(M%,N%)
; MATCOP(A(),B(),M%,N%)		COPIES B(M%,N%) TO A(M%,N%)
; MATMUL(A(),B(),C(),M%,P%,N%)	A(M%,N%) = B(M%,P%) * C(P%,N%)
; MULTRA(A(),B(),C(),M%,P%,N%)	A(M%,N%) = B(M%,P%) * ~C(N%,P%)
; TRAMUL(A(),B(),C(),M%,P%,N%)	A(M%,N%) = ~B(P%,M%) * C(N%,P%)
; MATINV(A(),M%)		INVERTS A(M%,M%) IN PLACE -- DET IN A(0%,0%)
; INPROD(R,I%,A(),J%,B(),N%)	INNER PRODUCT
; SUM(R,I%,A(),N%)		SUM
;
; NOTE THAT THE ARGUMENTS MUST BE OF THE CORRECT TYPE, EVEN WHEN LITERALS
; AS:
; 100 DIM A(5,5),B(3,7),C(1,5),D(6,1)
; 110 CALL "MATZER"(A(),5%,5%)
; 120 CALL "MATCON"(C(),1%,5%)
;
; MATRICES MUST BE DIMENSIONED BY A PREVIOUS DIM STATEMENT AND, IF THE
; ELEMENTS ARE TO BE ACCESSED MEANINGFULLY, MUST BE SIZED THE SAME AS
; IN THE MATPAK ROUTINE CALLS. THE ZERO ROW AND COLUMN ARE NOT USED.
; (EXCEPT THAT THE DETERMINANT IS RETURNED IN A(0,0) ON INVERSION.)
; NOTE THEREFORE THAT COLUMN OR ROW VECTORS MUST BE DIMENSIONED (M,1)
; OR (1,N) RESPECTIVELY (AS C AND D IN THE ABOVE EXAMPLE).
;
; THE FIRST DIMENSION REPRESENTS ROWS, AND THE SECOND COLUMNS I.E.
; A(M%,N%) WOULD APPEAR AS:
;
;	A11	A12	A13	...	A1N
;	A21	A22	A23	...	A2N
;	 .	 .	 .	...	 .
;	 .	 .	 .	...	 .
;	AM1	AM2	AM3	...	AMN
;
; ASSEMBLE AS:
; >MAC MATPAK=BMAC,ASSEM,MATPAK
;
; WHERE BMAC.MAC IS SUPPLIED ON THE BASIC SOURCE DISK.
;
; MODIFY AND RE-ASSEMBLE CALLI AS DESCRIBED IN THE BASIC USERS' GUIDE (CHAP 4)
; AND REBUILD AS DESCRIBED IN THE BASIC INSTALLATION GUIDE.
 
F0=%0
F1=%1
F2=%2
F3=%3
F4=%4
 
$LONGE=1		; LONG ERROR MESSAGES
LF=12
CR=15
 
	.GLOBL	ERRARG		; ?ARGUMENT ERROR SUBROUTINE
 
; ARGUMENT OFFSETS IN R5 CALL BLOCKS
 
A=2		; A(%5) -> A(0,0)
B=4		; B(%5) -> B(0,0)
C=6		; C(%5) -> C(0,0)
 
M1=4		; @M1(%5) = M%		} WHEN ONLY
N1=6		; @N1(%5) = N%		} 1 ARRAY
 
M2=6		; @M2(%5) = M%		} WHEN
N2=8.		; @N2(%5) = N%		} 2 ARRAYS
 
M3A=8.		; @M3A(%5) = M%		} ADD &
N3A=10.		; @N3A(%5) = N%		} SUBTRACT
 
M3=8.		; @M3(%5) = M%		}
P3=10.		; @P3(%5) = P% 		} MULTIPLIES
N3=12.		; @N3(%5) = N%		}
 
S=10.		; @S(%5) = S		  MATSCALE
; CLTOCL IS THE NO OF BYTES BETWEEN ADJACENT ELEMENTS (COLUMNS).
CLTOCL=4.
SCALE=2		; 2^SCALE=CLTOCL
; RWTORW IS THE INCREMENT TO SKIP COLUMN 0 WHEN ADVANCING FROM END OF
; ONE ROW TO START OF NEXT.
RWTORW=CLTOCL*2
; DEFINE DEBUG IF DEBUGGING CODE IS TO BE GENERATED
;DEBUG=0
.IIF DF DEBUG $$$=.
 
.MACRO	SAVE			; SAVE REGISTERS, NOTE SPECIAL ORDER
	MOV	%0,-(SP)
	MOV	%1,-(SP)
	MOV	%2,-(SP)
	MOV	%3,-(SP)
	MOV	%4,-(SP)	; %4 LAST
.ENDM
 
.MACRO	RESTORE			; RESTORE REGISTERS SAVED BY ABOVE
	MOV	(SP)+,%4	; %4 FIRST
	MOV	(SP)+,%3
	MOV	(SP)+,%2
	MOV	(SP)+,%1
	MOV	(SP)+,%0
.ENDM
 
; THIS MACRO GENERATES THE ROUTINE HEADER REQUIRED BY BASIC:
;
;	.GLOBL	NAME
; NAME:	.BYTE	1$-0$
; 0$:	.ASCII	"NAME"
; 1$:	.EVEN
;	.WORD	<START ADDRESS>
 
.MACRO	HEADER	NAME,START,?A,?B
	.GLOBL	NAME
NAME:	.BYTE	B-A
A:	.ASCII	"NAME"
B:	.EVEN
.IIF DF DEBUG	.GLOBL	START
	.WORD	START
.ENDM
 
; THIS MACRO CHECKS THAT THE NUMBER OF CALL ARGUMENTS IS RIGHT.
; JUMPS TO ERRARG IF NOT. (IN STANDARD BASIC ERROR HANDLER).
 
.MACRO	ARGCHK	N
	CMPB	@%5,#N
	BEQ	.+6
	JMP	ERRARG
.ENDM
 
; THIS MACRO GETS THE SCALED STRIDE = (N+1)*CLTOCL INTO R.
 
.MACRO	GETSTR	N,R
	MOV	@N(%5),R
	INC	R
	ASH	#SCALE,R
.ENDM
 
	.NLIST	TOC,TTM
	.LIST	ME
 
.PAGE
.SBTTL	MATZER, MATIDN & MATCON
;
; MATOP(A(),M%,N%)
;
; CALLED WITH:
;
;	____________
; %5 ->	| 202 |  3 |
;	+----------+
;	| ->A(0,0) |
;	+----------+
;	|   ->M%   |
;	+----------+
;	|   ->N%   |
;	+----------+
 
	HEADER	MATZER,ZERST
 
ZERST:	CLRF	F0		; STORE ZEROES
	BR	CONA		; PROCESS LIKE CON
 
	HEADER	MATIDN,IDNST
 
IDNST:	ARGCHK	2		; EXPECT 2 ARGUMENTS
	CLRF	F0		; 0'S IN ALL ELEMENTS
	LDF	#^F1.,F1	; NOT ON DIAGONAL, WHEN ONES
	MOV	@M1(%5),%3	; FETCH M%
	MOV	%3,%4		; N%=M% (SQUARE MATRIX)
	MOV	#PIDN,%0	; LOAD SUBROUTINE POINTER
	JMP	MATLP0		; GO TO COMMON LOOP START
 
	HEADER	MATCON,CONST
 
CONST:	LDF	#^F1.,F0	; STORE ONES
CONA:	ARGCHK	3		; EXPECT 3 ARGUMENTS
	MOV	#STF0,%0	; LOAD SUBROUTINE TO DO IT
	JMP	MATLP1		; EXIT ON COMPLETION OF LOOP
 
; ENTRY FOR IDN, TEST I=J AND SET 1 IF SO, 0 OTHERWISE
PIDN:	CMP	%3,%4		; I=J?
	BEQ	STF1		; THEN STORE 1
 
; ENTRY FOR ZER & CON -- STORE CONTENTS OF F0
STF0:	STF	F0,(%1)+	; STORE AND ADVANCE PTR
	RTS	PC		; EXIT
 
; ENTRY FOR IDN WITH I=J -- STORE CONTENTS OF F1
STF1:	STF	F1,(%1)+
	RTS	PC
 
.PAGE
.SBTTL	TRANSPOSE
 
; TRANSPOSE ROUTINE FOR BASIC MATRIX PACKAGE
;
; TRANS(A(),B(),M%,N%)
;
;	+----------+
; %5->  | 202 |  4 |
;	|----------|
;	| ->A(0,0) |
;	|----------|
;	| ->B(0,0) |
;	|----------|
;       |   ->M%   |
;	|----------|
;	|   ->N%   |
;	+----------+
 
	HEADER	MATTRN,TRNST
 
TRNST:	ARGCHK	4		; SHOULD BE 4 ARGS
	MOV	B(%5),%3	; GET ADR OF B(0,0)
	CMP	A(%5),%3	; SAME AS A?
	BEQ	TRANSI		; THEN TRANSPOSE IN PLACE
	GETSTR	M2,%4		; GET STRIDE(B) TO %4
	ADD	%4,%3		; %3->B(1,0)
	TSTF	(%3)+		; %3->B(1,1)
	MOV	%3,-(SP)	; WILL BE B COLUMN PTR
	MOV	%4,-(SP)	; PUSH STRIDE OF B
	MOV	2(SP),-(SP)	; SAVE PTR TO NEXT ELEMENT OF B
 
; WHEN PTRN IS CALLED FROM MATLOOP, STACK LOOKS LIKE:
;
;			+------------+
;		SP->	| PTRN R.A.  |
;			|------------|
;			| WORKSPACE  |
;			|------------|
;			| MATLP R.A. |
;			|------------|
COLPTR=6.;		| B COL PTR  |
;			|------------|
STRIDE=8.;		|  B STRIDE  |
;			|------------|
POINTR=10.;		|  B POINTER |
;			+------------+
 
	MOV	#PTRN,%0	; LOAD SUBROUTINE ADDRESS
	JSR	PC,MATLP2	; CALL LOOP ROUTINE
	ADD	#6.,SP		; PURGE STACK
	RTS	PC		; AND EXIT
 
; SUBROUTINE TO TRANSPOSE ELEMENTS
PTRN:	LDF	@POINTR(SP),F0	; GET NEXT ELEMENT OF B
	STF	F0,(%1)+	; STORE IN A
	CMP	%4,#1		; END OF ROW?
	BEQ	10$		; NEXT ROW IF SO
	ADD	STRIDE(SP),POINTR(SP) ; ELSE UPDATE B PTR 1 COLUMN
	RTS	PC
 
10$:	ADD	#CLTOCL,COLPTR(SP)    ; AT END OF ROW OF A, ADDRESS
	MOV	COLPTR(SP),POINTR(SP) ; NEXT B COLUMN
	RTS	PC
 
; TRANSPOSE IN PLACE, ASSUMING A SQUARE MATRIX
 
TRANSI:	MOV	#PTRNI,%0	; ADDRESS SUBROUTINE
	MOV	@M2(%5),%3	; GET M%
	MOV	@N2(%5),%4	; GET N%
	CMP	%3,%4		; MAKE SURE IT EQUALS M%
	BEQ	.+6
	JMP	ERRARG		; ERROR IF NOT
	CMP	-(SP),-(SP)	; MAKE WORKSPACE ON STACK
 
; STACK IS FILLED ON FIRST ENTRY TO PTRNI, LOCATIONS ARE:
;
;			+------------+
;		SP->	|  @%0 R.A.  |
;			|------------|
;			| WORKSPACE  |
;			|------------|
;			| MATLP R.A. |
;			|------------|
COLPTR=6.;		| B COL PTR  |
;			|------------|
POINTR=8.;		| B POINTER  |
;			+------------+
 
	JSR	PC,MATLP0	; CALL LOOP ROUTINE
	CMP	(SP)+,(SP)+	; PURGE STACK
	RTS	PC		; AND EXIT
 
PTRNI:	CMP	%3,%4		; COMPARE INDICES
	BLT	5$		; SKIP TRANSFER IF BELOW DIAGONAL
	BNE	1$		; GO SWAP IF NOT ON DIAGONAL
	CMP	%4,@N2(%5)	; FIRST TIME?
	BNE	5$		; NO, SKIP AGAIN, NO NEED TO SWAP DIAGONAL
; ON FIRST ENTRY, SET SOURCE POINTERS FROM MATLOOP INITIALISATION
	MOV	%1,COLPTR(SP)	; COLUMN POINTER
	MOV	%1,POINTR(SP)	; ENTRY POINTER
	BR	5$		; UPDATE FOR 1ST ENTRY
 
; SWAP A(I,J) AND A(J,I)
1$:	LDF	(%1),F1		; GET A(I,J)
	LDF	@POINTR(SP),F0	; AND A(J,I)
	STF	F0,(%1)		; SWAP
	STF	F1,@POINTR(SP)
 
; UPDATE POINTERS
5$:	TSTF	(%1)+		; ADVANCE DEST PTR
	CMP	%4,#1		; END-OF-ROW OF DEST?
	BEQ	10$		; BRANCH IF SO
	ADD	%2,POINTR(SP)	; ELSE UPDATE SOURCE PTR (%2=STRIDE)
	RTS	PC		; EXIT PTRNI
 
; AT END OF ROW OF DEST, ADDRESS NEXT SOURCE COLUMN
10$:	ADD	#CLTOCL,COLPTR(SP) ; ADVANCE TO NEXT COLUMN
	MOV	COLPTR(SP),POINTR(SP) ; SET FOR NEXT COLUMN
	RTS	PC		; EXIT PTRNI
.PAGE
.SBTTL MATCOP & MATSCALE
;
; MATOP(A(),B(),M%,N%,{S})
;
;		+----------+
;	%5->	| 202 | 4/5|
;		|----------|
;		| ->A(0,0) |
;		|----------|
;		| ->B(0,0) |
;		|----------|
;		|   ->M%   |
;		|----------|
;		|   ->N%   |
;		|----------|
;		|   ->S    |	{ MATSCA ONLY }
;		+----------+
 
	HEADER	MATSCA,SCAST
 
SCAST:	LDF	@S(%5),F0	; FETCH S
	MOV	#PSCA,%0	; LOAD SCALE SUBROUTINE POINTER
	CMPB	@%5,#5		; SHOULD BE 5 ARGS
	BR	COPA		; COMMON CODE FOR REST
 
	HEADER	MATCOP,COPST
 
COPST:	MOV	#PCOP,%0	; LOAD COPY SUBROUTINE POINTER
	CMPB	@%5,#4		; 4 ARGS
COPA:	BEQ	.+6		; OK IF ARG COUNT RIGHT
	JMP	ERRARG		; ELSE ERROR
	MOV	@N2(%5),-(SP)	; GET (N
	ADD	#2,@SP		;	+2
	ASL	@SP		;	  )*CLTOCL
	ASL	@SP
.IIF EQ CLTOCL-8. 	ASL	@SP
	ADD	B(%5),@SP	; ADDRESS B(1,1)
	JSR	PC,MATLP2	; PERFORM MATRIX LOOP
	TST	(SP)+		; PURGE STACK
	RTS	PC		; AND EXIT
 
; STACK WILL LOOK LIKE:
;
;			+------------+
;		SP->	|  @%0 R.A.  |
;			|------------|
;			| WORKSPACE  |
;			|------------|
;			| MATLP R.A. |
;			|------------|
POINTR=6.;		| B POINTER  |
;			+------------+
 
; SUBROUTINE TO BE PERFORMED EACH STEP OF LOOP
 
; ENTRY POINT FOR MATCOP ROUTINE
PCOP:	LDF	@POINTR(SP),F1	; GET B(I,J)
	BR	STFF1		; GO STORE IT
 
; ENTRY FOR MATSCALE ROUTINE
PSCA:	LDF	@POINTR(SP),F1	; GET B(I,J)
	MULF	F0,F1		; MULTIPLY BY F0=S
 
; COMMON CODE TO STORE RESULT
STFF1:	STF	F1,(%1)+	; STORE AND ADVANCE %1 POINTER
	CMP	%4,#1		; AT THE END OF A ROW?
	BEQ	NEWRW1		; YES, GO ADJUST FOR THAT
NEWCOL:	ADD	#CLTOCL,POINTR(SP) ; NO, NEXT COLUMN=NEXT ELEMENT
	RTS	PC		; EXIT
 
NEWRW1:	ADD	#RWTORW,POINTR(SP) ; SKIP COL 0
	RTS	PC		; AND EXIT
.PAGE
.SBTTL	MATADD & MATSUB
;
; MATOP(A(),B(),C(),M%,N%)
;
;		+----------+
;	%5->	| 202 | 5  |
;		|----------|
;		| ->A(0,0) |
;		|----------|
;		| ->B(0,0) |
;		|----------|
;		| ->C(0,0) |
;		|----------|
;		|   ->M%   |
;		|----------|
;		|   ->N%   |
;		+----------+
 
	HEADER	MATADD,ADDST
 
ADDST:	MOV	#PADD,%0	; ADD PROCEDURE
	BR	SUBA		; CONTINUE LIKE SUB
 
	HEADER	MATSUB,SUBST
 
SUBST:	MOV	#PSUB,%0	; SUBTRACT PROCEDURE
SUBA:	ARGCHK	5		; 5 ARGS
	MOV	@M3A(%5),%3	; GET M
	MOV	@N3A(%5),%4	; AND N
	MOV	%4,-(SP)	; PUSH (N
	ADD	#2,@SP		;	 +2)
	ASL	@SP		;	    *CLTOCL
	ASL	@SP
.IIF EQ CLTOCL-8. 	ASL	@SP
	MOV	@SP,-(SP)	; MAKE 2 COPIES ON STACK
	ADD	B(%5),@SP	; TOP POINTS TO B(1,1)
	ADD	C(%5),2(SP)	; SECOND TO C(1,1)
	JSR	PC,MATLP0	; PERFORM LOOP
	CMP	(SP)+,(SP)+	; PURGE STACK
	RTS	PC		; AND EXIT
 
; STACK CONTENTS IN SUBROUTINES (FILLED BY INIT2):
;
;			+------------+
;		SP->	|  @%0 R.A.  |
;			|------------|
;			| WORKSPACE  |
;			|------------|
;			| MATLP R.A. |
;			|------------|
BPTR=6.;		| B POINTER  |
;			|------------|
CPTR=8.;		| C POINTER  |
;			+------------+
 
; SUBROUTINE FOR ADDITION
PADD:	LDF	@BPTR(SP),F0	; GET ADDEND
	LDF	@CPTR(SP),F1
	ADDF	F1,F0		; ADD
	BR	UPPTRS		; UPDATE POINTERS
 
; SUBROUTINE FOR SUBRACTION
PSUB:	LDF	@BPTR(SP),F0	; GET MINUEND
	LDF	@CPTR(SP),F1	; GET SUBTRAHEND
	SUBF	F1,F0		; SUBTRACT
 
; STORE AND UPDATE POINTERS
UPPTRS:	STF	F0,(%1)+	; STORE RESULT
	CMP	%4,#1		; END OF ROW?
	BEQ	NEWRW2		; GO START NEW ONE
	ADD	#CLTOCL,BPTR(SP) ; ELSE ADVANCE TO NEXT ELEMENT
	ADD	#CLTOCL,CPTR(SP) ; IN BOTH ARRAYS
	RTS	PC
 
; TO START NEW ROW, SKIP COL 0
NEWRW2:	ADD	#RWTORW,BPTR(SP) ; UPDATE B PTR
	ADD	#RWTORW,CPTR(SP) ; AND C PTR
	RTS	PC
.PAGE
.SBTTL	MATMUL, MULTRA & TRAMUL
;
; MATOP(A(),B(),C(),M%,P%,N%)
;
;		+----------+
;	%5->	| 202 | 6  |
;		|----------|
;		| ->A(0,0) |
;		|----------|
;		| ->B(0,0) |
;		|----------|
;		| ->C(0,0) |
;		|----------|
;		|   ->M%   |
;		|----------|
;		|   ->P%   |
;		|----------|
;		|   ->N%   |
;		+----------+
 
; POINTERS:
; PRB	START OF CURRENT ROW/COLUMN OF B
; PRC	  "	"	"	"   OF C
; PC11	-> C(1,1)
; AT END OF EACH ROW OF A
; DPRB	 IS ADDED TO PRB
 
; AFTER EACH ELEMENT OF DESTINATION
; DPRC	IS ADDED TO PRC
 
; INNER LOOP PERFORMED BY INPROD:
; DPB	INCREMENT OF B
; DPC	INCREMENT OF C
; P	NO OF PRODUCTS
;
; INITIAL VALUES REQUIRED:
;
;		FOR MATMUL	    FOR MULTRA		FOR TRAMUL
;		----------	    ----------		----------
; DPC		(N+1)*CLTOCL	    CLTOCL		(N+1)*CLTOCL
; PC11		->C+(N+2)*CLTOCL    ->C+(P+2)*CLTOCL	->C+(N+2)*CLTOCL
; PRC			"		    "			"
; DPB		CLTOCL		    CLTOCL		(M+1)*CLTOCL
; PRB		->B+(P+2)*CLTOCL    ->B+(P+2)*CLTOCL	->B+(M+2)*CLTOCL
; DPRB		(P+1)*CLTOCL	    (P+1)*CLTOCL	CLTOCL
; DPRC		CLTOCL		    (P+1)*CLTOCL	CLTOCL
;
; INIT SETS UP THE STACK FOR MATMUL, AND LEAVES %0=CLTOCL, %3->B & %4->C
; THE STACK WILL NEED TO BE CORRECTED FOR THE OTHER ROUTINES.
.PAGE
; STACK CONTENTS IN PMUL SUBROUTINE:
;
;			+------------+
;		SP->	|     %4     |
;			|------------|
;			|     %3     |
;			|------------|
;			|     %2     |
;			|------------|
;			|     %1     |
;			|------------|
;			|     %0     |
;			|------------|
;			|  @%0 R.A.  |
;			|------------|
;			| WORKSPACE  |
;			|------------|
;			| MATLP R.A. |
; INITIAL ADDRESSES:	|------------|
DPRC=0.;		|    DPRC    |
;			|------------|
DPRB=2.;		|    DPRB    |
;
PRB=4.;			|  ROW B PTR |
;			|------------|
DPB=6.;			|    DPB     |
;			|------------|
PRC=8.;			| C ROW PTR  |
;			|------------|
PC11=10.;		|  ->C(1,1)  |
;			|------------|
DPC=12.;		|    DPC     |
;			+------------+
 
; MATMUL:    A    := B    * C
;             M,N     M,P    P,N
;
 
	HEADER	MATMUL,MULST
 
MULST:	JSR	%1,INIT		; INITIALISE
	BR	COMMON		; JOIN COMMON CODE
 
; MULTRA:    A    := B    * ~C
;             M,N     P,M     P,N
 
	HEADER	MULTRA,MTRAST
 
MTRAST:	JSR	%1,INIT		; INITIALISE
				; CORRECT STACK ENTRIES:
	MOV	DPRB(SP),(SP) 	; DPRC=SCALED STRIDE OF C
	MOV	%0,DPC(SP)	; DPC=CLTOCL
	MOV	PRB(SP),%1	; GET ->B+(P+2)*CLTOCL
	SUB	%3,%1		; AND REPLACE ->B BY
	ADD	%4,%1		; ->C, FORMING ->C+(P+2)*CLTOCL
	MOV	%1,PC11(SP)	; FOR PC11
	MOV	%1,PRC(SP)	; AND PRC
	BR	COMMON		; JOIN COMMON CODE
 
; TRAMUL:    A    := ~B    * C
;             M,N      P,M    P,N
 
	HEADER	TRAMUL,TRAMST
 
TRAMST:	JSR	%1,INIT		; COMMON ENTRY
	MOV	DPRB(SP),DPB(SP) ; CORRECT DPB TO (P+1)*CLTOCL
	MOV	%0,DPRB(SP)	; AND DPRB TO CLTOCL
	GETSTR	M3,%1		; COMPUTE STRIDE OF B
	MOV	%1,DPB(SP)	; FOR DPB
	ADD	#CLTOCL,%1	; ALSO FIND ADDRESS OF B(1,1)
	ADD	%3,%1		; = ->B+(M+1)*CLTOCL+CLTOCL
	MOV	%1,PRB(SP)	; INITIAL B ROW POINTER
				; FALL THROUGH TO COMMON CODE
COMMON:	MOV	#PMUL,%0	 ; LOAD PROC POINTER
	JSR	PC,MATLP3	 ; CALL LOOP ROUTINE
	ADD	#DPC+2,SP	; PURGE STACK
	RTS	PC		; AND EXIT
 
; THIS ROUTINE IS CALLED FOR EVERY ELEMENT OF A.
; IT USES INPROD TO COMPUTE THE INNER PRODUCT, THEN
; UPDATES POINTERS TO B & C ACCORDING TO DPB, DPRB, DPC, & DPRC.
;
; WE MUST INCREASE THE STACK OFFSETS DEFINED ABOVE TO ALLOW FOR SAVED
; REGISTERS AND NORMAL STACK USAGE OF MATLOOP.
OFFSET=16.
DPC=DPC+OFFSET
PC11=PC11+OFFSET
PRC=PRC+OFFSET
DPB=DPB+OFFSET
PRB=PRB+OFFSET
DPRB=DPRB+OFFSET
DPRC=DPRC+OFFSET
 
PMUL:	SAVE			; SAVE REGISTERS
	MOV	DPB(SP),%1	; INCREMENT OF B
	MOV	PRB(SP),%2	; START ADDRESS IN B
	MOV	DPC(SP),%3	; INCREMENT OF C
	MOV	PRC(SP),%4	; START ADDRESS IN C
	MOV	@P3(%5),%0	; SUM P TIMES
 
; THE FOLLOWING SPECIAL ENTRY TO INPROD SUMS INTO F0:
; %0		PRODUCTS ADDRESSED BY:
; %1 & %3	WITH INCREMENTS OF
; %2 & %4	RESPECTIVELY AT EACH STEP
	JSR	PC,INPRDA	; INPROD
	CMP	@SP,#1		; END OF ROW OF A?
	BEQ	10$		; BRANCH IF SO
	ADD	DPRC(SP),PRC(SP) ; ELSE JUST UPDATE C POINTER
	BR	20$		; GO STORE RESULT
; END OF ROW
10$:	ADD	DPRB(SP),PRB(SP) ; UPDATE B POINTER
	MOV	PC11(SP),PRC(SP) ; RESET C POINTER
20$:	RESTORE			; REGISTERS
	STF	F0,(%1)+	; STORE RESULT
	RTS	PC
 
; COMMON ENTRY CODE. SETS UP POINTERS FOR MATMUL (SEE ABOVE).
; THESE MAY NEED TO BE ALTERED ACCORDING TO ROUTINE ENTERED.
;
; ALSO RETURNS: %0=CLTOCL, %3->B(0,0), %4->C(0,0)
 
INIT:	CMPB	@%5,#6		; 6 ARGS?
	BNE	ERR		; NO, ERROR
	MOV	B(%5),%3	; FETCH ADDRESS OF B(0,0)
	CMP	A(%5),%3	; SAME AS A?
	BEQ	ERR		; ARRAYS MUST BE DIFFERENT
	MOV	C(%5),%4	; GET ADDRESS OF C(0,0) TOO
	CMP	A(%5),%4	; AND CHECK C<>A
	BEQ	ERR
	MOV	#CLTOCL,%0	; LOAD %0
	MOV	@N3(%5),@SP	; GET STRIDE OF C
	INC	@SP		; = (N%+1%)
	ASL	@SP		;	   *CLTOCL
	ASL	@SP		; [DPC]
.IIF EQ CLTOCL-8. 	ASL	@SP
	MOV	@SP,-(SP)	; NEXT ITEM ON STACK
	ADD	%4,@SP		; IS POINTER TO
	ADD	%0,@SP		;   C(1,1) [PC11]
	MOV	@SP,-(SP)	; NEXT IS [PRC], = PC11 INITIALLY
	MOV	%0,-(SP)	; NEXT IS [DPB], CLTOCL FOR MATMUL & MULTRA
	MOV	@P3(%5),-(SP)	; ALSO GET (P
	INC	@SP		;	     +1)
	ASL	@SP		;		*CLTOCL
	ASL	@SP		; [DPB]
.IIF EQ CLTOCL-8. 	ASL	@SP
	MOV	@SP,-(SP)	; MAKE ANOTHER COPY
	ADD	%3,2(SP)	; TO PUSH POINTER TO B(1,1)
	ADD	%0,2(SP)	; [PRB]
				; LEAVE [DPRB] INITIALISED TO STRIDE OF B
	MOV	%0,-(SP)	; FINALLY, [DPRC]
	JMP	@%1		; RETURN
 
ERR:	JMP	ERRARG		; ERROR JUMP
.PAGE
.SBTTL	MATRIX INVERSION IN PLACE.
; CACM ALGORITHM 140.
 
; MATINV(A(),M%)
;
;		+----------+
;	%5->	| 202 | 2  |
;		|----------|
;		| ->A(0,0) |
;		|----------|
;		|   ->M%   |
;		+----------+
;
; SET UP POINTERS ON STACK:
;
;			+------------+
;		SP->	|     I%     |
;			|------------|
A1I=2.;			|  ->A(1,I)  |
;			|------------|
AI1=4.;			|  ->A(I,1)  |
;			|------------|
AII=6.;			|  ->A(I,I)  |
;			|------------|
A11=8.;			|  ->A(1,1)  |
;			|------------|
STRIDE=10.;		|   STRIDE   |
;			+------------+
 
.IIF DF DEBUG .=$$$+2404
	HEADER	MATINV,INVST
 
INVST:	ARGCHK	2		; CHECK 2 ARGUMENTS
	MOV	@M1(%5),-(SP)	; FETCH M%
	INC	@SP		; FIND STRIDE OF A
	ASL	@SP
	ASL	@SP
.IIF EQ CLTOCL-8. 	ASL	@SP
	MOV	@SP,-(SP)	; COPY IT
	ADD	#CLTOCL,@SP	; ADD ANOTHER ELEMENT STEP TO
	ADD	A(%5),@SP	; INITIALISE A11 -> A(1,1)
	MOV	@SP,-(SP)	; AND AII -> A(I,I)
	MOV	@SP,-(SP)	; AND AI1 -> A(I,1)
	MOV	@SP,-(SP)	; AND A1I -> A(1,I)
 
	LDF	#^F1.,F0	; DET=1.0
	STF	F0,F4		; F4= CONST 1.0
	MOV	@M1(%5),-(SP)	; I=M ON STACK
; J = %4
; K = %0
 
ILOOP:	MOV	AII(SP),%4	; GET PTR TO A(I,I)
	LDF	@%4,F1		; GET A(I,I)
	STF	F1,F2		; COPY IT
	ABSF	F2		; TO INSPECT ABS VALUE
	CMPF	#^F1E-7,F2	; C.F. 1E-7
	CFCC
	BGT	FAILUR		; FAIL IF LESS
	MULF	F1,F0		; DET:=DET*Q
	LDF	F4,F2		; LOAD 1.0
	STF	F2,(%4)+	; A(I,I):=1.0
	ADD	STRIDE(SP),%4	; UPDATE AII POINTER
	MOV	%4,AII(SP)	; FOR NEXT TIME
	CMPF	F4,F1		; Q=1.0?
	CFCC
	BEQ	20$		; SKIP NEXT LOOP IF SO
 
; FOR J:=1 STEP 1 UNTIL N DO A(I,J):=A(I,J)/Q
	MOV	@M1(%5),%4	; FETCH M
	MOV	AI1(SP),%1	; GET ADDRESS TO A(I,J)
10$:	LDF	@%1,F2		; GET A(I,K)
	DIVF	F1,F2		; FORM A(I,K)/Q
	STF	F2,(%1)+	; STORE BACK & ADV PTR
	SOB	%4,10$		; REPEAT FOR ALL N
 
; FOR J:=1 STEP 1 UNTIL M DO
20$:	MOV	@M1(%5),%4	; LOAD M
	MOV	A11(SP),%1	; %1->A(J,1)  (J=1)
	MOV	A1I(SP),%3	; %3->A(J,I)
21$:	MOV	AI1(SP),%2	; %2->A(I,1)
	CMP	%4,@SP		; I=J?
	BEQ	30$		; SKIP NEXT SECTION IF I=J
	LDF	@%3,F1		; Q:=A(J,I)
	CLRF	@%3		; A(J,I):=0
 
; FOR K:=1 STEP 1 UNTIL M DO A(J,K):=A(J,K)-Q*A(I,K)
	MOV	%1,-(SP)	; SAVE A(J,1) POINTER
	MOV	@M1(%5),%0	; GET M
25$:	LDF	@%1,F3		; GET A(J,K)
	LDF	(%2)+,F2	; GET A(I,K)
	MULF	F1,F2		; FORM Q*A(I,K)
	SUBF	F2,F3		; A(J,K)-Q*A(I,K)
	STF	F3,(%1)+	; STORE OVER A(J,K)
	SOB	%0,25$		; FOR K=1 TO N
	MOV	(SP)+,%1	; RESTORE A(J,1)
 
; UPDATE POINTERS FOR NEXT J
30$:	ADD	STRIDE(SP),%1	; A(J,1)
	ADD	STRIDE(SP),%3	; A(J,I)
	SOB	%4,21$		; FOR J=1 TO N
 
; END OF I LOOP
	DEC	@SP		; MORE?
	BLE	SUCCES		; NO, FINISH
; UPDATE POINTERS FOR NEXT I
	ADD	STRIDE(SP),AI1(SP) ; A(I,1)
	ADD	#CLTOCL,A1I(SP)    ; A(1,I)
				   ; A(I,I) DONE ALREADY
	BR	ILOOP		   ; LOOP FOR MORE
 
 
FAILUR:	ERCALL	NIM,<NON-INVERTABLE MATRIX>
	CLRF	F0		; DET:=0
 
SUCCES:	STF	F0,@A(%5)	; STORE DET IN A(0,0)
	ADD	#STRIDE+2,SP	; PURGE STACK
	RTS	PC		; AND EXIT
.IF DF DEBUG
 
MSGERR:	TSTB	(%1)+
	BNE	MSGERR
	INC	%1
	BIC	#1,%1
	RTS	%1
.ENDC
.PAGE
.SBTTL INNER PRODUCT AND SUM
 
; INNER PRODUCT (INPROD)
;
; INPROD(R,I%,A(),J%,B(),N%)
;
; RETURNS IN R THE SUM OF [A]*[B] FOR N TERMS, WITH THE A POINTER
; INCREMENTED I% ELEMENTS EACH TIME, AND THE B POINTER J% ELEMENTS.
;
;		+----------+
;	%5->	| 202 | 6  |
;		|----------|
; 2.		|  -> R    |
;		|----------|
; 4.		|  -> I%   |
;		|----------|
; 6.		| ->A(X,Y) |
;		|----------|
; 8.		|  -> J%   |
;		|----------|
; 10.		| ->B(U,V) |
;		|----------|
; 12.		|   ->N%   |
;		+----------+
 
	HEADER	INPROD,INPRST
 
INPRST: ARGCHK	6		; MAKE SURE 6 ARGS
	CMP	(%5)+,(%5)+	; ADVANCE PTR
	MOV	@(%5)+,%1	; %1 = I%
	ASH	#SCALE,%1	; SCALE
	MOV	(%5)+,%2	; %2 -> A
	MOV	@(%5)+,%3	; %3 = J%
	ASH	#SCALE,%3	; SCALE
	MOV	(%5)+,%4	; %4 -> B
	MOV	@(%5)+,%0	; %0 = N%
	SUB	#14.,%5		; RESET POINTER
	JSR	PC,INPRDA	; CALL LOCAL ROUTINE
	STF	F0,@2(%5)	; STORE RESULT IN R
	RTS	PC		; AND EXIT
 
INPRDA:	SETD			; SUM IN DOUBLE-PRECISION
	CLRD	F0		; CLEAR RESULT
	TST	%0		; TEST N
	BLE	EXIT1		; EXIT IF <=0
10$:	LDCFD	@%2,F1		; LOAD 1ST HALF OF TERM
	LDCFD	@%4,F2		; AND 2ND, EXTENDING PRECISION
	MULD	F2,F1		; MULTIPLY IN D.P.
	ADDD	F1,F0		; DO SUM IN D.P.
	ADD	%1,%2		; UPDATE POINTERS
	ADD	%3,%4
	SOB	%0,10$		; DO ALL N TERMS
EXIT1:	SETF			; BACK TO S.P.
	RTS	PC		; AND RETURN
 
; SUM FUNCTION
;
; SUM(R,I%,A(),N%)
;
; SUMS [A] N TIMES, WITH THE A POINTER INCREMENTED BY I% ELEMENTS EACH TIME.
;
;		+----------+
;	%5->	| 202 | 4  |
;		|----------|
; 2.		|  -> R    |
;		|----------|
; 4.		|  -> I%   |
;		|----------|
; 6.		| ->A(X,Y) |
;		|----------|
; 8.		|   ->N%   |
;		+----------+
 
	HEADER	SUM,SUMST
 
SUMST:	ARGCHK	4		; 4 ARGS
	MOV	@4(%5),%1	; %1 = I%
	ASH	#SCALE,%1	; SCALE IT
	MOV	6(%5),%2	; %2 -> A
	MOV	@8.(%5),%3	; %3 = N%
	SETD			; DO SUMS IN D.P.
	CLRD	F0		; CLEAR RESULT
	BLE	EXIT2		; EXIT IF N%<=0
10$:	LDCFD	@%2,F1		; LOAD WORKING REGISTER -- S.P.->D.P.
	ADDD	F1,F0		; ADD TO RESULT
	ADD	%1,%2		; UPDATE POINTER
	SOB	%3,10$		; DO N TIMES
EXIT2:	SETF			; BACK TO S.P.
	STF	F0,@2(%5)	; STORE RESULT
	RTS	PC		; AND EXIT
.PAGE
.SBTTL	MATLOOP
 
; BASIC LOOP ROUTINE FOR MATRIX PACKAGE
; LOOPS THROUGH THE ARRAY ELEMENTS CALLING THE PROCEDURE ADDRESSED BY %0
; EACH TIME. THIS PROCEDURE SHOULD COMPUTE A FLOATING RESULT AND STORE
; IT AT @%1, ADVANCING %1. OTHER REGISTERS MUST BE PRESERVED.
;
; MATLOOP(A,M%,N%,F)
 
; COMPUTES ADDRESS OF A(1,1) IN %1, STRIDE OF A IN %2.
; OPTIONALLY LOADS %3 AND %4 WITH M% AND N% RESPECTIVELY.
; ENTER AT VARIOUS POINTS, DEPENDING ON CALLER:
 
; ENTRY FROM ROUTINES WITH 1 ARRAY
MATLP1:	MOV	@M1(%5),%3	; GET M%
	MOV	@N1(%5),%4	; AND N%
	BR	MATLP0		; BRANCH
 
; ROUTINES WITH 2 ARRAYS
MATLP2:	MOV	@M2(%5),%3
	MOV	@N2(%5),%4
	BR	MATLP0
 
; MULTIPLY ROUTINES WITH 3 ARRAYS
MATLP3:	MOV	@M3(%5),%3
	MOV	@N3(%5),%4
 
; ENTRY IF M% & N% HAVE ALREADY BEEN LOADED, JUST FIND PTR AND STRIDE.
MATLP0:	TST	%3		; DON'T DO ANYTHING
	BLE	EXIT3		; IF M%
	TST	%4		; OR N%
	BLE	EXIT3		; IS <= 0
	MOV	A(%5),%1	; GET -> A(0,0)
	MOV	%4,%2		; STRIDE IS
	INC	%2		; (N%+1%)
	ASH	#SCALE,%2	; *CLTOCL
	ADD	%2,%1		; %1 -> A(1,0)
	TSTF	(%1)+		; %1 -> A(1,1)
 
JLOOP:	MOV	%4,-(SP)	; SAVE N
10$:	JSR	PC,@%0		; CALL THE PROCEDURE
	SOB	%4,10$		; FOR ALL N
	TSTF	(%1)+		; UPDATE PTR TO NEXT ROW
				; (IGNORE COLUMN 0)
	MOV	(SP)+,%4	; GET N BACK
	SOB	%3,JLOOP	; DO ALL ROWS
EXIT3:	RTS	PC		; AND EXIT
 
	.END