C     GLSMOD
C
C     IBM PC - D.L. SMITH  21 DECEMBER 1992
C
      DIMENSION Y(40),EY(40),VY(40,40),CY(40,40),Q(40,40),QVY(40,40),
     1QVYI(40,40),YA(40),P(10),EP(10),VP(10,10),CP(10,10),PA(10),EA(10),
     2VA(10,10),CA(10,10),A(40,10),QN(40),WN(40,41)
C
C     INITIALIZATION AND CONTROL
C
      WRITE(*,1)
    1 FORMAT(' GLSMOD'/)
      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)
      READ(4,12) (YA(I),I=1,N)
      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)
      READ(4,12) (PA(I),I=1,M)
      DO 19 I=1,M
   19 READ(4,12) (VA(I,J),J=1,I)
      DO 20 I=1,M
      DO 20 J=1,I
   20 VA(J,I)=VA(I,J)
      DO 30 I=1,N
   30 READ(4,12) (A(I,J),J=1,M)
C
C     GENERAL LEAST-SQUARES ANALYSIS
C
      DO 1040 I=1,N
 1040 EY(I)=SQRT(VY(I,I))
      DO 1050 I=1,N
      DO 1050 J=1,N
 1050 CY(I,J)=VY(I,J)/EY(I)/EY(J)
      DO 1051 I=1,M
 1051 EA(I)=SQRT(VA(I,I))
      DO 1052 I=1,M
      DO 1052 J=1,M
 1052 CA(I,J)=VA(I,J)/EA(I)/EA(J)
      DO 1053 I=1,N
      DO 1053 J=1,N
      Q(I,J)=0.0
      DO 1053 K2=1,M
      DO 1053 K1=1,M
 1053 Q(I,J)=Q(I,J)+A(I,K2)*VA(K2,K1)*A(J,K1)
      DO 1054 I=1,N
      DO 1054 J=1,N
 1054 QVY(I,J)=Q(I,J)+VY(I,J)
      CALL MATINV(QN,QVY,QVYI,WN,NTEST,N,40,41)
      IF(NTEST.EQ.1) GO TO 1090
      WRITE(*,1080)
 1080 FORMAT(' NO INV'/)
      STOP
 1090 CONTINUE
      DO 1100 I=1,M
      DO 1100 J=1,M
      VP(I,J)=VA(I,J)
      DO 1100 K4=1,M
      DO 1100 K3=1,N
      DO 1100 K2=1,N
      DO 1100 K1=1,M
 1100 VP(I,J)=VP(I,J)-VA(I,K4)*A(K3,K4)*QVYI(K3,K2)*A(K2,K1)*VA(K1,J)
      DO 1110 I=1,N
 1110 QN(I)=Y(I)-YA(I)
      DO 1200 I=1,M
      P(I)=PA(I)
      DO 1200 K3=1,M
      DO 1200 K2=1,N
      DO 1200 K1=1,N
 1200 P(I)=P(I)+VA(I,K3)*A(K2,K3)*QVYI(K2,K1)*QN(K1)
      DO 1300 I=1,M
 1300 EP(I)=SQRT(VP(I,I))
      DO 1400 I=1,M
      DO 1400 J=1,M
 1400 CP(I,J)=VP(I,J)/EP(I)/EP(J)
      CHI2=0.0
      DO 1600 K2=1,N
      DO 1600 K1=1,N
 1600 CHI2=CHI2+QN(K2)*QVYI(K2,K1)*QN(K1)
      CHI2NM=CHI2/FLOAT(N)
C
C     PRINT OUTPUT TO FILE (UNIT 5)
C
      WRITE(5,1700)
 1700 FORMAT(' YA')
      WRITE(5,12) (YA(I),I=1,N)
      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,4000)
 4000 FORMAT(' PA')
      WRITE(5,12) (PA(I),I=1,M)
      WRITE(5,4100)
 4100 FORMAT(' EA')
      WRITE(5,12) (EA(I),I=1,M)
      WRITE(5,4200)
 4200 FORMAT(' VA')
      DO 4300 I=1,M
 4300 WRITE(5,12) (VA(I,J),J=1,I)
      WRITE(5,4400)
 4400 FORMAT(' CA')
      DO 4500 I=1,M
 4500 WRITE(5,12) (CA(I,J),J=1,I)
      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