.TITLE	MULM

;.MULM	X,I,J,K,Y,L,M,N,Z	;matrix multiplication: Z = X*Y

;X has dimension IxJxK, and Y has dimensions LxMxN. If a matrix has
;fewer than 3 dimensions, put in #1 for the missing dimension(s).
;One and only one of I, J or K must be a negative number, and
;one of L, M or N must be negative, in order to indicate
;that the sum of the products will be taken over these two indices.

;Example 1:	.MULM	#A,#2,#-3,#1,#B,#-3,#4,#1,#C
;		...
;	A:	.FLT2	1,2,3		;2X3 matrix
;		.FLT2	4,5,6
;	B:	.FLT2	7,8,9,10	;3X4 matrix
;		.FLT2	11,12,13,14
;		.FLT2	15,16,17,18
;	C:	.BLKB	4*2*4		;C must of course be 2*4

;The two #-3 numbers must be identical, of course. The computer calculates
;similarly to the following PASCAL program:

;FOR I:= 1 TO 2 DO
;	FOR K:= 1 TO 4 DO
;		BEGIN
;		C[I,K]:= 0;
;		FOR J:= 1 TO 3 DO  C[I,K]:= C[I,K] + A[I,J]*B[J,K];
;		END;

;Example 2:	.MULM	#X,#-2,#3,#4,#Y,#5,#-2,#1,#Z

;X has dimension 2x3x4.  Y has dimension 5x2.  Thus Z will be 3x4x5. The
;PASCAL equivalent is:

;FOR J:= 1 TO 3 DO
;	FOR K:= 1 TO 4 DO
;		FOR L:= 1 TO 5 DO
;			BEGIN
;			Z[J,K,L]:= 0;
;			FOR I:=1 TO 2 DO Z[J,K,L]:=Z[J,K,L] + X[I,J,K]*Y[L,I];
;			END;

;Example 3:	.MULM	#P,#2,#3,#-1,#Q,#-1,#1,#1,#R
;		...
;	P:	.FLT2	1,2,3	;2x3x1 matrix, that is, 2x3
;		.FLT2	4,5,6
;	Q:	.FLT2	2	;1 dimensional with one element (a scalar)
;	R:	.BLKB	4*2*3	;2x3 matrix

;This instruction multiplies matrix P by scalar 2.
	.GLOBL	.MULM,.CALL
	.MCALL	.TLQ
A0=R0
A1=R1
D1:	.BLKW	1	;ORDER IS DIFFERENT THAN USUAL
D2:	.BLKW	1
D3:	.BLKW	1
D4:	.BLKW	1
D5:	.BLKW	1
D6:	.BLKW	1
Z:	.BLKW	1	;ADDRESS OF Z
Y:	.BLKW	1	;ADDR OF Y
X:	.BLKW	1	;ADDR OF X
.MULM:	MOV	R0,-(SP)
	MOV	R1,-(SP)
	MOV	R2,-(SP)
	MOV	R3,-(SP)
	MOV	R4,-(SP)
	MOV	R5,-(SP)
	STF	A0,-(SP)
	STF	A1,-(SP)
	MOV	#D1,R0		;SCAN FOR ERRORS
	MOV	#2,R2
LOOP22:	CLR	R1
	MOV	#3,R3
LOOP33:	TST	(R0)+
	BEQ	ERROR
	BGT	END3
	INC	R1
END3:	SOB	R3,LOOP33
	CMP	R1,#1
	BNE	ERROR
	SOB	R2,LOOP22
	BR	ARRANG
ERROR:	.TLQ	<.MULM DIMENSION ERROR>
	JMP	.CALL

ARRANG:	MOV	#I1,R0
	MOV	#D1,R1
	MOV	#AIX,R2
	MOV	#ADX,R3
	MOV	#AI1,R4
	MOV	#AD1,R5
	CLR	DO6
LOOP6:	TST	(R1)
	BEQ	ERROR
	BLT	LT
	MOV	R0,(R4)+
	MOV	R1,(R5)+
	BR	NEXT
LT:	NEG	(R1)
	MOV	R0,(R2)+
	MOV	R1,(R3)+
NEXT:	ADD	#2,R0
	ADD	#2,R1
	INC	DO6
	CMP	DO6,#6
	BLT	LOOP6
	CMP	@ADX,@ADY
	BNE	ERROR

	MOV	Z,R0
	CLR	@AI1
LOOP1:	CLR	@AI2
LOOP2:	CLR	@AI3
LOOP3:	CLR	@AI4
LOOP4:	CLR	@AIX
	CLR	@AIY
	CLRF	A1
LOOPS:	MOV	I1,R5		;GET X(I1,I2,I3)
	MUL	D2,R5
	ADD	I2,R5
	MUL	D3,R5
	ADD	I3,R5
	ASH	#2,R5
	ADD	X,R5
	LDF	@R5,A0
	MOV	I4,R5		;GET Y(I4,I5,I6)	
	MUL	D5,R5
	ADD	I5,R5
	MUL	D6,R5
	ADD	I6,R5
	ASH	#2,R5
	ADD	Y,R5
	MULF	@R5,A0
	ADDF	A0,A1

	INC	@AIX
	INC	@AIY
	CMP	@AIX,@ADX
	BLT	LOOPS

	STF	A1,(R0)+
	INC	@AI4
	CMP	@AI4,@AD4
	BLT	LOOP4

	INC	@AI3
	CMP	@AI3,@AD3
	BLT	LOOP3

	INC	@AI2
	CMP	@AI2,@AD2
	BLT	LOOP2

	INC	@AI1
	CMP	@AI1,@AD1
	BLT	LOOP1
	LDF	(SP)+,A1
	LDF	(SP)+,A0
	MOV	(SP)+,R5
	MOV	(SP)+,R4
	MOV	(SP)+,R3
	MOV	(SP)+,R2
	MOV	(SP)+,R1
	MOV	(SP)+,R0
	RETURN
I1:	.BLKW	1
I2:	.BLKW	1
I3:	.BLKW	1
I4:	.BLKW	1
I5:	.BLKW	1
I6:	.BLKW	1
AD1:	.BLKW	1
AD2:	.BLKW	1
AD3:	.BLKW	1
AD4:	.BLKW	1

ADX:	.BLKW	1
ADY:	.BLKW	1

AI1:	.BLKW	1
AI2:	.BLKW	1
AI3:	.BLKW	1
AI4:	.BLKW	1

AIX:	.BLKW	1
AIY:	.BLKW	1

DO6:	.BLKW	1
	.END