      DIMENSION A(20,21),X(100,20),Y(100),YP(100)
      CHARACTER XS*10,YS*10,TN*14,MA*1,FN*14
      WRITE(*,*) '****************************************'
      WRITE(*,*) '***       LINEAR LEAST SQUARES       ***'
      WRITE(*,*) '****************************************'
      CALL FMA(M,N,X,Y,FN)
      DO 30 I=1,M
      DO 30 J=1,M
      A(I,J)=0.0
      DO 30 K=1,N
      A(I,J)=A(I,J)+X(K,I)*X(K,J)
   30 CONTINUE
      DO 40 I=1,M
      A(I,M+1)=0.0
      DO 40 J=1,N
      A(I,M+1)=A(I,M+1)+X(J,I)*Y(J)
   40 CONTINUE
      CALL CHORES(A,M)
      WRITE(*,*)'### ESTIMATED PARAMETERS ###'
      DO 50 I=1,M
      WRITE(*,'('' P('',I2,'')='',F10.5)')I,A(I,M+1)
   50 CONTINUE
      DO 60 I=1,N
      YP(I)=0.0
      DO 60 J=1,M
      YP(I)=YP(I)+A(J,M+1)*X(I,J)
   60 CONTINUE
      XB=0.0
      YB=0.0
      XX=0.0
      YY=0.0
      XY=0.0
      SS=0.0
      DO 65 I=1,N
      XB=XB+Y(I)
      YB=YB+YP(I)
      XX=XX+Y(I)*Y(I)
      YY=YY+YP(I)*YP(I)
      XY=XY+Y(I)*YP(I)
      SS=SS+(Y(I)-YP(I))*(Y(I)-YP(I))
   65 CONTINUE
      XB=XB/N
      YB=YB/N
      XX=XX/N
      YY=YY/N
      XY=XY/N
      R=(XY-XB*YB)/SQRT((XX-XB*XB)*(YY-YB*YB))
      AIC=N*ALOG(SS)+2*M
      WRITE(*,'('' R='',F10.5)') R
      WRITE(*,'('' AIC='',F10.5)') AIC
      CALL TMP
      X0=1E+10
      X9=-1E+10
      Y0=1E+10
      Y9=-1E+10
      DO 70 I=1,N
      IF(X0.GT.Y(I)) X0=Y(I)
      IF(X9.LT.Y(I)) X9=Y(I)
      IF(Y0.GT.YP(I)) Y0=YP(I)
      IF(Y9.LT.YP(I)) Y9=YP(I)
   70 CONTINUE
      XS='OBSERVED'
      YS='PRETICTED'
      CALL PLOT2(0,0,XT,YT,X0,X9,Y0,Y9,XS,YS,N,Y,YP)
      CALL TMP
      END

C ----  MODIFIED CHORESKY METHOD ----
      SUBROUTINE CHORES(A,NP)
      DIMENSION A(20,21),SL(20,21)
      IF(NP.EQ.1)THEN
      A(1,2)=A(1,2)/A(1,1)
      RETURN
      ENDIF
      DO 10 JS=1,NP
      SL(1,JS)=A(1,JS)
   10 CONTINUE
      DO 20 IS=2,NP
      DO 20 JS=IS,NP
      DS=0
      DO 40 KS=1,IS-1
      DS=DS+SL(KS,IS)*SL(KS,JS)/SL(KS,KS)
   40 CONTINUE
      SL(IS,JS)=A(IS,JS)-DS
   20 CONTINUE
      SL(1,NP+1)=A(1,NP+1)
      DO 50 IS=2,NP
      DS=0
      DO 60 KS=1,IS-1
      DS=DS+SL(KS,IS)*SL(KS,NP+1)/SL(KS,KS)
   60 CONTINUE
      SL(IS,NP+1)=A(IS,NP+1)-DS
   50 CONTINUE
      A(NP,NP+1)=SL(NP,NP+1)/SL(NP,NP)
      DO 70 IS=NP-1,1,-1
      DS=0
      DO 80 KS=IS+1,NP
      DS=DS+SL(IS,KS)*A(KS,NP+1)
   80 CONTINUE
      A(IS,NP+1)=(SL(IS,NP+1)-DS)/SL(IS,IS)
   70 CONTINUE
      END

      SUBROUTINE PLOT2(NS,NE,XT,YT,X0,X9,Y0,Y9,XS,YS,NP,TX,CY)
      DIMENSION XT(50),YT(50,5),XG(6),TX(250),CY(250)
      CHARACTER XS*10,YS*10,IP(51,17)
      XM=50.0
      YM=16.0
      DO 100 I=1,51
      DO 100 J=1,17
      IP(I,J)=' '
  100 CONTINUE
      DO 110 I=1,NS
      IX=INT(XM*(XT(I)-X0)/(X9-X0)+0.5)+1
      DO 110 J=1,NE
      IY=INT(YM*(Y9-YT(I,J))/(Y9-Y0)+0.5)+1
      IP(IX,IY)=CHAR(J+48)
  110 CONTINUE
      DO 107 I=1,NP
      IX=INT(XM*(TX(I)-X0)/(X9-X0)+0.5)+1
      IY=INT(YM*(Y9-CY(I))/(Y9-Y0)+0.5)+1
      IP(IX,IY)=CHAR(79)
  107 CONTINUE
      WRITE(*,'(5X,A10)') YS
      DO 130 K=1,17
      JJ=K-4*INT(K/4)
      IF(JJ.EQ.1) THEN
      YY=(Y9-Y0)*(16-K+1)/16.0+Y0
      IF((YY.GE.10000.0).OR.((YY.LT.0.1).AND.(YY.NE.0.0))) THEN
      WRITE(*,'(1X,E12.4,''+'',51A1)') YY,(IP(I,K),I=1,51)
      ELSE
      WRITE(*,'(5X,F8.3,''+'',51A1)') YY,(IP(I,K),I=1,51)
      END IF
      ELSE
      WRITE(*,'(13X,''I'',51A1)') (IP(I,K),I=1,51)
      END IF
  130 CONTINUE
      WRITE(*,'(13X,''++'',5A10)') ('---------+',I=1,5)
      DO 140 I=1,6
      XG(I)=(X9-X0)*(I-1)/5.0+X0
  140 CONTINUE
      WRITE(*,'(7X,6F10.2)') (XG(I),I=1,6)
      WRITE(*,'(35X,A10)') XS
      END

      SUBROUTINE TMP
      CHARACTER TP
      WRITE(*,*) 'DO YOU CONTINUE(Y/N)?'
      READ(*,'(A)') TP
      IF(TP.EQ.'N') STOP
      END
 
      SUBROUTINE FMA(M,N,X,Y,FN)
      DIMENSION X(100,20),Y(100)
      CHARACTER FN*14,IC*1
  100 WRITE(*,*) 'E)XIT, I)NPUT, R)EAD, W)RITE, S)TOP ?'
      READ(*,'(A)') IC
      IF(IC.EQ.'S') STOP
      IF(IC.EQ.'E'.AND.N.GT.0) RETURN
      IF(IC.EQ.'I') GOTO 1100
      IF(IC.EQ.'R') GOTO 1200
      IF(IC.EQ.'W') GOTO 1300
      GOTO 100
 1100 WRITE(*,'('' NUMBER OF DATA SETS ?'')')
      READ(*,*) N
      WRITE(*,'('' NUMBER OF PARAMETERS ?'')')
      READ(*,*) M
      DO 1110 I=1,N
      DO 1120 J=1,M
      WRITE(*,'('' X('',I2,'','',I2,'')?'')') I,J
      READ(*,*) X(I,J)
 1120 CONTINUE
      WRITE(*,'('' Y('',I2,'')?'')') I
      READ(*,*) Y(I)
 1110 CONTINUE
      GOTO 100
 1200 WRITE(*,*) 'FILE NAME TO READ ?'
      READ(*,'(A)') FN
      OPEN(UNIT=5,FILE=FN)
      READ(5,*) M,N
      WRITE(*,*) M,N
      DO 1210 I=1,N
      READ(5,*) (X(I,J),J=1,M)
      WRITE(*,*) (X(I,J),J=1,M)
      READ(5,*) Y(I)
      WRITE(*,*) Y(I)
 1210 CONTINUE
      CLOSE(UNIT=5)
      GOTO 100
 1300 WRITE(*,*) 'FILE NAME TO WRITE ?'
      READ(*,'(A)') FN
      OPEN(UNIT=6,FILE=FN,STATUS='NEW')
      WRITE(6,*) M,N
      DO 1310 I=1,N
      WRITE(6,*) (X(I,J),J=1,M)
      WRITE(6,*) Y(I)
 1310 CONTINUE
      CLOSE(UNIT=6)
      GOTO 100
      END
