C prefmap2.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. (1972). "Individual differences and multidimensional scaling". C In R.N.Shepard, A.K.Romney & S.B.Nerlove (Eds.), C "Multidimensional Scaling: Theory and Applications in the Behavioral C Sciences" (Vol.1, pp.105-155). New York and London: Seminar Press. C [Reprinted in P.Davies & A.P.M.Coxon (Eds.) (1984), "Key Texts in C Multidimensional Scaling", pp. 267-301, Exeter, NH: Heinemann.] C Carroll,J.D. (1980). "Models and methods for multidimensional C analysis of preferential choice (or other dominance) data". C In E.D.Lantermann & H.Feger (Eds.), "Similarity and Choice", pp. 234-289. C Bern: Hans Huber. C----+----@----+----@----+----@----+----@----+----@----+----@----+----@ $ FORTRAN C PREFMAP2 - PREFERENCE MAPPING VIA A GENERA- C LIZATION OF COOMBS UNFOLDING MODEL C BY CARROLL AND CHANG, NOVEMBER, 1972 C PARAMETER DEFINITION C N - NUMBER OF STIMULI (MAX. = 30) C K - NUMBER OF DIMENSIONS (MAX. = 5) C NSUB - NUMBER OF SUBJECTS (MAX. = 30) C ISV- 0, SMALLER VALUE IS THE MORE PREFERRED C ISV- 1, LARGER VALUE IS THE MORE PREFERRED C NORS - 1, NORMALIZE SCALE VALUE, MAKE LENGTH =1. C NORS - 0, DO NOT NORMALIZE SCALE VALUE. C IRX- 0, X IS PUNCHED N BY K C IRX-1, X IS PUNCHED K BY N C IPS - STARTING PHASE (MINIMUM = 1, MAXIMUM = 4). C IPE - ENDING PHASE (MINIMUM = 1, MAXIMUM = 4). C IRWT - 1, WHEN IPS=3, READ IN WEIGHTS, ONE FOR EACH DIMENSION. C IRWT - 0, DO NOT READ IN WEIGHTS. C LFITSW - 0, LINEAR FIT. C LFITSW - 1, MONOTONE FIT (PLAIN) C LFITSW - 2, MONOTONE FIT (BLOCK MONOTONE WITH ORDERING WITHIN). C LFITSW - 3, MONOTONE FIT (BLOCK MONOTONE WITH EQUALITY WITHIN). C IAV - 0, AVERAGE SUBJ SCALE VALUES ARE COMPUTED IN IPS PHASE. C IAV - 1, AVERAGE SUBJ SCALE VALUES ARE COMPUTED IN EACH PHASE. C MAXIT - (MONOTONE ONLY) MAXIMUM NUMBER OF ITERATIONS. C ISHAT - 0, AT THE BEGINNING OF A NEW PHASE, USE THE LAST SHAT FROM C THE PREVIOUS PHASE. C ISHAT - 1, AT THE BEGINNING OF A NEW PHASE, USE THE ORIGINAL S. C IPLOT - 0, PLOT AVERAGE SUBJECT ONLY. C IPLOT - 1, PLOT AVERAGE SUBJ, FOR EACH SUBJ DO FUNCTION PLOTS. C IPLOT - 2, PLOT AVERAGE S, FOR EACH S DO FUNCTION AND IDEAL PLOTS. C MDV - 1, DRAW ROTATION AND WEIGHTS FOR EACH IDEAL POINT. C MDV - 0, DRAW IDEAL POINTS WITHOUT VECTORS. C MLIP - 1, LABEL IDEAL POINTS WITH LETTERS. C MLIP - 0, LABEL IDEAL POINTS WITH NUMBERS. C MIDEN - 1, READ IN SYMBOLS TO LABEL STIMULI ON PLOTS. C MIDEN - 0, THE PROGRAM GENERATES NUMBERS TO LABEL STIMULI ON PLOTS. C MOPT2 - RELEVANT ONLY WHEN IPS=2 C =1, DO ROTATION AND WEIGHT C =0, NONE C MOPT3 - RELEVANT ONLY WHEN IPS=3 C =2, DO ROTATION AND WEIGHT C =1, DO WEIGHT ONLY C =0, NONE C IST - OPTION ON STIMULUS SPACE C =0, DOUBLE CENTER S - SUBTRACT ROW MEAN AND DIVIDE BY C SQRT(SSQ), THEN SUBTRACT COLUMN MEAN. C =1, SUBTRACT ROW MEAN ONLY. C =2, SUBTRACT ROW MEAN AND DIVIDE BY SQRT(SSQ) C =3, SUBTRACT ROW MEAN AND SUBTRACT COLUMN MEAN. C =4, READ IN X. C VARIABLES. C X - COORDINATES OF N STIMULI IN K DIMENSIONAL SPACE C S(I,J) - SCALE VALUE ON THE ITH STIMULUS FOR THE JTH SUBJECT DIMENSION U(5,5),ROOT(5) DIMENSION D(30),S(30,32),RSQ(32,4),Q(5,5),Z(5,5) DIMENSION X(30,5),P(5,32),T(5,5,32),SGN(5,32),XST(30,5),PST(5) DIMENSION PAL(32),BNUB(30),IX1(5),IY1(5) DIMENSION DHAT(30),SCR(30),SHAT(30,32),KC(30) DIMENSION W(30,21),WA(21,21),WI(975) DIMENSION WW(21,30),B(21),R(5,5),SS(30),LBLOCK(30) DIMENSION FAS(10),FVT(3),FHR(10) DIMENSION TIT(14),FPH(5),FOD(3),FS(3),FSA(3) DIMENSION FMTS(12),SIGN(5),FMTI(5),SDHS(30) COMMON N,K,XN,XK,NS,RMAX,IPH,NSUB,IS,NV,K1 COMMON D,S,RSQ,X,P,T,SGN,XST,PST,U,ROOT COMMON PAL,BNUB,NM,IX1,IY1 COMMON /FMIC/TIT,FPH,FOD,FUU,ZRO,FS,FSA DATA FVT(1)/18H(13HSQRD DISTANCE)/ DATA FHR(1)/54H(23HSCALE VALUE SUBJECTI3,3X,5HPHASEI2,3X,3HIT= 1I2)/ DATA FAS(1)/59H(28HSCALE VALUE AVERAGE SUBJI3,3X,5HPHASEI2,3X, 13HIT=I2)/ C SET UP MICROFILM PLOTTING IX1(1)=25 IX1(2)=1000 IX1(3)=1000 IX1(4)=25 IX1(5)=25 IY1(1)=25 IY1(2)=25 IY1(3)=1000 IY1(4)=1000 IY1(5)=25 NL=30 500 READ116,N,K,NSUB,ISV,NORS,IRX,IPS,IPE,IRWT,LFITSW,IAV,MAXIT,ISHAT, 1IPLOT,MDV,MLIP,MIDEN,MOPT2,MOPT3,IST IF(N.EQ.0) GOTO998 IF(MAXIT.EQ.0.AND.LFITSW.GT.0) MAXIT=15 PRINT124 IF(LFITSW.EQ.0) IAV=0 READ100,(TIT(I),I=2,13) XN=N XK=K C READ IN SCALE VALUES, EACH ROW REPRESENTS ONE SUBJECT READ100,FMTS DO211IM=1,NSUB READFMTS,(S(J,IM),J=1,N) IF(IPS.EQ.4) GOTO205 IF(ISV.EQ.0)GOTO135 209 DO137I=1,N 137 S(I,IM)=-S(I,IM) GOTO135 205 IF(ISV.EQ.0) GOTO209 C NORMALIZE S 135 IF(NORS)133,211,133 133 CALL NORMS(S(1,IM),N,XN) 211 CONTINUE C TEST IF READ IN COORDINATES OF STIMULUS POINTS IF(IST.NE.4) GOTO140 CALL READX(X,N,K,IRX,NL) 140 PRINT104,(TIT(I),I=2,13) PRINT101 PRINT114 PRINT117,N,K,NSUB,ISV,NORS,IRX,IPS,IPE,IRWT,LFITSW,IAV,MAXIT,ISHAT 1,IPLOT,MDV,MLIP,MIDEN,MOPT2,MOPT3,IST PRINT101 LFITSW=LFITSW-1 IF(IST.EQ.4)GOTO138 DO139I=1,N DO139J=1,NSUB 139 SHAT(I,J)=S(I,J) CALL STISP(SHAT,NL,N,NSUB,K,IST,X,T) DO299I=1,K 299 PRINT103,I,(X(J,I),J=1,N) C TEST IF IPS=3, IF SO SET OR READ IN WEIGHTS 138 IF(IPS.NE.3) GOTO42 IF(IRWT.EQ.0) GOTO43 READ100,(FMTS(I),I=1,12) READFMTS,(SIGN(I),I=1,K) GOTO42 43 DO44I=1,K 44 SIGN(I)=1. 42 IF(LFITSW.LT.0) GOTO45 READ118,CRIT C SET UP MICROFILM PLOTTING 45 IF(MIDEN)41,40,41 40 CALL LMIC(N,BNUB) GOTO46 41 READ100,(FMTI(I),I=1,5) READFMTI,(BNUB(I),I=1,N) 46 NS=NSUB&1 IF(MLIP)48,47,48 47 DO49I=1,NS 49 CALL PIBCI(I,PAL(I),4H(I2)) GOTO50 C USE ALPHABET TO LABEL IDEAL POINTS 48 CALL LET(PAL,NS) 50 NV=K*NS*2 XNS=NSUB K1=K-1 MK=((K&1)*(K&2))/2 IS=1 CALL ZEROD(Q,5,5) C TEST IF ENTER PHASE 2 AND OPTION ON ROTATION AND WEIGHT IF(IPS.EQ.2.AND.MOPT2.EQ.1) GOTO206 C TEST IF ENTER PHASE 3 IF(IPS.EQ.3) GOTO208 GOTO207 208 IF(MOPT3-1)207,210,206 210 M=2*K&1 CALL COMPW(W,X,WW,WA,WI,M,SIGN,1) CALL WGT(W,N,M,S,K,X,NSUB,SIGN) GOTO207 206 M=MK CALL COMPW(W,X,WW,WA,WI,M,SIGN,2) CALL ROWT(W,M,B,R,SIGN,WA) 207 IPH=IPS CALL ZEROD(RSQ,32,4) 95 GOTO(82,92,78,99),IPH 82 M=MK CALL COMPW(W,X,WW,WA,WI,MK,SIGN,2) CALL ROWT(W,M,B,R,SIGN,WA) GOTO79 99 M=K&1 GOTO79 78 M=K&2 GOTO79 92 M=2*K&1 79 CALL ZEROT(T,5,5,32) 81 PRINT106,IPH PRINT105 DO97I=1,K 97 PRINT103,I,(X(J,I),J=1,N) CALL COMPW(W,X,WW,WA,WI,M,SIGN,0) CALL ZEROS(SS,N) RMAX=0. DO90IS=1,NS IT=1 LAST=0 IF(IS.EQ.NS)GOTO141 C SUBJECTS PRINT101 PRINT113,IS IF(IPH.EQ.IPS) GOTO132 IF(ISHAT)132,131,132 C AVERAGE SUBJECT 141 PRINT115 IF(IPH.EQ.IPS)GOTO85 IF(IAV.EQ.0) GOTO131 C OPTION IAV=1 - COMPUTE AVERAGE SUBJ. FROM THE CURRENT PHASE 85 DO28I=1,N 28 S(I,IS)=SS(I)/XNS C PUT S IN SHAT 132 DO34I=1,N 34 SHAT(I,IS)=S(I,IS) 131 IF(IS.EQ.NS.OR.LFITSW.LT.0) GOTO134 C SORT S FOR MONOTONE FIT DO5I=1,N SCR(I)=S(I,IS) 5 KC(I)=I CALL SORTFL CALL SORT(SCR(1),SCR(2),N,KC) C SET UP LBLOCK FOR BLOCK MONOTONE IF(LFITSW.GT.0) CALL BMONO(SCR,LBLOCK,N) C PRINT SHAT 134 PRINT108 PRINT109,(SHAT(I,IS),I=1,N) C COMPUTE BETA IF(LFITSW.LT.0.OR.IS.EQ.NS) GOTO31 PRINT111 31 DO11I=1,M B(I)=0. DO11J=1,N 11 B(I)=B(I)&WW(I,J)*SHAT(J,IS) IF(IT.EQ.1) GOTO152 DO6I=1,N D(I)=0. DO6J=1,M 6 D(I)=D(I)&B(J)*W(I,J) GOTO155 C COMPUTE B 152 PRINT102,(B(I),I=1,M) DO150I=1,K 150 B(I)=B(I&1) IF(IPH.EQ.4)GOTO84 IF(IPH-2)33,32,35 35 II=K&2 GOTO32 C COMPUTE R 33 CALL COMPR(R,B,WA,K,XK) C COMPUTE WEIGHTS OF R*(RT) IF(IS.EQ.NS) GOTO177 SUMD=R(1,1)**2 SUMR=0. DO175I=2,K II=I-1 SUMD=SUMD&R(I,I)**2 DO175J=1,II 175 SUMR=SUMR&R(I,J)**2 SUMR=SUMD&2.*SUMR C COMPUTE R*(RT) DO179I=1,K DO179J=1,I Z(I,J)=0. DO179KK=1,K 179 Z(I,J)=Z(I,J)&R(I,KK)*R(J,KK) C COMPUTE P C INVERT R (R IN WA) 177 MM=1 CALL MTINV(WA,21,21,WI,21,21,K,MM) IF(MM.EQ.1) GOTO53 PRINT107 GOTO998 53 DO51I=1,K P(I,IS)=0. DO52J=1,K 52 P(I,IS)=P(I,IS)&B(J)*WA(J,I) 51 P(I,IS)=-.5*P(I,IS) KM=1 CALL MTEVV(R,5,5,U,5,5,K,KM) C COMPUTE T CALL COMPT(R) CALL COMPD(XST,PST,SGN(1,IS),SHAT(1,IS),RSQ(IS,IPH)) PRINT110,RSQ(IS,IPH) C COMPUTE Q=Q&W*R*(RT) IF(IPH.NE.1.OR.IS.EQ.NS) GOTO154 WEI=RSQ(IS,1)**2/SUMR DO178I=1,K DO178J=1,I 178 Q(I,J)=Q(I,J)&WEI*Z(I,J) 154 IF(IS.EQ.NS) GOTO71 IF(IPLOT.LT.2) GOTO155 71 IF(K1)153,155,153 153 CALL PTXP(IT) GOTO155 C PHASE 4 84 SUMB=0. DO136I=1,K 136 SUMB=SUMB&B(I)**2 SUMB=SQRT(SUMB) DO86I=1,K 86 P(I,IS)=B(I)/SUMB GOTO87 C PHASE 2 AND PHASE 3 32 DO55I=1,K IF(IPH-2)55,75,76 76 ROOT(I)=B(II)*SIGN(I) GOTO77 75 II=I&K&1 ROOT(I)=B(II) 77 IF(ROOT(I))56,57,57 56 SGN(I,IS)=-1. GOTO58 57 SGN(I,IS)=1. 58 T(I,I,IS)=ABS(ROOT(I)) IF(T(I,I,IS)-RMAX)55,55,59 59 RMAX=T(I,I,IS) 55 CONTINUE C COMPUTE P DO60I=1,K 60 P(I,IS)=-.5*B(I)/ROOT(I) 87 CALL COMPD(X,P(1,IS),ROOT,SHAT(1,IS),RSQ(IS,IPH)) PRINT110,RSQ(IS,IPH) 155 IF(IPH.NE.4) GOTO21 156 PRINT129 157 PRINT109,(D(I),I=1,N) 21 IF(IS.EQ.NS) GOTO73 IF(LFITSW.LT.0) GOTO173 C REARRANGE D IN DHAT, FIND MONOTONE FIT OF DHAT AND STORE IN D CALL MONO(DHAT,D,KC,N,LFITSW,SCR,LBLOCK) C REARRANGED AND UNNORMALIZED SHAT IN D AND SHAT IN SCR C NORMALIZE SHAT (IN SCR) IN ORIGINAL ORDER CALL NORMS(SCR,N,XN) IF(LAST.EQ.1) GOTO61 C TEST CRITERION ITEST=0 C SHAT HAS SHAT OF PREVIOUS IT. AND SCR HAS SHAT OF PRESENT IT. SUM=0. DO2I=1,N 2 SUM=SUM&(SHAT(I,IS)-SCR(I))**2 IF(SUM.GT.CRIT) ITEST=1 IF(IT.EQ.1) GOTO61 IF(ITEST.EQ.0) GOTO160 IF(IT.LT.MAXIT) GOTO203 PRINT121,MAXIT ITEST=0 GOTO161 160 PRINT122,IT 161 LAST=1 GOTO152 C PLOT S VS. DSQ C REARRANGE S IN SDHS 61 DO202I=1,N J=KC(I) 202 SDHS(I)=S(J,IS) IF(IPLOT.EQ.0) GOTO62 IF(IS.EQ.NS) GOTO204 CALL PSD(SDHS,DHAT,D,IT,TIT,FVT,FHR) GOTO62 204 CALL PSD(SDHS,DHAT,D,IT,TIT,FVT,FAS) C PRINT SUBJECT 62 IF(LAST.EQ.0.OR.IPH.EQ.4) GOTO203 64 CALL PRSUB PRINT101 C PUT SCR IN SHAT 203 DO63I=1,N 63 SHAT(I,IS)=SCR(I) IF(ITEST.EQ.0) GOTO172 IT=IT&1 GOTO31 173 IF(IPLOT.EQ.0) GOTO174 73 IT=0 CALL PLSD(S(1,IS),D,IT,TIT,FVT,FHR,FAS) 174 IF(IPH.EQ.4) GOTO172 PRINT101 CALL PRSUB PRINT101 172 IF(IPH.EQ.IPS) GOTO151 IF(LFITSW.LT.0) GOTO90 IF(IAV.EQ.0) GOTO90 151 DO27I=1,N 27 SS(I)=SS(I)&SHAT(I,IS) 90 CONTINUE IF(IPH.EQ.4)GOTO88 CALL PXPV(MDV) IPH=IPH&1 IF(IPH.GT.IPE) GOTO98 IF(IPH-3)91,94,12 C BEFORE PHASE 2 DO LINEAR TRANSFORMATION 91 DO180I=1,K DO180J=1,N 180 XST(J,I)=X(J,I) CALL TRANX(X,Q,Z,N ,K,XST,NL) GOTO95 94 DO181I=1,K DO181J=1,N 181 X(J,I)=XST(J,I) GOTO95 C FOR PHASE 4, REVERSE THE PREFERENCES, SO HIGHER VALUE MEANS C GREATER PREFERENCE. 12 DO22I=1,NS DO22J=1,N SHAT(J,I)=-SHAT(J,I) 22 S(J,I)=-S(J,I) GOTO95 88 CALL PJIDV(P,X,N,K,NS,TIT,BNUB,PAL) PRINT125 PRINT126,(I,I=1,K) DO201I=1,NSUB 201 PRINT127,I,(P(J,I),J=1,K) PRINT128,(P(J,NS),J=1,K) 98 CALL CPR(RSQ,W ,N,K,NS,XN,IPS,IPE) GOTO500 998 CALL CLEAN 100 FORMAT(12A6) 101 FORMAT(120H0****************************************************** 1*****************************************************************) 102 FORMAT(5H BETA/(2X,10E13.5)) 103 FORMAT(I4,10E12.4/(4X,10E12.4)) 104 FORMAT(1H012A6) 105 FORMAT(16H0STIMULUS MATRIX) 106 FORMAT(6H1PHASEI3) 107 FORMAT(21H0R CANNOT BE INVERTED) 108 FORMAT(27H0S (VECTOR OF SCALE VALUES)) 109 FORMAT(2X,10E13.5) 110 FORMAT(46H0R(CORRELATION BETWEEN DATA AND SQUARED DIST)=E13.5) 111 FORMAT(32H0BEGIN ITERATION ON MONOTONE FIT) 113 FORMAT(8H0SUBJECTI4) 114 FORMAT(115H0 N K NSUB ISV NORS IRX IPS IPE IRWT LFITSW I 1AV MAXIT ISHAT IPLOT MDV MLIP MIDEN MOPT2 MOPT3 IST) 115 FORMAT(16H0AVERAGE SUBJECT) 116 FORMAT(20I4) 117 FORMAT(2I3,3I5,I6,3I5,3I7,5I6,I8,2I7) 118 FORMAT(10F7.2) 121 FORMAT(28H0REACHED MAXIMUM ITERATION =I4) 122 FORMAT(30H0REACHED CRITERION, ITERATION=I4) 124 FORMAT(1H1//68H0MDSCALING VIA A GENERALIZATION OF COOMBS UNFOLDING 1 MODEL - PREFMAP2) 125 FORMAT(44H0DIRECTION COSINES OF FITTED SUBJECT VECTORS//21X,9HDIME 1NSION) 126 FORMAT(10H0 SUBJECTI5,9I9) 127 FORMAT(I7,3X,10F9.4) 128 FORMAT(4X,3HAVR,3X,10F9.4) 129 FORMAT(33H0PROJECTIONS ON THE FITTED VECTOR) 999 STOP END $ FORTRAN SUBROUTINE TRANX(X,Q,B,N,K,SP,NL) C DO NOT NORMALIZE DIMENSIONS DIMENSION X(NL,5),Q(5,5),B(5,5),SP(NL,5),P(5) DO6I=2,K II=I-1 DO6J=1,II 6 Q(J,I)=Q(I,J) PRINT100 DO10I=1,K 10 PRINT101,I,(Q(I,J),J=1,K) MX=1 CALL MTEVV(Q,5,5,B, 5, 5,K,MX) PRINT102 PRINT101,I,(Q(I,I),I=1,K) PRINT103 DO11J=1,K 11 PRINT101,J,(B(I,J),I=1,K) DO2J=1,N DO3I=1,K P(I)=0. DO3M=1,K 3 P(I)=P(I)&B(M,I)*X(J,M) DO4L=1,K 4 X(J,L)=P(L) 2 CONTINUE PRINT104 DO12I=1,K 12 PRINT101,I,(X(J,I),J=1,N) 100 FORMAT(2H0Q) 101 FORMAT(I3,2X,10E12.4/(5X,10E12.4)) 102 FORMAT(18H0EIGEN VALUES OF Q) 103 FORMAT(19H0EIGEN VECTORS OF Q) 104 FORMAT(29H0OPTIMAL GROUP STIMULUS SPACE) RETURN END $ FORTRAN SUBROUTINE ROWT(W,M,B,R,SIGN,WA) DIMENSION U(5,5),ROOT(5) DIMENSION D(30),S(30,32),RSQ(32,4) DIMENSION X(30,5),P(5,32),T(5,5,32),SGN(5,32),XST(30,5),PST(5) DIMENSION W(30,21),B(21),R(5,5),SIGN(5),WA(21,21),Z(21) COMMON N,K,XN,XK,NS,RMAX,IPH,NSUB,IS,NV,K1 COMMON D,S,RSQ,X,P,T,SGN,XST,PST,U,ROOT XK=K CALL OPTM(W,N,M,S,Z,K,NSUB) II=K&1 MM=M-II DO18I=1,MM II=II&1 18 B(II)=Z(I) CALL COMPR(R,B,WA,K,XK) KM=1 CALL MTEVV(R,5,5,U,5,5,K,KM) CALL COMPT(R) C KEEP SIGN FOR PHASE 3 DO5I=1,K 5 SIGN(I)=SGN(I,IS) DO93I=1,N DO93J=1,K 93 X(I,J)=XST(I,J) PRINT101 101 FORMAT(2H0X) DO200I=1,K PRINT103,I,(X(J,I),J=1,N) 200 PUNCH104,I,(XST(J,I),J=1,N) 103 FORMAT(I3,3X,10E12.4/(6X,10E12.4)) 104 FORMAT(I2,10F7.3/(2X,10F7.3)) CALL PROPT(1,DUMY) RETURN END $ FORTRAN SUBROUTINE WGT(W,N,M,S, K,X,NSUB,SIGN) DIMENSION W(30,21),S(30,32),T(5),X(30,5),SIGN(5) CALL OPTM(W,N,M,S,T,K,NSUB) DO2I=1,K IF(T(I))4,5,5 4 SIGN(I)=-1. GOTO6 5 SIGN(I)=1. 6 Q=ABS(T(I)) Q=SQRT(Q) DO2J=1,N 2 X(J,K)=Q*X(J,K) CALL PROPT(2,T) RETURN END $ FORTRAN SUBROUTINE OPTM(W,N,M,S,T,K,NSUB) DIMENSION W(30,21),S(30,32) DIMENSION Z(30,30),V(30,30),A(30,30),E(15) DIMENSION X1(15,30),X2(6,30),P(15,15),T(15) M2=K&1 DO2I=1,M2 DO2J=1,N 2 X2(I,J)=W(J,I) MM=M2&1 L=0 DO3I=MM,M L=L&1 DO3J=1,N 3 X1(L,J)=W(J,I) M1=L PRINT100 100 FORMAT(3H0X1) DO30I=1,M1 30 PRINT101,I,(X1(I,J),J=1,N) 101 FORMAT(I3,3X,10E12.4/(6X,10E12.4)) PRINT102 102 FORMAT(3H0X2) DO31I=1,M2 31 PRINT101,I,(X2(I,J),J=1,N) C COMPUTE A DO5I=1,M2 DO5J=1,M2 Z(I,J)=0. DO5L=1,N 5 Z(I,J)=Z(I,J)&X2(I,L)*X2(J,L) MX=1 CALL MTINV(Z,30,30,A,15,15,M2,MX) IF(MX.EQ.1) GOTO8 PRINT116 STOP 116 FORMAT(21H0Z CANNOT BE INVERTED) 8 DO12I=1,M2 DO12J=1,N V(I,J)=0. DO12L=1,M2 12 V(I,J)=V(I,J)&Z(I,L)*X2(L,J) DO4I=1,N DO13J=1,N A(I,J)=0. DO13L=1,M2 13 A(I,J)=A(I,J)-X2(L,I)*V(L,J) 4 A(I,I)=1.&A(I,I) PRINT103 103 FORMAT(2H0A) DO32I=1,N 32 PRINT101,I,(A(I,J),J=1,N) C COMPUTE P DO17I=1,N DO17J=1,NSUB V(I,J)=0. DO17L=1,N 17 V(I,J)=V(I,J)&A(I,L)*S(L,J) DO18I=1,M1 DO18J=1,NSUB Z(I,J)=0. DO18L=1,N 18 Z(I,J)=Z(I,J)&X1(I,L)*V(L,J) DO19I=1,M1 DO19J=1,M1 P(I,J)=0. DO19L=1,NSUB 19 P(I,J)=P(I,J)&Z(I,L)*Z(J,L) PRINT104 104 FORMAT(2H0P) DO33I=1,M1 33 PRINT101,I,(P(I,J),J=1,M1) C COMPUTE Q AND STORE IN A DO14I=1,N DO14J=1,M1 Z(I,J)=0. DO14L=1,N 14 Z(I,J)=Z(I,J)&A(I,L)*X1(J,L) DO15I=1,M1 DO15J=1,M1 A(I,J)=0. DO15L=1,N 15 A(I,J)=A(I,J)&X1(I,L)*Z(L,J) PRINT105 105 FORMAT(2H0Q) DO34I=1,M1 34 PRINT101,I,(A(I,J),J=1,M1) MX=1 CALL MTEVV(A,30,30,V,30,30,M1,MX) PRINT106 106 FORMAT(18H0EIGEN VALUES OF Q) PRINT107,(A(I,I),I=1,M1) 107 FORMAT(2X,10E12.4) PRINT108 108 FORMAT(19H0EIGEN VECTORS OF Q) DO35I=1,M1 35 PRINT101,I,(V(J,I),J=1,M1) 39 DO16I=1,M1 16 T(I)=1./SQRT(A(I,I)) C COMPUTE (VT)P DO20I=1,M1 DO20J=1,M1 Z(I,J)=0. DO20L=1,M1 20 Z(I,J)=Z(I,J)&V(L,I)*P(L,J) C COMPUTE (VT)PV DO21I=1,M1 DO21J=1,M1 A(I,J)=0. DO21L=1,M1 21 A(I,J)=A(I,J)&Z(I,L)*V(L,J) C COMPUTE PST AND STORE IN P DO22I=1,M1 DO22J=1,M1 22 A(I,J)=A(I,J)*T(I)*T(J) PRINT109 109 FORMAT(4H0PST) DO36I=1,M1 36 PRINT101,I,(A(I,J),J=1,M1) DO23I=1,M1 DO23J=1,M1 23 Z(I,J)=T(I)*V(J,I) C FIND ROOTS AND VECTORS OF PST MX=1 CALL MTEVV(A,30,30,V,30,30,M1,MX) PRINT110 110 FORMAT(20H0EIGEN VALUES OF PST) PRINT107,(A(I,I),I=1,M1) PRINT111 111 FORMAT(21H0EIGEN VECTORS OF PST) DO37I=1,M1 37 PRINT101,I,(V(J,I),J=1,M1) C COMPUTE T FROM TSTAR DO24J=1,M1 T(J)=0. DO24L=1,M1 24 T(J)=T(J)&V(L,1)*Z(L,J) PRINT112 112 FORMAT(15H0UNNORMALIZED T) PRINT107,(T(I),I=1,M1) C NORMALIZE T C COMPUTE C C=0. DO26I=1,M1 Z(I,1)=0. DO25J=1,M1 25 Z(I,1)=Z(I,1)&T(J)*P(J,I) 26 C=C&Z(I,1)*T(I) C=SQRT(C) DN=FLOAT(N*K)/C DO27I=1,M1 27 T(I)=DN*T(I) PRINT113,C 113 FORMAT(3H0C=F10.5) PRINT114,DN 114 FORMAT(24H0NORMALIZATION CONSTANT=F10.5) PRINT115 115 FORMAT(13H0NORMALIZED T) PRINT107,(T(I),I=1,M1) RETURN END $ FORTRAN SUBROUTINE COMPW(W,X,WW,WA,WI,M,SIGN,IW) DIMENSION W(30,21),X(30,5),WI(975),WW(21,30),WA(21,21),SIGN(5) COMMON N,K,XN,XK,NS,RMAX,IPH,NSUB,IS,NV,K1 C IW - PARAMETER FOR OPTIMIZATION C =0, THE NORMAL PROCEDURE FOR EACH PHASE. C =1, WEIGHT ONLY C =2, ROTATION AND WEIGHT C COMPUTE W MATRIX DO5I=1,N 5 W(I,1)=1. DO6J=1,K JJ=J&1 JK=JJ&K DO6I=1,N W(I,JJ)=X(I,J) 6 W(I,JK)=X(I,J)**2 IF(IW-1)73,50,70 73 GOTO(70,71,72,71),IPH 70 DO7I=1,N K2=2*K&1 DO7J=1,K1 JJ=J&1 DO7JW=JJ,K K2=K2&1 7 W(I,K2)=X(I,J)*X(I,JW) IF(IW.NE.0) GOTO50 GOTO71 72 DO20I=1,N SUM=0. DO21J=1,K 21 SUM=SUM&SIGN(J)*X(I,J)**2 20 W(I,M)=SUM 71 DO8I=1,M DO8J=1,M WA(I,J)=0. DO8JM=1,N 8 WA(I,J)=WA(I,J)&W(JM,I)*W(JM,J) MX=1 CALL MTINV(WA,21,21,WI,21,21,M,MX) IF(MX.EQ.1) GOTO9 PRINT100 100 FORMAT(22H0WA CANNOT BE INVERTED) STOP 9 DO10I=1,M DO10J=1,N WW(I,J)=0. DO10JM=1,M 10 WW(I,J)=WW(I,J)&WA(I,JM)*W(J,JM) 50 RETURN END $ FORTRAN SUBROUTINE PRSUB DIMENSION U(5,5),ROOT(5) DIMENSION TIT(14),FPH(5),FOD(3),FS(3),FSA(3) DIMENSION D(30),S(30,32),RSQ(32,4) DIMENSION X(30,5),P(5,32),T(5,5,32),SGN(5,32),XST(30,5),PST(5) DIMENSION PAL(52),BNUB(50),IX1(5),IY1(5) COMMON N,K,XN,XK,NS,RMAX,IPH,NSUB,IS,NV,K1 COMMON D,S,RSQ,X,P,T,SGN,XST,PST,U,ROOT COMMON PAL,BNUB,NM,IX1,IY1 COMMON /FMIC/TIT,FPH,FOD,FUU,ZRO,FS,FSA PRINT 1,IS IF(IPH-1) 100,100,200 100 PRINT 10 PRINT 11,(I,I=1,K) PRINT 12 DO 5 I=1,K 5 PRINT 13,I,(U(J,I),J=1,K) PRINT 30 DO 6 I=1,K 6 PRINT32,I,(T(J,I,IS),J=1,K) PRINT80 DO81I=1,K 81 PRINT82,I,(XST(J,I),J=1,N) PRINT 50 PRINT21,(PST(I),I=1,K) 200 PRINT 40 PRINT21,(P(I,IS),I=1,K) PRINT 20 PRINT 21,(ROOT(I),I=1,K) 1 FORMAT(8H0SUBJECTI6) 10 FORMAT(49H0DIRECTION COSINE OF NEW AXES WITH RESPECT TO OLD) 11 FORMAT(1H0,5X,3HOLDI8,4I15) 12 FORMAT(4H0NEW) 13 FORMAT(I8,1X,5E15.5) 20 FORMAT(24H0IMPORTANCES OF NEW AXES) 21 FORMAT(5X,5E15.5) 30 FORMAT(32H0COMPOSITE TRANSFORMATION MATRIX/) 32 FORMAT(I5,5E15.5) 40 FORMAT(37H0IDEAL POINT WITH RESPECT TO OLD AXES) 50 FORMAT(37H0IDEAL POINT WITH RESPECT TO NEW AXES) 80 FORMAT(14H0TRANSFORMED X) 82 FORMAT(I4,10E12.4/(4X,10E12.4)) RETURN END $ FORTRAN SUBROUTINE PROPT(IRW,WGHT) DIMENSION U(5,5),ROOT(5) DIMENSION D(30),S(30,32),RSQ(32,4) DIMENSION X(30,5),P(5,32),T(5,5,32),SGN(5,32),XST(30,5),PST(5) DIMENSION PAL(52),BNUB(50),IX1(5),IY1(5) DIMENSION WGHT(1) COMMON N,K,XN,XK,NS,RMAX,IPH,NSUB,IS,NV,K1 COMMON D,S,RSQ,X,P,T,SGN,XST,PST,U,ROOT COMMON PAL,BNUB,NM,IX1,IY1 GOTO(2,3),IRW 2 PRINT 101 101 FORMAT(28H0OPTIMAL ROTATION AND WEIGHT) 100 PRINT 10 PRINT 11,(I,I=1,K) PRINT 12 DO 5 I=1,K 5 PRINT 13,I,(U(J,I),J=1,K) PRINT 30 DO 6 I=1,K 6 PRINT32,I,(T(J,I,IS),J=1,K) PRINT80 DO81I=1,K 81 PRINT82,I,(XST(J,I),J=1,N) PRINT 20 PRINT 21,(ROOT(I),I=1,K) RETURN 3 PRINT102 102 FORMAT(15H0OPTIMAL WEIGHT) PRINT21,(WGHT(I),I=1,K) 1 FORMAT(8H0SUBJECTI6) 10 FORMAT(49H0DIRECTION COSINE OF NEW AXES WITH RESPECT TO OLD) 11 FORMAT(1H0,5X,3HOLDI8,4I15) 12 FORMAT(4H0NEW) 13 FORMAT(I8,1X,5E15.5) 20 FORMAT(24H0IMPORTANCES OF NEW AXES) 21 FORMAT(5X,5E15.5) 30 FORMAT(32H0COMPOSITE TRANSFORMATION MATRIX/) 32 FORMAT(I5,5E15.5) 80 FORMAT(14H0TRANSFORMED X) 82 FORMAT(I4,10E12.4/(4X,10E12.4)) RETURN END $ FORTRAN SUBROUTINE STISP(S,NL,N,NSUB,K,IST,X,R) DIMENSION S(NL,NL),X(NL,K),R(NL,NL) IST=IST&1 XN=N XNS=NSUB DO50IS=1,NSUB SUM=0. SSQ=0. DO10I=1,N 10 SUM=SUM&S(I,IS) XMN=SUM/XN DO15I=1,N S(I,IS)=S(I,IS)-XMN 15 SSQ=SSQ&S(I,IS)**2 GOTO(16,50,16,50),IST 16 SV=1./SQRT(SSQ) DO21I=1,N 21 S(I,IS)=S(I,IS)*SV 50 CONTINUE C TEST IF SUBTRACTING COLUMN MEAN GOTO(52,51,51,52),IST 52 DO22I=1,N SUM=0. DO23J=1,NSUB 23 SUM=SUM&S(I,J) XMC=SUM/XNS DO24J=1,NSUB 24 S(I,J)=S(I,J)-XMC 22 CONTINUE 51 DO30I=1,N DO30J=1,N R(I,J)=0. DO30M=1,NSUB 30 R(I,J)=R(I,J)&S(I,M)*S(J,M) M=1 CALL MTEVV(R,NL,NL,S,NL,NL,N,M) DO31I=1,K DO31J=1,N 31 X(J,I)=S(J,I) 60 RETURN END $ FORTRAN SUBROUTINE COMPT(R) DIMENSION R(5,5),R2(5) DIMENSION U(5,5),ROOT(5) DIMENSION X(30,5),P(5,32),T(5,5,32),SGN(5,32),XST(30,5),PST(5) DIMENSION D(30),S(30,32),RSQ(32,4) COMMON N,K,XN,XK,NS,RMAX,IPH,NSUB,IS,NV,K1 COMMON D,S,RSQ,X,P,T,SGN,XST,PST,U,ROOT DO16I=1,K ROOT(I)=R(I,I) IF(ROOT(I))17,18,18 17 SGN(I,IS)=-1. GOTO19 18 SGN(I,IS)=1. 19 R2(I)=ABS(ROOT(I)) IF(R2(I).GT.RMAX) RMAX=R2(I) DO82J=1,K 82 T(J,I,IS)=U(J,I)*SQRT(R2(I)) 16 CONTINUE C COMPUTE PST DO20J=1,K PST(J)=0. DO20I=1,K 20 PST(J)=PST(J)&P(I,IS)*T(I,J,IS) C COMPUTE X STAR DO21I=1,N DO21J=1,K XST(I,J)=0. DO21JJ=1,K 21 XST(I,J)=XST(I,J)&X(I,JJ)*T(JJ,J,IS) RETURN END $ FORTRAN SUBROUTINE COMPD(XST,P,SGN,S,R) DIMENSION XST(30,5),P(5),SGN(5),D(30),S(30) COMMON N,K,XN,XK,NS,RMAX,IPH,NSUB,IS,NV,K1 COMMON D SUMD=0. SUMS=0. SDS=0. SSD=0. SSS=0. DO25I=1,N D(I)=0. DO26J=1,K IF(IPH.EQ.4)GOTO20 D(I)=D(I)&SGN(J)*(XST(I,J)-P(J))**2 GOTO26 20 D(I)=D(I)&XST(I,J)*P(J) 26 CONTINUE SUMD=SUMD&D(I) SSD=SSD&D(I)**2 SUMS=SUMS&S(I) SSS=SSS&S(I)**2 25 SDS=SDS&D(I)*S(I) C COMPUTE CORRELATION DEN=(XN*SSD-SUMD**2)*(XN*SSS-SUMS**2) R=(XN*SDS-SUMD*SUMS)/SQRT(DEN) RETURN END $ FORTRAN SUBROUTINE COMPR(R,B,WA,K,XK) DIMENSION R(5,5),B(21),WA(21,21) 33 DO12I=1,K DO12J=1,I IF(I-J)14,15,14 15 II=I&K&1 R(I,I)=B(II) GOTO122 14 XI=J II=(XI&1.)*(XK-XI/2.) II=II&I&1 R(I,J)=.5*B(II) 122 WA(I,J)=R(I,J) WA(J,I)=WA(I,J) R(J,I)=R(I,J) 12 CONTINUE RETURN END $ FORTRAN SUBROUTINE CPR(RSQ,F,N,K,NS,XN,IPS,IPE) DIMENSION RSQ(32,4),F(32,10),DF(10),R(32,4) DIMENSION GF(10),KG(10),KD(10),XM(4) DO3I=1,NS DO3J=1,4 3 R(I,J)=RSQ(I,J)**2 C COMPUTE AND PRINT CORRELATIONS AND F RATIO DO19I=1,10 DF(I)=0. 19 GF(I)=0. CALL ZEROD(F,32,10) IPO=IPE-1 XM(1)=((K&1)*(K&2))/2 XM(2)=2*K&1 XM(3)=K&2 XM(4)=K&1 DO1I=IPS,IPE GF(I)=XM(I)-1. 1 DF(I)=XN-XM(I) DO6I=IPS,IPO II=I&1 DO6J=II,IPE IF(I.EQ.1) GOTO31 L=I&J&3 GOTO7 31 L=J-I&4 7 DF(L)=DF(I) 6 GF(L)=XM(I)-XM(J) DO4I=1,NS DO4J=1,4 4 F(I,J)=R(I,J)/(1.-R(I,J))*DF(J)/GF(J) DO2I=1,NS DO2I1=IPS,IPO II=I1&1 DO2I2=II,IPE IF(I1.EQ.1) GOTO11 L=I2&I1&3 GOTO2 11 L=I2-I1&4 2 F(I,L)=(R(I,I1)-R(I,I2))/(1.-R(I,I1))*DF(L)/GF(L) PRINT100 PRINT101 DO10I=1,10 KG(I)=GF(I) 10 KD(I)=DF(I) PRINT102,(KG(I),KD(I),I=1,4) NSUB=NS-1 XNSUB=NSUB DO20I=1,NSUB 20 PRINT103,I,(RSQ(I,J),J=1,4),(F(I,J1),J1=1,4) PRINT104,(RSQ(NS,J),J=1,4),(F(NS,J1),J1=1,4) PRINT106 PRINT107,(KG(I),KD(I),I=5,10) DO21I=1,NSUB 21 PRINT103,I,(F(I,J1),J1=5,10) PRINT104,(F(NS,J1),J1=5,10) C COMPUTE AND PRINT ROOT MEAN SQUARE PRINT105 DO25I=1,4 SUM=0. DO24J=1,NSUB 24 SUM=SUM&RSQ(J,I)**2 RMSQ=SQRT(SUM/XNSUB) 25 PRINT103,I,RMSQ 100 FORMAT(1H1/4X,19HCORRELATION (PHASE),36X,15HF RATIO (PHASE)) 101 FORMAT(10X,2HR1,11X,2HR2,11X,2HR3,11X,2HR4,11X,2HF1,11X,2HF2,11X,2 1HF3,11X,2HF4) 102 FORMAT(3H DF,50X,4(I10,I3)/5H SUBJ) 103 FORMAT(I3,2X,8E13.4) 104 FORMAT(4H AVE1X,8E13.4) 105 FORMAT(1H0/17H0ROOT MEAN SQUARE/7H PHASE) 106 FORMAT(1H0//4X,23HF RATIO (BETWEEN PHASE)) 107 FORMAT(10X,3HF12,10X,3HF13,10X,3HF14,10X,3HF23,10X,3HF24,10X,3HF34 1/3H DF,I8,I3,5(I10,I3)/5H SUBJ) RETURN END $ FORTRAN SUBROUTINE READX(X,N,K,IRX,NL) DIMENSION X(NL,K),FMT(12) READ100,(FMT(I),I=1,12) 100 FORMAT(12A6) 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 $ FORTRAN SUBROUTINE PXPV(MDV) DIMENSION U(5,5),ROOT(5) DIMENSION D(30),S(30,32),RSQ(32,4) DIMENSION X(30,5),P(5,32),T(5,5,32),SGN(5,32),XST(30,5),PST(5) DIMENSION PAL(32),BNUB(30),IX1(5),IY1(5) DIMENSION TIT(14),FPH(5),FOD(3),FS(3),FSA(3) DIMENSION TX(250),TY(250),ITX(2),ITY(2),LVEC(5) DIMENSION IP(3),WORD(2),XM(2),KPO(50),FFPH(4) COMMON N,K,XN,XK,NS,RMAX,IPH,NSUB,IS,NV,K1 COMMON D,S,RSQ,X,P,T,SGN,XST,PST,U,ROOT COMMON PAL,BNUB,NM,IX1,IY1 COMMON /FMIC/TIT,FPH,FOD,FUU,ZRO,FS,FSA DO47I=1,K 47 CALL PIBCI(I,LVEC(I),4H(I1)) DO260I=1,K1 JJ=I&1 DO260J=JJ,K KPO(1)=0 CALL ROLL(1) CALL REFRNC(1,1024,1024,1024) CALL TEXT1(45,1024,TIT(2),0,12) IF(K1)41,42,41 41 CALL MAMI(X(1,I),X(1,J),N,XM(1),XM(2)) RANGE=XM(1)-XM(2) XSL=(RANGE* 80.)/924. ALF=XSL/SQRT(RMAX) C LABEL PLANE AND PHASE IP(1)=IPH IP(2)=I IP(3)=J CALL FORCNV(IP,3,FFPH,FPH,DUM1,DUM2) CALL TEXT1(25,1,FFPH,0,4) C COMPUTE VECTORS FOR EACH IDEAL POINT C TEST RANGE OF IDEAL POINTS SRANG=XM(2)-RANGE GRANG=XM(1)&RANGE LL=0 DO10IM=1,NS L=0 II=I 12 L=L&1 IF(P(II,IM).LT.SRANG)GOTO13 IF(P(II,IM).LE.GRANG)GOTO16 P(II,IM)=GRANG GOTO15 13 P(II,IM)=SRANG 15 LL=LL&1 KPO(LL)=IM 16 IF(L.EQ.2)GOTO10 II=J GOTO12 10 CONTINUE PRINT100,IPH,I,J,(KPO(II),II=1,LL) IF(MDV)21,83,21 83 DO81IK=1,NS TX(IK)=P(I,IK) 81 TY(IK)=P(J,IK) CALL MAMI(TX,TY,NS,TMAX,TMIN) GOTO82 21 L2=0 DO35I1=1,NS DO35I2=1,K L1=L2&1 L2=L1&1 AI=ALF*T(I,I2,I1) TX(L1)=P(I,I1)&AI AJ=ALF*T(J,I2,I1) TY(L1)=P(J,I1)&AJ TX(L2)=P(I,I1)-AI 35 TY(L2)=P(J,I1)-AJ C CONVERT TO SCALE CALL MAMI(TX,TY,NV,TMAX,TMIN) 82 IF(XM(1)-TMAX)38,39,39 38 XM(1)=TMAX 39 IF(XM(2)-TMIN)40,40,43 43 XM(2)=TMIN GOTO40 C FOR ONE DIM. 42 CALL MM(X(1,1),N,XM(1),XM(2)) DO20L=1,NS IF(P(1,L)-XM(1))26,20,27 27 XM(1)=P(1,L) 26 IF(P(1,L)-XM(2))28,20,20 28 XM(2)=P(1,L) 20 CONTINUE CALL TEXT1(45,15,FOD,0,3) C FIND ZERO POINT 40 RANGE=XM(1)-XM(2) SCALE=924./RANGE IZR=-XM(2)*SCALE&50. C DRAW SQUARE CALL DVR2(IX1,IY1,5) C LABEL X AXIS DO5II=1,2 5 CALL PIBCI(XM(II),WORD(II),6H(F6.3)) XMIN=XM(2) CALL TEXT1(50,15,WORD(2),0,1) CALL TEXT1(IZR,15,ZRO,0,1) CALL TEXT1(15,IZR,ZRO,0,1) CALL TEXT1(979,15,WORD(1),0,1) C DRAW STIMULUS POINTS IF(K1)50,51,50 51 IAY=512 50 DO22KI=1,N IAX=(X(KI,I)-XMIN)*SCALE&50. IF(K1)53,22,53 53 IAY=(X(KI,J)-XMIN)*SCALE&50. 22 CALL TEXT1(IAX,IAY,BNUB(KI),0,1) IF(K1)61,60,61 C LABEL IDEAL POINTS 61 DO25KI=1,NSUB IAX=(P(I,KI)-XMIN)*SCALE&50. IAY=(P(J,KI)-XMIN)*SCALE&50. CALL TSP1(IAX,IAY,1H*,1) IAY=IAY&7 25 CALL TEXT1(IAX,IAY,PAL(KI),0,1) GOTO130 60 NAY=505 DO125KI=1,NSUB IAX=(P(I,KI)-XMIN)*SCALE&50. CALL TSP1(IAX,IAY,1H.,1) CALL TSP1(IAX,IAY,1H0,1) NAX=IAX-7 125 CALL TEXT1(NAX,NAY,PAL(KI),0,1) C DRAW AVERAGE IDEAL POINT 130 IAX=(P(I,NS)-XMIN)*SCALE&50. IF(K1)62,63,62 62 IAY=(P(J,NS)-XMIN)*SCALE&50. 63 CALL TSP1(IAX,IAY,1HO,1) CALL TSP1(IAX,IAY,1H&,1) IF(MDV)75,260,75 75 IF(K1)49,46,49 C DRAW VECTORS AND LABEL 49 KI=0 DO70IK=1,NS 71 DO70II=1,K 45 DO48I1=1,2 KI=KI&1 ITX(I1)=(TX(KI)-XMIN)*SCALE&50. 48 ITY(I1)=(TY(KI)-XMIN)*SCALE&50. CALL DVR2(ITX,ITY,2) CALL TEXT1(ITX(1),ITY(1),LVEC(II),0,1) IF(SGN(II,IK).LT.0.) CALL TEXT1(ITX,ITY,1HO,0,1) 70 CONTINUE 260 CONTINUE 100 FORMAT(6H0PHASEI3,3X,5HPLANE2I2,3X,19HPOINTS OUT OF RANGE/(5X,30I4 1)) 46 RETURN END $ FORTRAN SUBROUTINE PTXP(IT) DIMENSION LX(2),LY(2) DIMENSION IP(3),WORD(2),XLB(2),FSS(3),FFPH(5) DIMENSION TIT(14),FPH(5),FOD(3),FS(3),FSA(3) DIMENSION U(5,5),ROOT(5) DIMENSION D(30),S(30,32),RSQ(32,4) DIMENSION X(30,5),P(5,32),T(5,5,32),SGN(5,32),XST(30,5),PST(5) DIMENSION PAL(32),BNUB(30),IX1(5),IY1(5) COMMON N,K,XN,XK,NS,RMAX,IPH,NSUB,IS,NV,K1 COMMON D,S,RSQ,X,P,T,SGN,XST,PST,U,ROOT COMMON PAL,BNUB,NM,IX1,IY1 COMMON /FMIC/TIT,FPH,FOD,FUU,ZRO,FS,FSA DO260I=1,K1 JJ=I&1 DO260J=JJ,K IP(1)=IPH IP(2)=I IP(3)=J CALL MAMI(XST(1,I),XST(1,J),N,XLB(1),XLB(2)) IF(PST(I)-PST(J))15,15,17 17 PMAX=PST(I) PMIN=PST(J) GOTO18 15 PMAX=PST(J) PMIN=PST(I) 18 IF(PMAX-XLB(1))11,11,12 12 XLB(1)=PMAX 11 IF(PMIN-XLB(2))14,13,13 14 XLB(2)=PMIN 13 RANGE=XLB(1)-XLB(2) XSCALE=924./RANGE IZR=-XLB(2)*XSCALE&50. CALL ROLL(1) CALL REFRNC(1,1024,1024,1024) CALL TEXT1(45,1024,TIT(2),0,12) C DRAW SQUARE CALL DVR2(IX1,IY1,5) C LABEL X AXIS DO5II=1,2 5 CALL PIBCI(XLB(II),WORD(II),6H(F6.3)) CALL FORCNV(IP,3,FFPH,FPH,DUM1,DUM2) CALL TEXT1(25,1,FFPH,0,5) CALL TEXT1(50,15,WORD(2),0,1) CALL TEXT1(IZR,15,ZRO,0,1) CALL TEXT1(15,IZR,ZRO,0,1) CALL TEXT1(979,15,WORD(1),0,1) IF(IS.EQ.NS)GOTO23 IP(1)=IS IP(2)=IT CALL FORCNV(IP,2,FSS,FS,DUM1,DUM2) CALL TEXT1(300,1,FSS,0,3) GOTO26 23 CALL TEXT1(450,1,FSA,0,3) C TEST NEGATIVE ROOTS 26 IF(SGN(I,IS))4,6,6 4 CALL TEXT1(1007,25,FUU,0,1) 6 IF(SGN(J,IS))7,8,8 7 CALL TEXT1(25,1007,FUU,0,1) C DRAW STIMULUS POINTS 8 DO22KI=1,N IAX=(XST(KI,I)-XLB(2))*XSCALE&50. IAY=(XST(KI,J)-XLB(2))*XSCALE&50. 22 CALL TEXT1(IAX,IAY,BNUB(KI),0,1) C DRAW IDEAL POINT 25 IPX=(PST(I)-XLB(2))*XSCALE&50. IPY=(PST(J)-XLB(2))*XSCALE&50. CALL TSP1(IPX,IPY,1HO,1) C DRAW CROSS FOR IDEAL POINT LX(1)=IPX&13 LX(2)=IPX-13 LY(1)=IPY LY(2)=IPY CALL DVR2(LX,LY,2) LX(1)=IPX LX(2)=IPX LY(1)=IPY&13 LY(2)=IPY-13 CALL DVR2(LX,LY,2) 260 CONTINUE RETURN END $ FORTRAN SUBROUTINE NORMS(S,N,XN) DIMENSION S(100) SUM=0. SSQ=0. DO10I=1,N 10 SUM=SUM&S(I) XMN=SUM/XN DO15I=1,N S(I)=S(I)-XMN 15 SSQ=SSQ&S(I)**2 SV=1./SQRT(SSQ) DO16I=1,N 16 S(I)=S(I)*SV RETURN END $ FORTRAN BLOCK DATA COMMON/FMIC/TIT,FPH,FOD,FUU,ZRO,FS,FSA DIMENSION TIT(14),FPH(5),FOD(3),FS(3),FSA(3) DATA TIT(1)/6H(79H /,TIT(14)/6H )/ DATA FPH(1)/28H(5HPHASEI2,5X,5HPLANE2X,2I1)/ DATA FOD(1)/13HONE DIMENSION/ DATA FS(1)/18H(3HSBJI3,4H IT=I3)/ DATA FSA(1)/15HAVERAGE SUBJECT/ DATA FUU/1HM/ DATA ZRO/1H0/ END $ FORTRAN SUBROUTINE MONO(DHAT,D,KC,N,LFITSW,P,LBLOCK) DIMENSION DHAT(1),D(1),KC(1),P(1) DIMENSION LBLOCK(100),LA(100),XLA(100) EQUIVALENCE (LA,XLA) C REARRANGE D IN DHAT IF(LFITSW.GE.1) GOTO14 NPASS=0 GOTO15 14 NPASS=1 DO5I=1,N LA(I)=LBLOCK(I) 5 P(I)=I 15 DO8I=1,N M=KC(I) 8 DHAT(I)=D(M) CALL MFIT(DHAT,D,LA ,N,LFITSW,NPASS,P,P2) IF(LFITSW.NE.1) GOTO20 DO21I=1,N 21 XLA(I)=D(I) CALL SORTFL CALL SORT(P(1),P(2),N,XLA,DHAT) DO12I=1,N M=KC(I) 12 P(M)=XLA(I) RETURN 20 DO10I=1,N M=KC(I) 10 P(M)=D(I) RETURN END $ FORTRAN SUBROUTINE MAMI(X,Y,N,TMAX,TMIN) DIMENSION X(1),Y(1) CALL MM(X,N,XMAX,XMIN) CALL MM(Y,N,YMAX,YMIN) IF(YMAX.GT.XMAX) GOTO10 TMAX=XMAX GOTO13 10 TMAX=YMAX 13 IF(YMIN.LT.XMIN) GOTO15 TMIN=XMIN GOTO20 15 TMIN=YMIN 20 RETURN END $ FORTRAN SUBROUTINE BMONO(SCR,LBLOCK,N2) DIMENSION SCR(1),LBLOCK(1) L=1 LL=1 L1=1 DO17I=2,N2 II=I-1 IF(SCR(I).EQ.SCR(II))GOTO19 DO18J=L1,LL 18 LBLOCK(J)=L L1=L1&L L=1 GOTO16 19 L=L&1 16 LL=LL&1 17 CONTINUE DO20J=L1,LL 20 LBLOCK(J)=L RETURN END $ FORTRAN 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 $ FORTRAN CMTINV MTINV SUBROUTINE MTINV (A,I1,I2,B,I3,I4,N,M) MTPK1115 DIMENSION A(I1,I2), B(I3,I4) MTPK1116 C MTPK1117 C THIS SUBROUTINE CALCULATES THE INVERSE OF MATRIX A BY THE GAUSS- MTPK1118 C JORDAN ELIMINATION SCHEME. ALL CALCULATIONS ARE DONE IN DOUBLE- MTPK1119 C PRECISION ARITHMETIC MTPK1120 C MTPK1121 DOUBLE PRECISION B,BMAX,BMULT MTPK1122 EQUIVALENCE (BMAX,BMULT) MTPK1123 C MTPK1124 C SINGLE TO DOUBLE PRECISION TRANSFER MTPK1125 C MTPK1126 101 DO 103 I = 1,N MTPK1127 102 DO 103 J = 1,N MTPK1128 103 B(I,J) = A(I,J) MTPK1129 C MTPK1130 C FIND THE LARGEST ELEMENT IN THE N-K BY N-K LOWER RIGHT SUBMATRIX. MTPK1133 C MTPK1134 200 DO 430 K = 1,N MTPK1131 A(K,1) = FLOAT(K) MTPK1135 A(K,2) = A(K,1) 203 BMAX = DABS(B(K,K)) MTPK1137 210 DO 219 I = K,N MTPK1138 211 DO 219 J = K,N MTPK1139 IF (BMAX - DABS(B(I,J)) ) 213,219,219 213 BMAX = DABS ( B(I,J)) MTPK1141 214 A(K,1) = FLOAT(I) MTPK1142 215 A(K,2) = FLOAT(J) MTPK1143 219 CONTINUE MTPK1145 216 IF (BMAX - 1.D-38) 801,301,301 C MTPK1146 C EXCHANGE ROWS AND COLUMNS TO PUT B(I,J) ON DIAGONAL. MTPK1147 C MTPK1148 301 I = IFIX(A(K,1)) MTPK1149 302 J = IFIX(A(K,2)) MTPK1150 313 DO 314 M = 1,N MTPK1151 BMULT = B(I,M) MTPK1152 B(I,M) = B(K,M) MTPK1153 314 B(K,M) = BMULT MTPK1155 303 DO 304 M = 1,N MTPK1156 BMULT = B(M,J) MTPK1157 B(M,J) = B(M,K) MTPK1158 304 B(M,K) = BMULT MTPK1159 C MTPK1160 C ACTUAL INVERSION STEP. MTPK1161 C BEGIN ROW ITERATION MTPK1162 C MTPK1163 401 DO 425 I = 1,N MTPK1164 C IGNORE ROW K. MTPK1165 402 IF (I - K) 405,425,405 C SET ROW MULTIPLIER. MTPK1167 405 BMULT = B(I,K) / B(K,K) C MODIFY ROW ELEMENTS. MTPK1169 410 DO 420 J = 1,N MTPK1170 C IGNORE COLUMN K MTPK1171 411 IF (J - K) 414,412,414 412 B(I,J) = -BMULT MTPK1173 413 GO TO 420 MTPK1174 414 B(I,J) = B(I,J) - B(K,J) * BMULT MTPK1175 420 CONTINUE MTPK1176 425 CONTINUE MTPK1177 C MTPK1178 C DIVIDE PIVOT ROW BY PIVOT ELEMENT MTPK1179 C MTPK1180 426 BMULT = 1.D0/ B(K,K) 427 B(K,K) = 1.D0 C ACTUAL DIVISION MTPK1183 428 DO 429 J = 1,N MTPK1184 429 B(K,J) = B(K,J) * BMULT MTPK1185 430 CONTINUE MTPK1186 C MTPK1187 C REARRANGE AND REASSEMBLE INVERSE MATRIX. MTPK1188 C MTPK1189 501 DO 512 IJK = 1,N MTPK1190 502 I = N-IJK&1 MTPK1191 503 L = IFIX(A(I,2)) MTPK1192 504 J = IFIX(A(I,1)) MTPK1193 505 DO 508 M = 1,N MTPK1194 506 BMULT = B(I,M) MTPK1195 507 B(I,M) = B(L,M) MTPK1196 508 B(L,M) = BMULT MTPK1197 509 DO 512 M = 1,N MTPK1198 510 BMULT = B(M,I) MTPK1199 511 B(M,I) = B(M,J) MTPK1200 512 B(M,J) = BMULT MTPK1201 C MTPK1202 C MOVE A INVERSE BACK INTO A MTPK1203 C MTPK1204 601 DO 603 I = 1,N MTPK1205 602 DO 603 J = 1,N MTPK1206 603 A(I,J) = B(I,J) 701 M = 1 MTPK1208 702 RETURN MTPK1209 C MTPK1210 C ERROR RETURN AT SOME TIME BMAX WAS .LT. 0.1D-38 C MTPK1212 801 M = 2 MTPK1213 802 RETURN MTPK1214 999 END MTPK1215 $ FORTRAN *EIGENJ EIGENVALUE AND EIGENVECTOR SUBROUTINE EIGJ0020 * CD600D4.009 DATE 05/05/65 EIGJ0030 SUBROUTINE EIGENJ (A,V,NN,IV,E,I3,I4) DIMENSION A(5050),V( I3, I4),E(NN) LOGICAL IV EIGJ0060 * EIGENJ FINDS THE EIGENVALUES AND EIGENVECTORS OF EIGJ0070 * A SYMMETRIC MATRIX (A) USING A MODIFIED THRESHOLD JACOBI METHOD EIGJ0080 * ONLY THE UPPER TRIANGLE IS STORED EIGJ0090 * DIMENSION A(K=M(M&1)/2),V(M,M),E(M) WHERE M=MAXIMUM ORDER EIGJ0100 * DIMENSION A(K=M(M&1)/2),V(M,M),E(M) WHERE M=MAXIMUM ORDER EIGJ0110 N=NN EIGJ0120 IND=1 EIGJ0130 NM2=N-2 EIGJ0140 AMAX=0.0 EIGJ0150 NM=N-1 EIGJ0160 IF(.NOT.IV)GO TO 5 EIGJ0170 * SET UP VECTOR MATRIX EIGJ0180 1 DO 3 I=1,N EIGJ0190 DO 2 J=1,N EIGJ0200 2 V(I,J)=0.0 EIGJ0210 3 V(I,I)=1.0 EIGJ0220 5 K=2 EIGJ0230 DO 28 I=1,NM EIGJ0240 IP=I&1 EIGJ0250 DO 27 J=IP,N EIGJ0260 Y=A(K) EIGJ0270 X=ABS(Y) EIGJ0280 IF(X-AMAX)27,27,6 EIGJ0290 6 GO TO (7,8),IND EIGJ0300 7 AMAX=X EIGJ0310 GO TO 27 EIGJ0320 * TRANSFORMATION SETUP FOLLOWS EIGJ0330 8 II=K&I-J EIGJ0340 JJ=(J*(N&N-J&3))/2-N EIGJ0350 ITEST=1 EIGJ0360 X=A(II)-A(JJ) EIGJ0370 IF(X)10,11,10 EIGJ0380 10 X=Y/X EIGJ0390 Y=.5*X EIGJ0400 IF(ABS(Y)-.414213562)14,11,11 EIGJ0410 11 C=.707106781 EIGJ0420 IF(Y)12,12,13 EIGJ0430 12 S=-C EIGJ0440 GO TO 15 EIGJ0450 13 S=C EIGJ0460 GO TO 15 EIGJ0470 14 X=Y*Y EIGJ0480 C=1.&X EIGJ0490 S=2.*Y/C EIGJ0500 C=(1.-X)/C EIGJ0510 15 X=S*S EIGJ0520 Y=C*C EIGJ0530 XY=S*C EIGJ0540 AXY=2.*A(K)*XY EIGJ0550 Z=A(II)*Y&AXY&A(JJ)*X EIGJ0560 W=A(II)*X-AXY&A(JJ)*Y EIGJ0570 A(K)=A(K)*(Y-X)&XY*(A(JJ)-A(II)) EIGJ0580 A(II)=Z EIGJ0590 A(JJ)=W EIGJ0600 IF(NM2)24,24,16 EIGJ0610 16 IT=I EIGJ0620 JT=J EIGJ0630 * TRANSFORM A EIGJ0640 DO 23 M=1,NM2 EIGJ0650 NP=N-M EIGJ0660 IF(JT-K)20,17,18 EIGJ0670 17 IT=IT&1 EIGJ0680 JT=JT&NP EIGJ0690 ITEST=2 EIGJ0700 18 IF(IT-K)20,19,19 EIGJ0710 19 IT=IT&1 EIGJ0720 JT=JT&1 EIGJ0730 ITEST=3 EIGJ0740 20 X=A(IT)*C&A(JT)*S EIGJ0750 A(JT)=A(JT)*C-A(IT)*S EIGJ0760 A(IT)=X EIGJ0770 GO TO (21,22,23),ITEST EIGJ0780 21 IT=IT&NP EIGJ0790 JT=JT&NP EIGJ0800 GO TO 23 EIGJ0810 22 IT=IT&1 EIGJ0820 JT=JT&NP-1 EIGJ0830 23 CONTINUE EIGJ0840 24 IF(.NOT.IV)GO TO 27 EIGJ0850 * TRANSFORM V EIGJ0860 25 DO 26 M=1,N EIGJ0870 X=V(M,I)*C&V(M,J)*S EIGJ0880 V(M,J)=V(M,J)*C-V(M,I)*S EIGJ0890 26 V(M,I)=X EIGJ0900 27 K=K&1 EIGJ0910 28 K=K&1 EIGJ0920 * SEQUENCE TO ADJUST THRESHOLD FOR EACH SWEEP EIGJ0930 GO TO (29,31),IND EIGJ0940 29 IF(AMAX)30,35,30 EIGJ0950 30 AT=AMAX EIGJ0960 IND=2 EIGJ0970 F=.25 EIGJ0980 AMAX=F*AT EIGJ0990 MAX=6 EIGJ1000 GO TO 5 EIGJ1010 31 MAX=MAX-1 EIGJ1020 IF(MAX)34,33,32 EIGJ1030 32 F=2.0*F*F EIGJ1040 AMAX=AT*F EIGJ1050 33 C=0.0 EIGJ1060 GO TO 5 EIGJ1070 34 IF(C)33,35,33 EIGJ1080 * COPY EIGENVALUES INTO E EIGJ1090 35 MM=1 EIGJ1100 NM=N&1 EIGJ1110 DO 36 M=1,N EIGJ1120 E(M)=A(MM) EIGJ1130 36 MM=MM&NM-M EIGJ1140 RETURN EIGJ1150 END $ FORTRAN SUBROUTINE MTEVV(X,I1,I2,D,I3,I4,N,M) DIMENSION A(5050),X(I1,I2),D(I3,I4),E(100),NS(100),MS(100) K=1 DO1J=1,N DO1I=J,N A(K)=X(I,J) 1 K=K&1 CALL EIGENJ(A,D,N,M,E,I3,I4) CALL SORTFL(-1) IF(M)4,3,4 3 CALL SORT(E(1),E(2),N) DO 50 I=1,N 50 X(I,I)=E(I) RETURN 4 DO 5 I=1,N MS(I)=I 5 NS(I)=I CALL SORT(E(1),E(2),N,MS(1)) CALL SORTFX(1) CALL SORT(MS(1),MS(2),N,NS(1)) CALL ARREV(N,NS,D,I3,I4) DO 6 I=1,N 6 X(I,I)=E(I) RETURN END $ FORTRAN SUBROUTINE ARREV(N,NS,D,I3,I4) DIMENSION NS(N),D(I3,I4) DO 10 I=1,N 1 IF(NS(I)-I)2,10,2 2 CALL XCH(N,I,NS(I),D,NS,I3,I4) GOTO1 10 CONTINUE RETURN END $ FORTRAN SUBROUTINE XCH(N,L,M,D,NS,I3,I4) DIMENSION NS(N),D(I3,I4) I=L J=M DO 1 K=1,N T=D(K,I) D(K,I)=D(K,J) 1 D(K,J)=T NS(I)=NS(J) NS(J)=J RETURN END $ FORTRAN SUBROUTINE MFIT(DIST,DHAT,LBLOCK,MM,LFITSW,NPASS,PAS1,PAS2) 403 C FIT PERFORMS WEIGHTED-LEAST-SQUARES MONOTONE REGRESSION C THIS ROUTINE FINDS THE VALUES FOR DHAT WHICH ARE MONOTONIC C AND WHICH BEST APPROXIMATE THE VALUES OF DIST, C IN THE SENSE THAT THE SUM OF THE SQUARED DEVIATIONS, C EACH ONE WEIGHTED BY THE CORRESPONDING VALUE IN WW , C IS A MINIMUM. C IT DIRECTLY USES SORT. 407 408 DIMENSION PAS1(1),PAS2(1),DIST(1),DHAT(1),LBLOCK(1) 411 C FORM FIRST APPROXIMATION TO CORRECT PARTITON 412 C IF LFITSW=0, DO PLAIN REGRESSION C IF LFITSW GT 0, DO BLOCK REGRESSION C IF LFITSW=1, USE PRIMARY APPROACH. THUS WE SORT EACH 414 C BLOCK OF EQUAL VALUES OF DATA ACCORDING TO DIST VALUES. 415 C THEN WE CREATE PARTITON INTO BLOCKS OF SIZE 1 TO START. 416 417 C IF LFITSW=2, USE SECONDARY APPROACH. THUS WE START WITH 418 C PARTITION INTO BLOCKS OF EQUAL DATA VALUES. 419 420 C WITHIN EACH BLOCK OF WHATEVER SIZE, THE FIRST DHAT VALUE 421 C GIVES THE WEIGHTED TOTAL OF THE DIST VALUES IN THAT BLOCK, 423 C SIMILARLY, WITHIN EACH BLOCK, THE FIRST LBLOCK 424 C VALUE AND THE LAST LBLOCK VALUE BOTH GIVE THE SIZE OF THE 425 C BLOCK. C BLOCKS OF SIZE 1 FORM A PARTIAL EXCEPTION. EVERYTHIN IS THE SAME C EXCEPT THAT THE SUM OF THE W VALUES IS NOT FOUND IN THE C LAST DHAT VALUE IN THE BLOCK. 427 IF(LFITSW.GT.0)GOTO100 DO90M=1,MM DHAT(M)=DIST(M) 90 LBLOCK(M)=1 GOTO500 100 MA=1 430 110 K=LBLOCK(MA) MB=MA&K-1 432 GO TO ( 200, 300 ), LFITSW 433 434 C PRIMARY APPROACH 435 436 200 IF(K-1) 9999, 220, 210 437 C IF K=1, SAVE SORTING TIME 438 210 CALL SORTFL(&1) IF(NPASS.EQ.0)CALL SORT(DIST(MA),DIST(MA&1),K) IF(NPASS.EQ.1)CALL SORT(DIST(MA),DIST(MA&1),K,PAS1(MA)) IF(NPASS.EQ.2)CALL SORT(DIST(MA),DIST(MA&1),K,PAS1(MA),PAS2(MA)) 440 C 'SORT' WILL SORT THE K ELEMENTS OF 'DIST' IN ALGEBRAIC 441 C ORDER. 442 C BECAUSE THE FINAL ARGUMENT IS ZERO, SORTING WILL BE 443 C ASCENDING OR DESCENDING ACCORDING TO THE VALUE OF SDSWIT 444 C SUPPLIED DURING THE CALL FOR SORT IN MDSCAL. 445 C THE ELEMENTS OF'IJ' AND 'DATA' WILL BE REARRANGED IN 446 C EXACTLY THE SAME ORDER AS THOSE OF 'DIST'. 447 C ALSO THE ELEMENTS OF 'WW'. 448 220 DO 230 M=MA,MB 449 DHAT(M)=DIST(M) LBLOCK(M)=1 230 CONTINUE GO TO 400 452 453 C SECONDARY APPROACH 454 455 300 CONTINUE TEMP1=0.0 458 DO 310 M=MA,MB 459 TEMP1=TEMP1&DIST(M) 310 CONTINUE DHAT(MA)=TEMP1 461 464 C PROCEED TO NEXT BOCK. QUERY. IS IT THE LAST 465 466 400 MA=MA&K 467 IF(MA-MM-1) 110, 500, 9999 468 469 C START MAIN COMPUTATION ************************************* 470 471 500 MA=1 472 510 LUD=2 473 NSATIS=0 474 520 K=LBLOCK(MA) 475 MB=MA&K-1 476 WT=FLOAT(K) 550 DAV = DHAT(MA) / WT GO TO ( 600, 700 ), LUD 478 479 C IS BLOCK DOWN-SATISFIED. IF NOT, MERGE 480 481 600 IF(MA-1) 9999, 630, 610 482 610 MBD=MA-1 483 KD=LBLOCK(MBD) 484 MAD=MBD-KD&1 485 WTD=FLOAT(KD) 617 DAVD = DHAT(MAD) / WTD IF( DAVD-DAV ) 630, 620, 620 487 620 KNEW=K&KD 488 LBLOCK(MAD)=KNEW 489 LBLOCK(MB)=KNEW 490 DTONEW = DHAT(MAD) & DHAT(MA) DHAT(MAD) = DTONEW NSATIS=0 494 MA=MAD 495 GO TO 800 496 630 NSATIS=NSATIS&1 497 GO TO 800 498 499 C IS BLOCK UP-SATISFIED. IF NOT, MERGE 500 501 700 IF(MB-MM) 710, 730, 9999 502 710 MAU=MB&1 503 KU=LBLOCK(MAU) 504 MBU=MAU&KU-1 505 WTU=FLOAT(KU) 717 DAVU = DHAT(MAU) / WTU IF( DAV-DAVU ) 730, 720, 720 507 720 KNEW=K&KU 508 LBLOCK(MA)=KNEW 509 LBLOCK(MBU)=KNEW 510 DTONEW = DHAT(MA) & DHAT(MAU) DHAT(MA) = DTONEW NSATIS=0 514 GO TO 800 515 730 NSATIS=NSATIS&1 516 GO TO 800 517 518 C PROCEED TO NEXT BLOCK IF READY. 519 520 800 LUD = 3-LUD 521 C QUERY. IS BLOCK BOTH UP AND DOWN SATISFIED. IF NOT, RETUR 522 IF(NSATIS-1) 520, 520, 810 523 C QUERY. IS THIS LAST BLOCK. IF NOT, GO ON TO NEXT BLOCK. 524 810 IF(MB-MM) 820, 900, 9999 525 820 MA=MB&1 526 GO TO 510 527 528 C MAIN COMPUTATION COMPLETE. PLACE ANSWERS IN DHAT. 529 530 900 MA=1 531 910 K=LBLOCK(MA) 532 MB=MA&K-1 533 IF(K-1) 9999, 940, 920 534 920 TEMP1=DHAT(MA)/FLOAT(K) DO 930 M=MA,MB 536 930 DHAT(M)=TEMP1 537 GO TO 945 940 DHAT(MA) = DIST(MA) 945 MA = MB & 1 IF(MA-MM-1) 910, 950, 9999 539 950 RETURN 540 541 C TROUBLE EXIT 542 543 9999 PRINT99 99 FORMAT(58H0MFIT DIAGNOSTIC. IMPOSSIBLE BRANCH TAKEN ON IF STATEMEN 1T.) STOP 546 547 END 584 $ GMAP C SET DIVIDE CHECK TO 0 SYMREF ...... SETDC SAVE LDA A,DU STCA 7,70 TRA ...... A FLD =0.,DU RET 6 END $ FORTRAN C ************************************************************************** C THE FOLLOWING ROUTINES ARE FOR MICROFILM PLOTTING C ************************************************************************** $ FORTRAN SUBROUTINE PSD(X,Y,Z,IT,TIT,FVT,FHO) DIMENSION X(1),Y(1),Z(1),TIT(1),FVT(1),FHO(1) COMMON N,K,XN,XK,NS,RMAX,IPH,NSUB,IS,NV,K1 CALL MM(X,N,XMAX,XMIN) CALL MM(Y,N,YMAX,YMIN) CALL FRAME(TIT) CALL HORIZ(FHO,IS,IPH,IT) CALL VERT(FVT) CALL LABEL(XMAX,XMIN,YMAX,YMIN) CALL DRAW(X,Z,N,1) CALL DRAW(X,Y,N,0,1H*) RETURN END $ FORTRAN SUBROUTINE PLSD(X,Y,IT,TIT,FVT,FHO,FAS) DIMENSION X(1),Y(1),TIT(1),FVT(1),FHO(1),FAS(1) COMMON N,K,XN,XK,NS,RMAX,IPH,NSUB,IS,NV,K1 CALL MM(X,N,XMAX,XMIN) CALL MM(Y,N,YMAX,YMIN) CALL FRAME(TIT) IF(IS.EQ.NS) GOTO10 CALL HORIZ(FHO,IS,IPH,IT) GOTO12 10 CALL HORIZ(FAS,IS,IPH,IT) 12 CALL VERT(FVT) CALL LABEL(XMAX,XMIN,YMAX,YMIN) CALL DRAW(X,Y,N,0,1H*) RETURN END $ FORTRAN 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 $ FORTRAN SUBROUTINE MM(X,N,XMAX,XMIN) C FINDS THE MAXIMUM AND MINIMUM OF ARRAY X 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 $ FORTRAN SUBROUTINE PJIDV(P,X,N,K,M,TIT,BNUB,PAL) DIMENSION P(5,32),X(30,5),TIT(14),BNUB(50),PAL(52) DIMENSION FFM(2),FI(2),FMXS(3),FXS(3),FPH(2) DATA FFM(1)/10H(4HDIM.I3)/ DATA FMXS(1)/14H(6HSCALE=F7.3)/ DATA FPH(1)/11H(7HPHASE 4)/ PSCALE=460. DO5I=2,K II=I-1 DO5J=1,II CALL MAXI(X(1,I),X(1,J),XMAX,N) XSCALE=460./XMAX CALL ROLL CALL REFRNC(-512,512,1024,1024) CALL TEXT1(-460,512,TIT(2),0,12) C LABEL PHASE CALL TEXT1(-460,490,FPH,0,2) C DRAW AXES CALL DVR1(0,-460,0,460,1) CALL DVR1(-460,0,460,0,1) C DRAW AND LABEL POINTS DO28IK=1,N IX=X(IK,J)*XSCALE IY=X(IK,I)*XSCALE 28 CALL TEXT1(IX,IY,BNUB(IK),0,1) C DRAW AND LABEL IDEAL VECTORS DO22IK=1,M IPX=P(J,IK)*PSCALE IPY=P(I,IK)*PSCALE CALL DVR1(0,0,IPX,IPY,1) IPY=IPY&7 IF(IK.EQ.M) GOTO25 CALL TEXT1(IPX,IPY,PAL(IK),0,1) GOTO22 25 CALL TSP1(IPX,IPY,1HO,1) CALL TSP1(IPX,IPY,1H&,1) 22 CONTINUE C LABEL DIMENSIONS CALL FORCNV(J,1,FI,FFM,DUM1,DUM2) CALL TEXT1(440,7,FI,0,2) CALL FORCNV(I,1,FI,FFM,DUM1,DUM2) CALL TEXT1(-10,467,FI,0,2) C LABEL SCALE SCALE=460./XSCALE CALL FORCNV(SCALE,1,FXS,FMXS,DUM1,DUM2) CALL TEXT1(-460,-510,FXS,0,3) 5 CONTINUE RETURN END $ GMAP * * SUBROUTINE LET FOR J.-J. CHANG 5/12/67 * LET STORES N LETTERS (A TO Z AND AA TO ZZ) LEFT ADJUSTED IN ITS * FIRST ARGUMENT. * * CALLING SEQUENCE# * * CALL LET (NAME,N) * * NAME# THE ARRAY IN WHICH THE LETTERS WILL BE STORED. * N# FULL-WORD INTEGER - SIZE OF ARRAY * LET SAVE 2,3 LDX2 0,DU LDX3 0,DU LDA 2,1 STA ADD LDA 3,1* ALS 18 STA N LD1 LDQ A,3 STOR1 ADQ =O000100000000 STQ ADD,*2 ADX2 1,DU CMPX2 N TRC RTN CMPQ TEST,3 TNC STOR1 ADX3 1,DU CMPX3 3,DU TNC LD1 LD2 LDQ A,3 STOR2 ADQ =O000101000000 STQ ADD,*2 ADX2 1,DU CMPX2 N TRC RTN CMPQ TEST,3 TNC STOR2 ADX3 1,DU CMPX3 6,DU TNC LD2 RTN RETURN LET ADD ARG ** N BSS 1 A OCT 202020202020 OCT 204020202020 OCT 206120202020 OCT 202020202020 OCT 204040202020 OCT 206161202020 TEST OCT 203120202020 OCT 205120202020 OCT 207120202020 OCT 203131202020 OCT 205151202020 OCT 207171202020 END $ FORTRAN SUBROUTINE LMIC(N,BNUB) DIMENSION BNUB(1) 40 DO130I=1,N IF(I.GE.10)GOTO43 CALL PIBCI(I,BNUB(I),4H(I1)) GOTO130 43 IF(I.GT.99)GOTO45 CALL PIBCI(I,BNUB(I),4H(I2)) GOTO130 45 CALL PIBCI(I,BNUB(I),4H(I3)) 130 CONTINUE RETURN END *