C mdpref.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----+----@----+----@----+----@----+----@----+----@----+----@----+----@ C MDPREF - MULTIDIMENSIONAL ANALYSIS OF PREFERENCE DATA MDPF 10 C BY J D CARROLL AND J-J CHANG MDPF 20 C MAY 17,1968 MDPF 30 REAL*8 SCR(3600) MDPF 40 DIMENSION W(60,60) MDPF 50 DIMENSION TIT(18) , BLK(6) MDPF 60 DIMENSION D(60,60) MDPF 70 DIMENSION WT(60,60) MDPF 80 DIMENSION VI(60,60) MDPF 90 REAL*4 A(120),B(120) INTEGER*4 SP DIMENSION X(60,60) MDPF 100 DIMENSION SPT(3600) MDPF 110 EQUIVALENCE (VI(1,1),SPT(1)) MDPF 120 COMMON W MDPF 130 DATA BBK/' '/ MDPF 140 DO 5 I=1,6 MDPF 150 5 BLK(I) = BBK MDPF 151 60 READ303,(TIT(I),I=1,18) MDPF 160 DO 10 I=1,6 MDPF 170 IF(TIT(I)-BLK(I))11,10,11 MDPF 180 10 CONTINUE MDPF 190 GO TO 50 MDPF 200 11 READ 300,NP,NS,NF,NFP,IREAD,MDATA,NS1,NPUNCH,IPUNF MDPF 210 C PRINT TITLE AND OPTIONS MDPF 220 26 PRINT 301,(TIT(I),I=1,18) MDPF 230 PRINT 302,NP,NS,NF,NFP,IREAD,MDATA,NS1,NPUNCH ,IPUNF MDPF 240 IF(IREAD.GT.0) GO TO 42 MDPF 250 C COMPUTE FIRST SCORE MATRIX MDPF 260 41 CALL COMPW(NP,NS,D,W,WT,MDATA,IPUNF,NS1,SCR,VI) MDPF 270 GO TO 51 MDPF 280 C READ IN FIRST SCORE MATRIX MDPF 290 42 CALL RFSM(W,NP,NS,IREAD,IPUNF) MDPF 300 C COMPUTE SECOND SCORE MATRIX MDPF 310 51 CALL COMPR(NP,NS,NF,W,D,X,WT,NPUNCH) MDPF 320 C MDPF 330 C HERE IS WHERE PLOTTING ROUTINE WILL GO MDPF 340 C PEOPLE PLOT MDPF 350 WRITE(6,400) MDPF 360 PTWO=+1.5 MDPF 380 AMTWO=-1.5 MDPF 390 PONE = +1 MDPF 400 YMAX = +1. MDPF 410 YMIN = -1. MDPF 420 DO 410 I=1,NP A(I) = X(I,1) 410 B(I) = X(I,2) CALL PLOTX(A,B,NP) 400 FORMAT('1',10X,'PLOT OF SUBJECTS IN FIRST TWO DIMENSIONS ') WRITE(6,401) 401 FORMAT('1',10X,'PLOT OF STIMULUS POINTS IN FIRST TWO DIMENSIONS') DO 403 I=1,NS A(I) = WT(I,1) 403 B(I) = WT(I,2) CALL PLOTX(A,B,NS) WRITE(6,402) 402 FORMAT('1',10X,'PLOT OF STIMULI AND SUBJECTS IN FIRST TWO DIMENSIO XNS ') SP = NP+NS DO 417 I=1,NS A(I) = WT(I,1) 417 B(I) = WT(I,2) DO 404 I=1,NP J=I+NS A(J) = X(I,1) 404 B(J) = X(I,2) CALL PLOTX(A,B,SP) WRITE(6,405) WRITE(6,405) NS,NP 405 FORMAT(//3X,'FIRST',1X,I3,1X,'POINTS ARE STIMULI AND NEXT',1X,I3,1 XX,'POINTS ARE SUBJECTS '/) GO TO 60 MDPF 590 300 FORMAT(18I4) MDPF 600 301 FORMAT(1H1,18A4) MDPF 610 302 FORMAT(66H0 NP NS NF NFP IREAD MDATA NS1 NPUNCH IPUNF MIDMDPF 620 1EN LP KV/4I4,I6,I7,I6,3I7,I6,I4) MDPF 630 303 FORMAT(18A4) MDPF 640 50 CALL EXIT MDPF 650 STOP MDPF 660 END MDPF 670 SUBROUTINE COMPW(NP,NS,D,W,WT,MDATA,IPUNF,NS1,SCR,VI) MDPF 680 REAL*8 SCR(3600) MDPF 690 DIMENSION D(60,60) MDPF 700 DIMENSION WT(60,60) MDPF 710 DIMENSION W(60,60) MDPF 720 DIMENSION FMT(18),CODE(4) MDPF 730 DIMENSION Y(3600) MDPF 740 DIMENSION VI(60,60) MDPF 750 READ 100,(FMT(I),I=1,18) MDPF 760 XM=NS MDPF 770 READFMT,(CODE(I),I=1,4) MDPF 780 DO 18 L=1,NP MDPF 790 C READ IN MATRIX D AND CONVERT MDPF 800 IF(MDATA.GT.0) GO TO 11 MDPF 810 CALL READS(D,FMT,CODE,NS) MDPF 820 GO TO 406 MDPF 830 C READ IN DATA AND WEIGHT MATRIX MDPF 840 11 GO TO (12,13,12),MDATA MDPF 850 13 CALL READW(D,WT,FMT,CODE,NS) MDPF 860 GO TO 17 MDPF 870 C READ IN DATA AND GENERATE WEIGHT MATRIX FOR MISSING DATA MDPF 880 12 CALL READMW(D,WT,FMT,CODE,NS) MDPF 890 IF(MDATA.EQ.3) GO TO 19 MDPF 900 17 CALL LSQMD(D,WT,Y,NS,NP,W,L,NS1,SCR,VI) MDPF 910 GO TO 18 MDPF 920 C MISSING DATA SPECIAL BLOCK CONSTRUCTION MDPF 930 19 CALL MDN(Y,WT,NS) MDPF 940 C COMPUTE FIRST SCORE MATRIX MDPF 950 406 SWS=0. MDPF 960 DO 10 I=1,NS MDPF 970 SR=0. MDPF 980 SC=0. MDPF 990 9 DO 7 J=1,NS MDPF1000 SR=SR+D(I,J) MDPF1010 7 SC=SC+D(J,I) MDPF1020 5 IF(MDATA.EQ.3) GO TO 8 MDPF1030 W(L,I)=SR-SC MDPF1040 GO TO 10 MDPF1050 8 W(L,I)=(SR-SC)/Y(I) MDPF1060 10 SWS=SWS+W(L,I) MDPF1070 WMEAN=SWS/XM MDPF1080 SUM=0. MDPF1090 SUMN=0. MDPF1100 DO 15 I=1,NS MDPF1110 W(L,I)=W(L,I)-WMEAN MDPF1120 IF(MDATA.EQ.0) GO TO 15 MDPF1130 WSQ=W(L,I)**2 MDPF1140 SUMN=SUMN+WSQ*Y(I) MDPF1150 SUM=SUM+WSQ MDPF1160 15 CONTINUE MDPF1170 IF(MDATA.EQ.0) GO TO 18 MDPF1180 C=SQRT(SUMN/SUM) MDPF1190 DO 16 I=1,NS MDPF1200 16 W(L,I)=C*W(L,I) MDPF1210 18 CONTINUE MDPF1220 C PRINT W MDPF1230 PRINT 101 MDPF1240 DO 20 I=1,NP MDPF1250 20 PRINT 102,I,(W(I,J),J=1,NS) MDPF1260 IF(IPUNF)22,25,22 MDPF1270 C PUNCH FIRST SCORE MATRIX MDPF1280 22 PUNCH 106 MDPF1290 DO 21 I=1,NP MDPF1300 21 PUNCH 107,I,(W(I,J),J=1,NS) MDPF1310 100 FORMAT(18A4) MDPF1320 101 FORMAT(19H0FIRST SCORE MATRIX/8H0SUBJECT5X,8HSTIMULUS) MDPF1330 102 FORMAT(1H I4,10F12.4/(5X,10F12.4)) MDPF1340 103 FORMAT(3H C=3F13.4) MDPF1350 106 FORMAT(18HFIRST SCORE MATRIX/11H(2X,5F13.5)) MDPF1360 107 FORMAT(I2,5F13.5/(2X,5F13.5)) MDPF1370 25 RETURN MDPF1380 END MDPF1390 SUBROUTINE COMPR(NP,NS,NF,W,R,X,Y,NPUNCH) MDPF1400 DIMENSION Y(60,60) MDPF1410 DIMENSION R(60,60) MDPF1420 DIMENSION W(60,60) MDPF1430 DIMENSION X(60,60) MDPF1440 C COMPUTE CROSS PRODUCT MATRIX OF SUBJECTS Y=W(WT) MDPF1450 DO 41 I=1,NP MDPF1460 DO 41 J=I,NP MDPF1470 Y(I,J)=0. MDPF1480 DO 43 L=1,NS MDPF1490 43 Y(I,J)=Y(I,J)+W(I,L)*W(J,L) MDPF1500 41 Y(J,I)=Y(I,J) MDPF1510 PRINT 109 MDPF1520 DO 42 I=1,NP MDPF1530 42 PRINT 104,I,(Y(I,J),J=1,NP) MDPF1540 C COMPUTE CORRELATION MATRIX OF SUBJECTS AND PRINT OUT MDPF1550 DO 44 I=2,NP MDPF1560 II=I-1 MDPF1570 DO 44 J=1,II MDPF1580 R(I,J)=Y(I,J)/SQRT(Y(I,I)*Y(J,J)) MDPF1590 R(J,I)=R(I,J) MDPF1600 44 CONTINUE MDPF1610 DO 45 I=1,NP MDPF1620 45 R(I,I)=1. MDPF1630 PRINT 111 MDPF1640 DO 46 I=1,NP MDPF1650 46 PRINT 104,I,(R(I,J),J=1,NP) MDPF1660 C COMPUTE CROSS PRODUCT MATRIX OF STIMULI R=WT(W) MDPF1670 DO 20 I=1,NS MDPF1680 DO 20 J=I,NS MDPF1690 R(I,J)=0. MDPF1700 DO 21 K=1,NP MDPF1710 21 R(I,J)=R(I,J)+W(K,I)*W(K,J) MDPF1720 20 R(J,I)=R(I,J) MDPF1730 PRINT 110 MDPF1740 DO 40 I=1,NS MDPF1750 40 PRINT 104,I,(R(I,J),J=1,NS) MDPF1760 C FIND ROOTS AND VECTORS MDPF1770 M=1 MDPF1780 PRINT 311 MDPF1790 IF(NS.LE.NP) GO TO 61 MDPF1800 C THE PURPOSE OF FINDING ROOTS AND VECTORS OF THE LESSER ORDER IS TOMDPF1810 C SAVE COMPUTOR TIME MDPF1820 C FIND ROOTS AND VECTORS OF Y, IF NP.LT.NS. MDPF1830 C CALL MTEVV(Y,100,100,R,100,100,NP,M) MDPF1840 CALL MTEVV(Y,60,60,R,60,60,NP,M) MDPF1850 C CALL MTEVV(Y,20,20,R,20,20,NP,M) MDPF1860 PRINT 312,(Y(I,I),I=1,NP) MDPF1870 DO 47 I=1,NF MDPF1880 47 Y(I,I)=SQRT(Y(I,I)) MDPF1890 C COMPUTE POPULATION MATRIX X MDPF1900 DO 62 J=1,NP MDPF1910 SUMQ=0. MDPF1920 DO 67 I=1,NF MDPF1930 X(J,I)=R(J,I)*Y(I,I) MDPF1940 SUMQ=SUMQ+X(J,I)**2 MDPF1950 67 R(J,I)=R(J,I)/Y(I,I) MDPF1960 SUMQ=SQRT(SUMQ) MDPF1970 DO 66 I=1,NF MDPF1980 66 X(J,I)=X(J,I)/SUMQ MDPF1990 62 CONTINUE MDPF2000 C COMPUTE STIMULUS MATRIX Y MDPF2010 DO 63 J=1,NF MDPF2020 SUM=0. MDPF2030 DO 68 I=1,NS MDPF2040 Y(I,J)=0. MDPF2050 DO 64 M=1,NP MDPF2060 64 Y(I,J)=Y(I,J)+W(M,I)*R(M,J) MDPF2070 68 SUM=SUM+Y(I,J)**2 MDPF2080 SUM=SQRT(SUM) MDPF2090 C NORMALIZE Y MDPF2100 DO 65 L=1,NS MDPF2110 65 Y(L,J)=Y(L,J)/SUM MDPF2120 63 CONTINUE MDPF2130 GO TO 28 MDPF2140 C FIND ROOTS AND VECTORS OF R, IF NS.LE.NP. MDPF2150 C61 CALL MTEVV(R,100,100,Y,100,100,NS,M) MDPF2160 61 CALL MTEVV(R,60,60,Y,60,60,NS,M) MDPF2170 C COMPUTE POPULATION MATRIX X MDPF2180 PRINT 312,(R(I,I),I=1,NS) MDPF2190 DO 30 I=1,NP MDPF2200 SUMQ=0. MDPF2210 DO 25 K=1,NF MDPF2220 X(I,K)=0. MDPF2230 DO 22 J=1,NS MDPF2240 22 X(I,K)=X(I,K)+W(I,J)*Y(J,K) MDPF2250 25 SUMQ=SUMQ+X(I,K)**2 MDPF2260 SUMQ=SQRT(SUMQ) MDPF2270 DO 29 K=1,NF MDPF2280 29 X(I,K)=X(I,K)/SUMQ MDPF2290 30 CONTINUE MDPF2300 C COMPUTE SECOND SCORE MATRIX AND STORE IN W MDPF2310 28 DO 329 I=1,NP MDPF2320 DO 329 J=1,NS MDPF2330 W(I,J)=0. MDPF2340 DO 328 K=1,NF MDPF2350 328 W(I,J)=W(I,J)+X(I,K)*Y(J,K) MDPF2360 329 CONTINUE MDPF2370 PRINT 300 MDPF2380 DO 35 I=1,NP MDPF2390 35 PRINT 104,I,(W(I,J),J=1,NS) MDPF2400 C PRINT AND PUNCH X AND Y MDPF2410 PRINT 103 MDPF2420 DO 53 I=1,NF MDPF2430 53 PRINT 104,I,(X(J,I),J=1,NP) MDPF2440 PRINT 105 MDPF2450 DO 54 I=1,NF MDPF2460 54 PRINT 104,I,(Y(J,I),J=1,NS) MDPF2470 IF(NPUNCH)51,50,51 MDPF2480 51 PUNCH 106 MDPF2490 DO 57 I=1,NF MDPF2500 57 PUNCH 107,I,(Y(J,I),J=1,NS) MDPF2510 PUNCH 108 MDPF2520 DO 58 I=1,NF MDPF2530 58 PUNCH 107,I,(X(J,I),J=1,NP) MDPF2540 103 FORMAT(18H0POPULATION MATRIX/7H0FACTOR) MDPF2550 104 FORMAT(1H I4,10F12.4/(5X,10F12.4)) MDPF2560 105 FORMAT(16H0STIMULUS MATRIX/7H0FACTOR) MDPF2570 106 FORMAT(15HSTIMULUS MATRIX/11H(2X,5F13.5)) MDPF2580 107 FORMAT(I2,5F13.5/(2X,5F13.5)) MDPF2590 108 FORMAT(17HPOPULATION MATRIX/11H(2X,5F13.5)) MDPF2600 109 FORMAT(33H0CROSS PRODUCT MATRIX OF SUBJECTS) MDPF2610 110 FORMAT(32H0CROSS PRODUCT MATRIX OF STIMULI) MDPF2620 111 FORMAT(31H0CORRELATION MATRIX OF SUBJECTS) MDPF2630 300 FORMAT(20H0SECOND SCORE MATRIX/8H0SUBJECT5X,8HSTIMULUS) MDPF2640 311 FORMAT(32H0ROOTS OF THE FIRST SCORE MATRIX) MDPF2650 312 FORMAT(3X,10F12.4) MDPF2660 50 RETURN MDPF2670 END MDPF2680 SUBROUTINE LSQMD(D,WT,Y,NS,NP,W,L,NS1,SCR,VI) MDPF2690 REAL*8 SCR(3600) MDPF2700 DIMENSION Y(3600) MDPF2710 DIMENSION D(60,60) MDPF2720 DIMENSION WT(60,60) MDPF2730 DIMENSION VI(60,60) MDPF2740 DIMENSION W(60,60) MDPF2750 C BEFORE MAKING ANY CHANGES IN THIS ROUTINE, CHECK THE EQUIVALENCE MDPF2760 C STATEMENT IN MAIN PROGRAM INVOLVING SCR,D AND VI. MDPF2770 C DATA IN D AND WEIGHT IN WT MDPF2780 PRINT 100,L MDPF2790 BIG=10.**8 MDPF2800 C COMPUTE TOTAL VARIANCE MDPF2810 TVAR=0. MDPF2820 DO 2 I=1,NS MDPF2830 Y(I)=0. MDPF2840 DO 3 J=1,NS MDPF2850 TVAR=TVAR+WT(I,J)*D(I,J)**2 MDPF2860 3 Y(I)=Y(I)+WT(I,J)*D(I,J)-WT(J,I)*D(J,I) MDPF2870 2 CONTINUE MDPF2880 IF(L.EQ.1.OR.L.GT.NS1) GO TO 50 MDPF2890 GO TO 52 MDPF2900 C COMPUTE U AND PUT IN D MDPF2910 50 DO 5 I=1,NS MDPF2920 D(I,I)=0. MDPF2930 DO 5 J=1,NS MDPF2940 IF(I.EQ.J) GO TO 5 MDPF2950 D(I,J)=-(WT(I,J)+WT(J,I)) MDPF2960 D(I,I)=D(I,I)-D(I,J) MDPF2970 5 CONTINUE MDPF2980 DO 8 I=1,NS MDPF2990 DO 8 J=1,I MDPF3000 WT(I,J)=ABS(D(I,J)) MDPF3010 IF(D(I,J).EQ.0.) WT(I,J)=BIG MDPF3020 8 WT(J,I)=WT(I,J) MDPF3030 C PARTITION STIMULI IN BLOCKS MDPF3040 CALL BLOCK(WT,NS) MDPF3050 20 DO 21 I=1,NS MDPF3060 DO 21 J=1,NS MDPF3070 IF(I.EQ.J) GO TO 22 MDPF3080 IF(WT(I,J).GE.BIG) GO TO 24 MDPF3090 22 WT(I,J)=1. MDPF3100 GO TO 21 MDPF3110 24 WT(I,J)=0. MDPF3120 21 CONTINUE MDPF3130 C COMPUTE V MDPF3140 DO 25 I=1,NS MDPF3150 DO 25 J=1,NS MDPF3160 25 WT(I,J)=WT(I,J)+D(I,J) MDPF3170 MM=1 MDPF3180 C ON RETURN THE INVERSE IS STORED IN WT (ORDER NS). SCR IS USED AS MDPF3190 C SCRATCH MDPF3200 CALL MTINV(WT,60,60,SCR,60,60,NS,NN) MDPF3210 C CALL MTINV(WT,20,20,SCR,20,20,NS,NM) MDPF3220 C CALL MTINV(WT,100,100,SCR,100,100,NS,MM) MDPF3230 IF(MM.EQ.1) GO TO 53 MDPF3240 C ON RETURN IF MM IS NOT 1 INDICATES WT IS ILL CONDITIONED. MDPF3250 C PROGRAM PRINTS OUT MESSAGE AND STOPS. MDPF3260 PRINT 111,MM MDPF3270 STOP MDPF3280 53 DO 51 I=1,NS MDPF3290 DO 51 J=1,NS MDPF3300 51 VI(I,J)=WT(I,J) MDPF3310 52 DO 26 I=1,NS MDPF3320 W(L,I)=0. MDPF3330 DO 26 J=1,NS MDPF3340 26 W(L,I)=W(L,I)+Y(J)*VI(J,I) MDPF3350 PRINT 109 MDPF3360 PRINT 102,I,(Y(I),I=1,NS) MDPF3370 C COMPUTE VARIANCE ACCOUNTED FOR MDPF3380 VAF=0. MDPF3390 DEN=0. MDPF3400 DO 27 I=1,NS MDPF3410 VAF=VAF+W(L,I)*Y(I) MDPF3420 27 DEN=DEN+W(L,I)**2 MDPF3430 PVAF=VAF/TVAR MDPF3440 C NORMALIZE MDPF3450 C=SQRT(VAF/DEN) MDPF3460 DO 28 I=1,NS MDPF3470 28 W(L,I)=W(L,I)*C MDPF3480 PRINT 107,TVAR,VAF,PVAF,C MDPF3490 100 FORMAT(8H0SUBJECTI4) MDPF3500 102 FORMAT(I4,10F12.4/(4X,10F12.4)) MDPF3510 107 FORMAT(6H0TVAR=F12.4,3X,4HVAF=F12.4,3X,5HPVAF=F12.4,3X,2HC=F12.4) MDPF3520 109 FORMAT(9H0Y VECTOR) MDPF3530 111 FORMAT(4H0MM=I1,21HWT CANNOT BE INVERTED) MDPF3540 RETURN MDPF3550 END MDPF3560 SUBROUTINE RFSM(W,NP,NS,IREAD,IPUNF) MDPF3570 DIMENSION W(60,60) MDPF3580 DIMENSION WMEAN(100),SD(100),FMT(18) MDPF3590 XM=NS MDPF3600 DEN=1./XM MDPF3610 READ 100,(FMT(I),I=1,18) MDPF3620 DO 5 I=1,NP MDPF3630 5 READFMT,(W(I,J),J=1,NS) MDPF3640 C COMPUTE MEAN FOR EACH SUBJECT MDPF3650 DO 6 I=1,NP MDPF3660 SUM=0. MDPF3670 DO 4 J=1,NS MDPF3680 4 SUM=SUM+W(I,J) MDPF3690 WMEAN(I)=SUM/XM MDPF3700 C NORMALIZE W - SUBTRACT THE ROW MEAN MDPF3710 SUM=0. MDPF3720 DO 7 J=1,NS MDPF3730 W(I,J)=W(I,J)-WMEAN(I) MDPF3740 IF(IREAD.EQ.2)SUM=SUM+W(I,J)**2 MDPF3750 7 CONTINUE MDPF3760 IF(IREAD.EQ.1) GO TO 6 MDPF3770 C IREAD=2, NORMALIZED W - DIVIDE BY THE ROW SD MDPF3780 SD(I)=SQRT(SUM*DEN) MDPF3790 DO 8 J=1,NS MDPF3800 8 W(I,J)=W(I,J)/SD(I) MDPF3810 6 CONTINUE MDPF3820 C PRINT MEAN MDPF3830 PRINT 103 MDPF3840 PRINT 104,(WMEAN(I),I=1,NP) MDPF3850 IF(IREAD.EQ.1) GO TO 10 MDPF3860 PRINT 105 MDPF3870 PRINT 104,(SD(I),I=1,NP) MDPF3880 C PRINT W MDPF3890 10 PRINT 101 MDPF3900 DO 20 I=1,NP MDPF3910 20 PRINT 102,I,(W(I,J),J=1,NS) MDPF3920 IF(IPUNF)22,25,22 MDPF3930 C PUNCH FIRST SCORE MATRIX MDPF3940 22 PUNCH 106 MDPF3950 DO 21 I=1,NP MDPF3960 21 PUNCH 107,I,(W(I,J),J=1,NS) MDPF3970 100 FORMAT(18A4) MDPF3980 101 FORMAT(19H0FIRST SCORE MATRIX/8H0SUBJECT5X,8HSTIMULUS) MDPF3990 102 FORMAT(1H I4,10F12.4/(5X,10F12.4)) MDPF4000 103 FORMAT(37H0MEAN OF THE RAW SCORES (BY SUBJECT)) MDPF4010 104 FORMAT(5X,10F12.4) MDPF4020 105 FORMAT(35H0SD OF THE RAW SCORES (BY SUBJECT)) MDPF4030 106 FORMAT(18HFIRST SCORE MATRIX/11H(2X,5F13.5)) MDPF4040 107 FORMAT(I2,5F13.5/(2X,5F13.5)) MDPF4050 25 RETURN MDPF4060 END MDPF4070 SUBROUTINE MDN(Y,WT,NS) MDPF4080 DIMENSION Y(3600) MDPF4090 DIMENSION WT(60,60) MDPF4100 DIMENSION W(60,60) MDPF4110 COMMON W MDPF4120 DO 6 I=1,NS MDPF4130 SUM=0. MDPF4140 SSQ=0. MDPF4150 SDSQ=0. MDPF4160 DO 5 J=1,NS MDPF4170 IF(I.EQ.J) GO TO 5 MDPF4180 SUM=SUM+WT(I,J)+WT(J,I) MDPF4190 SDSQ=SDSQ+(WT(J,I)-WT(I,J))**2 MDPF4200 SSQ=SSQ+WT(J,I)**2+WT(I,J)**2 MDPF4210 5 CONTINUE MDPF4220 IF(SSQ)4,7,4 MDPF4230 7 Y(I)=2. MDPF4240 GO TO 6 MDPF4250 4 Y(I)=2.+SUM-SQRT(SDSQ/SSQ) MDPF4260 6 CONTINUE MDPF4270 PRINT 100 MDPF4280 PRINT 101,(Y(I),I=1,NS) MDPF4290 100 FORMAT(9H0VECTOR N) MDPF4300 101 FORMAT(2X,10F12.4) MDPF4310 RETURN MDPF4320 END MDPF4330 SUBROUTINE READS(D,FMT,CODE,NS,WT) MDPF4340 C READ IN DATA MATRIX WITHOUT MISSING DATA MDPF4350 DIMENSION D(60,60),WT(60,60) MDPF4360 DIMENSION FMT(1), CODE(1) MDPF4370 N=1 MDPF4380 C READ IN DATA MATRIX MDPF4390 70 DO 15 I=1,NS MDPF4400 15 READFMT,(D(I,J),J=1,NS) MDPF4410 GO TO(1,2,3),N MDPF4420 C NO MISSING DATA MDPF4430 1 DO 50 I=1,NS MDPF4440 DO 50 J=1,NS MDPF4450 IF(I.EQ.J) GO TO 45 MDPF4460 DO 32 K=1,4 MDPF4470 IF(D(I,J)-CODE(K))32,41,32 MDPF4480 41 GO TO(10,20,45,45),K MDPF4490 10 D(I,J)=1. MDPF4500 GO TO 50 MDPF4510 20 D(I,J)=-1. MDPF4520 GO TO 50 MDPF4530 32 CONTINUE MDPF4540 45 D(I,J)=0. MDPF4550 50 CONTINUE MDPF4560 RETURN MDPF4570 ENTRY READMW(D,WT,FMT,CODE,NS) MDPF4580 C MISSING DATA, WEIGHT IS EITHER 0 OR 1 MDPF4590 N=2 MDPF4600 GO TO 70 MDPF4610 2 DO 55 I=1,NS MDPF4620 DO 55 J=1,NS MDPF4630 IF(I.EQ.J) GO TO 55 MDPF4640 IF(D(I,J).EQ.CODE(4)) GO TO 53 MDPF4650 WT(I,J)=1. MDPF4660 GO TO 55 MDPF4670 53 WT(I,J)=0. MDPF4680 55 CONTINUE MDPF4690 GO TO 1 MDPF4700 ENTRY READW(D,WT,FMT,CODE,NS) MDPF4710 C READ IN WEIGHT MATRIX MDPF4720 N=3 MDPF4730 GO TO 70 MDPF4740 3 DO 62 I=1,NS MDPF4750 62 READ 100,(WT(I,J),J=1,NS) MDPF4760 100 FORMAT(10F7.2) MDPF4770 DO 65 I=1,NS MDPF4780 DO 65 J=1,NS MDPF4790 IF(I.EQ.J) GO TO 65 MDPF4800 IF(D(I,J).EQ.CODE(4)) WT(I,J)=0. MDPF4810 65 CONTINUE MDPF4820 GO TO 1 MDPF4830 END MDPF4840 SUBROUTINE BLOCK(D,N) MDPF4850 DIMENSION D(60,60) MDPF4860 C COMPUTE MINIMUM PATH MDPF4870 DO 5 I=1,N MDPF4880 DO 5 J=2,N MDPF4890 IF(I.EQ.J) GO TO 5 MDPF4900 JJ=J-1 MDPF4910 DO 4 K=1,JJ MDPF4920 IF(I.EQ.J) GO TO 4 MDPF4930 DSUM=D(I,J)+D(I,K) MDPF4940 IF(D(J,K).GT.DSUM) D(J,K)=DSUM MDPF4950 D(K,J)=D(J,K) MDPF4960 4 CONTINUE MDPF4970 5 CONTINUE MDPF4980 RETURN MDPF4990 END MDPF5000 SUBROUTINE SORT (A, N, B, C, D, K, SWITCH ) MDPF5010 CSORT MDPF5020 C MDPF5030 C THIS ROUTINE SORTS INPUT ARRAY 'A' AND REARRANGES, OPTIONALLY, MDPF5040 C ARRAYS 'B', 'C', AND 'D', IN ORDER CORRESPONDING TO 'A'. MDPF5050 C N = NUMBER OF ITEMS IN 'A' (AND 'B', 'C', 'D', IF USED) MDPF5060 C K = 0--SORT 'A' ONLY, 1--REARRANGE 'B', 2--REARRANGE 'B' AND 'C', MDPF5070 C 3--REARRANGE 'B', 'C', AND 'D'. MDPF5080 C IF 'SWITCH' IS POSITIVE, SORT WILL BE IN ASCENDING ORDER, MDPF5090 C IF ZERO OR NEGATIVE, IN DESCENDING ORDER. MDPF5100 C ALGORITHM FROM CACM, JULY 1959, PAGE 30 BY D. L. SHELL MDPF5110 C MDPF5120 DIMENSION A(N),B(N),C(N),D(N) MDPF5130 INTEGER SWITCH MDPF5140 105 KP1=K+1 MDPF5150 IF(N.LE.1) GO TO 999 MDPF5160 M = 1 MDPF5170 106 M = M + M MDPF5180 IF( M .LT. N ) GO TO 106 MDPF5190 M = M - 1 MDPF5200 994 M = M/2 MDPF5210 IF(M.EQ.0) GO TO 999 MDPF5220 KK = N-M MDPF5230 J = 1 MDPF5240 992 I = J MDPF5250 996 IM = I + M MDPF5260 IF(SWITCH) 810,810,800 MDPF5270 800 IF(A(I).GT.A(IM)) GO TO 110 MDPF5280 GO TO 995 MDPF5290 810 IF(A(I).LT.A(IM)) GO TO 110 MDPF5300 995 J = J+1 MDPF5310 IF(J.GT.KK) GO TO 994 MDPF5320 GO TO 992 MDPF5330 110 TEMP=A(I) MDPF5340 A(I) = A(IM) MDPF5350 A(IM) = TEMP MDPF5360 GO TO ( 140, 130, 120, 115), KP1 MDPF5370 115 TEMP = D(I) MDPF5380 D(I) = D(IM) MDPF5390 D(IM) = TEMP MDPF5400 120 TEMP=C(I) MDPF5410 C(I) = C(IM) MDPF5420 C(IM) = TEMP MDPF5430 130 TEMP=B(I) MDPF5440 B(I) = B(IM) MDPF5450 B(IM) = TEMP MDPF5460 140 I = I-M MDPF5470 IF(I.LT.1) GO TO 995 MDPF5480 GO TO 996 MDPF5490 999 RETURN MDPF5500 END MDPF5510 SUBROUTINE MTINV (A,I1,I2,B,I3,I4,N,M) MDPF5520 C INCODE IBMF MDPF5530 C FORTRAN DECK,LSTOU,STAB MDPF5540 C MTINV MDPF5550 DOUBLE PRECISION DABS MDPF5560 REAL*8 B(I3,I4),BMAX,BMULT MDPF5570 DIMENSION A(I1,I2) MDPF5580 C MDPF5590 C THIS SUBROUTINE CALCULATES THE INVERSE OF MATRIX A BY THE GAUSS- MDPF5600 C JORDAN ELIMINATION SCHEME. ALL CALCULATIONS ARE DONE IN DOUBLE- MDPF5610 C PRECISION ARITHMETIC MDPF5620 C MDPF5630 EQUIVALENCE (BMAX,BMULT) MDPF5640 C MDPF5650 C SINGLE TO DOUBLE PRECISION TRANSFER MDPF5660 C MDPF5670 101 DO 103 I = 1,N MDPF5680 102 DO 103 J = 1,N MDPF5690 103 B(I,J) = A(I,J) MDPF5700 C MDPF5710 C FIND THE LARGEST ELEMENT IN THE N-K BY N-K LOWER RIGHT SUBMARTIX. MDPF5720 C MDPF5730 200 DO 430 K = 1,N MDPF5740 A(K,1) = FLOAT(K) MDPF5750 A(K,2) = A(K,1) MDPF5760 203 BMAX = DABS(B(K,K)) MDPF5770 210 DO 219 I = K,N MDPF5780 211 DO 219 J = K,N MDPF5790 IF (BMAX - DABS(B(I,J)) ) 213,219,219 MDPF5800 213 BMAX = DABS ( B(I,J)) MDPF5810 214 A(K,1) = FLOAT(I) MDPF5820 215 A(K,2) = FLOAT(J) MDPF5830 219 CONTINUE MDPF5840 216 IF (BMAX - 1.D-38) 801,301,301 MDPF5850 C MDPF5860 C EXCHANGE ROWS AND COLUMNS TO PUT B(I,J) ON DIAGONAL. MDPF5870 C MDPF5880 301 I = IFIX(A(K,1)) MDPF5890 302 J = IFIX(A(K,2)) MDPF5900 313 DO 314 M = 1,N MDPF5910 BMULT = B(I,M) MDPF5920 B(I,M) = B(K,M) MDPF5930 314 B(K,M) = BMULT MDPF5940 303 DO 304 M = 1,N MDPF5950 BMULT = B(M,J) MDPF5960 B(M,J) = B(M,K) MDPF5970 304 B(M,K) = BMULT MDPF5980 C MDPF5990 C ACTUAL INVERSION STEP. MDPF6000 C BEGIN ROW ITERATION MDPF6010 C IGNORE ROW K MDPF6020 401 DO 425 I = 1,N MDPF6030 C MDPF6040 402 IF (I - K) 405,425,405 MDPF6050 C SET ROW MULTIPLIER. MDPF6060 405 BMULT = B(I,K) / B(K,K) MDPF6070 C MODIFY ROW ELEMENTS. MDPF6080 410 DO 420 J = 1,N MDPF6090 C IGNORE COLUMN K MDPF6100 411 IF (J - K) 414,412,414 MDPF6110 412 B(I,J) = -BMULT MDPF6120 413 GO TO 420 MDPF6130 414 B(I,J) = B(I,J) - B(K,J) * BMULT MDPF6140 420 CONTINUE MDPF6150 425 CONTINUE MDPF6160 C MDPF6170 C DIVIDE PIVOT ROW BY PIVOT ELEMENT MDPF6180 C MDPF6190 426 BMULT = 1.D0/B(K,K) MDPF6200 427 B(K,K) = 1.D0 MDPF6210 C ACTUAL DIVISION MDPF6220 428 DO 429 J = 1,N MDPF6230 429 B(K,J) = B(K,J) * BMULT MDPF6240 430 CONTINUE MDPF6250 C MDPF6260 C REARRANGE AND REASSEMBLE INVERSE MATRIX. MDPF6270 C MDPF6280 501 DO 512 IJK = 1,N MDPF6290 502 I = N-IJK+1 MDPF6300 503 L = IFIX(A(I,2)) MDPF6310 504 J = IFIX(A(I,1)) MDPF6320 505 DO 508 M = 1,N MDPF6330 506 BMULT = B(I,M) MDPF6340 507 B(I,M) = B(L,M) MDPF6350 508 B(L,M) = BMULT MDPF6360 509 DO 512 M = 1,N MDPF6370 510 BMULT = B(M,I) MDPF6380 511 B(M,I) = B(M,J) MDPF6390 512 B(M,J) = BMULT MDPF6400 C MDPF6410 C MOVE A INVERSE BACK INTO A MDPF6420 C MDPF6430 601 DO 603 I = 1,N MDPF6440 602 DO 603 J = 1,N MDPF6450 603 A(I,J) = B(I,J) MDPF6460 701 M = 1 MDPF6470 702 RETURN MDPF6480 C MDPF6490 C ERROR RETURN AT SOME TIME BMAX WAS .LT. 0.1D-38 MDPF6500 C MDPF6510 801 M = 2 MDPF6520 802 RETURN MDPF6530 END MDPF6540 SUBROUTINE ARREV(N,NS,D,I3,I4) MDPF6550 DIMENSION NS(N),D(I3,I4) MDPF6560 DO 10 I=1,N MDPF6570 1 CONTINUE MDPF6580 IF(NS(I)-I) 2,10,2 MDPF6590 2 CONTINUE MDPF6600 KK=I MDPF6610 JJ=NS(I) MDPF6620 CALL XCH(N,KK,JJ,D,NS,I3,I4) MDPF6630 GO TO 1 MDPF6640 10 CONTINUE MDPF6650 RETURN MDPF6660 END MDPF6670 SUBROUTINE XCH(N,L,M,D,NS,I3,I4) MDPF6680 DIMENSION NS(N),D(I3,I4) MDPF6690 I=L MDPF6700 J=M MDPF6710 DO 1 K=1,N MDPF6720 T=D(K,I) MDPF6730 D(K,I)=D(K,J) MDPF6740 1 D(K,J)=T MDPF6750 NS(I)=NS(J) MDPF6760 NS(J)=J MDPF6770 RETURN MDPF6780 END MDPF6790 SUBROUTINE SORT1(A, N, B, C, D, K, SWITCH ) MDPF6800 DIMENSION A(N),B(N),C(N),D(N) MDPF6810 INTEGER SWITCH, SW, A, B, C, D, TEMP MDPF6820 IF(SWITCH) 100, 105 , 100 MDPF6830 100 SW = SWITCH MDPF6840 105 KP1=K+1 MDPF6850 NM1 = N-1 MDPF6860 IF(NM1) 300,300,107 MDPF6870 107 DO 150 I=1,NM1 MDPF6880 IP1 = I+1 MDPF6890 DO 140 J=IP1,N MDPF6900 IF(A(I)-A(J)) 140 , 140 , 110 MDPF6910 110 TEMP=A(I) MDPF6920 A(I) = TEMP MDPF6930 A(J) = TEMP MDPF6940 GO TO(140,130, 120,115), KP1 MDPF6950 115 TEMP=D(I) MDPF6960 D(I) = D(J) MDPF6970 D(J) = TEMP MDPF6980 120 TEMP= C(I) MDPF6990 C(I) = C(J) MDPF7000 C(J) = TEMP MDPF7010 130 TEMP= B(I) MDPF7020 B(I) = B(J) MDPF7030 B(J) = TEMP MDPF7040 140 CONTINUE MDPF7050 150 CONTINUE MDPF7060 IF(SW ) 200,300,300 MDPF7070 200 NHALF=N/2 MDPF7080 DO 230 I=1,NHALF MDPF7090 IC=N+1-I MDPF7100 TEMP= A(I) MDPF7110 A(I) = A(IC) MDPF7120 A(IC) = TEMP MDPF7130 GO TO ( 230,220,210,205 ) ,KP1 MDPF7140 205 TEMP= D(I) MDPF7150 D(I) = D(J) MDPF7160 D(IC) = TEMP MDPF7170 210 TEMP= C(I) MDPF7180 C(I) = C(IC) MDPF7190 C(IC) = TEMP MDPF7200 220 TEMP = B(I) MDPF7210 B(I) = B(IC) MDPF7220 B(IC) = TEMP MDPF7230 230 CONTINUE MDPF7240 300 CONTINUE MDPF7250 RETURN MDPF7260 END MDPF7270 SUBROUTINE MTEVV(X,I1,I2,D,I3,I4,N,M) MDPF7280 DIMENSION A(1830),X(60,60),D(60,60),E(1830),S(220),F(60,60) MDPF7290 K=1 MDPF7300 DO 1 J=1,N MDPF7310 DO 1 I=1,J MDPF7320 A(K)=X(I,J) MDPF7330 1 K=K+1 MDPF7340 CALL EIGEN(A,E,N,0) MDPF7350 II = 1 MDPF7360 DO 974 I=1,N MDPF7370 DO 974 J=1,N MDPF7380 D(J,I) = E(II) MDPF7390 II = II + 1 MDPF7400 974 CONTINUE MDPF7410 II = 0 MDPF7420 DO 973 J=1,N MDPF7430 II = II + J MDPF7440 973 X(J,J) = A(II) MDPF7450 RETURN MDPF7460 END MDPF7470 SUBROUTINE EIGEN(A,R,N,MV) MDPF7480 C MDPF7490 C SUBROUTINE EIGEN MDPF7500 C MDPF7510 C PURPOSE MDPF7520 C COMPUTE EIGENVALUES AND EIGENVECTORS OF A REAL SYMMETRIC MDPF7530 C MATRIX MDPF7540 C MDPF7550 C USAGE MDPF7560 C CALL EIGEN(A,R,N,MV) MDPF7570 C MDPF7580 C DESCRIPTION OF PARAMETERS MDPF7590 C A - ORIGINAL MATRIX (SYMMETRIC), DESTROYED IN COMPUTATION. MDPF7600 C RESULTANT EIGENVALUES ARE DEVELOPED IN DIAGONAL OF MDPF7610 C MATRIX A IN DECENDING ORDER. MDPF7620 C R - RESULTANT MATRIX OF EIGENVECTORS (STORED COLUMNWISE, MDPF7630 C IN SAME SEQUENCE AS EIGENVALUES) MDPF7640 C N - ORDER OF MATRICES A AND R MDPF7650 C MV- INPUT CODE MDPF7660 C 0 COMPUTE EIGENVALUES AND EIGENVECTORS MDPF7670 C 1 COMPUTE EIGENVALUES ONLY (R NEED NOT BE MDPF7680 C DIMENSIONED BUT MUST STILL APPEAR IN CALLING MDPF7690 C SEQUENCE) MDPF7700 C MDPF7710 C REMARKS MDPF7720 C ORIGINAL MATRIX A MUST BE REAL SYMMETRIC (STORAGE MODE=1) MDPF7730 C MATRIX A CANNOT BE IN THE SAME LOCATION AS MATRIX R MDPF7740 C MDPF7750 C SUBROUTINES USED BY THIS SUBROUTINE MDPF7760 C NONE MDPF7770 C MDPF7780 C METHOD MDPF7790 C DIAGONALIZATION METHOD ORIGINATED BY JACOBI AND ADAPTED MDPF7800 C BY VON NEUMANN FOR LARGE COMPUTERS AS FOUND IN MATHEMATICAL MDPF7810 C METHODS FOR DIGITAL COMPUTERS, EDITED BY A. RALSTON AND MDPF7820 C H.S. WILF, JOHN WILEY AND SONS, NEW YORK, 1962, CHAPTER 7 MDPF7830 C MDPF7840 C ..................................................................MDPF7850 C MDPF7860 C MDPF7870 DIMENSION A(1830), R(1830) MDPF7880 C GENERATE IDENTITY MATRIX MDPF7890 C MDPF7900 IF(MV-1) 22,14,22 MDPF7910 22 IQ=-N MDPF7920 DO 1 J=1,N MDPF7930 IQ=IQ+N MDPF7940 DO 1 I=1,N MDPF7950 IJ=IQ+I MDPF7960 R(IJ)=0.0 MDPF7970 IF(I-J) 1,31,1 MDPF7980 31 R(IJ)=1.0 MDPF7990 1 CONTINUE MDPF8000 C MDPF8010 C COMPUTE INITIAL AND FINAL NORMS (ANORM AND ANORMX) MDPF8020 C MDPF8030 14 ANORM=0.0 MDPF8040 DO 2 I=1,N MDPF8050 DO 2 J=I,N MDPF8060 IF(I-J) 15,2,15 MDPF8070 15 IA=I+(J*J-J)/2 MDPF8080 ANORM=ANORM+A(IA)*A(IA) MDPF8090 2 CONTINUE MDPF8100 IF(ANORM) 39, 39, 50 MDPF8110 50 ANORM=2.0*SQRT(ANORM) MDPF8120 ANRMX =ANORM*1.0E-8 MDPF8130 C MDPF8140 C INITIALIZE INDICATORS AND COMPUTE THRESHOLD, THR MDPF8150 C MDPF8160 IND=0 MDPF8170 THR=ANORM MDPF8180 3 THR=THR/FLOAT(N) MDPF8190 13 L=1 MDPF8200 4 M=L+1 MDPF8210 C MDPF8220 C COMPUTE SIN AND COS MDPF8230 C MDPF8240 5 MQ=(M*M-M)/2 MDPF8250 LQ=(L*L-L)/2 MDPF8260 LM=L+MQ MDPF8270 IF(ABS(A(LM))-THR) 7,32,32 MDPF8280 32 IND=1 MDPF8290 LL=L+LQ MDPF8300 MM=M+MQ MDPF8310 X=0.5*(A(LL)-A(MM)) MDPF8320 Y=-A(LM)/SQRT(A(LM)*A(LM)+X*X) MDPF8330 IF(X) 33,34,34 MDPF8340 33 Y=-Y MDPF8350 34 SINX=Y/SQRT(2.0*(1.0+(SQRT(1.0-Y*Y)))) MDPF8360 SINX2=SINX*SINX MDPF8370 COSX=SQRT(1.0-SINX2) MDPF8380 COSX2=COSX*COSX MDPF8390 SINCS =SINX*COSX MDPF8400 C MDPF8410 C ROTATE L AND M COLUMNS MDPF8420 C MDPF8430 ILQ=N*(L-1) MDPF8440 IMQ=N*(M-1) MDPF8450 DO 6 I=1,N MDPF8460 IQ=(I*I-I)/2 MDPF8470 IF(I-L) 16,11,16 MDPF8480 16 IF(I-M) 17,11,18 MDPF8490 17 IM=I+MQ MDPF8500 GO TO 19 MDPF8510 18 IM=M+IQ MDPF8520 19 IF(I-L) 37,23,23 MDPF8530 37 IL=I+LQ MDPF8540 GO TO 24 MDPF8550 23 IL=L+IQ MDPF8560 24 X=A(IL)*COSX-A(IM)*SINX MDPF8570 A(IM)=A(IL)*SINX+A(IM)*COSX MDPF8580 A(IL)=X MDPF8590 11 IF(MV-1) 26,6,26 MDPF8600 26 ILR=ILQ+I MDPF8610 IMR=IMQ+I MDPF8620 X=R(ILR)*COSX-R(IMR)*SINX MDPF8630 R(IMR)=R(ILR)*SINX+R(IMR)*COSX MDPF8640 R(ILR)=X MDPF8650 6 CONTINUE MDPF8660 X=2.0*A(LM)*SINCS MDPF8670 Y=A(LL)*COSX2+A(MM)*SINX2-X MDPF8680 X=A(LL)*SINX2+A(MM)*COSX2+X MDPF8690 A(LM)=(A(LL)-A(MM))*SINCS+A(LM)*(COSX2-SINX2) MDPF8700 A(LL)=Y MDPF8710 A(MM)=X MDPF8720 C MDPF8730 C TESTS FOR COMPLETION MDPF8740 C MDPF8750 C TEST FOR M = LAST COLUMN MDPF8760 7 IF(M-N) 35,8,35 MDPF8770 35 M=M+1 MDPF8780 GO TO 5 MDPF8790 C TEST FOR L = SECOND FROM LAST COLUMN MDPF8800 8 IF(L-(N-1)) 36,9,36 MDPF8810 36 L=L+1 MDPF8820 GO TO 4 MDPF8830 9 IF(IND-1) 10,38,10 MDPF8840 38 IND=0 MDPF8850 GO TO 13 MDPF8860 C COMPARE THRESHOLD WITH FINAL NORM MDPF8870 10 IF(THR-ANRMX ) 39,39,3 MDPF8880 C MDPF8890 C SORT EIGENVALUES AND EIGENVECTORS MDPF8900 C MDPF8910 39 IQ=-N MDPF8920 DO 30 I=1,N MDPF8930 IQ=IQ+N MDPF8940 LL=I+(I*I-I)/2 MDPF8950 JQ=N*(I-2) MDPF8960 DO 30 J=I,N MDPF8970 JQ=JQ+N MDPF8980 MM=J+(J*J-J)/2 MDPF8990 IF(A(LL)-A(MM)) 40,30,30 MDPF9000 40 X=A(LL) MDPF9010 A(LL)=A(MM) MDPF9020 A(MM)=X MDPF9030 IF(MV-1) 28,30,28 MDPF9040 28 DO 20 K=1,N MDPF9050 ILR=IQ+K MDPF9060 IMR=JQ+K MDPF9070 X=R(ILR) MDPF9080 R(ILR)=R(IMR) MDPF9090 20 R(IMR)=X MDPF9100 30 CONTINUE MDPF9110 RETURN MDPF9120 END MDPF9130 SUBROUTINE PLOTX(A,B,N) DIMENSION A(1),B(1) CALL GRAPH(+1.2,-1.2,1.0,-1.0) CALL AXES CALL LOC2(A,B,N,1) CALL PLOT CALL IDENT RETURN END SUBROUTINE GRAPH(XMX,XMN,YMX,YMN) GRPH 10 DIMENSION Z(1),X(1),Y(1),SMALL(13),CHAR(50),HOLL(11) GRPH 20 REAL ITEM(55,121) GRPH 30 DATA CHAR /'1','2','3','4','5','6','7','8','9','A','B','C', GRPH 50 1 'D','E','F','G','H','I','J','K','L','M','N','O','P','Q','R','S', GRPH 60 2 'T','U','V','W','X','Y','Z','+','/','=','*','&','$','@','%','?', GRPH 70 3 '<','(',')','"',';',' '/ GRPH 80 DATA HOLL /' ','X','2','3','4','5','6','7','8','9','M'/ GRPH 90 DATA BLANK,DD,PLUS,AIE,AMINUS,ORIG,AMULT,GT50/' ','.','+','\', GRPH 100 1 '-','0','#','>'/ GRPH 110 DO 115 I = 1,55 GRPH 120 DO 115 J = 1,121 GRPH 130 115 ITEM(I,J) = BLANK GRPH 140 XMAX = XMX GRPH 150 XMIN = XMN GRPH 160 YMAX = YMX GRPH 170 YMIN = YMN GRPH 180 RETURN GRPH 190 ENTRY LIMIT(Z,NP,IS) GRPH 200 ZMAX = -1.0E9 GRPH 210 ZMIN = +1.0E9 GRPH 220 DO 10 I = 1,NP GRPH 230 IF(Z(I) .GT. ZMAX) ZMAX = Z(I) GRPH 240 10 IF(Z(I) .LT. ZMIN) ZMIN = Z(I) GRPH 250 IF(IS .EQ. 2) GO TO 11 GRPH 260 IF(ZMAX .GT. XMAX) XMAX = ZMAX GRPH 270 IF(ZMIN .LT. XMIN) XMIN = ZMIN GRPH 280 GO TO 12 GRPH 290 11 IF(ZMAX .GT. YMAX) YMAX = ZMAX GRPH 300 IF(ZMIN .LT. YMIN) YMIN = ZMIN GRPH 310 12 RETURN GRPH 320 ENTRY SCALE(IP) GRPH 330 XEXP = 0.05 * (XMAX - XMIN) GRPH 340 XMAX = XMAX + XEXP GRPH 350 XMIN = XMIN - XEXP GRPH 360 YEXP = 0.05 * (YMAX - YMIN) GRPH 370 YMAX = YMAX + YEXP GRPH 380 YMIN = YMIN - YEXP GRPH 390 IF(IP .NE. 1) RETURN GRPH 400 IF(YMAX - XMAX) 20,21,22 GRPH 410 20 YMAX = XMAX GRPH 420 GO TO 21 GRPH 430 22 XMAX = YMAX GRPH 440 21 CONTINUE GRPH 450 IF(YMIN - XMIN) 23,24,25 GRPH 460 23 XMIN = YMIN GRPH 470 GO TO 24 GRPH 480 25 YMIN = XMIN GRPH 490 24 EXP = 0.1667 * (XMAX - XMIN) GRPH 500 XMAX = XMAX + EXP GRPH 510 XMIN = XMIN - EXP GRPH 520 RETURN GRPH 530 ENTRY AXES GRPH 540 DELX = (XMAX - XMIN)/120.0 GRPH 550 DELY = (YMAX - YMIN)/54.0 GRPH 560 K = 0 GRPH 570 M = 0 GRPH 580 IF(XMIN .LT. 0.0 .AND. XMAX .GT. 0.0) GO TO 30 GRPH 590 GO TO 31 GRPH 600 30 K = (-XMIN/DELX) + 1.5 GRPH 610 DO 32 I = 1,55 GRPH 620 32 ITEM(I,K) = AIE GRPH 630 31 CONTINUE GRPH 640 IF(YMIN .LT. 0.0 .AND. YMAX .GT. 0.0) GO TO 33 GRPH 650 GO TO 34 GRPH 660 33 M = (YMAX/DELY) + 1.5 GRPH 670 DO 35 J = 1,121 GRPH 680 35 ITEM(M,J) = AMINUS GRPH 690 IF(K .EQ. 0 .OR. M .EQ. 0) GO TO 34 GRPH 700 ITEM(M,K) = ORIG GRPH 710 34 RETURN GRPH 720 ENTRY LOC1(X,Y,NPOI) GRPH 730 DELX = (XMAX - XMIN)/120.0 GRPH 740 DELY = (YMAX - YMIN)/54.0 GRPH 750 DO 50 II = 1,NPOI GRPH 760 I = (YMAX - Y(II))/DELY + 1.5 GRPH 770 J = (X(II) - XMIN)/DELX + 1.5 GRPH 780 IF(I.GT.55.OR.I.LT.1.OR.J.GT.121.OR.J.LT.1) GO TO 50 GRPH 790 DO 52 JJ = 1,10 GRPH 800 IF(ITEM(I,J) .EQ. HOLL(JJ)) GO TO 51 GRPH 810 52 CONTINUE GRPH 820 IF(ITEM(I,J).EQ.AIE.OR.ITEM(I,J).EQ.AMINUS.OR.ITEM(I,J).EQ.ORIG) GRPH 830 1 ITEM(I,J) = HOLL(2) GRPH 840 GO TO 50 GRPH 850 51 ITEM(I,J) = HOLL(JJ+1) GRPH 860 50 CONTINUE GRPH 870 RETURN GRPH 880 ENTRY LOC2(X,Y,NPOI,ISTART) GRPH 890 DELX = (XMAX - XMIN)/120.0 GRPH 900 DELY = (YMAX - YMIN)/54.0 GRPH 910 ISPN = ISTART + NPOI - 1 GRPH 920 DO 53 II = ISTART,ISPN GRPH 930 I = (YMAX - Y(II))/DELY + 1.5 GRPH 940 J = (X(II) - XMIN)/DELX + 1.5 GRPH 950 IF(I.GT.55.OR.I.LT.1.OR.J.GT.121.OR.J.LT.1) GO TO 53 GRPH 960 IF(ITEM(I,J).EQ.BLANK.OR.ITEM(I,J).EQ.AIE.OR.ITEM(I,J).EQ.AMINUS GRPH 970 1 .OR.ITEM(I,J).EQ.ORIG) GO TO 54 GRPH 980 ITEM(I,J) = AMULT GRPH 990 GO TO 53 GRPH1000 54 CONTINUE GRPH1010 IF(II .GT. 50) GO TO 55 GRPH1020 ITEM(I,J) = CHAR(II) GRPH1030 GO TO 53 GRPH1040 55 ITEM(I,J) = GT50 GRPH1050 53 CONTINUE GRPH1060 RETURN GRPH1070 ENTRY LOC3(X,Y,NPOI,POINT) GRPH1080 DELX = (XMAX - XMIN)/120.0 GRPH1090 DELY = (YMAX - YMIN)/54.0 GRPH1100 DO 56 II = 1,NPOI GRPH1110 I = (YMAX - Y(II))/DELY + 1.5 GRPH1120 J = (X(II) - XMIN)/DELX + 1.5 GRPH1130 IF(I.GT.55.OR.I.LT.1.OR.J.GT.121.OR.J.LT.1) GO TO 56 GRPH1140 IF(ITEM(I,J).EQ.BLANK.OR.ITEM(I,J).EQ.AIE.OR.ITEM(I,J).EQ.AMINUS GRPH1150 1 .OR. ITEM(I,J) .EQ. ORIG .OR. ITEM(I,J) .EQ. POINT) GO TO 57 GRPH1160 ITEM(I,J) = AMULT GRPH1170 GO TO 56 GRPH1180 57 ITEM(I,J) = POINT GRPH1190 56 CONTINUE GRPH1200 RETURN GRPH1210 ENTRY PLOT GRPH1220 SMALL(1) = XMIN GRPH1230 DO 120 I = 1,12 GRPH1240 120 SMALL(I+1) = SMALL(I) + 10.0 * DELX GRPH1250 VALUE = YMAX + DELY GRPH1260 WRITE(6,301) GRPH1270 301 FORMAT(9X,' +.........+.........+.........+.........+.........+...GRPH1280 1......+.........+.........+.........+.........+.........+.........GRPH1290 2+') GRPH1300 DO 330 I = 1,55 GRPH1310 VALUE = VALUE - DELY GRPH1320 II = I - 1 GRPH1330 L = II + 6 GRPH1340 IF(L/6*6 - L) 320,310,320 GRPH1350 310 A = PLUS GRPH1360 WRITE(6,302) VALUE,A,(ITEM(I,J), J=1,121),A GRPH1370 302 FORMAT(1X,F8.3,123A1) GRPH1380 GO TO 330 GRPH1390 320 A = DD GRPH1400 WRITE(6,303) A,(ITEM(I,J), J=1,121),A GRPH1410 303 FORMAT(9X,123A1) GRPH1420 330 CONTINUE GRPH1430 WRITE(6,301) GRPH1440 WRITE(6,304) (SMALL(K), K=1,12) GRPH1450 304 FORMAT(4X,12F10.3) GRPH1460 WRITE(6,305) SMALL(13) GRPH1470 305 FORMAT(1H+,123X,F8.2/) GRPH1480 RETURN GRPH1490 ENTRY IDENT GRPH1500 WRITE(6,400) (J, J=1,25) GRPH1510 400 FORMAT(/35X,'*****IDENTIFICATION KEY FOR PLOTS WITH IDENTIFIED POIGRPH1520 1NTS*****'//' PT #',25I5) GRPH1530 WRITE(6,401) (CHAR(I), I=1,25) GRPH1540 401 FORMAT(' CHAR',25(4X,A1)/) GRPH1550 WRITE(6,402) (J, J=26,50) GRPH1560 402 FORMAT(' PT #', 25I5) GRPH1570 WRITE(6,403) (CHAR(I), I=26,50) GRPH1580 403 FORMAT(' CHAR', 25(4X,A1)//) GRPH1590 WRITE(6,404) GT50, AMULT GRPH1600 404 FORMAT(' POINT NUMBERS ABOVE 50 IDENTIFIED AS ',A1/' MULTIPLE POIGRPH1610 1NTS IDENTIFIED AS ',A1/) GRPH1620 RETURN GRPH1630 END GRPH1640