SUBROUTINE MATINV(ARRAY,N,DET)
C  
C  MATRIX INVERSION ROUTINE--BEVINGTON,P.302
C
C  ARRAY IS THE SQUARE,SYMMETRIC MATRIX TO BE INVERTED
C  N IS THE ORDER, DET IS THE CALCULATED DETERMINANT
C  THE INVERTED MATRIX REPLACES THE ORIGINAL ONE
C
	DIMENSION	ARRAY(N,N)
	DIMENSION IK(20),JK(20)
	DET	= 1.
	DO 100 K	= 1,N
C
C  FIND LARGEST ELEMENT,REORDER SO IT IS ON THE DIAGONAL
C  PARTIAL PIVOTING
C
	AMAX	= 0.0
21	DO 30 I	= K,N
	DO 30 J	= K,N
	IF(ABS(AMAX) .GT. ABS(ARRAY(I,J))) GO TO 30
	AMAX	= ARRAY(I,J)
	IK(K)	= I
	JK(K)	= J
30	CONTINUE
	IF(AMAX .NE. 0.) GO TO 41
	DET	= 0.
	GO TO 140
41	I	= IK(K)
	IF(I-K) 21,51,43
43	DO 50 J	= 1,N
	SAVE	= ARRAY(K,J)
	ARRAY(K,J)	= ARRAY(I,J)
50	ARRAY(I,J)	= -SAVE
51	J	= JK(K)
	IF(J-K)21,61,53
53	DO 60 I	= 1,N
	SAVE	= ARRAY(I,K)
	ARRAY(I,K)	= ARRAY(I,J)
60	ARRAY(I,J)	= -SAVE
C
C  INVERT MATRIX
C
61	DO 70 I	= 1,N
	IF(I .EQ. K) GO TO 70
	ARRAY(I,K)	= -ARRAY(I,K)/AMAX
70	CONTINUE
71	DO 80 I	= 1,N
	DO 80 J	= 1,N
	IF(I .EQ. K) GO TO 80
	IF(J .EQ. K) GO TO 80
	ARRAY(I,J)	= ARRAY(I,J)+ARRAY(I,K)*ARRAY(K,J)
80	CONTINUE
81	DO 90 J	= 1,N
	IF(J .EQ. K) GO TO 90
	ARRAY(K,J)	= ARRAY(K,J)/AMAX
90	CONTINUE
	ARRAY(K,K)	= 1./AMAX
100	DET	= DET*AMAX
C
C  RESTORE ORDERING OF MATRIX
C
101	DO 130 L	= 1,N
	K	= N-L+1
	J	= IK(K)
	IF(J .LE. K) GO TO 111
	DO 110 I	= 1,N
	SAVE	= ARRAY(I,K)
	ARRAY(I,K)	= -ARRAY(I,J)
110	ARRAY(I,J)	= SAVE
111	I	= JK(K)
	IF(I .LE. K) GO TO 130
	DO 120 J	= 1,N
	SAVE	= ARRAY(K,J)
	ARRAY(K,J)	= -ARRAY(I,J)
120	ARRAY(I,J)	= SAVE
130	CONTINUE
140	RETURN
	END