C profit.f
C The authors of this software are J D Carroll and Jih-Jie Chang.

C Copyright (c) 1993 by AT&T.
C Permission to use, copy, modify, and distribute this software for any
C purpose without fee is hereby granted, provided that this entire
C notice is included in all copies of any software which is or includes
C a copy or modification of this software and in all copies of the
C supporting documentation for such software.
C THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMP-
C LIED WARRANTY.  IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY
C REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
C OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.

C This software comes from the FIRST MDS Package of AT&T Bell Laboratories

C For explanation of the method this software implements, see
C Carroll,J.D. & Chang,J.J. (1964), "A general index of nonlinear
C correlation and its application to the problem of relating physical
C and psychological dimensions" (Abstract), in
C American Psychologist, 19, 540.
C----+----@----+----@----+----@----+----@----+----@----+----@----+----@

C
C     PROFIT - PERFORMS REGRESSION ON OPTIMIZATION OF LINEAR AND
C              NONLINEAR FIT.
C     2/69,  BY CARROLL AND CHANG
C
C     ADAPTED TO 360 SYSTEM BY VITHAL R. RAO
C
C     INPUT PARAMETERS
C     LANA - 1, DO LINEAR REGRESSION ONLY
C            2, DO NONLINEAR REGRESSION ONLY
C            3, DO BOTH LINEAR AND NONLINEAR REGRESSION
C     N - NUMBER OF STIMULI  (MAXIMUM=150)
C     K - NUMBER OF DIMENSIONS  (MAXIMUM=10)
C     M - NUMBER OF PROPERTIES  (MAXIMUM=20)
C     IRX - 0, X IS PUNCHED N BY K
C           1, X IS PUNCHED K BY N
C     IWGT - SEE WRITE UP
C     IPLOT - 0, PLOT PROJECTIONS OF PROPERTIES ONLY
C             1, PLOT PROJECTIONS OF PROPERTIES AND FUNCTION PLOTS
C             2, DO ALL PLOTS (SEE WRITE UP)
C     BCO - SEE WRITE UP
C     X - COORDINATES OF STIMULUS POINTS
C     P - PROPERTIES
C
C    CARD SETUP
C
C     CARD 1   INPUT PARAMETERS
C      COLS. 1-4 LANA
C             1-4    LANA
C             5-8    N
C             9-12   K
C            13-16   M
C            17-20   IRX
C            21-24   IWGT
C            25-28   IPLOT
C            29-35   BCO
C     CARD 2 . TITLE FOR THE RUN (20A4)
C     CARD 3. VARIABLE FORMAT FOR READING THE STIMULUS CONFIGURATION AS
C             PER THE OPTION CHOSEN FOR IRX (20A4)
C     CARD(S) 3. DATA CARDS OF STIMULUS CONFIGURATION
C     CARD 4. VARIABLE FORMAT FOR READING THE PROPERTY MATRIX
C     CARD(S) 5. DATA CARDS OF PROPERTY MATRIX
C
C     NAME OF PROPERTY ON A CARD AND ONE OR MORE CARDS FOR SCALE VALUES
C     AS PER FORMAT FOR EACH PROPERTY.
C
C     REPEAT THE SET OF CARDS 1-5 FOR SUBSEQUENT PROBLEMS IF ANY. IF NO
C     MORE PROBLEMS PLACE A BLANK CARD.
C
C
C
C
C
C
C
C
C
C
C
C
      REAL*8 ZU(10,10)
      DIMENSION X(150,10),P(150,20),PNAME(20,20),T(10,20)
      DIMENSION XMEAN(10),SINV(10),RR(10),BS(10,10),H(150),ROOT(10)
      DIMENSION NUMB(150),PNUM(20)
      DIMENSION A(22500),FMT(20),PX(170),PY(170)
      DIMENSION XMO(4),OLANA(2),TLANA(2)
      DIMENSION TIT(20),FMPS(5),FMXS(5),FORA(3),FORS(3),FFM(2)
C     EQUIVALENCE (A,PX),(A(21),PY)
      COMMON N,K,M,XN,CDEN,SUMML,ASIN,TANG,IWGT,BCO
      COMMON X,P,PNAME,T,XMEAN,SINV,RR,BS,H,ROOT
      COMMON XSCALE,PSCALE,NUMB,PNUM
      COMMON /FMIC/TIT,FMPS,FMXS,FORA,FORS,FFM
C     DATA FORV(1)/33H(28HPROJECTIONS ON FITTED VECTOR)/
C     DATA OLANA(1)/12H LINEAR-    /
C     DATA TLANA(1)/12H NONLINEAR- /
C     DATA FORH(1)/34H(27HORIGINAL VALUES ON PROPERTYI3)/
      NL=150
 5    READ102,LANA,N,K,M,IRX,      IWGT,IPLOT,BCO
      WRITE(6,121)
      IF(LANA.EQ.0) GOTO99
 95   READ101,(TIT(I),I=1,20)
C     READ IN CONFIGURATION
      CALL READX(X,N,K,IRX,NL)
C     READ IN PROPERTY NAMES AND PROPERTIES
      READ101,(FMT(I),I=1,20)
      DO4J=1,M
      READ101,(PNAME(JW,J),JW=1,20)
 4    READFMT,(P(I,J),I=1,N)
C     OPTION TO LABEL STIMULI ON MICROFILM
C     IF(MIDEN.NE.0)GOTO2
C     CALL LMIC(N,NUMB)
C     GOTO3
C2    READ101,(FMT(I),I=1,20)
C     READFMT,(NUMB(I),I=1,N)
C     GENERATE ALPHABET TO LABEL VECTORS ON MICROFILM
C3    CALL LET(PNUM,M)
    3 CONTINUE
      XN=N
      XN1=XN-1.
      DND=XN*XN1
      CONS=2.*BCO/DND
C     NORMALIZE X
      DO10J=1,K
      SUM=0.
      DO11I=1,N
 11   SUM=SUM+X(I,J)
      XMEAN(J)=SUM/XN
      DO12I=1,N
 12   X(I,J)=X(I,J)-XMEAN(J)
 10   CONTINUE
 15   PRINT104,(TIT(I),I=1,20)
      IF(LANA.NE.2)GOTO21
      CALL XSTAR(A,A(101),NL)
      PRINT105
      GOTO22
   21 CONTINUE
      PRINT106
 22   PRINT107
      DO13I=1,K
 13   PRINT108,I,(X(J,I),J=1,N)
      IF(LANA.EQ.2)GOTO23
C     LINEAR REGRESSION
      CALL TRANS(A,A(101),A(302),NL)
      GOTO25
 23   CALL CSD(X,SINV,N,K,XN,NL)
 25   DO50JM=1,M
      PRINT109,JM
      IF(LANA.EQ.2)GOTO26
C     COMPUTE DIRECTION COSINES FOR L.R.
      CALL DCOSL(A(302),JM,NL,XL)
      GOTO27
 26   SUM=0.
      DO24I=2,N
      II=I-1
      DO24J=1,II
 24   SUM=SUM+(P(I,JM)-P(J,JM))**2
      CDEN=SUM*CONS
      CALL MOMENT(A,P(1,JM),XMO,NL)
C     COMPUTE DIRECTION COSINES FOR N.L.R.
      CALL DCOSN(A,A(101),JM,XL,XMO)
 27   DO35J=1,K
 35   T(J,JM)=T(J,JM)/XL
C     COMPUTE PROJECTIONS OF STIMULI ON FITTED VECTORS-H
      DO28I=1,N
      H(I)=0.
      DO28J=1,K
 28   H(I)=H(I)+X(I,J)*T(J,JM)
      IF(LANA.NE.2) CALL CORR(P(1,JM),H,ROOT(JM),N,XN)
C     PRINT ORIGINAL (P) AND OBTAINED (H) PROPERTY
      PRINT110,JM
      PRINT111,(P(J,JM),J=1,N)
      PRINT120
      PRINT111,(H(J),J=1,N)
      IF(IPLOT.EQ.0) GOTO50
C     PLOT OBTAINED VERSUS GIVEN PROPERTY
      WRITE(6,170) JM
  170 FORMAT('1'10X,'PLOT OF ORIGINAL (X-AXIS) VERSUS OBTAINED (Y-AXIS)
     XFOR PROPERTY VECTOR NO.',I3/)
      IP=1
      DO 171 J=1,N
      PX(J) = P(J,JM)
  171 PY(J) = H(J)
      CALL PLOTX(PX,PY,N,IP)
      CALL CORREL(PX,PY,N,JM)
C50   CALL PH(P(1,JM),H,FORH,FORV,TIT,N,JM)
   50 CONTINUE
C     PRINT TABLE 1
 51   IF(LANA.NE.2)GOTO29
      PRINT112
      PRINT113
      GOTO32
 29   PRINT114
      PRINT 115
 32   DO30I=1,M
 30   PRINT116,I,ROOT(I),(PNAME(J,I),J=1,20)
C     PRINT TABLE 2 MATRIX OF DIRECTION COSINE
      PRINT117
      PRINT118,(I,I=1,K)
      DO31I=1,M
 31   PRINT119,I,(T(J,I),J=1,K)
C     PRINT TABLE 3 COSINE OF ANGLES BETWEEN VECTORS
C     COMPUTE COSINE OF ANGLES BETWEEN VECTORS AND PLOT PAIRS OF VECTORS
      CALL PVECT(B,T,M,K)
C     CALL PVECT(A,IPLOT)
      IF(LANA.EQ.2) GOTO33
      DO34J=1,K
      DO34I=1,N
 34   X(I,J)=X(I,J)+XMEAN(J)
      GOTO36
 33   CALL RSTORX(A,A(201))
C     PLOT PROJECTION OF VECTORS IN ORIGINAL SPACE
C     IF(IPLOT.EQ.0) GOTO51
C36   KK=K-1
C     DO40I=1,KK
C     JJ=I+1
C     DO40J=JJ,K
C     DO41JM=1,M
C     PX(JM)=T(I,JM)
C41   PY(JM)=T(J,JM)
C     CALL MAXI(X(1,I),X(1,J),XMAX,N)
C     XSCALE=460./XMAX
C     CALL MAXI(PX,PY,PMAX,M)
C     PSCALE=460./PMAX
C40   CALL PJPRO(PX,PY,X(1,I),X(1,J),I,J)
   36 CONTINUE
      NM=N+M
      WRITE(6,172)
  172 FORMAT('1',10X,'PLOT FOR FIRST TWO DIMENSIONS OF STIMULUS POINTS
     XAND DIRECTION COSINES OF FITTED PROPERT VECTORS'/)
      DO 173 I=1,N
      PX(I) = X(I,1)
  173 PY(I) = X(I,2)
      IN=N+1
      DO 174 I=IN,NM
      II=I-N
      PX(I) = T(1,II)
  174 PY(I) = T(2,II)
      IP=0
      CALL PLOTX(PX,PY,NM,IP)
   40 CONTINUE
      IF(LANA.LT.3)GOTO5
      LANA=2
      GOTO15
   99 CALL EXIT
  101 FORMAT(20A4)
  102 FORMAT(7I4,F7.3)
  121 FORMAT('1',10X,'PROFIT PROGRAM BY CHANG AND CARROLL'/10X,'ADAPTED
     XTO IBM 360 SYSTEM BY VITHAL R. RAO'//)
  104 FORMAT(//10X,20A4)
 105  FORMAT(21H0NONLINEAR REGRESSION)
 106  FORMAT(18H0LINEAR REGRESSION)
 107  FORMAT(25H0NORMALIZED CONFIGURATION)
 108  FORMAT(I4,3X,10E12.4/(7X,10E12.4))
  109 FORMAT('1',10X,'PROPERTY',I3/)
 110  FORMAT(31H0   ORIGINAL VALUES ON PROPERTYI3)
 111  FORMAT(6X,10E12.4)
 112  FORMAT(45H1TABLE 1.  THE SMALLEST EIGEN ROOT ASSOCIATED/11X,18HWIT
     1H EACH PROPERTY)
 113  FORMAT(1H0,12X,4HROOT,8X,8HPROPERTY)
 114  FORMAT(55H1TABLE 1.  THE MAXIMUM CORRELATION BETWEEN THE PROPERTY/
     111X,36HAND THE PROJECTIONS ON FITTED VECTOR)
 115  FORMAT(1H0,12X,3HRHO,9X,8HPROPERTY)
 116  FORMAT(I5,3X,E12.4,3X,20A4)
 117  FORMAT(46H0TABLE 2.  DIRECTION COSINES OF FITTED VECTORS/11X,20H I
     1N NORMALIZED SPACE//21X,9HDIMENSION)
 118  FORMAT(10H0   VECTORI5,9I9)
 119  FORMAT(I7,3X,10F9.4)
 120  FORMAT(33H0   PROJECTIONS ON FITTED VECTORS)
      STOP
      END
      SUBROUTINE DCOSL(C,JM,NL,XL)
      DIMENSION C(10,NL)
      DIMENSION X(150,10),P(150,20),PNAME(20,20),T(10,20)
      DIMENSION XMEAN(10),SINV(10),RR(10),BS(10,10),H(150),ROOT(10)
      COMMON N,K,M,XN,CDEN,SUMML,ASIN,TANG,IWGT,BCO
      COMMON X,P,PNAME,T,XMEAN,SINV,RR,BS,H,ROOT
      SSQ=0.
      DO24I=1,K
      T(I,JM)=0.
      DO23II=1,N
 23   T(I,JM)=T(I,JM)+C(I,II)*P(II,JM)
 24   SSQ=SSQ+T(I,JM)**2
      XL=SQRT(SSQ)
      DO5I=1,K
 5    PRINT100,(C(I,J),J=1,N)
 100  FORMAT(2X,10E12.4)
      PRINT100,(P(I,JM),I=1,N)
      PRINT100,(T(I,JM),I=1,K)
      PRINT101,SSQ,XL
 101  FORMAT(5H0SSQ=E13.5,3X,3HXL=E13.5)
      RETURN
      END
      SUBROUTINE DCOSN(A,B,JM,XL,XMO)
C     COMPUTE DIRECTION COSINS
      DIMENSION X(150,10),P(150,20),PNAME(20,20),T(10,20)
      DIMENSION A(10,10),B(10,10),XMO(4)
      DIMENSION XMEAN(10),SINV(10),RR(10),BS(10,10),H(150),ROOT(10)
      COMMON N,K,M,XN,CDEN,SUMML,ASIN,TANG,IWGT,BCO
      COMMON X,P,PNAME,T,XMEAN,SINV,RR,BS,H,ROOT
      XN1=XN-1.
      IF(IWGT.LT.2)GOTO5
      DO6I=1,N
      H(I)=1.
      DO6J=1,N
      IF(I.EQ.J)GOTO6
      PD=ABS(P(I,JM)-P(J,JM))
      IF(PD.GT.BCO)GOTO6
      H(I)=H(I)+1.
 6    CONTINUE
 5    DO25I=1,K
      DO25J=1,I
      Q=0.
      DO20IK=2,N
      KK=IK-1
      DO20L=1,KK
      PD=ABS(P(IK,JM)-P(L,JM))
      XX=(X(IK,I)-X(L,I))*(X(IK,J)-X(L,J))
      IF(IWGT.EQ.0)GOTO21
      IF(PD.GT.BCO)GOTO20
      GOTO(3,4),IWGT
 3    PP=1.
      GOTO22
 4    PP=1./SQRT(H(IK)*H(L))
      GOTO22
 21   PP=1./(PD**2+CDEN)
 22   Q=Q+XX*PP
 20   CONTINUE
      A(I,J)=SINV(I)*SINV(J)*Q
 25   A(J,I)=A(I,J)
C     FIND ROOTS AND VECTORS OF A
      MX=1
      CALL MTEVV(A,10,10,B,10,10,K,MX)
      ROOT(JM)=A(K,K)
C     COMPUTE ZSQ
      Z=(XN1*ROOT(JM)-XMO(1))/SQRT(XMO(2))
      ZSQ=Z**2
      PRINT102
 102  FORMAT(50H0  THE STANDARDIZED INDEX OF NONLINEAR CORRELATION)
      PRINT101,Z,ZSQ
 101  FORMAT(5H0  Z=E14.5,5X,4HZSQ=E14.5)
      SSQ=0.
      DO26I=1,K
      T(I,JM)=B(I,K)*SINV(I)
 26   SSQ=SSQ+T(I,JM)**2
      XL=SQRT(SSQ)
      RETURN
      END
      SUBROUTINE TRANS(A,V,C,NL)
      REAL*8 ZU(10,10)
      DIMENSION X(150,10),P(150,20),PNAME(20,20),T(10,20)
      DIMENSION XMEAN(10),SINV(10),RR(10),BS(10,10),H(150),ROOT(10)
      DIMENSION A(10,10),V(201),C(10,NL)
      COMMON N,K,M,XN,CDEN,SUMML,ASIN,TANG,IWGT,BCO
      COMMON X,P,PNAME,T,XMEAN,SINV,RR,BS,H,ROOT
      DO5I=1,K
      DO5J=1,K
      A(I,J)=0.
      DO5L=1,N
 5    A(I,J)=A(I,J)+X(L,J)*X(L,I)
      PRINT100
      DO6I=1,K
 6    PRINT101,I,(A(I,J),J=1,K)
C     COMPUTE (XT)X INVERSE BY MEANS OF ROOTS AND VECTORS,STORE IN A.
      MX=1
      CALL MTINV(A,10,10,ZU,10,10,K,MX)
      IF(MX.EQ.1) GOTO8
      PRINT103
      STOP
C     COMPUTE X*(XTX) INVERSE AND STORE IN C.
 8    DO12I=1,K
      DO12J=1,N
      C(I,J)=0.
      DO12L=1,K
 12   C(I,J)=C(I,J)+X(J,L)*A(L,I)
      PRINT102
      DO13I=1,K
 13   PRINT101,I,(C(I,J),J=1,N)
 100  FORMAT(18H0COVARIANCE MATRIX)
 101  FORMAT(I4,3X,10E12.4/(7X,10E12.4))
 102  FORMAT(16H0X*(XTX) INVERSE)
 103  FORMAT(4H0MX=I1,3X,20HA CANNOT BE INVERTED)
      RETURN
      END
      SUBROUTINE READX(X,N,K,IRX,NL)
      DIMENSION X(NL,K),FMT(20)
      READ100,(FMT(I),I=1,20)
  100 FORMAT(20A4)
      IF(IRX.EQ.0) GOTO10
      DO6I=1,K
 6    READFMT,(X(J,I),J=1,N)
      GOTO20
 10   DO8I=1,N
 8    READFMT,(X(I,J),J=1,K)
 20   RETURN
      END
      SUBROUTINE CORR(P,H,RHO,N,XN)
      DIMENSION P(1),H(1)
C     COMPUTE CORRELATION BETWEEN GIVEN AND OBTAINED VALUES ON ONE
C     PROPERTY
      SUMP=0.
      SUMPP=0.
      SQP=0.
      SQXL=0.
      DO5I=1,N
 5    SUMP=SUMP+P(I)
      PMEAN=SUMP/XN
      DO6I=1,N
      PDV=P(I)-PMEAN
      SQP=SQP+PDV**2
        SUMPP=SUMPP+PDV*H(I)
 6    SQXL=SQXL+H(I)**2
      RHO=SUMPP/SQRT(SQP*SQXL)
      RETURN
      END
      SUBROUTINE XSTAR(A,B,NL)
      DIMENSION X(150,10),P(150,20),PNAME(20,20),T(10,20)
      DIMENSION XMEAN(10),SINV(10),RR(10),BS(10,10),H(150),ROOT(10)
      DIMENSION A(10,10),B(10,10)
      COMMON N,K,M,XN,CDEN,SUMML,ASIN,TANG,IWGT,BCO
      COMMON X,P,PNAME,T,XMEAN,SINV,RR,BS,H,ROOT
C     ROTATE CONFIGURATION TO PRINCIPLE AXES WITH EQUAL VARIANCE
      DO10I=1,K
      DO10J=1,K
      A(I,J)=0.
      DO10JJ=1,N
 10   A(I,J)=A(I,J)+X(JJ,I)*X(JJ,J)
      MX=1
      CALL MTEVV(A,10,10,B,10,10,K,MX)
      DO15I=1,K
      RR(I)=A(I,I)
      DO15J=1,K
 15   BS(J,I)=B(J,I)/SQRT(RR(I))
      DO20I=1,N
      DO21J=1,K
      B(J,1)=0.
      DO21JJ=1,K
 21   B(J,1)=B(J,1)+X(I,JJ)*BS(JJ,J)
      DO22J1=1,K
 22   X(I,J1)=B(J1,1)
 20   CONTINUE
      RETURN
      END
      SUBROUTINE CSD(X,SINV,N,K,XN,NL)
      DIMENSION X(NL,K),SINV(10)
      DO10J=1,K
      SUMV=0.
      SUM=0.
      DO12I=1,N
 12   SUM=SUM+X(I,J)
      CMEAN=SUM/XN
      DO15I=1,N
 15   SUMV=SUMV+(X(I,J)-CMEAN)**2
 10   SINV(J)=1./SQRT(SUMV)
      RETURN
      END
      SUBROUTINE RSTORX(A,B)
      DIMENSION X(150,10),P(150,20),PNAME(20,20),T(10,20)
      DIMENSION XMEAN(10),SINV(10),RR(10),BS(10,10),H(150),ROOT(10)
      DIMENSION A(10,10),B(10,20)
      COMMON N,K,M,XN,CDEN,SUMML,ASIN,TANG,IWGT,BCO
      COMMON X,P,PNAME,T,XMEAN,SINV,RR,BS,H,ROOT
      DO5JM=1,M
      SSQ=0.
      DO6I=1,K
      B(I,JM)=0.
      DO7JJ=1,K
 7    B(I,JM)=B(I,JM)+BS(I,JJ)*T(JJ,JM)
 6    SSQ=SSQ+B(I,JM)**2
      DEN=SQRT(SSQ)
      DO8I=1,K
 8    T(I,JM)=B(I,JM)/DEN
 5    CONTINUE
      DO10I=1,K
      DO10J=1,K
 10   A(I,J)=BS(J,I)*RR(I)
      DO11I=1,N
      DO12J=1,K
      SINV(J)=0.
      DO12JJ=1,K
 12   SINV(J)=SINV(J)+X(I,JJ)*A(JJ,J)
      DO13JB=1,K
 13   X(I,JB)=SINV(JB)+XMEAN(JB)
 11   CONTINUE
C     PRINT DIRECTION COSINE WRT ORIGINAL SPACE
      PRINT101
      PRINT102,(I,I=1,K)
      DO21I=1,M
 21   PRINT103,I,(T(J,I),J=1,K)
 101  FORMAT(46H0TABLE 4.  DIRECTION COSINES OF FITTED VECTORS/11X,18H I
     1N ORIGINAL SPACE//21X,9HDIMENSION)
 102  FORMAT(11H0   VECTORSI5,9I9)
 103  FORMAT(I7,3X,10F9.4)
      RETURN
      END
      SUBROUTINE ZEROS(A,N)
      DIMENSION A(1)
      K=N
 10   DO5II=1,K
 5    A(II)=0.
      GOTO50
      ENTRY ZEROD(A,I,J)
      K=I*J
      GOTO10
      ENTRY ZEROT(A,I,J,M)
      K=I*J*M
      GOTO10
 50   RETURN
      END
      SUBROUTINE MOMENT(A,P,XMO,NL)
      DIMENSION A(NL,NL),P(NL),SBU(4),ALFA(3),XMO(4)
      COMMON N,K,M,XN,CDEN,SUMML,ASIN,TANG,IWGT,BCO
C     COMPUTE BETA MU UP TO THE FOURTH POWER
C     COMPUTE A
      DO3I=1,N
      A(I,I)=0.
      DO3J=1,N
      IF(I-J)2,3,4
 2    A(I,J)=-1./((P(I)-P(J))**2+CDEN)
      A(J,I)=A(I,J)
 4    A(I,I)=A(I,I)-A(I,J)
 3    CONTINUE
      DO1I=1,4
 1    SBU(I)=0.
      DO5I1=1,N
      SBU(1)=SBU(1)+A(I1,I1)
      DO5I2=1,N
      SBU(2)=SBU(2)+A(I1,I2)**2
      DO5I3=1,N
      AA=A(I1,I2)*A(I2,I3)
      SBU(3)=SBU(3)+AA*A(I3,I1)
      DO5I4=1,N
 5    SBU(4)=SBU(4)+AA*A(I3,I4)*A(I4,I1)
C     COMPUTE MOMENTS
      XMO(1)=SBU(1)
      XN1=XN-1.
      DO16I=1,4
 16   SBU(I)=SBU(I)/XN1
      BSQ=SBU(1)**2
      BTH=BSQ*SBU(1)
      BFO=BSQ**2
      SBSM2=XN1*(SBU(2)-BSQ)
      SBSM3=XN1*(SBU(3)-3.*SBU(1)*SBU(2)+2.*BTH)
      SBSM4=XN1*(SBU(4)-4.*SBU(1)*SBU(3)+6.*BSQ*SBU(2)-3.*BFO)
      XNN1=N+1
      XNIS=XN1**2
      XN3=N+3
      XN4=N+4
      Q=XNN1*XN3
      XMO(2)=(2.*XN1/XNN1)*SBSM2
      XMO(3)=(4.*XNIS/Q)*SBSM3
      XMO(4)=(24.*XNIS*XN1/(Q*XN4))*(2.*SBSM4+SBSM2**2)
      SKEW=XMO(3)/(XMO(2)**1.5)
      CURT=XMO(4)/(3.*XMO(2)**2)
      PRINT101,(J,J=1,4),(XMO(I),I=1,4)
 101  FORMAT(10H0  MOMENTS/1H I13,3I15/4X,4E15.5)
      PRINT102,SKEW,CURT
 102  FORMAT(12H0  SKEWNESS=E14.5,5X,9HKURTOSIS=E14.5)
      RETURN
      END
CMTINV             MTINV
      SUBROUTINE MTINV (A,I1,I2,B,I3,I4,N,M)
      DIMENSION A(I1,I2), B(I3,I4)
C
C     THIS SUBROUTINE CALCULATES THE INVERSE OF MATRIX A BY THE GAUSS-
C     JORDAN ELIMINATION SCHEME.  ALL CALCULATIONS ARE DONE IN DOUBLE-
C     PRECISION ARITHMETIC
C
      DOUBLE PRECISION B,BMAX,BMULT
      EQUIVALENCE (BMAX,BMULT)
C
C     SINGLE TO DOUBLE PRECISION TRANSFER
C
  101 DO 103 I = 1,N
  102 DO 103 J = 1,N
  103 B(I,J) = A(I,J)
C
C     FIND THE LARGEST ELEMENT IN THE N-K BY N-K LOWER RIGHT SUBMATRIX.
C
  200 DO 430 K = 1,N
      A(K,1) = FLOAT(K)
      A(K,2) = A(K,1)
  203 BMAX = DABS(B(K,K))
  210 DO 219 I = K,N
  211 DO 219 J = K,N
      IF (BMAX - DABS(B(I,J)) )  213,219,219
  213 BMAX = DABS ( B(I,J))
  214 A(K,1) = FLOAT(I)
  215 A(K,2) = FLOAT(J)
  219 CONTINUE
  216 IF (BMAX  -   1.D-38)   801,301,301
C
C     EXCHANGE ROWS AND COLUMNS TO PUT B(I,J) ON DIAGONAL.
C
  301 I = IFIX(A(K,1))
  302 J = IFIX(A(K,2))
  313 DO 314 M = 1,N
      BMULT = B(I,M)
      B(I,M) = B(K,M)
  314 B(K,M) = BMULT
  303 DO 304 M = 1,N
      BMULT = B(M,J)
      B(M,J) = B(M,K)
  304 B(M,K) = BMULT
C
C     ACTUAL INVERSION STEP.
C     BEGIN ROW ITERATION
C
  401 DO 425 I = 1,N
C                                      IGNORE ROW K.
  402 IF (I  -   K) 405,425,405
C                                      SET ROW MULTIPLIER.
 405  BMULT = B(I,K) / B(K,K)
C                                      MODIFY ROW ELEMENTS.
  410 DO 420 J = 1,N
C                                      IGNORE COLUMN K
  411 IF (J  -   K)  414,412,414
  412 B(I,J) = -BMULT
  413 GO TO 420
  414 B(I,J) = B(I,J) - B(K,J) * BMULT
  420 CONTINUE
  425 CONTINUE
C
C     DIVIDE PIVOT ROW BY PIVOT ELEMENT
C
  426 BMULT = 1.D0/ B(K,K)
  427 B(K,K) = 1.D0
C                                      ACTUAL DIVISION
  428 DO 429 J = 1,N
  429 B(K,J) = B(K,J) * BMULT
  430 CONTINUE
C
C     REARRANGE AND REASSEMBLE INVERSE MATRIX.
C
  501 DO 512 IJK = 1,N
  502 I = N-IJK+1
  503 L = IFIX(A(I,2))
  504 J = IFIX(A(I,1))
  505 DO 508 M = 1,N
  506 BMULT = B(I,M)
  507 B(I,M) = B(L,M)
  508 B(L,M) = BMULT
  509 DO 512 M = 1,N
  510 BMULT = B(M,I)
  511 B(M,I) = B(M,J)
  512 B(M,J) = BMULT
C
C     MOVE A INVERSE BACK INTO A
C
  601 DO 603 I = 1,N
  602 DO 603 J = 1,N
  603 A(I,J) = B(I,J)
  701 M = 1
  702 RETURN
C
C     ERROR RETURN AT SOME TIME BMAX WAS .LT. 0.1D-38
C
  801 M = 2
  802 RETURN
  999 CONTINUE
      END
      SUBROUTINE MTEVV(X,I1,I2,D,I3,I4,N,M)
      DIMENSION A(55),X(I1,I2),D(I1,I2),E(100) ,NS(10),MS(10)
      K=1
      DO 1 J=1,N
      DO 1 I=1,J
      A(K) = X(I,J)
    1 K=K+1
      CALL EIGEN(A,E,N,0)
      II=1
      DO 974 I=1,N
      DO 974 J=1,N
      D(J,I) = E(II)
      II=II+1
  974 CONTINUE
      II=0
      DO 973 J=1,N
      II=II+J
  973 X(J,J) = A(II)
      RETURN
      END
      SUBROUTINE MAXI(X,Y,Z,N)
C     MAXI FINDS THE MAGNITUDE OF ARRAYS X AND Y EACH OF LENGTH N.
C     ON RETURN ABSOLUTE MAXIMUM IS STORED IN Z.
      DIMENSION X(1),Y(1)
      Z=ABS(X(1))
      DO10I=1,N
      AXO=ABS(X(I))
      IF(AXO.GT.Z)  Z=AXO
      AYO=ABS(Y(I))
      IF(AYO.GT.Z)  Z=AYO
 10   CONTINUE
      RETURN
      END
      SUBROUTINE MM(X,N,XMAX,XMIN)
      DIMENSION X(100)
      XMAX=X(1)
      XMIN=XMAX
      DO20I=2,N
      IF(X(I)-XMAX)15,20,10
 10   XMAX=X(I)
      GOTO20
 15   IF(X(I)-XMIN)16,20,20
 16   XMIN=X(I)
 20   CONTINUE
      RETURN
      END
      SUBROUTINE PVECT(B,T,M,K)
      DIMENSION B(20,20),T(10,20)
C     SUBROUTINE PVECT(B,IPLOT)
C     DIMENSION X(150,10),P(150,20),PNAME(12,20),T(10,20)
C     DIMENSION XMEAN(10),SINV(10),RR(10),BS(10,10),H(150),ROOT(10)
C     DIMENSION B(20,20)
C     COMMON N,K,M,XN,CDEN,SUMML,ASIN,TANG,IWGT,BCO
C     COMMON X,P,PNAME,T,XMEAN,SINV,RR,BS,H,ROOT
      DO5I=2,M
      JJ=I-1
      DO5J=1,JJ
      SUMML=0.
      DO10KK=1,K
 10   SUMML=SUMML+T(KK,I)*T(KK,J)
      ASIN=1.-SUMML**2
      ASIN=SQRT(ASIN)
      B(I,J)=SUMML
C     TANG=ASIN/SUMML
      IF(IPLOT.LT.2) GOTO5
C     CALL PRJECT(T(1,J),T(1,I),PNAME(1,J),PNAME(1,I),X)
 5    CONTINUE
C     PRINT TABLE 3
      PRINT101
      PRINT102,(I,I=1,M)
      DO20I=2,M
      JJ=I-1
 20   PRINT103,I,(B(I,J),J=1,JJ)
 101  FORMAT(43H0TABLE 3.  COSINE OF ANGLES BETWEEN VECTORS)
 102  FORMAT(10H0   VECTORI5,9I9/6X,10I9)
 103  FORMAT(I7,3X,10F9.4/(10X,10F9.4))
      RETURN
      END
      SUBROUTINE EIGEN(A,R,N,MV)
C
C        SUBROUTINE EIGEN
C
C        PURPOSE
C           COMPUTE EIGENVALUES AND EIGENVECTORS OF A REAL SYMMETRIC
C           MATRIX
C
C        USAGE
C           CALL EIGEN(A,R,N,MV)
C
C        DESCRIPTION OF PARAMETERS
C           A - ORIGINAL MATRIX (SYMMETRIC), DESTROYED IN COMPUTATION.
C               RESULTANT EIGENVALUES ARE DEVELOPED IN DIAGONAL OF
C               MATRIX A IN DECENDING ORDER.
C           R - RESULTANT MATRIX OF EIGENVECTORS (STORED COLUMNWISE,
C               IN SAME SEQUENCE AS EIGENVALUES)
C           N - ORDER OF MATRICES A AND R
C           MV- INPUT CODE
C                       SEQUENCE)
C
C        REMARKS
C           ORIGINAL MATRIX A MUST BE REAL SYMMETRIC (STORAGE MODE=1)
C           MATRIX A CANNOT BE IN THE SAME LOCATION AS MATRIX R
C
C        SUBROUTINES USED BY THIS SUBROUTINE
C           NONE
C
C        METHOD
C           DIAGONALIZATION METHOD ORIGINATED BY JACOBI AND ADAPTED
C           BY VON NEUMANN FOR LARGE COMPUTERS AS FOUND IN MATHEMATICAL
C           METHODS FOR DIGITAL COMPUTERS, EDITED BY A. RALSTON AND
C           H.S. WILF, JOHN WILEY AND SONS, NEW YORK, 1962, CHAPTER 7
C
C     ..................................................................
C
C
      DIMENSION A(1),R(1)
C        GENERATE IDENTITY MATRIX
C
      IF(MV-1) 22,14,22
   22 IQ=-N
      DO 1 J=1,N
      IQ=IQ+N
      DO 1 I=1,N
      IJ=IQ+I
      R(IJ)=0.0
      IF(I-J) 1,31,1
   31 R(IJ)=1.0
    1 CONTINUE
C
C        COMPUTE INITIAL AND FINAL NORMS (ANORM AND ANORMX)
C
   14 ANORM=0.0
      DO 2 I=1,N
      DO 2 J=I,N
      IF(I-J) 15,2,15
   15 IA=I+(J*J-J)/2
      ANORM=ANORM+A(IA)*A(IA)
    2 CONTINUE
      IF(ANORM) 39, 39, 50
   50 ANORM=2.0*SQRT(ANORM)
      ANRMX =ANORM*1.0E-8
C
C        INITIALIZE INDICATORS AND COMPUTE THRESHOLD, THR
C
      IND=0
      THR=ANORM
    3 THR=THR/FLOAT(N)
   13 L=1
    4 M=L+1
C
C        COMPUTE SIN AND COS
C
    5 MQ=(M*M-M)/2
      LQ=(L*L-L)/2
      LM=L+MQ
      IF(ABS(A(LM))-THR) 7,32,32
   32 IND=1
      LL=L+LQ
      MM=M+MQ
      X=0.5*(A(LL)-A(MM))
      Y=-A(LM)/SQRT(A(LM)*A(LM)+X*X)
      IF(X) 33,34,34
   33 Y=-Y
   34 SINX=Y/SQRT(2.0*(1.0+(SQRT(1.0-Y*Y))))
      SINX2=SINX*SINX
      COSX=SQRT(1.0-SINX2)
      COSX2=COSX*COSX
      SINCS =SINX*COSX
C
C        ROTATE L AND M COLUMNS
C
      ILQ=N*(L-1)
      IMQ=N*(M-1)
      DO 6 I=1,N
      IQ=(I*I-I)/2
      IF(I-L) 16,11,16
   16 IF(I-M) 17,11,18
   17 IM=I+MQ
      GO TO 19
   18 IM=M+IQ
   19 IF(I-L) 37,23,23
   37 IL=I+LQ
      GO TO 24
   23 IL=L+IQ
   24 X=A(IL)*COSX-A(IM)*SINX
      A(IM)=A(IL)*SINX+A(IM)*COSX
      A(IL)=X
   11 IF(MV-1) 26,6,26
   26 ILR=ILQ+I
      IMR=IMQ+I
      X=R(ILR)*COSX-R(IMR)*SINX
      R(IMR)=R(ILR)*SINX+R(IMR)*COSX
      R(ILR)=X
    6 CONTINUE
      X=2.0*A(LM)*SINCS
      Y=A(LL)*COSX2+A(MM)*SINX2-X
      X=A(LL)*SINX2+A(MM)*COSX2+X
      A(LM)=(A(LL)-A(MM))*SINCS+A(LM)*(COSX2-SINX2)
      A(LL)=Y
      A(MM)=X
C
C        TESTS FOR COMPLETION
C
C        TEST FOR M = LAST COLUMN
    7 IF(M-N) 35,8,35
   35 M=M+1
      GO TO 5
C        TEST FOR L = SECOND FROM LAST COLUMN
    8 IF(L-(N-1)) 36,9,36
   36 L=L+1
      GO TO 4
    9 IF(IND-1) 10,38,10
   38 IND=0
      GO TO 13
C        COMPARE THRESHOLD WITH FINAL NORM
   10 IF(THR-ANRMX ) 39,39,3
C
C        SORT EIGENVALUES AND EIGENVECTORS
C
   39 IQ=-N
      DO 30 I=1,N
      IQ=IQ+N
      LL=I+(I*I-I)/2
      JQ=N*(I-2)
      DO 30 J=I,N
      JQ=JQ+N
      MM=J+(J*J-J)/2
      IF(A(LL)-A(MM)) 40,30,30
   40 X=A(LL)
      A(LL)=A(MM)
      A(MM)=X
      IF(MV-1) 28,30,28
   28 DO 20 K=1,N
      ILR=IQ+K
      IMR=JQ+K
      X=R(ILR)
      R(ILR)=R(IMR)
   20 R(IMR)=X
   30 CONTINUE
      RETURN
      END
      SUBROUTINE PLOTX(A,B,N,IP)
      DIMENSION A(1),B(1)
      IF(IP.EQ.0) GOTO10
      CALL GRAPH(4.0,0.0,4.0,0.0)
      CALL LIMIT(A,N,1)
      CALL LIMIT(B,N,2)
C     CALL SCALE(1)
      GO TO 11
   10 CALL GRAPH(+1.2,-1.2,1.0,-1.0)
      CALL AXES
   11 CONTINUE
      CALL LOC2(A,B,N,1)
      CALL PLOT
      IF(IP.EQ.1) GO TO 12
      CALL IDENT
   12 CONTINUE
      RETURN
      END
      SUBROUTINE GRAPH(XMX,XMN,YMX,YMN)
      DIMENSION Z(1),X(1),Y(1),SMALL(13),CHAR(50),HOLL(11)
      REAL ITEM(55,121)
      DATA CHAR     /'1','2','3','4','5','6','7','8','9','A','B','C',
     1 'D','E','F','G','H','I','J','K','L','M','N','O','P','Q','R','S',
     2 'T','U','V','W','X','Y','Z','+','/','=','*','&','$','@','%','?',
     3 '<','(',')','"',';',' '/
      DATA HOLL     /' ','X','2','3','4','5','6','7','8','9','M'/
      DATA BLANK,DD,PLUS,AIE,AMINUS,ORIG,AMULT,GT50/' ','.','+','\',
     1 '-','0','#','>'/
      DO 115 I = 1,55
      DO 115 J = 1,121
  115 ITEM(I,J) = BLANK
      XMAX = XMX
      XMIN = XMN
      YMAX = YMX
      YMIN = YMN
      RETURN
      ENTRY LIMIT(Z,NP,IS)
      ZMAX = -1.0E9
      ZMIN = +1.0E9
      DO 10 I = 1,NP
      IF(Z(I) .GT. ZMAX) ZMAX = Z(I)
   10 IF(Z(I) .LT. ZMIN) ZMIN = Z(I)
      IF(IS .EQ. 2) GO TO 11
      XMAX=ZMAX
      XMIN=ZMIN
      GO TO 12
   11 CONTINUE
      YMAX=ZMAX
      YMIN=ZMIN
   12 RETURN
      ENTRY SCALE(IP)
      XEXP = 0.05 * (XMAX - XMIN)
      XMAX = XMAX + XEXP
      XMIN = XMIN - XEXP
      YEXP = 0.05 * (YMAX - YMIN)
      YMAX = YMAX + YEXP
      YMIN = YMIN - YEXP
      IF(IP .NE. 1) RETURN
      IF(YMAX - XMAX) 20,21,22
   20 YMAX = XMAX
      GO TO 21
   22 XMAX = YMAX
   21 CONTINUE
      IF(YMIN - XMIN) 23,24,25
   23 XMIN = YMIN
      GO TO 24
   25 YMIN = XMIN
   24 EXP = 0.1667 * (XMAX - XMIN)
      XMAX = XMAX + EXP
      XMIN = XMIN - EXP
      RETURN
      ENTRY AXES
      DELX = (XMAX - XMIN)/120.0
      DELY = (YMAX - YMIN)/54.0
      K = 0
      M = 0
      IF(XMIN .LT. 0.0 .AND. XMAX .GT. 0.0) GO TO 30
      GO TO 31
   30 K = (-XMIN/DELX) + 1.5
      DO 32 I = 1,55
   32 ITEM(I,K) = AIE
   31 CONTINUE
      IF(YMIN .LT. 0.0 .AND. YMAX .GT. 0.0) GO TO 33
      GO TO 34
   33 M = (YMAX/DELY) + 1.5
      DO 35 J = 1,121
   35 ITEM(M,J) = AMINUS
      IF(K .EQ. 0 .OR. M .EQ. 0) GO TO 34
      ITEM(M,K) = ORIG
   34 RETURN
      ENTRY LOC1(X,Y,NPOI)
      DELX = (XMAX - XMIN)/120.0
      DELY = (YMAX - YMIN)/54.0
      DO 50 II = 1,NPOI
      I = (YMAX - Y(II))/DELY + 1.5
      J = (X(II) - XMIN)/DELX + 1.5
      IF(I.GT.55.OR.I.LT.1.OR.J.GT.121.OR.J.LT.1) GO TO 50
      DO 52 JJ = 1,10
      IF(ITEM(I,J) .EQ. HOLL(JJ)) GO TO 51
   52 CONTINUE
      IF(ITEM(I,J).EQ.AIE.OR.ITEM(I,J).EQ.AMINUS.OR.ITEM(I,J).EQ.ORIG)
     1  ITEM(I,J) = HOLL(2)
      GO TO 50
   51 ITEM(I,J) = HOLL(JJ+1)
   50 CONTINUE
      RETURN
      ENTRY LOC2(X,Y,NPOI,ISTART)
      DELX = (XMAX - XMIN)/120.0
      DELY = (YMAX - YMIN)/54.0
      ISPN = ISTART + NPOI - 1
      DO 53 II = ISTART,ISPN
      I = (YMAX - Y(II))/DELY + 1.5
      J = (X(II) - XMIN)/DELX + 1.5
      IF(I.GT.55.OR.I.LT.1.OR.J.GT.121.OR.J.LT.1) GO TO 53
      IF(ITEM(I,J).EQ.BLANK.OR.ITEM(I,J).EQ.AIE.OR.ITEM(I,J).EQ.AMINUS
     1 .OR.ITEM(I,J).EQ.ORIG) GO TO 54
      ITEM(I,J) = AMULT
      GO TO 53
   54 CONTINUE
      IF(II .GT. 50) GO TO 55
      ITEM(I,J) = CHAR(II)
      GO TO 53
   55 ITEM(I,J) = GT50
   53 CONTINUE
      RETURN
      ENTRY LOC3(X,Y,NPOI,POINT)
      DELX = (XMAX - XMIN)/120.0
      DELY = (YMAX - YMIN)/54.0
      DO 56 II = 1,NPOI
      I = (YMAX - Y(II))/DELY + 1.5
      J = (X(II) - XMIN)/DELX + 1.5
      IF(I.GT.55.OR.I.LT.1.OR.J.GT.121.OR.J.LT.1)  GO TO 56
      IF(ITEM(I,J).EQ.BLANK.OR.ITEM(I,J).EQ.AIE.OR.ITEM(I,J).EQ.AMINUS
     1 .OR. ITEM(I,J) .EQ. ORIG .OR. ITEM(I,J) .EQ. POINT) GO TO 57
      ITEM(I,J) = AMULT
      GO TO 56
   57 ITEM(I,J) = POINT
   56 CONTINUE
      RETURN
      ENTRY PLOT
      SMALL(1) = XMIN
      DO 120 I = 1,12
  120 SMALL(I+1) = SMALL(I) + 10.0 * DELX
      VALUE = YMAX + DELY
      WRITE(6,301)
  301 FORMAT(9X,' +.........+.........+.........+.........+.........+...
     1......+.........+.........+.........+.........+.........+.........
     2+')
      DO 330 I = 1,55
      VALUE = VALUE - DELY
      II = I - 1
      L = II + 6
      IF(L/6*6 - L) 320,310,320
  310 A = PLUS
      WRITE(6,302) VALUE,A,(ITEM(I,J), J=1,121),A
  302 FORMAT(1X,F8.3,123A1)
      GO TO 330
  320 A = DD
      WRITE(6,303)  A,(ITEM(I,J), J=1,121),A
  303 FORMAT(9X,123A1)
  330 CONTINUE
      WRITE(6,301)
      WRITE(6,304)  (SMALL(K), K=1,12)
  304 FORMAT(4X,12F10.3)
      WRITE(6,305) SMALL(13)
  305 FORMAT(1H+,123X,F8.2/)
      RETURN
      ENTRY IDENT
      WRITE(6,400) (J, J=1,25)
  400 FORMAT(/35X,'*****IDENTIFICATION KEY FOR PLOTS WITH IDENTIFIED POI
     1NTS*****'//' PT #',25I5)
      WRITE(6,401)  (CHAR(I), I=1,25)
  401 FORMAT(' CHAR',25(4X,A1)/)
      WRITE(6,402)  (J, J=26,50)
  402 FORMAT(' PT #', 25I5)
      WRITE(6,403) (CHAR(I), I=26,50)
  403 FORMAT(' CHAR', 25(4X,A1)//)
      WRITE(6,404) GT50, AMULT
  404 FORMAT(' POINT NUMBERS ABOVE 50 IDENTIFIED AS  ',A1/' MULTIPLE POI
     1NTS IDENTIFIED AS  ',A1/)
      RETURN
      END
      SUBROUTINE CORREL(X,Y,N,ID)
      REAL*4 X(1),Y(1)
      SX=0.0
      SY=0.0
      SXX=0.0
      SXY=0.0
      SYY=0.0
      FN=N
      DO 10 I=1,N
      SX=SX+X(I)/FN
      SY=SY+Y(I)/FN
      SXX=SXX+X(I)*X(I)/FN
      SYY=SYY+Y(I)*Y(I)/FN
      SXY=SXY+X(I)*Y(I)/FN
   10 CONTINUE
      A=SXY-SX*SY
      B=SXX-SX*SX
      C=SYY-SY*SY
      R=A/SQRT(B*C)
      RS=R*R
      WRITE(6,1) ID,R,RS
    1 FORMAT(//5X,'CORRELATION BETWEEN ORIGINAL AND FITTED VECTORS FOR
     XPROPERTY NO.',I3,1X,'IS ',F8.3,' , RSQ=',F8.3/)
      RETURN
      END