real*4 mat(10000,10),MTM(10,10),Y(10000),YTM(10)
c 	M is the number of observations, number of rows in mat()
c 	N is the number of variables, number of columns in mat()
c	form i,j th element of mtm by dotting ith into j'th column of mat.
c	form ytm by dotting y unto each column of mat
c	multiply ytm by inverse of mtm, gives best coefficients.

c	form i,j th element of mtm by dotting ith into j'th column of mat.
	do 100 i=1,n
	do 100 j=1,i
	x=0
	do 90 k=1,m
90	x=x+mat(k,i)*mat(k,j)
	mtm(i,j)=x
100	mtm(j,i)=x

c	form ytm by dotting y unto each column of mat
	do 200 i=1,n
	x=0
	do 190 k=1,m
190	x=x+mat(k,i)*y(k)
200	ytm(i)=x
	
c	multiply ytm by inverse of mtm, gives best coefficients.
	call matinv(mtm,n,x)
	do 300 i=1,n
	x=0
	do 290 j=1,n
290	x=x + mtm(i,j)*ytm(j)
300	p(i)=x
	end



	SUBROUTINE MATINV (ARRAY, NORDER, DET)
	DOUBLE PRECISION ARRAY, AMAX, SAVE
2	DIMENSION ARRAY(4,4), IK (10), JK (10)
10	DET = 1
11 	DO 100  K=1, NORDER
C
C	FIND LARGEST ELEMENT ARRAY(I,J) IN REST OF MATRIX
C
	AMAX=0
21 	DO 30 I=K, NORDER
     	DO 30 J=K, NORDER
23 	IF (DABS (AMAX) - DABS (ARRAY(I,J) ) ) 24, 24, 30
24 	AMAX=ARRAY(I,J)
     	IK(K) =I
     	JK(K) =J
30 	CONTINUE
C
C INTERCHANGE ROWS AND COLUMNS TO PUT AMAX IN ARRAY(K,K)
C
31 	IF (AMAX) 41, 32, 41
32 	DET=0.
     	GO TO 140
41 	I=IK (K)
     	IF (I-K) 21, 51, 43
43 	DO 50 J=1, NORDER
     	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, NORDER
     	SAVE= ARRAY(I,K)
     	ARRAY(I,K)= ARRAY(I,J)
60 	ARRAY(I,J)= -SAVE
C
C	ACCUMULATE ELEMENTS OF INVERSE MATRIX
C
61 	DO 70 I=1, NORDER
     	IF(I-K) 63, 70, 63
63 	ARRAY(I,K) = -ARRAY(I,K)  / AMAX
70 	CONTINUE
71 	DO 80 I=1, NORDER
     	DO 80 J=1, NORDER
     	IF (I-K) 74, 80, 74
74 	IF (J-K) 75, 80, 75
75 	ARRAY(I,J) = ARRAY(I,J) + ARRAY(I,K) *ARRAY(K,J)
80 	CONTINUE
81 	DO 90 J=1, NORDER
     	IF (J-K) 83, 90, 83
83 	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, NORDER
       	K= NORDER  - L +1
       	J= IK (K)
       	IF (J-K) 111, 111, 105
105 	DO 110  I=1, NORDER
       	SAVE= ARRAY(I,K)
        ARRAY(I,K)= -ARRAY(I,J)
110 	ARRAY(I,J)= SAVE
111 	I= JK(K)
       	IF (I-K) 130, 130, 113
113 	DO 120  J=1, NORDER
       	SAVE= ARRAY(K,J)
       	ARRAY(K,J) = -ARRAY(I,J)
120 	ARRAY(I,J)= SAVE
130 	CONTINUE
140  	RETURN
	END