C
C       INITIAL SETUP
C
        REAL A(25,25), B(25), C(25), Z(25)
        LOGICAL	SMALLL, BIGP, KN, JN
C
	IC=2
	IP=3
	IT=10
	IK=11
C	**********
	CALL ASSIGN (IP,'MATRIX.OUT')
C	CALL FOPEN (IP,'MATRIX.OUT')
C	**********
500     WRITE(IT,1000)
        READ(IK,1001) NOINV
	IF(NOINV)502,502,501
502	CALL EXIT
501     WRITE(IT,1002)
        READ(IK,1001) N
        IF((N .LE. 25) .AND. (N .GT. 2)) GO TO 510
        WRITE(IT,1004)
        GO TO 500
C
C       CREATE A MATRIX WHOSE INVERS IS KNOWN
C       CACM ALGORITHM 52
C
510     CC=	N*(N+1)*(N+N-5)/6
        D=	1./CC
        A(N,N)=	-D
        DO 521 I=1,N-1
        W=D*I
        A(I,N)=	W
        A(N,I)=	W
        A(I,I)=	D*(CC-I*I)
        IF (I .EQ. 1) GO TO 521
        DO 520 J=1,I-1
        W=-D*I*J
        A(I,J)=	W
520     A(J,I)=	W
521     CONTINUE
C
C       LIST THE ORIGINAL MARTIX
C
        WRITE(IP,1005)
        DO 530 I=1,N
530     WRITE(IP,1006) I, (A(I,J),J=1,N)
      WRITE(1,1009)
        EPSLON=	1.E-30
C
C
C       PERFORM THE REQUESTED NO OF INVERSIONS
C
        DO 600 NT=1,NOINV
C
C       MATRIX INVERSION BY GAUSS-JORDAN ELIMINATION
C       FORTRAN RECODING OF CACM ALGORITHM 120
C
C
        DELTA=	1.
C
C       INITIALIZE INDEX VECTOR
C
        DO 10 J=1,N
10      Z(J)=	J
C
C       MAJOR LOOP
C
        DO 100 I=1,N
        K=	I
        Y=	A(I,I)
        L=	I-1
        SMALLL=	.FALSE.
        IF(L .LT. 1) SMALLL=	.TRUE.
        P=	I+1
        BIGP=	.FALSE.
        IF(P .GT. N) BIGP=	.TRUE.
C
C       FIND PIVOT
C
        IF (BIGP)	GO TO 25
        DO 20 J=IFIX(P),N
        W=	A(I,J)
        IF(ABS(W) .LE. ABS(Y)) GO TO 20
        K=	J
        Y=	W
20      CONTINUE
C
C       DEVELOP DETERMINANT
C
25      DELTA=	DELTA*Y
C
C       CHECK FOR SINGULARITY
C
        IF(ABS(Y) .LT. EPSLON) GO TO 900
        Y=	1./Y
        DO 30 J=1,N
        C(J)=	A(J,K)
        A(J,K)=	A(J,I)
        A(J,I)=	-C(J)*Y
        B(J)=	A(I,J)*Y
30      A(I,J)=	B(J)
        A(I,I)=	Y
        J=	Z(I)
        Z(I)=	Z(K)
        Z(K)=	J
        IF (SMALLL)	GO TO 40
        K1=	1
        K2=	L
        KN=	.TRUE.
        GO TO 50
40      IF (BIGP)	GO TO 100
        K1=	P
        K2=	N
        KN=	.FALSE.
50      DO 90 K=K1,K2
        IF (SMALLL)	GO TO 60
        J1=	1
        J2=	L
        JN=	.TRUE.
        GO TO 70
60      IF (BIGP)	GO TO 90
        J1=	P
        J2=	N
        JN=	.FALSE.
70      DO 80 J=J1,J2
80      A(K,J)=	A(K,J)-B(J)*C(K)
        IF (JN) GO TO 60
90      CONTINUE
        IF (KN) GO TO 40
100     CONTINUE
C
C       RESTORE MATRIX ORDERING
C
        DO 130 I=1,N
110     K=	Z(I)
        IF(K .EQ. I) GO TO 130
        DO 120 J=1,N
        W=	A(I,J)
        A(I,J)=	A(K,J)
120     A(K,J)=	W
        P=	Z(I)
        Z(I)=	Z(K)
        Z(K)=	P
        DELTA=	-DELTA
        GO TO 110
130     CONTINUE
C
C       END OF MATRIX INVERSION
C
600     CONTINUE
C
C
C       LIST THE INVERTED MATRIX AND ITS DETERMINANT
C
      WRITE(IT,1010)
        WRITE(IP,1007) DELTA
        DO 610 I=1,N
610     WRITE(IP,1006) I, (A(I,J),J=1,N)
      WRITE(IT,1011)
      READ(IK,1012)TIME
      WRITE(IP,1013)NOINV,TIME
        GO TO 500
C
C       SINGULARITY
C
900     WRITE(IP,1008) EPSLON
        GO TO 500
C
C
1000    FORMAT(' MATRIX INVERSION'/
     1  ' HOW MANY TIMES?'/)
1001    FORMAT(I2)
1002    FORMAT(' SIZE OF MATRIX?'/)
1004    FORMAT(' MATRIX SIZE MUST BE BETWEEN 3 AND 25')
1005    FORMAT('1 ORIGINAL MATRIX'/)
1006    FORMAT(' ROW 'I2/(5G14.6))
1007    FORMAT('0 AFTER INVERSION, DET IS 'F14.6/)
1008    FORMAT(' MATRIX IS SINGULAR WITHIN EPSILON'G14.6)
 1009 FORMAT(' GO'/)
 1010 FORMAT(' DONE'/)
 1011 FORMAT(' TIME ?'/)
 1012 FORMAT(F6.2)
 1013 FORMAT(1H I2,' INVERSIONS TOOK 'F6.2,' SECONDS'/)
        END