C COPYRIGHT (C) 1983 GLENN EVERHART
C PERMISSION IS GIVEN TO ANYONE TO USE, DISTRIBUTE, OR COPY THIS
C PROGRAM FREELY BUT NOT TO SELL IT COMMERICALLY.
C       SOLVE THE MATRIX EQUATION AX=B, WHERE A AND B ARE MATRICES,
C       N BY N FOR A, N BY M FOR B
C RETURNS VALUE OF X IN LOCATION USED FOR B ON INPUT.
C        THE STORAGE CONVENTION FOR FORTRAN IS
C       A(COL, ROW), B(COL,ROW).
       SUBROUTINE MTXEQU(A1,A2,B1,B2,N,M,D)
C A1,A2 ARE DIMENSIONS OF A SUBMATRIX ORIGIN IN XVBLS
C B1,B2 ARE DIMS OF B SUBMATRIX
       INCLUDE 'VKLUGPRM.FTN'
C
C NOTE THIS PROGRAM MUST BE MODIFIED TO WORK WITHIN THE SPREAD
C SHEET MATRIX RATHER THAN JUST ASSUMING THAT THE N DIMENSION
C AND M DIMENSION GIVE THE STORAGE ADDRESSES... ALTERNATIVELY,
C THE PROGRAM MUST OPERATE ONLY ON COPIED, DENSELY STORED
C MATRICES.
C
C
C   ORIGINAL PROGRAM TEXT FOLLOWS:
C       DIMENSION A(1),B(1)
CC ALTER DECLARATIONS FOR USE WITH SPREAD SHEET
C	REAL*8 A,B
C       KMAX=N-1
C       DO 90 K=1,KMAX
C       AMAX=0.
C       J2=K
C       DO 20 J1=K,N
C       IK=(J1-1)*N+K
C       IF(ABS(AMAX)-ABS(A(IK)))10,20,20
C10       AMAX=A(IK)
C       J2=J1
C20       CONTINUE
CC       EXCHANGE ROW K,J2 IF NECESSARY
C       IF(J2-K)30,60,30
C30       DO 40 J=K,N
C       J3=(K-1)*N+J
C       J4=(J2-1)*N+J
C       SAVE=A(J3)
C       A(J3)=A(J4)
C       A(J4)=SAVE
C40       CONTINUE
C       DO 50 J=1,M
C       J3=(K-1)*M+J
C       J4=(J2-1)*M+J
C       SAVE=B(J3)
C       B(J3)=B(J4)
C50       B(J4)=SAVE
CC       REDUCTION
C60       K1=K+1
C       KK=(K-1)*N+K
C       DO 80 I=K1,N
C       IK=(I-1)*N+K
C       DO 70 J=K1,N
C       IJ=(I-1)*M+J
C       KJ=(K-1)*M+J
C70       A(IJ)=A(IJ)-A(KJ)*A(IK)/A(KK)
C       DO 80 J=1,M
C       IJ=(I-1)*M+J
C       KJ=(K-1)*N+J
C80       B(IJ)=B(IJ)-B(KJ)*A(IK)/A(KK)
C90       CONTINUE
CC       SUBSTITUTE BACK
CC       NN=(N-1)*N+N
C       NN=N*N
C       DO 110 J=1,M
C       NJ=(N-1)*M+J
C       B(NJ)=B(NJ)/A(NN)
C       I1MAX=N-1
C       IF(I1MAX)110,110,95
C95       DO 111 I1=1,I1MAX
C       I=N-I1
C       IJ=(I-1)*M+J
C       II=(I-1)*N+I
C       I2=I+1
C       DO 100 L=I2,N
C       IL=(I-1)*N+L
C       LJ=(L-1)*M+J
C100       B(IJ)=B(IJ)-A(IL)*B(LJ)
C       B(IJ)=B(IJ)/A(II)
C111       CONTINUE
C110       CONTINUE
C       RETURN
C       END
	INTEGER A1,A2,B1,B2
C       DIMENSION A(1),B(1)
C ALTER DECLARATIONS FOR USE WITH SPREAD SHEET
C NOTE THAT OUR COLUMN DIMENSION IS RRW, NOT N OR M HERE
C SUBSCRIPTS ARE (ROW-1)*COL-DIMENSION + COL
C  THEREFORE, CHANGE *N OR *M IN SUBSCRIPT COMPUTATIONS TO
C  *RRW
	REAL*8 A,B,AW1,AW2,BW1,BW2,AW3,AW4
	INTEGER ABASE,BBASE
	ABASE=(A2-1)*RRW+A1-1
	BBASE=(B2-1)*RRW+B1-1
	D=1.
       KMAX=N-1
       DO 90 K=1,KMAX
       AMAX=0.
       J2=K
       DO 20 J1=K,N
       IK=(J1-1)*RRW+K
	CALL XVBLGT(IK+ABASE,1,A)
       IF(ABS(AMAX)-ABS(A))10,20,20
10       AMAX=A
       J2=J1
20       CONTINUE
C       EXCHANGE ROW K,J2 IF NECESSARY
       IF(J2-K)30,60,30
30       DO 40 J=K,N
       J3=(K-1)*RRW+J
       J4=(J2-1)*RRW+J
	CALL XVBLGT(J3+ABASE,1,SAVE)
C       SAVE=A(J3)
	CALL XVBLGT(J4+ABASE,1,AW1)
	CALL XVBLST(J3+ABASE,1,AW1)
	CALL XVBLST(J4+ABASE,1,SAVE)
C       A(J3)=A(J4)
C       A(J4)=SAVE
40       CONTINUE
       DO 50 J=1,M
       J3=(K-1)*RRW+J
       J4=(J2-1)*RRW+J
C       SAVE=B(J3)
C       B(J3)=B(J4)
C50       B(J4)=SAVE
	CALL XVBLGT(J3+BBASE,1,SAVE)
	CALL XVBLGT(J4+BBASE,1,BW1)
	CALL XVBLST(J3+BBASE,1,BW1)
	CALL XVBLST(J4+BBASE,1,SAVE)
50	CONTINUE
C       REDUCTION
60       K1=K+1
       KK=(K-1)*RRW+K
	CALL XVBLGT(KK+ABASE,1,A)
	IF(A.EQ.0)GOTO 999
C	IF(A(KK).EQ.0.)GOTO 999
       DO 80 I=K1,N
       IK=(I-1)*RRW+K
       DO 70 J=K1,N
       IJ=(I-1)*RRW+J
       KJ=(K-1)*RRW+J
C70       A(IJ)=A(IJ)-A(KJ)*A(IK)/A(KK)
	CALL XVBLGT(IJ+ABASE,1,AW1)
	CALL XVBLGT(KJ+ABASE,1,AW2)
	CALL XVBLGT(IK+ABASE,1,AW3)
	CALL XVBLGT(KK+ABASE,1,AW4)
	AW1=AW1-AW2*AW3/AW4
	CALL XVBLST(IJ+ABASE,1,AW1)
70	CONTINUE
       DO 80 J=1,M
       IJ=(I-1)*RRW+J
       KJ=(K-1)*RRW+J
C80       B(IJ)=B(IJ)-B(KJ)*A(IK)/A(KK)
	CALL XVBLGT(IJ+BBASE,1,BW1)
	CALL XVBLGT(KJ+BBASE,1,BW2)
	BW1=BW1-BW2*AW3/AW4
	CALL XVBLST(IJ+BBASE,1,BW1)
80	CONTINUE
90       CONTINUE
C       SUBSTITUTE BACK
       NN=(N-1)*RRW+N
C       NN=N*N
	CALL XVBLGT(NN+ABASE,1,AW1)
	IF(AW1.EQ.0.)GOTO 999
       DO 110 J=1,M
       NJ=(N-1)*RRW+J
C       B(NJ)=B(NJ)/A(NN)
	CALL XVBLGT(NJ+BBASE,1,BW1)
	BW1=BW1/AW1
	CALL XVBLST(NJ+BBASE,1,BW1)
       I1MAX=N-1
       IF(I1MAX)110,110,95
95       DO 111 I1=1,I1MAX
       I=N-I1
       IJ=(I-1)*RRW+J
       II=(I-1)*RRW+I
       I2=I+1
	CALL XVBLGT(II+ABASE,1,AW1)
       DO 100 L=I2,N
       IL=(I-1)*RRW+L
       LJ=(L-1)*RRW+J
C100       B(IJ)=B(IJ)-A(IL)*B(LJ)
	CALL XVBLGT(IJ+BBASE,1,BW1)
	CALL XVBLGT(IL+ABASE,1,AW2)
	CALL XVBLGT(LJ+BBASE,1,BW2)
	BW1=BW1-AW2*BW2
	CALL XVBLST(IJ+BBASE,1,BW1)
100	CONTINUE
C       B(IJ)=B(IJ)/A(II)
	BW1=BW1/AW1
	CALL XVBLST(IJ+BBASE,1,BW1)
111       CONTINUE
110       CONTINUE
       RETURN
999	CONTINUE
	D=0.
	RETURN
       END
