C mdprefs.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 "Individual Differences and Multidimensional Scaling" by J Douglas C Carroll in "Multidimensional Scaling: Theory and Applications in the C Behavioral Sciences", vol. 1, Theory, Seminar Press, New York and C London, 1972 C----+----@----+----@----+----@----+----@----+----@----+----@----+----@ $ FORTRAN C MDPREFS - MULTIDIMENSIONAL ANALYSIS OF PREFERENCE DATA C BY J D CARROLL AND J-J CHANG C FEBRUARY 1972 DIMENSION D( 64, 64),WT( 64, 64),X( 64,10),W(4096),BLK(6) DIMENSION TIT(12),BNUM( 64) DIMENSION FMIC(12),BIDEN( 64),SCR( 8193),VI( 64, 64) EQUIVALENCE (SCR,D),(SCR( 4097),VI) C SUBROUTINES PIBCI AND CLEAN ARE SYSTEM ROUTINES UNDER THE C MICROFILM PACKAGE. C SUBROUTINE PPS IS FOR PLOTTING PURPOSE ONLY. C NP- NUMBER OF SUBJECTS (MAX.=64) C NS- NUMBER OF STIMULI (MAX.=64) C NF- NUMBER OF FACTORS (MAX=10) C NFP- NUMBER OF FACTORS TO BE PLOTTED. C IREAD- 0, READ IN PAIRED-COMPARISON DATA. C IREAD- 1, READ IN FSM. C IREAD- 2, READ IN FSM AND NORMALIZE EACH ROW. C MDATA- 0, NO MISSING DATA. C MDATA- 1, MISSING DATA. C MDATA- 2, MISSING DATA WITH WEIGHTS TO BE READ IN. C MDATA- 3, MISSING DATA AND NON-EMPTY CELLS ARE IN BLOCK PATTERN. C NS1- IF MDATA=1 OR 2, INDICATE THE NUMBER OF SUBJECTS HAVING THE C SAME MISSING PATTERN. C NPUNCH- 1, PUNCH SOLUTION ON CARDS. C NPUNCH- 0, NO PUNCH. C IPUNF- 1, PUNCH NORMALIZED FSM ON CARDS. C IPUNF- 0, NO PUNCH. C MIDEN- 1, READ IN STIMULUS LABELS ON PLOTS. C MIDEN- 0, DO NOT READ IN LABELS. C LP- 0, LABEL PERSON VECTORS ON PLOTS. C LP- 1, NO LABEL. C KV- 0, DRAW VECTORS TO REPRESENT SUBJECTS. C KV- 1, DRAW ASTERISKS. C NORP- 1, NORMALIZE SUBJECTS VECTORS. C NORP- 0, DO NOT NORMALIZE SUBJECTS VECTORS. DATA BLK/6*6H / 60 READ303,(TIT(I),I=1,12) DO10I=1,6 IF(TIT(I)-BLK)11,10,11 10 CONTINUE GOTO50 11 READ300,NP,NS,NF,NFP,IREAD,MDATA,NS1,NPUNCH,IPUNF,MIDEN,LP,KV,NORP C GENERATE NUMBERS IN BCI FORMAT FOR LABELLING MICROFILM PLOTS IF(LP.EQ.1) GOTO25 DO186IB=1,NP 186 CALL PIBCI(IB,BNUM(IB),4H(I2)) 25 IF(MIDEN.EQ.0)GOTO27 READ303,(FMIC(I),I=1,12) READFMIC,(BIDEN(I),I=1,NS) GOTO26 27 DO30I=1,NS IF(I.LT.10)GOTO31 32 CALL PIBCI(I,BIDEN(I),4H(I2)) GOTO30 31 CALL PIBCI(I,BIDEN(I),4H(I1)) 30 CONTINUE C PRINT TITLE AND OPTIONS 26 PRINT301,(TIT(I),I=1,12) PRINT302,NP,NS,NF,NFP,IREAD,MDATA,NS1,NPUNCH,IPUNF,MIDEN,LP,KV,NOR 1P IF(IREAD.GT.0) GOTO42 C COMPUTE FIRST SCORE MATRIX 41 CALL COMPW(NP,NS,D,W,WT,MDATA,IPUNF,NS1,SCR,VI) GOTO51 C READ IN FIRST SCORE MATRIX 42 CALL RFSM(W,NP,NS,IREAD,IPUNF) C COMPUTE SECOND SCORE MATRIX 51 CALL COMPR(NP,NS,NF,W,D,X,WT,NPUNCH,NORP) C PLOT SUBJECTS AND STIMULI CALL PPS(X,WT,NP,NS,NFP,TIT,BNUM,KV,LP,BIDEN) GOTO60 300 FORMAT(18I4) 301 FORMAT(1H112A6) 302 FORMAT(72H0 NP NS NF NFP IREAD MDATA NS1 NPUNCH IPUNF MID 1EN LP KV NORP/4I4,I6,I7,I6,3I7,I6,2I4) 303 FORMAT(12A6) 50 CALL CLEAN STOP END $ FORTRAN SUBROUTINE COMPR(NP,NS,NF,W,R,X,Y,NPUNCH,NORP) C 1/26/72 C CALLS SUBROUTINE CPV, WHICH COMPUTES PROPORTION OF VARIANCE. C COMPUTES CORRELATION BETWEEN FIRST AND SECOND SCORE MATRIX DIMENSION W(NP,NS),R( 64, 64),Y( 64, 64),X(NP,NF),COR(64) XNS=NS C COMPUTE CROSS PRODUCT MATRIX OF SUBJECTS Y=W(WT) NL=64 DO41I=1,NP DO41J=I,NP Y(I,J)=0. DO43L=1,NS 43 Y(I,J)=Y(I,J)&W(I,L)*W(J,L) 41 Y(J,I)=Y(I,J) C COMPUTE CORRELATION MATRIX OF SUBJECTS AND PRINT OUT. DO44I=2,NP II=I-1 DO44J=1,II R(I,J)=Y(I,J)/SQRT(Y(I,I)*Y(J,J)) R(J,I)=R(I,J) 44 CONTINUE DO45I=1,NP 45 R(I,I)=1. PRINT111 DO46I=1,NP 46 PRINT104,I,(R(I,J),J=1,NP) C FIND ROOTS AND VECTORS M=1 PRINT311 IF(NS.LE.NP) GOTO61 C THE PURPOSE OF FINDING ROOTS AND VECTORS OF THE LESSER ORDER IS TO C SAVE COMPUTER TIME. C FIND ROOTS AND VECTORS OF Y, IF NP.LT.NS. CALL MTEVV(Y, 64, 64,R, 64, 64,NP,M) PRINT312,(Y(I,I),I=1,NP) CALL CPV(Y,NL,NP) DO47I=1,NF 47 Y(I,I)=SQRT(Y(I,I)) C COMPUTE POPULATION MATRIX X DO62J=1,NP SUMQ=0. DO67I=1,NF X(J,I)=R(J,I)*Y(I,I) SUMQ=SUMQ&X(J,I)**2 67 R(J,I)=R(J,I)/Y(I,I) SUMQ=SQRT(SUMQ) IF(NORP.EQ.0) GOTO62 DO66I=1,NF 66 X(J,I)=X(J,I)/SUMQ 62 CONTINUE C COMPUTE STIMULUS MATRIX Y DO63J=1,NF SUM=0. DO68I=1,NS Y(I,J)=0. DO64M=1,NP 64 Y(I,J)=Y(I,J)&W(M,I)*R(M,J) 68 SUM=SUM&Y(I,J)**2 SUM=SQRT(SUM) C NORMALIZE Y DO65L=1,NS 65 Y(L,J)=Y(L,J)/SUM 63 CONTINUE GOTO28 C COMPUTE CROSS PRODUCT MATRIX OF STIMULI R=(WT)W 61 DO20I=1,NS DO20J=I,NS R(I,J)=0. DO21K=1,NP 21 R(I,J)=R(I,J)&W(K,I)*W(K,J) 20 R(J,I)=R(I,J) C FIND ROOTS AND VECTORS OF R, IF NS.LE.NP. CALL MTEVV(R, 64, 64,Y, 64, 64,NS,M) C COMPUTE POPULATION MATRIX X PRINT312,(R(I,I),I=1,NS) CALL CPV(R,NL,NS) DO30I=1,NP SUMQ=0. DO25K=1,NF X(I,K)=0. DO22J=1,NS 22 X(I,K)=X(I,K)&W(I,J)*Y(J,K) 25 SUMQ=SUMQ&X(I,K)**2 SUMQ=SQRT(SUMQ) IF(NORP.EQ.0) GOTO30 DO29K=1,NF 29 X(I,K)=X(I,K)/SUMQ 30 CONTINUE C COMPUTE SECOND SCORE MATRIX AND STORE IN W C COMPUTE CORRELATION BETWEEN FIRST AND SECOND SCORE MATRIX 28 DO329I=1,NP SX=0. SY=0. SQX=0. SQY=0. SXY=0. DO328J=1,NS SX=SX&W(I,J) SQX=SQX&W(I,J)**2 SUM=0. DO327K=1,NF 327 SUM=SUM&X(I,K)*Y(J,K) SXY=SXY&W(I,J)*SUM W(I,J) SUM SY=SY&W(I,J) 328 SQY=SQY&W(I,J)**2 329 COR(I)=(XNS*SXY-SX*SY)/SQRT((XNS*SQX-SX**2)*(XNS*SQY-SY**2)) PRINT300 DO35I=1,NP 35 PRINT104,I,(W(I,J),J=1,NS) PRINT101 DO330I=1,NP 330 PRINT102,I,COR(I) C PRINT AND PUNCH X AND Y PRINT103 DO53I=1,NF 53 PRINT104,I,(X(J,I),J=1,NP) PRINT105 DO54I=1,NF 54 PRINT104,I,(Y(J,I),J=1,NS) IF(NPUNCH)51,50,51 51 PUNCH106 DO57I=1,NF 57 PUNCH107,I,(Y(J,I),J=1,NS) PUNCH108 DO58I=1,NF 58 PUNCH107,I,(X(J,I),J=1,NP) 101 FORMAT(50H0CORRELATION BETWEEN FIRST AND SECOND SCORE MATRIX/8H SU 1BJECT,7X,1HR) 102 FORMAT(I5,5X,F10.6) 103 FORMAT(18H0POPULATION MATRIX/7H0FACTOR) 104 FORMAT(1H I4,10E12.4/(5X,10E12.4)) 105 FORMAT(16H0STIMULUS MATRIX/7H0FACTOR) 106 FORMAT(15HSTIMULUS MATRIX/11H(2X,5E13.5)) 107 FORMAT(I2,5E13.5/(2X,5E13.5)) 108 FORMAT(17HPOPULATION MATRIX/11H(2X,5E13.5)) 111 FORMAT(31H0CORRELATION MATRIX OF SUBJECTS) 300 FORMAT(20H0SECOND SCORE MATRIX/8H0SUBJECT5X,8HSTIMULUS) 311 FORMAT(32H0ROOTS OF THE FIRST SCORE MATRIX) 312 FORMAT(3X,10E12.4) 50 RETURN END $ FORTRAN SUBROUTINE COMPW(NP,NS,D,W,WT,MDATA,IPUNF,NS1,SCR,VI) DIMENSION D( 64, 64),W(NP,NS),FMT(12),CODE(4),Y( 64),WT( 64, 64) DIMENSION SCR( 8193),VI( 64, 64) XM=NS READ100,(FMT(I),I=1,12) READFMT,(CODE(I),I=1,4) DO18L=1,NP C READ IN MATRIX D AND CONVERT IF(MDATA.GT.0) GOTO11 CALL READS(D,FMT,CODE,NS) GOTO406 C READ IN DATA AND WEIGHT MATRIX 11 GOTO(12,13,12),MDATA GOTO17 13 CALL READW(D,WT,FMT,CODE,NS) C READ IN DATA AND GENERATE WEIGHT MATRIX FOR MISSING DATA. 12 CALL READMW(D,WT,FMT,CODE,NS) IF(MDATA.EQ.3) GOTO19 17 CALL LSQMD(D,WT,Y,NS,NP,W,L,NS1,SCR,VI) GOTO18 C MISSING DATA SPECIAL BLOCK CONSTRUCTION 19 CALL MDN(Y,WT,NS) C COMPUTE FIRST SCORE MATRIX 406 SWS=0. DO10I=1,NS SR=0. SC=0. 9 DO7J=1,NS SR=SR&D(I,J) 7 SC=SC&D(J,I) 5 IF(MDATA.EQ.3)GOTO8 W(L,I)=SR-SC GOTO10 8 W(L,I)=(SR-SC)/Y(I) 10 SWS=SWS&W(L,I) WMEAN=SWS/XM SUM=0. SUMN=0. DO15I=1,NS W(L,I)=W(L,I)-WMEAN IF(MDATA.EQ.0)GOTO15 WSQ=W(L,I)**2 SUMN=SUMN&WSQ*Y(I) SUM=SUM&WSQ 15 CONTINUE IF(MDATA.EQ.0)GOTO18 C=SQRT(SUMN/SUM) DO16I=1,NS 16 W(L,I)=C*W(L,I) 18 CONTINUE C PRINT W PRINT 101 DO20I=1,NP 20 PRINT102,I,(W(I,J),J=1,NS) IF(IPUNF)22,25,22 C PUNCH FIRST SCORE MATRIX 22 PUNCH 106 DO21I=1,NP 21 PUNCH107,I,(W(I,J),J=1,NS) 100 FORMAT(12A6) 101 FORMAT(19H0FIRST SCORE MATRIX/8H0SUBJECT5X,8HSTIMULUS) 102 FORMAT(1H I4,10E12.4/(5X,10E12.4)) 103 FORMAT(3H C=3E13.4) 106 FORMAT(18HFIRST SCORE MATRIX/11H(2X,5E13.5)) 107 FORMAT(I2,5E13.5/(2X,5E13.5)) 25 RETURN END $ FORTRAN SUBROUTINE LSQMD(D,WT,Y,NS,NP,W,L,NS1,SCR,VI) DIMENSION D( 64, 64),WT( 64, 64),W(NP,NS),Y( 64) DIMENSION VI( 64, 64),SCR( 8193) C BEFORE MAKING ANY CHANGES IN THIS ROUTINE, CHECK THE EQUIVALENCE C STATEMENT IN MAIN PROGRAM INVOLVING SCR, D AND VI. C DATA IN D AND WEIGHT IN WT PRINT100,L BIG=10.**8 C COMPUTE TOTAL VARIANCE TVAR=0. DO2I=1,NS Y(I)=0. DO3J=1,NS TVAR=TVAR&WT(I,J)*D(I,J)**2 3 Y(I)=Y(I)&WT(I,J)*D(I,J)-WT(J,I)*D(J,I) 2 CONTINUE IF(L.EQ.1.OR.L.GT.NS1) GOTO50 GOTO52 C COMPUTE U AND PUT IN D 50 DO5I=1,NS D(I,I)=0. DO5J=1,NS IF(I.EQ.J)GOTO5 D(I,J)=-(WT(I,J)&WT(J,I)) D(I,I)=D(I,I)-D(I,J) 5 CONTINUE DO8I=1,NS DO8J=1,I WT(I,J)=ABS(D(I,J)) IF(D(I,J).EQ.0.) WT(I,J)=BIG 8 WT(J,I)=WT(I,J) C PARTITION STIMULI IN BLOCKS CALL BLOCK(WT,NS) 20 DO21I=1,NS DO21J=1,NS IF(I.EQ.J) GOTO22 IF(WT(I,J).GE.BIG) GOTO24 22 WT(I,J)=1. GOTO21 24 WT(I,J)=0. 21 CONTINUE C COMPUTE V DO25I=1,NS DO25J=1,NS 25 WT(I,J)=WT(I,J)&D(I,J) MM=1 C ON RETURN THE INVERSE IS STORED IN WT (ORDER NS). SCR IS USED AS C SCRATCH. CALL MTINV(WT, 64, 64,SCR, 64, 64,NS,MM) IF(MM.EQ.1) GOTO53 C ON RETURN IF MM IS NOT 1 INDICATES WT IS ILL CONDITIONED. C PROGRAM PRINTS OUT A MESSAGE AND STOPS. PRINT111,MM STOP 53 DO51I=1,NS DO51J=1,NS 51 VI(I,J)=WT(I,J) 52 DO26I=1,NS W(L,I)=0. DO26J=1,NS 26 W(L,I)=W(L,I)&Y(J)*VI(J,I) PRINT109 PRINT102,I,(Y(I),I=1,NS) C COMPUTE VARIANCE ACCOUNTED FOR VAF=0. DEN=0. DO27I=1,NS VAF=VAF&W(L,I)*Y(I) 27 DEN=DEN&W(L,I)**2 PVAF=VAF/TVAR C NORMALIZE W C=SQRT(VAF/DEN) DO28I=1,NS 28 W(L,I)=W(L,I)*C PRINT107,TVAR,VAF,PVAF,C 100 FORMAT(8H0SUBJECTI4) 102 FORMAT(I4,10E12.4/(4X,10E12.4)) 107 FORMAT(6H0TVAR=E12.4,3X,4HVAF=E12.4,3X,5HPVAF=E12.4,3X,2HC=E12.4) 109 FORMAT(9H0Y VECTOR) 111 FORMAT(4H0MM=I1,21HWT CANNOT BE INVERTED) RETURN END $ FORTRAN SUBROUTINE CPV(R,NL,N) DIMENSION R(NL,NL),P(64) SUM=0. DO2I=1,N 2 SUM=SUM&R(I,I) DO3I=1,N 3 P(I)=R(I,I)/SUM PRINT100,(I,I=1,N) PRINT101,(P(I),I=1,N) SUM=P(1) DO4I=2,N SUM=SUM&P(I) 4 P(I)=SUM PRINT102,(I,I=1,N) PRINT101,(P(I),I=1,N) 100 FORMAT(52H0PROPORTION OF VARIANCE ACCOUNTED FOR BY EACH FACTOR/(I9 1,9I12)) 101 FORMAT(2X,10E12.4) 102 FORMAT(48H0CUMULATIVE PROPORTION OF VARIANCE ACCOUNTED FOR/(I9,9I1 12)) RETURN END $ FORTRAN SUBROUTINE MDN(Y,WT,NS) DIMENSION Y( 64),WT( 64, 64) DO6I=1,NS SUM=0. SSQ=0. SDSQ=0. DO5J=1,NS IF(I.EQ.J)GOTO5 SUM=SUM&WT(I,J)&WT(J,I) SDSQ=SDSQ&(WT(J,I)-WT(I,J))**2 SSQ=SSQ&WT(J,I)**2&WT(I,J)**2 5 CONTINUE IF(SSQ)4,7,4 7 Y(I)=2. GOTO6 4 Y(I)=2.&SUM-SQRT(SDSQ/SSQ) 6 CONTINUE PRINT100 PRINT101,(Y(I),I=1,NS) 100 FORMAT(9H0VECTOR N) 101 FORMAT(2X,10E12.4) RETURN END $ FORTRAN SUBROUTINE RFSM(W,NP,NS,IREAD,IPUNF) DIMENSION W(NP,NS),WMEAN( 64),SD( 64),FMT(12) XM=NS DEN=1./XM READ100,(FMT(I),I=1,12) DO5I=1,NP 5 READFMT,(W(I,J),J=1,NS) C COMPUTE MEAN FOR EACH SUBJECT DO6I=1,NP SUM=0. DO4J=1,NS 4 SUM=SUM&W(I,J) WMEAN(I)=SUM/XM C NORMALIZE W - SUBTRACT THE ROW MEAN. SUM=0. DO7J=1,NS W(I,J)=W(I,J)-WMEAN(I) IF(IREAD.EQ.2)SUM=SUM&W(I,J)**2 7 CONTINUE IF(IREAD.EQ.1)GOTO6 C IREAD=2, NORMALIZED W - DIVIDE BY THE ROW SD SD(I)=SQRT(SUM*DEN) DO8J=1,NS 8 W(I,J)=W(I,J)/SD(I) 6 CONTINUE C PRINT MEAN PRINT103 PRINT104,(WMEAN(I),I=1,NP) IF(IREAD.EQ.1)GOTO10 PRINT105 PRINT104,(SD(I),I=1,NP) C PRINT W 10 PRINT101 DO20I=1,NP 20 PRINT102,I,(W(I,J),J=1,NS) IF(IPUNF)22,25,22 C PUNCH FIRST SCORE MATRIX 22 PUNCH 106 DO21I=1,NP 21 PUNCH107,I,(W(I,J),J=1,NS) 100 FORMAT(12A6) 101 FORMAT(19H0FIRST SCORE MATRIX/8H0SUBJECT5X,8HSTIMULUS) 102 FORMAT(1H I4,10E12.4/(5X,10E12.4)) 103 FORMAT(37H0MEAN OF THE RAW SCORES (BY SUBJECT)) 104 FORMAT(5X,10E12.4) 105 FORMAT(35H0SD OF THE RAW SCORES (BY SUBJECT)) 106 FORMAT(18HFIRST SCORE MATRIX/11H(2X,5E13.5)) 107 FORMAT(I2,5E13.5/(2X,5E13.5)) 25 RETURN END $ FORTRAN SUBROUTINE BLOCK(D,N) DIMENSION D( 64, 64) C COMPUTE MINIMUM PATH DO5I=1,N DO5J=2,N IF(I.EQ.J)GOTO5 JJ=J-1 DO4K=1,JJ IF(I.EQ.K)GOTO4 DSUM=D(I,J)&D(I,K) IF(D(J,K).GT.DSUM) D(J,K)=DSUM D(K,J)=D(J,K) 4 CONTINUE 5 CONTINUE RETURN END $ FORTRAN SUBROUTINE READS(D,FMT,CODE,NS) C READ IN DATA MATRIX WITHOUT MISSING DATA DIMENSION D( 64, 64),FMT(1),CODE(1),Y( 64),WT( 64, 64) N=1 C READ IN DATA MATRIX 70 DO15I=1,NS 15 READFMT,(D(I,J),J=1,NS) GOTO(1,2,3),N C NO MISSING DATA 1 DO50I=1,NS DO50J=1,NS IF(I.EQ.J)GOTO45 DO32K=1,4 IF(D(I,J)-CODE(K))32,41,32 41 GOTO(10,20,45,45),K 10 D(I,J)=1. GOTO50 20 D(I,J)=-1. GOTO50 32 CONTINUE 45 D(I,J)=0. 50 CONTINUE RETURN ENTRY READMW(D,WT,FMT,CODE,NS) C MISSING DATA , WEIGHT IS EITHER 0 OR 1 N=2 GOTO70 2 DO55I=1,NS DO55J=1,NS IF(I.EQ.J)GOTO55 IF(D(I,J).EQ.CODE(4))GOTO53 WT(I,J)=1. GOTO55 53 WT(I,J)=0. 55 CONTINUE GOTO1 ENTRY READW(D,WT,FMT,CODE,NS) C READ IN WEIGHT MATRIX N=3 GOTO70 3 DO62I=1,NS 62 READ100,(WT(I,J),J=1,NS) 100 FORMAT(10F7.2) DO65I=1,NS DO65J=1,NS IF(I.EQ.J) GOTO65 IF(D(I,J).EQ.CODE(4)) WT(I,J)=0. 65 CONTINUE GOTO1 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) C THE FOLLOWING SYSTEM SORT PACKAGE ROUTINES SORT THE FLOATING POINT C ARRAY E OF N ELEMENTS IN DESCENDING ORDER. 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 C THE FOLLOWING ROUTINE SORTS THE FLOATING POINT ARRAY E OF N C ELEMENTS IN DESCENDING ORDER AND ARRAY MS IS THE PASSIVE LIST C WHICH WILL BE REARRANGED IN THE SAME ORDER AS E. CALL SORT(E(1),E(2),N,MS(1)) C THE FOLLOWING ROUTINES SORT THE FIXED POINT ARRAY MS OF N ELEMENTS C IN ASCENDING ORDER WITH ARRAY NS AS THE PASSIVE LIST. 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 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 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 PPS(X,Y,NP,NS,NF,TIT,BNUM,KV,LP,BIDEN) C PROGRAM TO PLOT NS STIMULI(Y) AND NP (X) PEOPLE IN NF DIMENSIONS. C ALL SUBROUTINES EXCEPT MAXI AND ZEROS ARE SYSTEM SUBROUTINES UNDER C MICROFILM PACKAGE. C SUBROUTINE ZEROS(A,N) SETS THE SINGLE ARRAYED VARIABLE A OF SIZE N C TO ZERO. C SUBROUTINE MAXI(X,Y,XMAX,N) FINDS THE MAGNITUDE OF ARRAYS X AND Y C IN ORDER TO SCALE THE PLOT OF PROJECTIONS OF POINTS ON XY PLANE. DIMENSION FDM(4),FFDM(3),FFM(2),FI(2),FJ(2) DIMENSION X(NP,NF),Y( 64, 64),TIT(12),FXS(4) DIMENSION FMXS(5),FMYS(5),FYS(5),FLXY(4),BNUM( 64),BIDEN( 64) DATA FMXS(1)/26H(17HSCALE FOR PEOPLE=F7.2)/ DATA FMYS(1)/28H(19HSCALE FOR STIMULUS=F7.2)/ DATA FLXY(1)/22H* OR VECTOR FOR PEOPLE/ DATA FDM(1)/20H(I2,12H DIMENSIONAL)/ DATA FFM(1)/10H(4HDIM.I3)/ C PLOTTING AREA IS 460 BY 460, ORIGIN AT CENTER OF SQUARE C LABEL NF DIMENSIONAL SOLUTION CALL FORCNV(NF,1,FFDM,FDM,DUM1,DUM2) NNF=NF-1 DO260I=1,NNF JJ=I&1 DO260J=JJ,NF IF(NF.GT.1)GOTO35 CALL ZEROS(X(1,2),NP) CALL ZEROS(Y(1,2),NS) 35 CALL MAXI(X(1,I),X(1,J),XMAX,NP) XSCALE=460./XMAX CALL MAXI(Y(1,I),Y(1,J),YMAX,NS) YSCALE=460./YMAX CALL ROLL CALL REFRNC(-512,512,1024,1024) C PRINT TITLE CALL TEXT(-460,512,TIT,0,12) CALL TEXT(-460,460,FFDM,0,3) C DRAW PEOPLE DO31K=1,NP IAX=X(K,I)*XSCALE IAY=X(K,J)*XSCALE IF(KV.EQ.1)GOTO34 CALL DVR1(0,0,IAX,IAY,1) GOTO33 34 CALL TSP(IAX,IAY,1H*,1) C LABEL PEOPLE 33 IF(LP.EQ.1)GOTO31 NAX=IAX-7 NAY=IAY&7 CALL TEXT1(NAX,NAY,BNUM(K),0,1) 31 CONTINUE C DRAW POINTS DO22K=1,NS IAX=Y(K,I)*YSCALE IAY=Y(K,J)*YSCALE CALL TEXT1(IAX,IAY,BIDEN(K),0,1) 22 CONTINUE C LABEL SCALE SCALE=460./XSCALE CALL FORCNV(SCALE,1,FXS,FMXS,DUM1,DUM2) CALL TEXT(-460,-490,FXS,0,4) SCALE=460./YSCALE CALL FORCNV(SCALE,1,FYS,FMYS,DUM1,DUM2) CALL TEXT(-460,-510,FYS,0,5) CALL TEXT(0,-490,FLXY,0,4) IF(NF.EQ.1)GOTO50 C DRAW AXES CALL DVR1(-460,0,460,0,1) CALL DVR1(0,-460,0,460,1) C LABEL AXES CALL FORCNV(J,1,FJ,FFM,DUM1,DUM2) CALL TEXT(-10,467,FJ,0,2) CALL FORCNV(I,1,FI,FFM,DUM1,DUM2) CALL TEXT(440,7,FI,0,2) 260 CONTINUE 50 RETURN END