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(A,B,N,M,D)
       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
       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
	D=1.
       KMAX=N-1
       DO 90 K=1,KMAX
       AMAX=0.
       J2=K
       DO 20 J1=K,N
       IK=(J1-1)*RRW+K
       IF(ABS(AMAX)-ABS(A(IK)))10,20,20
10       AMAX=A(IK)
       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
       SAVE=A(J3)
       A(J3)=A(J4)
       A(J4)=SAVE
40       CONTINUE
       DO 50 J=1,M
       J3=(K-1)*RRW+J
       J4=(J2-1)*RRW+J
       SAVE=B(J3)
       B(J3)=B(J4)
50       B(J4)=SAVE
C       REDUCTION
60       K1=K+1
       KK=(K-1)*RRW+K
	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
70       A(IJ)=A(IJ)-A(KJ)*A(IK)/A(KK)
       DO 80 J=1,M
       IJ=(I-1)*RRW+J
       KJ=(K-1)*RRW+J
80       B(IJ)=B(IJ)-B(KJ)*A(IK)/A(KK)
90       CONTINUE
C       SUBSTITUTE BACK
       NN=(N-1)*RRW+N
C       NN=N*N
	IF(A(NN).EQ.0.)GOTO 999
       DO 110 J=1,M
       NJ=(N-1)*RRW+J
       B(NJ)=B(NJ)/A(NN)
       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
       DO 100 L=I2,N
       IL=(I-1)*RRW+L
       LJ=(L-1)*RRW+J
100       B(IJ)=B(IJ)-A(IL)*B(LJ)
       B(IJ)=B(IJ)/A(II)
111       CONTINUE
110       CONTINUE
       RETURN
999	CONTINUE
	D=0.
	RETURN
       END