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