C     LSMOD
C
C     IBM PC - D.L. SMITH AND R.T. MAINARDI - 17 JULY 1988
C
C     REVISED 12/18/92. ALTERED I/O AND LEAST-SQUARES ROUTINE.
C
      DIMENSION Y(50),EY(50),CY(50,50),VY(50,50),VYI(50,50),
     1P(10),EP(10),CP(10,10),VP(10,10),VPI(10,10),A(50,10),
     2QN(50),WN(50,51),QM(10),WM(10,11)
C
C     INITIALIZATION AND CONTROL
C
      WRITE(*,1)
    1 FORMAT(' LSMOD'/)
      WRITE(*,2)
    2 FORMAT(' ENTER INPUT FILE NAME (UNIT 4)'/)
      OPEN(4,FILE=' ',STATUS='OLD ')
      WRITE(*,3)
    3 FORMAT(' ENTER OUTPUT FILE NAME (UNIT 5)'/)
      OPEN(5,FILE=' ',STATUS='OLD ')
C
C     READ INPUT FROM FILE (UNIT 4)
C
      READ(4,10) N,M
   10 FORMAT(16I5)
      READ(4,12) (Y(I),I=1,N)
   12 FORMAT(6E12.6)
      DO 17 I=1,N
   17 READ(4,12) (VY(I,J),J=1,I)
      DO 18 I=1,N
      DO 18 J=1,I
   18 VY(J,I)=VY(I,J)
      DO 19 I=1,N
   19 READ(4,12) (A(I,J),J=1,M)
C
C     ORDINARY LEAST-SQUARES ANALYSIS
C
      DO 1004 I=1,N
 1004 EY(I)=SQRT(VY(I,I))
      DO 1005 I=1,N
      DO 1005 J=1,N
 1005 CY(I,J)=VY(I,J)/EY(I)/EY(J)
      CALL MATINV(QN,VY,VYI,WN,NTEST,N,50,51)
      IF(NTEST.EQ.1) GO TO 1009
 1007 WRITE(*,1008)
 1008 FORMAT(' NO INV'/)
      STOP
 1009 DO 1010 I=1,M
      DO 1010 J=1,M
      VPI(I,J)=0.0
      DO 1010 K2=1,N
      DO 1010 K1=1,N
 1010 VPI(I,J)=VPI(I,J)+A(K2,I)*VYI(K2,K1)*A(K1,J)
      CALL MATINV(QM,VPI,VP,WM,NTEST,M,10,11)
      IF(NTEST.EQ.0) GO TO 1007
      DO 1011 I=1,M
      P(I)=0.0
      DO 1011 K3=1,M
      DO 1011 K2=1,N
      DO 1011 K1=1,N
 1011 P(I)=P(I)+VP(I,K3)*A(K2,K3)*VYI(K2,K1)*Y(K1)
      DO 1012 I=1,M
 1012 EP(I)=SQRT(VP(I,I))
      DO 1013 I=1,M
      DO 1013 J=1,M
 1013 CP(I,J)=VP(I,J)/EP(I)/EP(J)
      DO 1014 I=1,N
      QN(I)=Y(I)
      DO 1014 K1=1,M
 1014 QN(I)=QN(I)-A(I,K1)*P(K1)
      CHI2=0.0
      DO 1015 K2=1,N
      DO 1015 K1=1,N
 1015 CHI2=CHI2+QN(K2)*VYI(K2,K1)*QN(K1)
      CHI2NM=CHI2/FLOAT(N-M)
C
C     PRINT OUTPUT TO FILE (UNIT 5)
C
      WRITE(5,40)
   40 FORMAT(' Y')
   41 WRITE(5,12) (Y(I),I=1,N)
      WRITE(5,410)
  410 FORMAT(' EY')
      WRITE(5,12) (EY(I),I=1,N)
      WRITE(5,411)
  411 FORMAT(' VY')
      DO 412 I=1,N
  412 WRITE(5,12) (VY(I,J),J=1,I)
      WRITE(5,43)
   43 FORMAT(' CY')
      DO 44 I=1,N
   44 WRITE(5,12) (CY(I,J),J=1,I)
      WRITE(5,441)
  441 FORMAT(' A')
      DO 442 I=1,N
  442 WRITE(5,12) (A(I,J),J=1,M)
      WRITE (5,45)
   45 FORMAT(' P')
      WRITE(5,12) (P(I),I=1,M)
      WRITE(5,451)
  451 FORMAT(' EP')
      WRITE(5,12) (EP(I),I=1,M)
      WRITE(5,452)
  452 FORMAT(' VP')
      DO 453 I=1,M
  453 WRITE(5,12) (VP(I,J),J=1,I)
      WRITE(5,47)
   47 FORMAT(' CP')
      DO 48 I=1,M
   48 WRITE(5,12) (CP(I,J),J=1,I)
      WRITE(5,49)
   49 FORMAT(' CHI2,CHI2NM')
      WRITE(5,12) CHI2,CHI2NM
      STOP
      END
      SUBROUTINE MATINV(B,D,Q,E,NTEST,NS,NARA,NMAX)
      DIMENSION B(NARA),D(NARA,NARA),Q(NARA,NARA),E(NARA,NMAX)
      IP=NS+1
      BIG=0.0
      DO 555 I=1,NS
      DO 555 J=1,NS
      ABD=ABS(D(I,J))
      IF(ABD-BIG) 555,555,554
  554 BIG=ABD
  555 CONTINUE
      FACT=SQRT(BIG)
      I=1
    1 IF(I-NS) 2,2,20
    2 J=1
    3 IF(J-NS) 4,4,8
    4 K=1
    5 IF(K-NS) 6,6,7
    6 E(J,K)=D(K,J)/FACT
      K=K+1
      GO TO 5
    7 J=J+1
      GO TO 3
    8 L=1
    9 IF(L-NS) 10,10,14
   10 IF(L-I) 11,13,11
   11 E(L,IP)=0.0
   12 L=L+1
      GO TO 9
   13 E(L,IP)=1.0
      GO TO 12
   14 CALL JORDAN(B,E,NTEST,NS,NARA,NMAX)
      IF(NTEST) 15,15,16
   15 RETURN
   16 M=1
   17 IF(M-NS) 18,18,19
   18 Q(I,M)=E(M,IP)/FACT
      M=M+1
      GO TO 17
   19 I=I+1
      GO TO 1
   20 RETURN
      END
      SUBROUTINE JORDAN(B,C,INDEX,N,NARA,NMAX)
C
C     SUBROUTINE JORDAN SOLVES A SYSTEM OF LINEAR NONHOMOGENEOUS
C     EQUATIONS BY THE METHOD OF GAUSS-JORDAN REDUCTION. IF THE SYSTEM
C     FAILS TO HAVE A SOLUTION, A FLAG IS SET WHICH SIGNALS THE MAIN
C     PROGRAM.
C
      DIMENSION B(NARA),C(NARA,NMAX)
      K=1
    1 IF(K-N) 2,2,22
    2 IF(C(K,K)) 10,3,10
    3 L=K+1
    4 IF(L-N) 5,5,21
    5 IF(C(L,K)) 7,6,7
    6 L=L+1
      GO TO 4
    7 M=1
    8 IF(M-N-1) 9,9,2
    9 B(M)=C(K,M)
      C(K,M)=C(L,M)
      C(L,M)=B(M)
      M=M+1
      GO TO 8
   10 J=N+1
   11 IF(J-K) 13,12,12
   12 C(K,J)=C(K,J)/C(K,K)
      J=J-1
      GO TO 11
   13 I=1
   14 IF(I-N) 16,16,15
   15 K=K+1
      GO TO 1
   16 IF(I-K) 18,17,18
   17 I=I+1
      GO TO 14
   18 II=N+1
   19 IF(II-K) 17,20,20
   20 C(I,II)=C(I,II)-C(I,K)*C(K,II)
      II=II-1
      GO TO 19
   21 INDEX=0
      GO TO 23
   22 INDEX=1
   23 RETURN
      END