C	LESQ
C	MEAN SQUARE PROGRAM
C	MANUAL IN LESQ.CO
C	REQUIRES DET3 AND LLSQ FROM SSP
C	ALSO USER FUNCTION
	COMMON IFEL,IX,NK,NX,X(30),Y(30),X2(30)
	DIMENSION AUX(12),XM(6),IPIV(6),A(180)
	DIMENSION KOL(5)
	DIMENSION B(36),B2(6,6),YD(7),ARG(10)
	EQUIVALENCE(B2(1,1),B(1))
1300	FORMAT(I5,' OBSERVATIONS')
1500	FORMAT(I3,/(A6))
1201	FORMAT(I3,7(E11.4))
1600	FORMAT(' ERRORS' ,2I4)
1810	FORMAT(14X,1HX,19X,1HY,18X,2HBY,18X,2HDB,/,(4G20.5))
1707	FORMAT(6E13.3)
1705	FORMAT(' VAR/COV:' )
1710	FORMAT('0STANDARD ERRORS:',/,1H ,6G13.3)
1701	FORMAT(' VAR.:' ,/,(6G13.3))
1702	FORMAT(/,'1 CODE=',I5)
1700	FORMAT('0 CONST:' ,/,(6G13.3))
1115	FORMAT ('  OPTION NO ?',$)
1116	FORMAT(' MEAN CHANGE =',G13.3,'  MAX=',G13.3/
     1' CHANGE %' ,/,(6G13.3))
1002	FORMAT(' INPUT ',I5,' CONSTANTS')
1011	FORMAT(7I3)
1010	FORMAT(2I1,5I3)
1000	FORMAT(////' KTGGG...' )
	DATA TMIN,TMAX/0.1,50./
3	REWIND 6
	READ(6,1500)NG,A(1),(A(I),I=1,NG)
31	WRITE(0,1000)
	READ(4,1010)NK,ITYP,KOL
	IF(NK.NE.0)GOTO 5
	CALL FUNCT(ARG,X,Y)
	GOTO 31
5	IF(NK.GT.6)GOTO 31
	IF (KOL(1).EQ.0)GOTO 7
	IX=1
4	READ(6,1201)I,X(IX),Y(IX),X2(IX)
	IF(I.EQ.0)GOTO 7
	DO 6 J=1,5
	IF(KOL(J).EQ.I)GOTO 8
6	CONTINUE
	GOTO 4
8	IX=IX+1
	IF(IX.LT.31)GOTO 4
7	NX=IX-1
	WRITE(0,1300)NX
2	NVARV=5
	WRITE(0,1002)NK
	READ(4,1200)ARG
1200	FORMAT(10F13.0)
	SMAX=200.
	IVAR=1
1	IFEL=0
	DO 100 IX=1,NX
	S=Y(IX)
	IF(ABS(S).GT.1.0E-30)GOTO 10
	H=1.0E-4
	GOTO 15
10	H=S/1000.
15	CONTINUE
	DO 20 J=1,7
	FJ=J
	YC=S+H*(FJ-4.)
20	YD(J)=FUNCT(ARG,X(IX),YC)
	IF(IFEL.NE.0)GOTO 790
	CALL DET3(H,YD,YD,7,IFEL)
	DFDY=YD(4)
	IF(DFDY.EQ.0)GOTO 760
	DO 100 J=1,NK
	A1=ARG(J)
	IF(ABS(A1).GT.1.0E-30)GOTO 30
	H=1.0E-4
	GOTO 35
30	H=A1/1000.
35	CONTINUE
	DO 40 K=1,7
	FJ=K
	ARG(J)=A1+H*(FJ-4.)
40	YD(K)=FUNCT(ARG,X(IX),S)
	IF(IFEL.NE.0)GOTO 790
	CALL DET3(H,YD,YD,7,IFEL)
	DFDK=YD(4)
	I3=IX+(J-1)*NX
	A(I3)=DFDK/(DFDY*S**ITYP)
80	ARG(J)=A1
	B(IX)=-FUNCT(ARG,X(IX),S)/(DFDY*S**ITYP)
100	CONTINUE
	GOTO (110,500),IVAR
110	CALL LLSQ(A,B,NX,NK,1,XM,IPIV,1.0E-7,IFEL,AUX)
	SS=AUX(1)
	IF(IFEL.NE.0)GOTO 780
	A1=0
	H=0
	DO 116 I=1,NK
	YD(I)=XM(I)/ARG(I)*100
	YC=ABS(YD(I))
	A1=AMAX1(A1,YC)
	H=H+YC
116	XM(I)=XM(I)+ARG(I)
	WRITE(0,1700)(XM(I),I=1,NK)
	H=H/NK
	WRITE(0,1116)H,A1,(YD(I),I=1,NK)
	IF(A1.LE.TMIN)GOTO 159
	IF(A1.GE.TMAX.OR.H.GT.1.1*SMAX)GOTO 160
	SMAX=H
	NVARV=NVARV-1
	IF(NVARV.GT.0)GOTO 220
160	WRITE(0,1115)
	READ(4,1011)IVAR,NVARV
C	IOPT A=
C	0	QUIT HERE
C	1	CARRY ON ITERATIONS
C	2	GO TO ENDRUTINE
C	3	READ NEW ARGUMENTS
C	4	READ NEW X/Y
C	5	AS FOR 1 BUT NVARV SET T0 5
	IF (IVAR.LE.0)CALL EXIT
	IF(IVAR.GT.5)GOTO 160
	GOTO (220,1,2,3,260),IVAR
159	IVAR=2
	GOTO 1
260	NVARV=5
	IVAR=1
220	DO 250 I=1,NK
250	ARG(I)=XM(I)
	SMAX=H
	GOTO 1
500	DO 600 I=1,NK
	DO 600 J=1,NK
	S=0.
	DO 520 K=1,NX
	I3=K+(I-1)*NX
	I4=K+(J-1)*NX
520	S=S+A(I3)*A(I4)
	B2(I,J)=S
600	B2(J,I)=S
C	S MINV FROM PERSSON
	I3=NK-1
	N2=0
	DO 611 L=1,NK
	AUX(NK)=1.0/B2(1,1)
	N2=N2+1
	DO 612 K=1,I3
	KP1=K+1
612	AUX(K)=B2(1,KP1)*AUX(NK)
	DO 613 I=1,I3
	IP1=I+1
	B2(I,NK)=-AUX(NK)*B2(IP1,1)
	DO 613 J=1,I3
	JP1=J+1
613	B2(I,J)=B2(IP1,JP1)-AUX(J)*B2(IP1,1)
	DO 611 J=1,NK
611	B2(NK,J)=AUX(J)
	DO 610 I=1,NK
610	ARG(I)=XM(I)
	S=0
	DO 700 IX=1,NX
	YB=FUNCT(ARG,X(IX),Y(IX))
	A(IX)=(Y(IX)-YB)
700	S=S+YB**2/Y(IX)**ITYP

	WRITE(3,1702)ITYP
	WRITE(3,1700)(ARG(I),I=1,NK)
	DFDK=NX-NK-1
	SSD=SS/DFDK
	WRITE(3,1701)SS,S,SSD
	WRITE(0,1701)SS,S,SSD
705	WRITE(3,1705)
	DO 706 I=1,NK
	DO 709 J=1,I
709	YD(J)=B2(I,J)*SS/DFDK
706	WRITE(3,1707)(YD(J),J=1,I)
	DO 710 I=1,NK
710	YD(I)=SQRT(B2(I,I)*SS/DFDK)
	WRITE(3,1710)(YD(I),I=1,NK)
	DO 811 I=1,NX
811	B(I)=(Y(I)-A(I))
	WRITE(3,1810)(X(I),Y(I),A(I),B(I),I=1,NX)
	GOTO 160
C	DFDY=0
760	IFEL2=-1
	IFEL=0
C	ERROR CODES ARE
C	IFEL2		IFEL
C	-1		 DFDY=0
C	1		ERRORS IN LLSQ SEE THIS(P162 SSP)
C	2		ERROR IN FUNCT
999	WRITE(0,1600)IFEL2,IFEL
	GOTO 160
C	ERROR IN LLSQ
780	IFEL2=1
	GOTO 999
C	ERROR IN FUNCT
790	IFEL2=2
	GOTO 999
	END
NS
C	2	GO TO ENDRUTINE
C	3	READ NEW ARGUMENTS
C	4	READ NEW X/Y
C	5	AS FOR 1 BUT NVARV SET T0 5
	IF (IVAR.LE.0)CALL EXIT
	IF(IVAR.GT.5)GOTO 160
	GOTO (220,1,2,3,260),IVAR