C prefmap.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----+----@----+----@----+----@----+----@----+----@----+----@----+----@ C MAPREF - PREFERENCE MAPPING VIA A GENERALIZATION OF COOMBS CCNM 100 C UNFOLDING MODEL. CCNM 110 C BY CARROLL AND CHANG, FEBRUARY, 1969 CCNM 120 REAL*8 WZ(21,21) DIMENSION TIT(20) DIMENSION PAL(50),BNUB(100),NUMB(5),IX1(5),IY1(5) DIMENSION D(100),S(100,50),RSQ(50,4) DIMENSION X(100,5),P(5,50),T(5,5,50),SGN(5,50),XST(100,5) ,PST(5) CCNM 180 DIMENSION U(5,5),ROOT(5) , ROOT1(5,50), D1(50,100) DIMENSION FMTS(20),SIGN(5),FMTI(5),SDHS(100) CCNM 190 DIMENSION DHAT(100),SCR(100),SHAT(100,50),KC(100) CCNM 200 DIMENSION W(100,21),WA(21,21),WI(975) CCNM 210 DIMENSION WW(21,100),B(21),R(5,5),SS(100),LBLOCK(100) CCNM 220 COMMON WZ,N,K,XN,XK,NS,RMAX,IPH,NSUB,IS,NV,K1 COMMON D,S,RSQ,X,P,T,SGN,XST,PST,U,ROOT CCNM 240 COMMON PAL,NUMB,BNUB,NM COMMON ROOT1,D1 CCNM 260 COMMON /WORKK/ LWORK(21), MWORK(21) C PARAMETER DEFINITION CCNM 270 C N - NUMBER OF STIMULI (MAX. = 100) CCNM 280 C K - NUMBER OF DIMENSIONS (MAX. = 5) CCNM 290 C NSUB - NUMBER OF SUBJECTS (MAX. = 49) CCNM 300 C ISV- 0, SMALLER VALUE IS THE MORE PREFERRED CCNM 310 C ISV- 1, LARGER VALUE IS THE MORE PREFERRED CCNM 320 C NORS - 1, NORMALIZE SCALE VALUE, MAKE LENGTH =1. CCNM 330 C NORS - 0, DO NOT NORMALIZE SCALE VALUE. CCNM 340 C IPS - STARTING PHASE (MINIMUM = 1, MAXIMUM = 4). CCNM 350 C IPE - ENDING PHASE (MINIMUM = 1, MAXIMUM = 4). CCNM 360 C IPLOT - 0, PLOT AVERAGE SUBJECT ONLY. CCNM 370 C IPLOT - 1, PLOT AVERAGE SUBJ, FOR EACH SUBJ DO FUNCTION PLOTS. CCNM 380 C IPLOT - 2, PLOT AVERAGE S, FOR EACH S DO FUNCTION AND IDEAL PLOTS.CCNM 390 C MDV - 1, DRAW ROTATION AND WEIGHTS FOR EACH IDEAL POINT. CCNM 400 C MDV - 0, DRAW IDEAL POINTS WITHOUT VECTORS. CCNM 410 C MLIP - 1, LABEL IDEAL POINTS WITH LETTERS. CCNM 420 C MLIP - 0, LABEL IDEAL POINTS WITH NUMBERS. CCNM 430 C MIDEN - 1, READ IN SYMBOLS TO LABEL STIMULI ON PLOTS. CCNM 440 C MIDEN - 0, THE PROGRAM GENERATES NUMBERS TO LABEL STIMULI ON PLOTSCCNM 450 C LFITSW - 0, LINEAR FIT. CCNM 460 C LFITSW - 1, MONOTONE FIT (REGULAR). CCNM 470 C LFITSW - 2, MONOTONE FIT (BLOCK MONOTONE WITH ORDERING WITHIN). CCNM 480 C LFITSW - 3, MONOTONE FIT (BLOCK MONOTONE WITH EQUALITY WITHIN). CCNM 490 C IAV - 0, AVERAGE SUBJ SCALE VALUES ARE COMPUTED IN IPS PHASE. CCNM 500 C IAV - 1, AVERAGE SUBJ SCALE VALUES ARE COMPUTED IN EACH PHASE. CCNM 510 C IRX- 0, X IS PUNCHED N BY K CCNM 520 C IRX-1, X IS PUNCHED K BY N CCNM 530 C MAXIT - (MONOTONE ONLY) MAXIMUM NUMBER OF ITERATIONS. CCNM 540 C IRWT - 1, WHEN IPS=3, READ IN WEIGHTS, ONE FOR EACH DIMENSION. CCNM 550 C IRWT - 0, DO NOT READ IN WEIGHTS. CCNM 560 C ISHAT - 0, AT THE BEGINNING OF A NEW PHASE, USE THE LAST SHAT FROMCCNM 570 C THE PREVIOUS PHASE. CCNM 580 C ISHAT - 1, AT THE BEGINNING OF A NEW PHASE, USE THE ORIGINAL S. CCNM 590 C VARIABLES. CCNM 600 C X - COORDINATES OF N STIMULI IN K DIMENSIONAL SPACE CCNM 610 C S(I,J) - SCALE VALUE ON THE ITH STIMULUS FOR THE JTH SUBJECT CCNM 620 C NL=100 CCNM 630 NL=100 500 READ116,N,K,NSUB,ISV,NORS,IPS,IPE,IPLOT,MDV,MLIP,MIDEN,LFITSW,IAV,CCNM 640 1IRX,MAXIT,IRWT,ISHAT CCNM 650 IF(N.EQ.0) GOTO998 CCNM 660 IF(MAXIT.EQ.0.AND.LFITSW.GT.0) MAXIT=15 CCNM 670 PRINT124 CCNM 680 IF(LFITSW.EQ.0) IAV=0 CCNM 690 READ118,CRIT CCNM 700 PRINT101 CCNM 710 PRINT114 CCNM 720 PRINT117,N,K,NSUB,ISV,NORS,IPS,IPE,IPLOT,MDV,MLIP,MIDEN,LFITSW,IAVCCNM 730 1,IRX,MAXIT,IRWT,ISHAT,CRIT CCNM 740 PRINT101 CCNM 750 READ(5,100) TIT CCNM 760 CALL READX(X,N,K,IRX,NL) CCNM 770 WRITE(6,104) TIT CCNM 780 LFITSW=LFITSW-1 CCNM 790 DO13I=1,K CCNM 800 13 NUMB(I)=I CCNM 810 C TEST IF IPS=3, IF SO SET OR READ IN WEIGHTS CCNM 820 IF(IPS.NE.3) GOTO42 CCNM 830 IF(IRWT.EQ.0) GOTO43 CCNM 840 READ(5,100) FMTS CCNM 850 READFMTS,(SIGN(I),I=1,K) CCNM 860 GOTO42 CCNM 870 43 DO44I=1,K CCNM 880 44 SIGN(I)=1. CCNM 890 42 CONTINUE CCNM 900 NS = NSUB + 1 CCNM 910 READ(5,100) FMTS CCNM 920 50 NV=K*NS*2 CCNM 930 XNS=NSUB CCNM 940 XN=N CCNM 950 XK=K CCNM 960 K1=K-1 CCNM 970 IPH=IPS CCNM 980 C CALL ZEROD(RSQ,50,4) CCNM 990 CALL ZEROD(RSQ,50,4,200) C CALL ZEROD(RSQ,4,4,16) 95 IF(IPH.GT.IPE) GOTO98 CCNM1000 GOTO(82,92,78,99),IPH CCNM1010 82 M=((K+1)*(K+2))/2 CCNM1020 GOTO81 CCNM1030 99 M=K+1 CCNM1040 GOTO81 CCNM1050 78 M=K+2 CCNM1060 GOTO79 CCNM1070 92 M=2*K+1 CCNM1080 C79 CALL ZEROT(T,5,5,10) CCNM1090 C79 CALL ZEROT(T,3,3,4,36) 79 CALL ZEROT(T,5,5,50,1250) 81 PRINT106,IPH CCNM1100 PRINT105 CCNM1110 DO97I=1,K CCNM1120 97 PRINT103,I,(X(J,I),J=1,N) CCNM1130 CALL COMPW(W,X,WW,WA,WI,M,SIGN) CCNM1140 CALL ZEROS(SS,N) CCNM1150 RMAX=0. CCNM1160 PRINT106,IPH CCNM1170 DO90IS=1,NS CCNM1180 IT=1 CCNM1190 LAST=0 CCNM1200 IF(IS.EQ.NS)GOTO141 CCNM1210 C SUBJECTS CCNM1220 PRINT113,IS CCNM1230 IF(IPH.EQ.IPS) GOTO30 CCNM1240 IF(ISHAT)132,131,132 CCNM1250 30 READFMTS,(S(JS,IS),JS=1,N) CCNM1260 IF(ISV.EQ.0)GOTO135 CCNM1270 DO137I=1,N CCNM1280 137 S(I,IS)=-S(I,IS) CCNM1290 C NORMALIZE S CCNM1300 135 IF(NORS)133,132,133 CCNM1310 133 CALL NORMS(S(1,IS),N,XN) CCNM1320 GOTO132 CCNM1330 C AVERAGE SUBJECT CCNM1340 141 PRINT115 CCNM1350 IF(IPH.EQ.IPS)GOTO85 CCNM1360 IF(IAV.EQ.0) GOTO131 CCNM1370 C OPTION IAV=1 - COMPUTE AVERAGE SUBJ. FROM THE CURRENT PHASE CCNM1380 85 DO28I=1,N CCNM1390 28 S(I,IS)=SS(I)/XNS CCNM1400 C PUT S IN SHAT CCNM1410 132 DO34I=1,N CCNM1420 34 SHAT(I,IS)=S(I,IS) CCNM1430 131 IF(IS.EQ.NS.OR.LFITSW.LT.0) GOTO134 CCNM1440 C SORT S FOR MONOTONE FIT CCNM1450 DO5I=1,N CCNM1460 SCR(I)=S(I,IS) CCNM1470 5 KC(I)=I CCNM1480 CALL SORT(SCR(1), N, KC, KC, KC, 1, 1) CCNM1490 C SET UP LBLOCK FOR BLOCK MONOTONE CCNM1500 IF(LFITSW.GT.0) CALL BMONO(SCR,LBLOCK,N) CCNM1510 C PRINT SHAT CCNM1520 134 PRINT108 CCNM1530 PRINT109,(SHAT(I,IS),I=1,N) CCNM1540 C COMPUTE BETA CCNM1550 IF(LFITSW.LT.0.OR.IS.EQ.NS) GOTO31 CCNM1560 PRINT111 CCNM1570 36 PRINT119,IT CCNM1580 31 DO11I=1,M CCNM1590 B(I)=0. CCNM1600 DO11J=1,N CCNM1610 11 B(I)=B(I)+WW(I,J)*SHAT(J,IS) CCNM1620 PRINT102,(B(I),I=1,M) CCNM1630 IF(IT.EQ.1) GOTO152 CCNM1640 DO6I=1,N CCNM1650 D(I)=0. CCNM1660 DO6J=1,M CCNM1670 6 D(I)=D(I)+B(J)*W(I,J) CCNM1680 GOTO155 CCNM1690 C COMPUTE B CCNM1700 152 DO150I=1,K CCNM1710 150 B(I)=B(I+1) CCNM1720 IF(IPH.EQ.4)GOTO84 CCNM1730 IF(IPH-2)33,32,35 CCNM1740 35 II=K+2 CCNM1750 GOTO32 CCNM1760 C COMPUTE R CCNM1770 33 CALL COMPR(R,B,WA,K,XK) CCNM1780 C COMPUTE P CCNM1790 C INVERT R (R IN WA) CCNM1800 C CALL MTINV(WA,21,21,WI,21,21,K,MM) CCNM1820 MM = 36 CALL ARRAY (2, K, K,21,21, WA, WA) CALL MINV( WA, K, DETR, LWORK, MWORK) CALL ARRAY (1, K, K,21,21, WA, WA) 53 DO51I=1,K CCNM1860 P(I,IS)=0. CCNM1870 DO52J=1,K CCNM1880 52 P(I,IS)=P(I,IS)+B(J)*WA(J,I) CCNM1890 51 P(I,IS)=-.5*P(I,IS) CCNM1900 KM=1 CCNM1910 CALL PAYA(3, R, 5, 5, WI, 975, K) KM = 100 KU = 25 CALL EIGEN(WI,U,K,0, KM, KU ) CALL ARRAY(1,K,K,5,5,U,U, KU) CALL PAYA(2, R, 5, 5, WI, 975, K) C COMPUTE T CCNM1930 CALL COMPT(R) CCNM1940 CALL COMPD(XST,PST,SGN(1,IS),SHAT(1,IS),RSQ(IS,IPH)) CCNM1950 PRINT110,RSQ(IS,IPH) CCNM1960 154 IF(IS.EQ.NS) GOTO71 CCNM1970 IF(IPLOT.LT.2) GOTO155 CCNM1980 71 IF(K1)153,155,153 CCNM1990 153 CONTINUE CCNM2000 GOTO155 CCNM2010 C PHASE 4 CCNM2020 84 SUMB=0. CCNM2030 DO136I=1,K CCNM2040 136 SUMB=SUMB+B(I)**2 CCNM2050 SUMB=SQRT(SUMB) CCNM2060 DO86I=1,K CCNM2070 86 P(I,IS)=B(I)/SUMB CCNM2080 GOTO87 CCNM2090 C PHASE 2 AND PHASE 3 CCNM2100 32 DO55I=1,K CCNM2110 IF(IPH-2)55,75,76 CCNM2120 76 ROOT(I)=B(II)*SIGN(I) CCNM2130 GOTO77 CCNM2140 75 II=I+K+1 CCNM2150 ROOT(I)=B(II) CCNM2160 77 IF(ROOT(I))56,57,57 CCNM2170 56 SGN(I,IS)=-1. CCNM2180 GOTO58 CCNM2190 57 SGN(I,IS)=1. CCNM2200 58 T(I,I,IS)=ABS(ROOT(I)) CCNM2210 IF(T(I,I,IS)-RMAX)55,55,59 CCNM2220 59 RMAX=T(I,I,IS) CCNM2230 55 CONTINUE CCNM2240 C COMPUTE P CCNM2250 DO60I=1,K CCNM2260 60 P(I,IS)=-.5*B(I)/ROOT(I) CCNM2270 87 CALL COMPD(X,P(1,IS),ROOT,SHAT(1,IS),RSQ(IS,IPH)) CCNM2280 PRINT110,RSQ(IS,IPH) CCNM2290 155 IF(IPH.EQ.4) GOTO156 CCNM2300 PRINT112 CCNM2310 GOTO157 CCNM2320 156 PRINT129 CCNM2330 157 PRINT109,(D(I),I=1,N) CCNM2340 21 IF(IS.EQ.NS) GOTO73 CCNM2350 IF(LFITSW.LT.0) GOTO173 CCNM2360 C REARRANGE D IN DHAT, FIND MONOTONE FIT OF DHAT AND STORE IN D CCNM2370 CALL MONO(DHAT,D,KC,N,LFITSW,SCR,LBLOCK) CCNM2380 C REARRANGED AND UNNORMALIZED SHAT IN D AND SHAT IN SCR CCNM2390 C NORMALIZE SHAT (IN SCR) IN ORIGINAL ORDER CCNM2400 CALL NORMS(SCR,N,XN) CCNM2410 IF(LAST.EQ.1) GOTO61 CCNM2420 C TEST CRITERION CCNM2430 ITEST=0 CCNM2440 C SHAT HAS SHAT OF PREVIOUS IT. AND SCR HAS SHAT OF PRESENT IT. CCNM2450 SUM=0. CCNM2460 DO2I=1,N CCNM2470 2 SUM=SUM+(SHAT(I,IS)-SCR(I))**2 CCNM2480 IF(SUM.GT.CRIT) ITEST=1 CCNM2490 IF(IT.EQ.1) GOTO61 CCNM2500 IF(ITEST.EQ.0) GOTO160 CCNM2510 IF(IT.LT.MAXIT) GOTO203 CCNM2520 PRINT121,MAXIT CCNM2530 ITEST=0 CCNM2540 GOTO161 CCNM2550 160 PRINT122 CCNM2560 161 LAST=1 CCNM2570 GOTO152 CCNM2580 C PLOT S VS. DSQ CCNM2590 C REARRANGE S IN SDHS CCNM2600 61 DO202I=1,N CCNM2610 J=KC(I) CCNM2620 202 SDHS(I)=S(J,IS) CCNM2630 IF(IPLOT.EQ.0) GOTO62 CCNM2640 FC = 0.0 CCNM2650 WRITE(6,725) 725 FORMAT(1H1) CALL PLOTXY(SDHS, DHAT, N, -1, -1, FC, FC, FC, FC ) CCNM2660 C REMOVED SOME CODING FOR PLOTTING CCNM2670 C PRINT SUBJECT CCNM2680 62 IF(LAST.EQ.0.OR.IPH.EQ.4) GOTO203 CCNM2690 64 PRINT101 CCNM2700 CALL PRSUB CCNM2710 PRINT101 CCNM2720 C PUT SCR IN SHAT CCNM2730 203 DO63I=1,N CCNM2740 63 SHAT(I,IS)=SCR(I) CCNM2750 IF(ITEST.EQ.0) GOTO172 CCNM2760 IT=IT+1 CCNM2770 GOTO36 CCNM2780 173 IF(IPLOT.EQ.0) GOTO174 CCNM2790 73 IT=0 CCNM2800 C REMOVED SOME CODING FOR PLOTTING CCNM2810 174 IF(IPH.EQ.4) GOTO172 CCNM2820 PRINT101 CCNM2830 CALL PRSUB CCNM2840 PRINT101 CCNM2850 172 CONTINUE CCNM2860 DO 162 I=1,K CCNM2870 162 ROOT1(I,IS) = ROOT(I) CCNM2880 DO 163 I=1,N CCNM2890 163 D1(IS,I) = D(I) CCNM2900 IF(IPH.EQ.IPS) GOTO151 CCNM2910 IF(LFITSW.LT.0) GOTO90 CCNM2920 IF(IAV.EQ.0) GOTO90 CCNM2930 151 DO27I=1,N CCNM2940 27 SS(I)=SS(I)+SHAT(I,IS) CCNM2950 90 CONTINUE CCNM2960 WRITE(6,725) IF(IPH.EQ.3.OR.IPH.EQ.4) CALL PR1SUB(IPLOT) IF(IPH.EQ.4)GOTO88 CCNM2980 IPH=IPH+1 CCNM2990 IF(IPH-3)91,94,12 CCNM3000 91 DO93I=1,N CCNM3010 DO93J=1,K CCNM3020 93 X(I,J)=XST(I,J) CCNM3030 94 DO96I=1,K CCNM3040 96 SIGN(I)=SGN(I,NS) CCNM3050 GOTO95 CCNM3060 C FOR PHASE 4, REVERSE THE PREFERENCES, SO HIGHER VALUE MEANS CCNM3070 C GREATER PREFERENCE. CCNM3080 12 DO22I=1,NS CCNM3090 DO22J=1,N CCNM3100 SHAT(J,I)=-SHAT(J,I) CCNM3110 22 S(J,I)=-S(J,I) CCNM3120 GOTO95 CCNM3130 88 CONTINUE CCNM3140 PRINT125 CCNM3150 PRINT126,(I,I=1,K) CCNM3160 DO201I=1,NSUB CCNM3170 201 PRINT127,I,(P(J,I),J=1,K) CCNM3180 PRINT128,(P(J,NS),J=1,K) CCNM3190 98 CALL CPR(RSQ,XST,N,K,NS,XN,IPS,IPE) CCNM3200 GOTO500 CCNM3210 998 CONTINUE CCNM3220 100 FORMAT(20A4) CCNM3230 101 FORMAT(120H0******************************************************CCNM3240 1*****************************************************************)CCNM3250 102 FORMAT(5H BETA/(2X,10F13.5)) CCNM3260 103 FORMAT(I4,10F12.4/(4X,10F12.4)) CCNM3270 104 FORMAT(1H020A4) CCNM3280 105 FORMAT(9H0X MATRIX) CCNM3290 106 FORMAT(6H1PHASEI3) CCNM3300 107 FORMAT(21H0R CANNOT BE INVERTED) CCNM3310 108 FORMAT(27H0S (VECTOR OF SCALE VALUES)) CCNM3320 109 FORMAT(2X,10F13.5) CCNM3330 110 FORMAT(17H0R (CORRELATION)=F13.5) CCNM3340 111 FORMAT(32H0BEGIN ITERATION ON MONOTONE FIT) CCNM3350 112 FORMAT(4H DSQ) CCNM3360 113 FORMAT(8H0SUBJECTI4) CCNM3370 114 FORMAT(106H0 N K NSUB ISV NORS IPS IPE IPLOT MDV MLIPCCNM3380 1 MIDEN LFITSW IAV IRX MAXIT IRWT ISHAT CRIT) CCNM3390 115 FORMAT(16H0AVERAGE SUBJECT) CCNM3400 116 FORMAT(19I4) CCNM3410 117 FORMAT(4X,2I3,3I5,3I6,I5,2I6,2I7,I5,I6,2I7,F10.5) CCNM3420 118 FORMAT(10F7.2) CCNM3430 119 FORMAT(10H0ITERATIONI3) CCNM3440 121 FORMAT(28H0REACHED MAXIMUM ITERATION =I4) CCNM3450 122 FORMAT(36H0END OF ITERATION, REACHED CRITERION) CCNM3460 124 FORMAT(1H1//57H0MDSCALING VIA A GENERALIZATION OF COOMBS UNFOLDINGCCNM3470 1 MODEL) CCNM3480 125 FORMAT(44H0DIRECTION COSINES OF FITTED SUBJECT VECTORS//21X,9HDIMECCNM3490 1NSION) CCNM3500 126 FORMAT(10H0 SUBJECTI5,9I9) CCNM3510 127 FORMAT(I7,3X,10F9.4) CCNM3520 128 FORMAT(4X,3HAVR,3X,10F9.4) CCNM3530 129 FORMAT(33H0PROJECTIONS ON THE FITTED VECTOR) CCNM3540 999 STOP CCNM3550 END CCNM3560 SUBROUTINE COMPW(W,X,WW,WA,WI,M,SIGN) COPW 100 REAL*8 WZ(21,21) DIMENSION W(100,21),X(100,5),WI( 975),WW(21,100),WA(21,21),SIGN(5)COPW 110 COMMON WZ,N,K,XN,XK,NS,RMAX,IPH,NSUB,IS,NV,K1 COMMON /WORKK/ LWORK(21), MWORK(21) C COMPUTE W MATRIX COPW 130 DO5I=1,N COPW 140 5 W(I,1)=1. COPW 150 DO6J=1,K COPW 160 JJ=J+1 COPW 170 JK=JJ+K COPW 180 DO6I=1,N COPW 190 W(I,JJ)=X(I,J) COPW 200 6 W(I,JK)=X(I,J)**2 COPW 210 GOTO(70,71,72,71),IPH COPW 220 70 DO7I=1,N COPW 230 K2=2*K+1 COPW 240 DO7J=1,K1 COPW 250 JJ=J+1 COPW 260 DO7JW=JJ,K COPW 270 K2=K2+1 COPW 280 7 W(I,K2)=X(I,J)*X(I,JW) COPW 290 GOTO71 COPW 300 72 DO20I=1,N COPW 310 SUM=0. COPW 320 DO21J=1,K COPW 330 21 SUM=SUM+SIGN(J)*X(I,J)**2 COPW 340 20 W(I,M)=SUM COPW 350 71 DO8I=1,M COPW 360 DO8J=1,M COPW 370 WA(I,J)=0. COPW 380 DO8JM=1,N COPW 390 8 WA(I,J)=WA(I,J)+W(JM,I)*W(JM,J) COPW 400 C CALL MTINV(WA,21,21,WI,21,21,M,MX) COPW 420 CALL ARRAY (2, M, M,21,21, WA, WA) CALL MINV( WA, M, DETR, LWORK, MWORK) CALL ARRAY (1, M, M,21,21, WA, WA) 9 DO10I=1,M COPW 470 DO10J=1,N COPW 480 WW(I,J)=0. COPW 490 DO10JM=1,M COPW 500 10 WW(I,J)=WW(I,J)+WA(I,JM)*W(J,JM) COPW 510 RETURN COPW 520 END COPW 530 SUBROUTINE COMPR(R,B,WA,K,XK) COPR 100 DIMENSION R(5,5),B(21),WA(21,21) COPR 110 33 DO12I=1,K COPR 120 DO12J=1,I COPR 130 IF(I-J)14,15,14 COPR 140 15 II=I+K+1 COPR 150 R(I,I)=B(II) COPR 160 GOTO122 COPR 170 14 XI=J COPR 180 II=(XI+1.)*(XK-XI/2.) COPR 190 II=II+I+1 COPR 200 R(I,J)=.5*B(II) COPR 210 122 WA(I,J)=R(I,J) COPR 220 WA(J,I)=WA(I,J) COPR 230 R(J,I)=R(I,J) COPR 240 12 CONTINUE COPR 250 RETURN COPR 260 END COPR 270 SUBROUTINE PRSUB PRSU 100 REAL*8 WZ(21,21) DIMENSION D(100),S(100,50),RSQ(50,4) DIMENSION TIT(20) DIMENSION PAL(50),BNUB(100),U(5,5),ROOT(5),IX1(5),IY1(5),NUMB(5) PRSU 130 DIMENSION X(100,5),P(5,50),T(5,5,50),SGN(5,50),XST(100,5) ,PST(5) CCNM 180 COMMON WZ,N,K,XN,XK,NS,RMAX,IPH,NSUB,IS,NV,K1 COMMON D,S,RSQ,X,P,T,SGN,XST,PST,U,ROOT PRSU 160 COMMON PAL,NUMB,BNUB,NM PRINT 1,IS PRSU 180 IF(IPH-1) 100,100,200 PRSU 190 100 PRINT 10 PRSU 200 PRINT 11,(NUMB(I),I=1,K) PRSU 210 PRINT 12 PRSU 220 DO 5 I=1,K PRSU 230 5 PRINT 13,I,(U(J,I),J=1,K) PRSU 240 PRINT 30 PRSU 250 DO 6 I=1,K PRSU 260 6 PRINT32,I,(T(J,I,IS),J=1,K) PRSU 270 PRINT80 PRSU 280 DO81I=1,K PRSU 290 81 PRINT82,I,(XST(J,I),J=1,N) PRSU 300 PRINT 50 PRSU 310 PRINT21,(PST(I),I=1,K) PRSU 320 200 PRINT 40 PRSU 330 PRINT21,(P(I,IS),I=1,K) PRSU 340 PRINT 20 PRSU 350 PRINT 21,(ROOT(I),I=1,K) PRSU 360 1 FORMAT(8H0SUBJECTI6) PRSU 370 10 FORMAT(49H0DIRECTION COSINE OF NEW AXES WITH RESPECT TO OLD) PRSU 380 11 FORMAT(1H0,5X,3HOLDI8,4I15) PRSU 390 12 FORMAT(4H0NEW) PRSU 400 13 FORMAT(I8,1X,5F15.5) PRSU 410 20 FORMAT(24H0IMPORTANCES OF NEW AXES) PRSU 420 21 FORMAT(5X,5F15.5) PRSU 430 30 FORMAT(32H0COMPOSITE TRANSFORMATION MATRIX/) PRSU 440 32 FORMAT(I5,5F15.5) PRSU 450 40 FORMAT(37H0IDEAL POINT WITH RESPECT TO OLD AXES) PRSU 460 50 FORMAT(37H0IDEAL POINT WITH RESPECT TO NEW AXES) PRSU 470 80 FORMAT(14H0TRANSFORMED X) PRSU 480 82 FORMAT(I4,10F12.4/(4X,10F12.4)) PRSU 490 RETURN PRSU 500 END PRSU 510 SUBROUTINE CPR(RSQ,F,N,K,NS,XN,IPS,IPE) CPR 100 DIMENSION RSQ(50,4),F(100,5),DF(10),R(50,4) DIMENSION GF(10),KG(10),KD(10),XM(4) CPR 120 SMALL=0.001 GREAT=1000.0 DO3I=1,NS CPR 130 DO 3 J=1,4 3 R(I,J)=RSQ(I,J)**2 CPR 150 C COMPUTE AND PRINT CORRELATIONS AND F RATIO CPR 160 DO19I=1,10 CPR 170 DF(I)=0. CPR 180 19 GF(I)=0. CPR 190 CALL ZEROD(F,50,10,500) C CALL ZEROD(F,50,10) CPR 200 IPO=IPE-1 CPR 210 XM(1)=((K+1)*(K+2))/2 CPR 220 XM(2)=2*K+1 CPR 230 XM(3)=K+2 CPR 240 XM(4)=K+1 CPR 250 DO 1 I=1,4 GF(I)=XM(I)-1. CPR 270 1 DF(I)=XN-XM(I) CPR 280 DO 6 I=1,3 II=I+1 CPR 300 DO 6 J=II,4 IF(I.EQ.1) GOTO31 CPR 320 L=I+J+3 CPR 330 GOTO7 CPR 340 31 L=J-I+4 CPR 350 7 DF(L)=DF(I) CPR 360 6 GF(L)=XM(I)-XM(J) CPR 370 DO4I=1,NS CPR 380 DO 4 J=1,4 IF(1.0-R(I,J).GT.SMALL) GO TO 5 F(I,J) = GREAT GO TO 4 5 F(I,J) = R(I,J)/(1.-R(I,J))*DF(J)/GF(J) 4 CONTINUE DO2I=1,NS CPR 410 DO 2 I1=1,3 II=I1+1 CPR 430 DO 2 I2=II,4 IF(I1.EQ.1) GOTO11 CPR 450 L=I2+I1+3 CPR 460 GO TO 201 11 L=I2-I1+4 CPR 480 201 CONTINUE IF(1.-R(I,I1).GT.SMALL) GO TO 200 F(I,L) = GREAT GO TO 2 200 F(I,L)=(R(I,I1)-R(I,I2))/(1.-R(I,I1))*DF(L)/GF(L) CPR 490 2 CONTINUE PRINT100 CPR 500 PRINT101 CPR 510 DO10I=1,10 CPR 520 KG(I)=GF(I) CPR 530 10 KD(I)=DF(I) CPR 540 PRINT102,(KG(I),KD(I),I=1,4) CPR 550 NSUB=NS-1 CPR 560 XNSUB=NSUB CPR 570 DO20I=1,NSUB CPR 580 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 CPR 610 PRINT107,(KG(I),KD(I),I=5,10) CPR 620 DO21I=1,NSUB CPR 630 21 PRINT103,I,(F(I,J1),J1=5,10) CPR 640 PRINT104,(F(NS,J1),J1=5,10) CPR 650 C COMPUTE AND PRINT ROOT MEAN SQUARE CPR 660 PRINT105 CPR 670 DO 25 I=1,4 SUM=0. CPR 690 DO24J=1,NSUB CPR 700 24 SUM=SUM+RSQ(J,I)**2 CPR 710 RMSQ=SQRT(SUM/XNSUB) CPR 720 25 PRINT103,I,RMSQ CPR 730 WRITE(6,300) 300 FORMAT(//5X,' A F - VALUE OF 1000.0 IN THE ABOVE TABLE INDICATES ' X/5X,'A POSSIBLE DIVISION BY ZERO. I.E. R IS VERY CLOSE TO 1.00'//) 100 FORMAT(1H1/4X,19HCORRELATION (PHASE),36X,15HF RATIO (PHASE)) CPR 740 101 FORMAT(10X,2HR1,11X,2HR2,11X,2HR3,11X,2HR4,11X,2HF1,11X,2HF2,11X,2CPR 750 1HF3,11X,2HF4) CPR 760 102 FORMAT(3H DF,50X,4(I10,I3)/5H SUBJ) CPR 770 103 FORMAT(I3,2X,8F13.4) CPR 780 104 FORMAT(4H AVE1X,8F13.4) CPR 790 105 FORMAT(1H0/17H0ROOT MEAN SQUARE/7H PHASE) CPR 800 106 FORMAT(1H0//4X,23HF RATIO (BETWEEN PHASE)) CPR 810 107 FORMAT(10X,3HF12,10X,3HF13,10X,3HF14,10X,3HF23,10X,3HF24,10X,3HF34CPR 820 1/3H DF,I8,I3,5(I10,I3)/5H SUBJ) CPR 830 RETURN CPR 840 END CPR 850 SUBROUTINE COMPT(R) COPT 100 REAL*8 WZ(21,21) DIMENSION D(100),S(100,50),RSQ(50,4) DIMENSION R(5,5),R2(5) COPT 120 DIMENSION U(5,5),ROOT(5) COPT 130 DIMENSION X(100,5),P(5,50),T(5,5,50),SGN(5,50),XST(100,5) ,PST(5) CCNM 180 COMMON WZ,N,K,XN,XK,NS,RMAX,IPH,NSUB,IS,NV,K1 COMMON D,S,RSQ,X,P,T,SGN,XST,PST,U,ROOT COPT 160 DO16I=1,K COPT 170 ROOT(I)=R(I,I) COPT 180 IF(ROOT(I))17,18,18 COPT 190 17 SGN(I,IS)=-1. COPT 200 GOTO19 COPT 210 18 SGN(I,IS)=1. COPT 220 19 R2(I)=ABS(ROOT(I)) COPT 230 IF(R2(I).GT.RMAX) RMAX=R2(I) COPT 240 DO82J=1,K COPT 250 82 T(J,I,IS)=U(J,I)*SQRT(R2(I)) COPT 260 16 CONTINUE COPT 270 C COMPUTE PST COPT 280 DO20J=1,K COPT 290 PST(J)=0. COPT 300 DO20I=1,K COPT 310 20 PST(J)=PST(J)+P(I,IS)*T(I,J,IS) COPT 320 C COMPUTE X STAR COPT 330 DO21I=1,N COPT 340 DO21J=1,K COPT 350 XST(I,J)=0. COPT 360 DO21JJ=1,K COPT 370 21 XST(I,J)=XST(I,J)+X(I,JJ)*T(JJ,J,IS) COPT 380 RETURN COPT 390 END COPT 400 SUBROUTINE COMPD(XST,P,SGN,S,R) COPD 100 REAL*8 WZ(21,21) DIMENSION XST(100,5),P(5),SGN(5),D(100),S(100) COMMON WZ,N,K,XN,XK,NS,RMAX,IPH,NSUB,IS,NV,K1 COMMON D COPD 130 SUMD=0. COPD 140 SUMS=0. COPD 150 SDS=0. COPD 160 SSD=0. COPD 170 SSS=0. COPD 180 DO25I=1,N COPD 190 D(I)=0. COPD 200 DO26J=1,K COPD 210 IF(IPH.EQ.4)GOTO20 COPD 220 D(I)=D(I)+SGN(J)*(XST(I,J)-P(J))**2 COPD 230 GOTO26 COPD 240 20 D(I)=D(I)+XST(I,J)*P(J) COPD 250 26 CONTINUE COPD 260 SUMD=SUMD+D(I) COPD 270 SSD=SSD+D(I)**2 COPD 280 SUMS=SUMS+S(I) COPD 290 SSS=SSS+S(I)**2 COPD 300 25 SDS=SDS+D(I)*S(I) COPD 310 C COMPUTE CORRELATION COPD 320 DEN=(XN*SSD-SUMD**2)*(XN*SSS-SUMS**2) COPD 330 R=(XN*SDS-SUMD*SUMS)/SQRT(DEN) COPD 340 RETURN COPD 350 END COPD 360 SUBROUTINE NORMS(S,N,XN) NORM 100 DIMENSION S(100) NORM 110 SUM=0. NORM 120 SSQ=0. NORM 130 DO10I=1,N NORM 140 10 SUM=SUM+S(I) NORM 150 XMN=SUM/XN NORM 160 DO15I=1,N NORM 170 S(I)=S(I)-XMN NORM 180 15 SSQ=SSQ+S(I)**2 NORM 190 SV=1./SQRT(SSQ) NORM 200 DO16I=1,N NORM 210 16 S(I)=S(I)*SV NORM 220 RETURN NORM 230 END NORM 240 SUBROUTINE MONO(DHAT,D,KC,N,LFITSW,P,LBLOCK) MONO 100 DIMENSION DHAT(N),D(N),KC(N),P(N) ,P2(100) DIMENSION LBLOCK(100),LA(100),XLA(100) MONO 120 EQUIVALENCE(LA(1),XLA(1)) C REARRANGE D IN DHAT MONO 140 IF(LFITSW.EQ.1) GOTO14 MONO 150 NPASS=0 MONO 160 GOTO15 MONO 170 14 NPASS=1 MONO 180 DO5I=1,N MONO 190 LA(I)=LBLOCK(I) MONO 200 5 P(I)=I MONO 210 15 DO8I=1,N MONO 220 M=KC(I) MONO 230 8 DHAT(I)=D(M) MONO 240 CALL MFIT(DHAT,D,LA ,N,LFITSW,NPASS,P,P2) MONO 250 IF(LFITSW.NE.1) GOTO20 MONO 260 DO21I=1,N MONO 270 21 XLA(I)=D(I) MONO 280 CALL SORT(P(1), N, XLA, DHAT, DHAT, 2, 1 ) MONO 290 DO12I=1,N MONO 300 M=KC(I) MONO 310 12 P(M)=XLA(I) MONO 320 RETURN MONO 330 20 DO10I=1,N MONO 340 M=KC(I) MONO 350 10 P(M)=D(I) MONO 360 RETURN MONO 370 END MONO 380 SUBROUTINE MAMI(X,Y,N,TMAX,TMIN) MAMI 100 DIMENSION X(N),Y(N) CALL MM(X,N,XMAX,XMIN) MAMI 120 CALL MM(Y,N,YMAX,YMIN) MAMI 130 IF(YMAX.GT.XMAX) GOTO10 MAMI 140 TMAX=XMAX MAMI 150 GOTO13 MAMI 160 10 TMAX=YMAX MAMI 170 13 IF(YMIN.LT.XMIN) GOTO15 MAMI 180 TMIN=XMIN MAMI 190 GOTO20 MAMI 200 15 TMIN=YMIN MAMI 210 20 RETURN MAMI 220 END MAMI 230 SUBROUTINE BMONO(SCR,LBLOCK,N2) BMON 100 DIMENSION SCR(N2),LBLOCK(N2) BMON 110 L=1 BMON 120 LL=1 BMON 130 L1=1 BMON 140 DO17I=2,N2 BMON 150 II=I-1 BMON 160 IF(SCR(I).EQ.SCR(II))GOTO19 BMON 170 DO18J=L1,LL BMON 180 18 LBLOCK(J)=L BMON 190 L1=L1+L BMON 200 L=1 BMON 210 GOTO16 BMON 220 19 L=L+1 BMON 230 16 LL=LL+1 BMON 240 17 CONTINUE BMON 250 DO20J=L1,LL BMON 260 20 LBLOCK(J)=L BMON 270 RETURN BMON 280 END BMON 290 SUBROUTINE READX(X,N,K,IRX,NL) READ 100 DIMENSION X(NL,K),FMT(20) READ 110 READ(5,100) FMT READ 120 100 FORMAT(20A4) READ 130 IF(IRX.EQ.0) GOTO10 READ 140 DO6I=1,K READ 150 6 READFMT,(X(J,I),J=1,N) READ 160 GOTO20 READ 170 10 DO8I=1,N READ 180 8 READFMT,(X(I,J),J=1,K) READ 190 20 RETURN READ 200 END READ 210 SUBROUTINE MAXI(X,Y,Z,N) MAXI 100 C MAXI FINDS THE MAGNITUDE OF ARRAYS X AND Y EACH OF LENGTH N. MAXI 110 C ON RETURN ABSOLUTE MAXIMUM IS STORED IN Z. MAXI 120 DIMENSION X(N),Y(N) MAXI 130 Z=ABS(X(1)) MAXI 140 DO10I=1,N MAXI 150 AXO=ABS(X(I)) MAXI 160 IF(AXO.GT.Z) Z=AXO MAXI 170 AYO=ABS(Y(I)) MAXI 180 IF(AYO.GT.Z) Z=AYO MAXI 190 10 CONTINUE MAXI 200 RETURN MAXI 210 END MAXI 220 SUBROUTINE MM(X,N,XMAX,XMIN) MM 100 C FINDS THE MAXIMUM AND MINIMUM OF ARRAY X MM 110 DIMENSION X(100) MM 120 XMAX=X(1) MM 130 XMIN=XMAX MM 140 DO20I=2,N MM 150 IF(X(I)-XMAX)15,20,10 MM 160 10 XMAX=X(I) MM 170 GOTO20 MM 180 15 IF(X(I)-XMIN)16,20,20 MM 190 16 XMIN=X(I) MM 200 20 CONTINUE MM 210 RETURN MM 220 END MM 230 SUBROUTINE ZEROS(A,N) ZERO 100 DIMENSION A(N) ZERO 110 K=N ZERO 120 10 DO5II=1,K ZERO 130 5 A(II)=0. ZERO 140 GOTO50 ZERO 150 ENTRY ZEROD(A,I,J, N ) ZERO 160 K=I*J ZERO 170 GOTO10 ZERO 180 ENTRY ZEROT(A,I,J,M, N ) ZERO 190 K=I*J*M ZERO 200 GOTO10 ZERO 210 50 RETURN ZERO 220 END ZERO 230 SUBROUTINE MFIT(DIST,DHAT,LBLOCK,MM,LFITSW,NPASS,PAS1,PAS2) MFIT 100 C MFIT 110 C FIT PERFORMS WEIGHTED-LEAST-SQUARES MONOTONE REGRESSION MFIT 120 C THIS ROUTINE FINDS THE VALUES FOR DHAT WHICH ARE MONOTONIC MFIT 130 C AND WHICH BEST APPROXIMATE THE VALUES OF DIST, MFIT 140 C IN THE SENSE THAT THE SUM OF THE SQUARED DEVIATIONS, MFIT 150 C EACH ONE WEIGHTED BY THE CORRESPONDING VALUE IN WW , MFIT 160 C IS A MINIMUM. MFIT 170 C IT DIRECTLY USES SORT. MFIT 180 C MFIT 190 DIMENSION PAS1(MM),PAS2(MM),DIST(MM),DHAT(MM),LBLOCK(MM) MFIT 200 C MFIT 210 C FORM FIRST APPROXIMATION TO CORRECT PARTITON MFIT 220 C IF LFITSW=0, DO PLAIN REGRESSION MFIT 230 C IF LFITSW GT 0, DO BLOCK REGRESSION MFIT 240 C IF LFITSW=1, USE PRIMARY APPROACH. THUS WE SORT EACH MFIT 250 C BLOCK OF EQUAL VALUES OF DATA ACCORDING TO DIST VALUES. MFIT 260 C THEN WE CREATE PARTITON INTO BLOCKS OF SIZE 1 TO START. MFIT 270 MFIT 280 C IF LFITSW=2, USE SECONDARY APPROACH. THUS WE START WITHMFIT 290 C PARTITION INTO BLOCKS OF EQUAL DATA VALUES. MFIT 300 MFIT 310 C WITHIN EACH BLOCK OF WHATEVER SIZE, THE FIRST DHAT VALUEMFIT 320 C GIVES THE WEIGHTED TOTAL OF THE DIST VALUES IN THAT BLOCK, MFIT 330 C SIMILARLY, WITHIN EACH BLOCK, THE FIRST LBLOCK MFIT 340 C VALUE AND THE LAST LBLOCK VALUE BOTH GIVE THE SIZE OF THEMFIT 350 C BLOCK. MFIT 360 C BLOCKS OF SIZE 1 FORM A PARTIAL EXCEPTION. EVERYTHIN IS THE SAME MFIT 370 C EXCEPT THAT THE SUM OF THE W VALUES IS NOT FOUND IN THE MFIT 380 C LAST DHAT VALUE IN THE BLOCK. MFIT 390 C MFIT 400 IF(LFITSW.GT.0)GOTO100 MFIT 410 DO90M=1,MM MFIT 420 DHAT(M)=DIST(M) MFIT 430 90 LBLOCK(M)=1 MFIT 440 GOTO500 MFIT 450 100 MA=1 MFIT 460 110 K=LBLOCK(MA) MFIT 470 MB=MA+K-1 MFIT 480 GO TO ( 200, 300 ), LFITSW MFIT 490 C PRIMARY APPROACH MFIT 510 MFIT 520 200 IF(K-1) 9999, 220, 210 MFIT 530 C IF K=1, SAVE SORTING TIME MFIT 540 210 CONTINUE MFIT 550 IF(NPASS.EQ.0) CALL SORT(DIST(MA),K,PAS1(MA),PAS2(MA), MFIT 560 1 PAS2(MA), 0 , 1) MFIT 570 IF(NPASS.EQ.1) CALL SORT(DIST(MA),K,PAS1(MA),PAS2(MA), MFIT 580 1 PAS2(MA), 1 , 1) MFIT 590 IF(NPASS.EQ.2) CALL SORT(DIST(MA),K,PAS1(MA),PAS2(MA), MFIT 600 1 PAS2(MA), 2 , 1) MFIT 610 MFIT 620 C 'SORT' WILL SORT THE K ELEMENTS OF 'DIST' IN ALGEBRAIC MFIT 630 C ORDER. MFIT 640 C BECAUSE THE FINAL ARGUMENT IS ZERO, SORTING WILL BE MFIT 650 C ASCENDING OR DESCENDING ACCORDING TO THE VALUE OF SDSWIT MFIT 660 C SUPPLIED DURING THE CALL FOR SORT IN MDSCAL. MFIT 670 C THE ELEMENTS OF'IJ' AND 'DATA' WILL BE REARRANGED IN MFIT 680 C EXACTLY THE SAME ORDER AS THOSE OF 'DIST'. MFIT 690 C ALSO THE ELEMENTS OF 'WW'. MFIT 700 MFIT 710 220 DO 230 M=MA,MB MFIT 720 DHAT(M)=DIST(M) MFIT 730 LBLOCK(M)=1 MFIT 740 230 CONTINUE MFIT 750 GO TO 400 MFIT 760 MFIT 770 C SECONDARY APPROACH MFIT 780 MFIT 790 300 CONTINUE MFIT 800 TEMP1=0.0 MFIT 810 DO 310 M=MA,MB MFIT 820 TEMP1=TEMP1+DIST(M) MFIT 830 310 CONTINUE MFIT 840 DHAT(MA)=TEMP1 MFIT 850 C PROCEED TO NEXT BOCK. QUERY. IS IT THE LAST MFIT 870 400 MA=MA+K MFIT 890 IF(MA-MM-1) 110, 500, 9999 MFIT 900 C START MAIN COMPUTATION *************************************MFIT 920 500 MA=1 MFIT 940 510 LUD=2 MFIT 950 NSATIS=0 MFIT 960 520 K=LBLOCK(MA) MFIT 970 MB=MA+K-1 MFIT 980 WT=FLOAT(K) MFIT 990 550 DAV = DHAT(MA) / WT MFIT1000 GO TO ( 600, 700 ), LUD MFIT1010 C IS BLOCK DOWN-SATISFIED. IF NOT, MERGE MFIT1030 600 IF(MA-1) 9999, 630, 610 MFIT1050 610 MBD=MA-1 MFIT1060 KD=LBLOCK(MBD) MFIT1070 MAD=MBD-KD+1 MFIT1080 WTD=FLOAT(KD) MFIT1090 617 DAVD = DHAT(MAD) / WTD MFIT1100 IF( DAVD-DAV ) 630, 620, 620 MFIT1110 620 KNEW=K+KD MFIT1120 LBLOCK(MAD)=KNEW MFIT1130 LBLOCK(MB)=KNEW MFIT1140 DTONEW = DHAT(MAD) + DHAT(MA) MFIT1150 DHAT(MAD) = DTONEW MFIT1160 NSATIS=0 MFIT1170 MA=MAD MFIT1180 GO TO 800 MFIT1190 630 NSATIS=NSATIS+1 MFIT1200 GO TO 800 MFIT1210 C IS BLOCK UP-SATISFIED. IF NOT, MERGE MFIT1230 700 IF(MB-MM) 710, 730, 9999 MFIT1250 710 MAU=MB+1 MFIT1260 KU=LBLOCK(MAU) MFIT1270 MBU=MAU+KU-1 MFIT1280 WTU=FLOAT(KU) MFIT1290 717 DAVU = DHAT(MAU) / WTU MFIT1300 IF( DAV-DAVU ) 730, 720, 720 MFIT1310 720 KNEW=K+KU MFIT1320 LBLOCK(MA)=KNEW MFIT1330 LBLOCK(MBU)=KNEW MFIT1340 DTONEW = DHAT(MA) + DHAT(MAU) MFIT1350 DHAT(MA) = DTONEW MFIT1360 NSATIS=0 MFIT1370 GO TO 800 MFIT1380 730 NSATIS=NSATIS+1 MFIT1390 GO TO 800 MFIT1400 C PROCEED TO NEXT BLOCK IF READY. MFIT1420 800 LUD = 3-LUD MFIT1440 C QUERY. IS BLOCK BOTH UP AND DOWN SATISFIED. IF NOT, RETURMFIT1450 IF(NSATIS-1) 520, 520, 810 MFIT1460 C QUERY. IS THIS LAST BLOCK. IF NOT, GO ON TO NEXT BLOCK. MFIT1470 810 IF(MB-MM) 820, 900, 9999 MFIT1480 820 MA=MB+1 MFIT1490 GO TO 510 MFIT1500 C MAIN COMPUTATION COMPLETE. PLACE ANSWERS IN DHAT. MFIT1520 900 MA=1 MFIT1540 910 K=LBLOCK(MA) MFIT1550 MB=MA+K-1 MFIT1560 IF(K-1) 9999, 940, 920 MFIT1570 920 TEMP1=DHAT(MA)/FLOAT(K) MFIT1580 DO 930 M=MA,MB MFIT1590 930 DHAT(M)=TEMP1 MFIT1600 GO TO 945 MFIT1610 940 DHAT(MA) = DIST(MA) MFIT1620 945 MA = MB + 1 MFIT1630 IF(MA-MM-1) 910, 950, 9999 MFIT1640 950 RETURN MFIT1650 C TROUBLE EXIT MFIT1670 9999 PRINT99 MFIT1690 99 FORMAT(58H0MFIT DIAGNOSTIC. IMPOSSIBLE BRANCH TAKEN ON IF STATEMENMFIT1700 1T.) MFIT1710 STOP MFIT1720 END MFIT1740 SUBROUTINE PLOTXY( Z , Y ,NPOI,AXES,IDENT,XMAX,XMIN,YMAX,YMIN) PLOT 100 C PLOT 110 C ROUTINE TO GENERATE A SQUARE PLOT OF ARRAY X VS. Y. PLOT 120 C PLOT 130 C THE INPUT VECTORS HAVE NPOI ENTRIES (A NUMBER LESS THAN OR EQUAL TO PLOT 140 C THE LENGTH OF THE ARGUMENT ARRAYS IN THE CALLING PROGRAM. PLOT 150 C PLOT 160 C THE PLOTTING IS DONE ON TAPE -OUT-. PLOT 170 C PLOT 180 C IF -AXES- IS POSITIVE, AXES WILL APPEAR ON THE PLOT. PLOT 190 C IF -AXES- IS NEGATIVE, NO AXES WILL APPEAR. PLOT 200 C PLOT 210 C IF -XMAX- = -XMIN- THE PROGRAM WILL DETERMINE ITS OWN SCALE FOR PLOT 220 C THE X AXIS, AND SIMILARLY FOR THE Y AXIS IF -YMAX- = -YMIN-. PLOT 230 C IF -XMAX- IS NOT EQUAL TO -XMIN-, THE PROGRAM WILL USE THESE VALUES PLOT 240 C FOR THE UPPER AND LOWER BOUNDS OF THE X AXIS, AND SIMILARLY FOR THE PLOT 250 C Y AXIS. PLOT 260 C PLOT 270 C IF -IDENT- IS POSITIVE, THE PLOTTED POINTS WILL BE IDENTIFIED BY PLOT 280 C NUMBER. THERE ARE FIFTY (50) IDENTIFICATION CHARACTERS AVAILABLE.PLOT 290 C IF -IDENT- IS NEGATIVE, THE PLOTTED POINTS WILL APPEAR AS "X", AND PLOT 300 C MULTIPLE POINTS WILL BE IDENTIFIED BY THE NUMBER OF SUPERIMPOSED PLOT 310 C POINTS AT THAT LOCATION. ***NOTE*** THIS OPTION SHOULD BE USED WHEN PLOT 320 C NPOI SIGNIFICANTLY EXCEEDS 50. PLOT 330 C PLOT 340 C PLOT 350 C WRITTEN BY FORREST W. YOUNG, NOVEMBER 1965. PLOT 360 C VERSION 3, APRIL 1967 PLOT 370 C PLOT 380 C----------ADAPTED BY F.J. CARMONE, JR., 5/1/67. PLOT 390 C----------MODIFIED BY R.F. MCCRACKEN 11/15/67. PLOT 400 C PLOT 410 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCPLOT 420 C PLOT 430 REAL*8 WZ(21,21) DIMENSION Z(NPOI), Y(NPOI) DIMENSION SMALL(25),CHAR(50),HOLL(11) DIMENSION VALY(55),VECT(30) PLOT 450 REAL ITEM(55,121),N1,N2,N3 PLOT 470 INTEGER OUT,AXES PLOT 480 LOGICAL FLAG/.FALSE./,WIPE/.FALSE./ PLOT 490 COMMON WZ,N,KK1,XN,XK,NS,RMAX,IPH DATA CHAR/'1','2','3','4','5','6','7','8','9','A','B','C','D', PLOT 500 1'E','F','G','H','I','J','K','L','M','N','O','P','Q','R','S','T', PLOT 510 2'U','V','W','X','Y','Z','+','/','=','*','&','$','@','%','?','<', PLOT 520 3'(',')','"','#',' '/ PLOT 530 DATA HOLL/' ','X','2','3','4','5','6','7','8','9','M'/ PLOT 540 DATA BLANK/' '/,DD/'.'/,PLUS/'*'/,AIE/'\'/,AMINUS/'-'/,ORIG/'0'/, PLOT 550 1 AMULT/';'/,GT50/'>'/ PLOT 560 OUT = 6 PLOT 570 IF(FLAG) GO TO 400 PLOT 580 DO 115 I = 1,55 PLOT 590 DO 115 J=1,121 PLOT 600 115 ITEM(I,J) = BLANK PLOT 610 FLAG = .TRUE. PLOT 620 400 CONTINUE PLOT 630 IF(XMAX .NE. XMIN) GO TO 222 PLOT 640 XMAX = -1.0E7 PLOT 650 XMIN = +1.0E7 PLOT 660 DO 221 I = 1,NPOI PLOT 670 IF(Z(I) .GT. XMAX) XMAX = Z(I) PLOT 680 IF(Z(I) .LT. XMIN) XMIN = Z(I) PLOT 690 221 CONTINUE PLOT 700 222 IF(YMAX .NE. YMIN) GO TO 224 PLOT 710 YMAX = -1.0E7 PLOT 720 YMIN = +1.0E7 PLOT 730 DO 223 I = 1,NPOI PLOT 740 IF(Y(I) .GT. YMAX) YMAX = Y(I) PLOT 750 IF(Y(I) .LT. YMIN) YMIN = Y(I) PLOT 760 223 CONTINUE PLOT 770 224 CONTINUE PLOT 780 SPANX = XMAX - XMIN PLOT 790 SPANY = YMAX - YMIN PLOT 800 DELX = SPANX/119.0 DELY = SPANY/54.0 PLOT 820 IF(AXES .LE. 0) GO TO 119 PLOT 830 K = 0 PLOT 840 M = 0 PLOT 850 IF(XMIN .LT. 0 .AND. XMAX .GT. 0) GO TO 200 PLOT 860 GO TO 202 PLOT 870 200 K = (-XMIN/DELX) + 1.5 PLOT 880 DO 201 I = 1,55 PLOT 890 201 ITEM(I,K) = AIE PLOT 900 202 CONTINUE PLOT 910 IF(YMIN .LT. 0 .AND. YMAX .GT. 0) GO TO 203 PLOT 920 GO TO 119 PLOT 930 203 M = (-YMIN/DELY) + 1.5 PLOT 940 M = 56 - M PLOT 950 DO 204 J = 1,119 204 ITEM(M,J) = AMINUS PLOT 970 IF(K .EQ. 0 .OR. M .EQ. 0) GO TO 119 PLOT 980 ITEM(M,K) = ORIG PLOT 990 119 CONTINUE PLOT1000 VALUE = YMAX + DELY PLOT1010 SMALL(1) = XMIN PLOT1020 DO 120 I=1,24 PLOT1030 120 SMALL(I+1) = SMALL(I) + 5.0*DELX PLOT1040 IF(IDENT .GT. 0) GO TO 140 PLOT1050 DO 135 II=1,NPOI PLOT1060 I = (YMAX - Y(II))/DELY + 1.5 PLOT1070 J = (Z(II) - XMIN)/DELX + 1.5 PLOT1080 IF(I.GT.55.OR.I.LT.1.OR.J.GT.119.OR.J.LT.1) GO TO 135 P DO 134 JJ=1,10 PLOT1100 IF(ITEM(I,J) .EQ. HOLL(JJ)) GO TO 136 PLOT1110 134 CONTINUE PLOT1120 IF(ITEM(I,J) .EQ. AIE .OR. ITEM(I,J) .EQ. AMINUS .OR. ITEM(I,J) PLOT1130 1 .EQ. ORIG) ITEM(I,J) = HOLL(2) PLOT1140 GO TO 135 PLOT1150 136 ITEM(I,J) = HOLL(JJ+1) PLOT1160 135 CONTINUE PLOT1170 GOTO355 PLOT1180 140 CONTINUE PLOT1190 DO 137 II = 1,NPOI PLOT1200 I = (YMAX - Y(II))/DELY + 1.5 PLOT1210 J = (Z(II) - XMIN)/DELX + 1.5 PLOT1220 IF(I.GT.55.OR.I.LT.1.OR.J.GT.119.OR.J.LT.1) GO TO 137 IF(ITEM(I,J) .EQ. BLANK .OR. ITEM(I,J) .EQ. AIE .OR. ITEM(I,J) PLOT1240 1 .EQ. AMINUS .OR. ITEM(I,J) .EQ. ORIG) GO TO 138 PLOT1250 ITEM(I,J) = AMULT PLOT1260 GO TO 137 PLOT1270 138 CONTINUE PLOT1280 IF(II .GT. 50) GO TO 139 PLOT1290 ITEM(I,J) = CHAR(II) PLOT1300 GO TO 137 PLOT1310 139 ITEM(I,J) = GT50 PLOT1320 137 CONTINUE PLOT1330 141 CONTINUE PLOT1340 355 WRITE(OUT,3001) PLOT1350 3001 FORMAT(10X,122H.*....*....*....*....*....*....*....*....*....*....PLOT1360 1*....*....*....*....*....*....*....*....*....*....*....*....*....*PLOT1370 2....*) PLOT1380 DO 330 I = 1,55 PLOT1390 VALUE = VALUE - DELY PLOT1400 A = BLANK PLOT1410 B = DD PLOT1420 L = I + 2 PLOT1430 IF(L/5*5-L) 330,310,330 PLOT1440 310 B = PLUS PLOT1450 IF(L/2*2-L) 330,320,330 PLOT1460 320 A = PLUS PLOT1470 330 WRITE(OUT,3300) VALUE,A,B,(ITEM(I,J), J=1,119), B, A 3300 FORMAT(1X,F8.3,123A1) PLOT1490 WRITE(OUT,3001) PLOT1500 WRITE(OUT,3301) DD,(SMALL(I),DD, I=2,24,2) PLOT1510 3301 FORMAT(11X,A1,12(F9.4,A1)) PLOT1520 WRITE(OUT,3302) (SMALL(I), I=1,23,2) PLOT1530 3302 FORMAT(6X,12F10.4) PLOT1540 IF(IDENT .LE. 0) GO TO 150 PLOT1550 WRITE(OUT,151) (J, J=1,25) PLOT1560 151 FORMAT(//40X,'POINT IDENTIFICATION KEY'//1X,'PT #',25I5/) PLOT1570 WRITE(OUT,152) (CHAR(I), I=1,25) PLOT1580 152 FORMAT(1X,'CHAR',25(4X,A1)///) PLOT1590 WRITE(OUT,153) (J, J=26,50) PLOT1600 153 FORMAT(1X,'PT #',25I5/) PLOT1610 WRITE(OUT,154) (CHAR(I), I=26,50) PLOT1620 154 FORMAT(1X,'CHAR',25(4X,A1)) PLOT1630 WRITE(OUT,155) PLOT1640 155 FORMAT(////1X,'ALL POINT NUMBERS ABOVE 50 ARE DESIGNATED BY THE SYPLOT1650 1MBOL >'//1X,'MULTIPLE POINTS AT THE SAME LOCATION ARE DESIGNATEDPLOT1660 2 BY THE SYMBOL ;') PLOT1670 N1= GT50 N2= GT50 N3= GT50 IF(N.GT.50) GO TO 3303 N1= CHAR(N) IF(N.EQ.50) GO TO 3303 N2= CHAR(N+1) IF((N+NS).GE.50) GO TO 3303 N3= CHAR(N+NS) 3303 WRITE(6,156)N1,N2,N3 156 FORMAT(//5X,'CHARACTERS 1 TO ',A4,'ARE STIMULI AND CHARACTERS 'PLOT1730 1,A4,'TO ',A4,'ARE IDEAL POINTS'//) PLOT1740 150 CONTINUE PLOT1750 IF(AXES .LE. 0) GO TO 302 PLOT1760 K = (-XMIN/DELX) + 1.5 PLOT1770 DO 300 I = 1,55 PLOT1780 300 ITEM(I,K) = BLANK PLOT1790 M = (-YMIN/DELY) + 1.5 PLOT1800 M = 56 - M PLOT1810 DO 301 J = 1,121 PLOT1820 301 ITEM(M,J) = BLANK PLOT1830 302 CONTINUE PLOT1840 DO 303 II = 1,NPOI PLOT1850 I = (YMAX - Y(II))/DELY + 1.5 PLOT1860 J = (Z(II) - XMIN)/DELX + 1.5 PLOT1870 IF(I.GT.55.OR.I.LT.1.OR.J.GT.121.OR.J.LT.1) GO TO 303 PLOT1880 ITEM(I,J) = BLANK PLOT1890 303 CONTINUE PLOT1900 RETURN PLOT1910 END PLOT1920 SUBROUTINE PR1SUB(MSPLOT) PR1S 100 REAL*8 WZ(21,21) DIMENSION D(100),S(100,50),RSQ(50,4) PR1S 110 DIMENSION PAL(50),BNUB(100),NUMB(5 ),IX1(5),IY1(5) PR1S 120 DIMENSION U(5,5),ROOT(5),ROOT1(5,50),D1(50,100) DIMENSION X(100,5),P(5,50),T(5,5,50),SGN(5,50),XST(100,5),PST(5) PR1S 140 DIMENSION GRAPHX(150),GRAPHY(150),LWORK(150) COMMON WZ,N,K,XN,XK,NS,RMAX,IPH,NSUB,IS,NV,K1 COMMON D,S,RSQ,X,P,T,SGN,XST,PST,U,ROOT PR1S 170 COMMON PAL,NUMB,BNUB,NM COMMON ROOT1,D1 PR1S 190 98 WRITE(6,304) PR1S 200 304 FORMAT(1H1,40X,'STIMULI COORDINATES'//) PR1S 210 WRITE(6,305)(NUMB(I),I=1,K) PR1S 220 305 FORMAT(1H0,9X,'DIMENSION',I8,4I15) PR1S 230 WRITE(6,306) PR1S 240 306 FORMAT(1H0,'STIMULI') PR1S 250 DO307I=1,N PR1S 260 307 WRITE(6,308)I,(X(I,J),J=1,K) PR1S 270 308 FORMAT(1H0,I8,8X,5F15.5) PR1S 280 IF(IPH.EQ.4)GOTO321 PR1S 290 WRITE(6,309) PR1S 300 309 FORMAT(1H0,//40X,'COORDINATES OF IDEAL POINTS'//) PR1S 310 WRITE(6,305)(NUMB(I),I=1,K) PR1S 320 WRITE(6,310) PR1S 330 310 FORMAT(1H0,'SUBJECTS') PR1S 340 DO311J=1,NS PR1S 350 311 WRITE(6,308)J,(P(I,J),I=1,K) PR1S 360 WRITE(6,312)NS PR1S 370 312 FORMAT(1H0,'SUBJECT ',I3,' IS THE AVERAGE SUBJECT') PR1S 380 WRITE(6,314) 314 FORMAT(1H0,40X,'WEIGHTS OF AXES'//) PR1S 400 WRITE(6,305)(NUMB(I),I=1,K) PR1S 410 WRITE(6,310) PR1S 420 DO315J=1,NS PR1S 430 315 WRITE(6,308)J,(ROOT1(I,J),I=1,K) PR1S 440 WRITE(6,312)NS PR1S 450 321 WRITE(6,3007) 3007 FORMAT(1H1) GO TO 3000 IF(MSPLOT.NE.1.OR.IPH.EQ.4)GOTO325 PR1S 460 DO322IS=1,NS PR1S 470 DO323I=1,N PR1S 480 GRAPHX(I)=S(I,IS) PR1S 490 LWORK(I)=GRAPHX(I) PR1S 500 323 GRAPHY(I)=ABS(D1(IS,I)) PR1S 510 WRITE(6,300) PR1S 520 NPOI=N PR1S 530 XMAX=XMIN PR1S 540 YMAX=YMIN PR1S 550 CALLPLOTXY(GRAPHX,GRAPHY,NPOI,-1,-1,XMAX,XMIN,YMAX,YMIN) PR1S 560 322 WRITE(6,324)(LWORK(I),I=1,N) PR1S 570 324 FORMAT(1H0,15X,'RANK ORDER OF STIMULI IS ',20(2X,I2)) PR1S 580 325 IF(IPH.EQ.4)GOTO326 PR1S 590 WRITE(6,316) PR1S 600 316 FORMAT(1H0,40X,'SQUARED DISTANCES BETWEEN IDEAL POINT AND ', PR1S 610 1'STIMULI'//) PR1S 620 DO318J=1,NS PR1S 630 DO318I=1,N,7 PR1S 640 N2=I+6 PR1S 650 IF(N2.GT.N)N2=N PR1S 660 WRITE(6,320)(NUMB(N1),N1=I,N2) PR1S 670 IF(I.EQ.1.AND.J.EQ.1)WRITE(6,319) PR1S 680 319 FORMAT(1H0,'SUBJECT IDEAL') PR1S 690 WRITE(6,317)J,(D1(J,N1),N1=I,N2) PR1S 700 320 FORMAT(1H0,14X,'STIMULUS',I10,6I15) PR1S 710 317 FORMAT(1H0,I8,13X,7F15.5) PR1S 720 318 CONTINUE PR1S 730 WRITE(6,312)NS PR1S 740 326 IF(IPH.LT.3.OR.K.GT.2)RETURN PR1S 750 WRITE(6,300) PR1S 760 300 FORMAT(1H1) PR1S 770 3000 CONTINUE CALL ZEROS(GRAPHY, 150) CALL ZEROS(GRAPHX, 150) NPOI=N+NS PR1S 800 IRANGE=1 PR1S 810 DO301I=1,N PR1S 820 GRAPHX(I)=X(I,1) PR1S 830 GRAPHY(I)=X(I,2) PR1S 840 301 IF(GRAPHX(I).GT.3..OR.GRAPHY(I).GT.3.)IRANGE=2 PR1S 850 DO303I=1,NS PR1S 860 GRAPHX(N+I)=P(1,I) PR1S 870 GRAPHY(N+I)=P(2,I) PR1S 880 303 IF(GRAPHX(I).GT.3..OR.GRAPHY(I).GT.3.)IRANGE=2 PR1S 890 IF(IPH.EQ.3.OR.IPH.EQ.4) GO TO 3001 XMAX=+4.01**IRANGE PR1S 910 XMIN=-4.02**IRANGE PR1S 920 YMAX=+3.01**IRANGE PR1S 930 YMIN=-3.02**IRANGE PR1S 940 GO TO 3002 3001 XMAX=2.0 XMIN=-2.0 YMAX=1.5 YMIN=-1.5 3002 CALLPLOTXY(GRAPHX,GRAPHY,NPOI,+1,+1,XMAX,XMIN,YMAX,YMIN) RETURN PR1S 970 END PR1S 980 C MINV 001 C ..................................................................MINV 002 C MINV 003 C SUBROUTINE MINV MINV 004 C MINV 005 C PURPOSE MINV 006 C INVERT A MATRIX MINV 007 C MINV 008 C USAGE MINV 009 C CALL MINV(A,N,D,L,M) MINV 010 C MINV 011 C DESCRIPTION OF PARAMETERS MINV 012 C A - INPUT MATRIX, DESTROYED IN COMPUTATION AND REPLACED BY MINV 013 C RESULTANT INVERSE. MINV 014 C N - NUMBER OF ROWS AND COLUMNS MINV 015 C D - RESULTANT DETERMINANT MINV 016 C L - WORK VECTOR OF LENGTH N MINV 017 C M - WORK VECTOR OF LENGTH N MINV 018 C MINV 019 C REMARKS MINV 020 C MATRIX A MUST BE A GENERAL MATRIX MINV 021 C MINV 022 C SUBROUTINES USED BY THIS SUBROUTINE MINV 023 C NONE MINV 024 C MINV 025 C METHOD MINV 026 C THE STANDARD GAUSS-JORDAN METHOD IS USED. THE DETERMINANT MINV 027 C IS ALSO CALCULATED. A DETERMINANT OF ZERO INDICATES THAT MINV 028 C THE MATRIX IS SINGULAR. MINV 029 C MINV 030 C ..................................................................MINV 031 C MINV 032 SUBROUTINE MINV(A,N,D,L,M, MM ) DIMENSION A(MM), L(21), M(21) C MINV 035 C SEARCH FOR LARGEST ELEMENT MINV 036 C MINV 037 D=1.0 MINV 038 NK=-N MINV 039 DO 80 K=1,N MINV 040 NK=NK+N MINV 041 L(K)=K MINV 042 M(K)=K MINV 043 KK=NK+K MINV 044 BIGA=A(KK) MINV 045 DO 20 J=K,N MINV 046 IZ=N*(J-1) MINV 047 DO 20 I=K,N MINV 048 IJ=IZ+I MINV 049 IF(ABS(BIGA)-ABS(A(IJ))) 15,20,20 MINV 050 15 BIGA=A(IJ) MINV 051 L(K)=I MINV 052 M(K)=J MINV 053 20 CONTINUE MINV 054 C MINV 055 C INTERCHANGE ROWS MINV 056 C MINV 057 J=L(K) MINV 058 IF(L(K)-K) 35,35,25 MINV 059 25 KI=K-N MINV 060 DO 30 I=1,N MINV 061 KI=KI+N MINV 062 HOLD=-A(KI) MINV 063 JI=KI-K+J MINV 064 A(KI)=A(JI) MINV 065 30 A(JI) =HOLD MINV 066 C MINV 067 C INTERCHANGE COLUMNS MINV 068 C MINV 069 35 I=M(K) MINV 070 IF(M(K)-K) 45,45,38 MINV 071 38 JP=N*(I-1) MINV 072 DO 40 J=1,N MINV 073 JK=NK+J MINV 074 JI=JP+J MINV 075 HOLD=-A(JK) MINV 076 A(JK)=A(JI) MINV 077 40 A(JI) =HOLD MINV 078 C MINV 079 C DIVIDE COLUMN BY MINUS PIVOT MINV 080 C MINV 081 45 DO 55 I=1,N MINV 082 IF(I-K) 50,55,50 MINV 083 50 IK=NK+I MINV 084 A(IK)=A(IK)/(-A(KK)) MINV 085 55 CONTINUE MINV 086 C MINV 087 C REDUCE MATRIX MINV 088 C MINV 089 DO 65 I=1,N MINV 090 IK=NK+I MINV 091 IJ=I-N MINV 092 DO 65 J=1,N MINV 093 IJ=IJ+N MINV 094 IF(I-K) 60,65,60 MINV 095 60 IF(J-K) 62,65,62 MINV 096 62 KJ=IJ-I+K MINV 097 A(IJ)=A(IK)*A(KJ)+A(IJ) MINV 098 65 CONTINUE MINV 099 C MINV 100 C DIVIDE ROW BY PIVOT MINV 101 C MINV 102 KJ=K-N MINV 103 DO 75 J=1,N MINV 104 KJ=KJ+N MINV 105 IF(J-K) 70,75,70 MINV 106 70 A(KJ)=A(KJ)/A(KK) MINV 107 C MINV 109 75 CONTINUE MINV 108 C PRODUCT OF PIVOTS MINV 110 C MINV 111 D=D*A(KK) MINV 112 C MINV 113 C REPLACE PIVOT BY RECIPROCAL MINV 114 C MINV 115 A(KK)=1.0/A(KK) MINV 116 80 CONTINUE MINV 117 C MINV 118 C FINAL ROW AND COLUMN INTERCHANGE MINV 119 C MINV 120 K=N MINV 121 100 K=(K-1) MINV 122 IF(K) 150,150,105 MINV 123 105 I=L(K) MINV 124 IF(I-K) 120,120,108 MINV 125 108 JQ=N*(K-1) MINV 126 JR=N*(I-1) MINV 127 DO 110 J=1,N MINV 128 JK=JQ+J MINV 129 HOLD=A(JK) MINV 130 JI=JR+J MINV 131 A(JK)=-A(JI) MINV 132 110 A(JI) =HOLD MINV 133 120 J=M(K) MINV 134 IF(J-K) 100,100,125 MINV 135 125 KI=K-N MINV 136 DO 130 I=1,N MINV 137 KI=KI+N MINV 138 HOLD=A(KI) MINV 139 JI=KI-K+J MINV 140 A(KI)=-A(JI) MINV 141 130 A(JI) =HOLD MINV 142 GO TO 100 MINV 143 150 RETURN MINV 144 END SUBROUTINE PAYA( KLA, R, N, NA,WI, M, K) DIMENSION R(N, N), WI(M) KK=K*K PAYA6200 II=0 PAYA6210 IF(KLA.EQ.2)GO TO 14 PAYA6220 DO 11 J=1,K PAYA6230 DO 11 I=1,J PAYA6240 II=II+1 PAYA6250 WI(II)=R(I,J) PAYA6260 11 CONTINUE PAYA6270 GO TO 15 PAYA6280 14 DO 13 J=1,K PAYA6290 II=II+J PAYA6300 13 R(J,J)=WI(II) PAYA6310 15 RETURN PAYA6320 END PAYA6330 SUBROUTINE EIGEN(A,R,N,MV, KM, KU ) C EIGEN001 C SUBROUTINE EIGEN EIGEN002 C EIGEN003 C PURPOSE EIGEN004 C COMPUTE EIGENVALUES AND EIGENVECTORS OF A REAL SYMMETRIC EIGEN005 C MATRIX EIGEN006 C EIGEN007 C USAGE EIGEN008 C CALL EIGEN(A,R,N,MV) EIGEN009 C EIGEN010 C DESCRIPTION OF PARAMETERS EIGEN011 C A - ORIGINAL MATRIX (SYMMETRIC), DESTROYED IN COMPUTATION. EIGEN012 C RESULTANT EIGENVALUES ARE DEVELOPED IN DIAGONAL OF EIGEN013 C MATRIX A IN DECENDING ORDER. EIGEN014 C R - RESULTANT MATRIX OF EIGENVECTORS (STORED COLUMNWISE, EIGEN015 C IN SAME SEQUENCE AS EIGENVALUES) EIGEN016 C N - ORDER OF MATRICES A AND R EIGEN017 C MV- INPUT CODE EIGEN018 C 0 COMPUTE EIGENVALUES AND EIGENVECTORS EIGEN019 C 1 COMPUTE EIGENVALUES ONLY (R NEED NOT BE EIGEN020 C DIMENSIONED BUT MUST STILL APPEAR IN CALLING EIGEN021 C SEQUENCE) EIGEN022 C EIGEN023 C REMARKS EIGEN024 C ORIGINAL MATRIX A MUST BE REAL SYMMETRIC (STORAGE MODE=1) EIGEN025 C MATRIX A CANNOT BE IN THE SAME LOCATION AS MATRIX R EIGEN026 C EIGEN027 C SUBROUTINES USED BY THIS SUBROUTINE EIGEN028 C NONE EIGEN029 C EIGEN030 C METHOD EIGEN031 C DIAGONALIZATION METHOD ORIGINATED BY JACOBI AND ADAPTED EIGEN032 C BY VON NEUMANN FOR LARGE COMPUTERS AS FOUND IN MATHEMATICAL EIGEN033 C METHODS FOR DIGITAL COMPUTERS, EDITED BY A. RALSTON AND EIGEN034 C H.S. WILF, JOHN WILEY AND SONS, NEW YORK, 1962, CHAPTER 7 EIGEN035 C EIGEN036 C ..................................................................EIGEN037 C EIGEN038 C EIGEN041 DIMENSION A(KM), R(KU) C GENERATE IDENTITY MATRIX EIGEN042 C EIGEN043 IF(MV-1) 22,14,22 EIGEN044 22 IQ=-N EIGEN045 DO 1 J=1,N EIGEN046 IQ=IQ+N EIGEN047 DO 1 I=1,N EIGEN048 IJ=IQ+I EIGEN049 R(IJ)=0.0 EIGEN050 IF(I-J) 1,31,1 EIGEN051 31 R(IJ)=1.0 EIGEN052 1 CONTINUE EIGEN053 C EIGEN054 C COMPUTE INITIAL AND FINAL NORMS (ANORM AND ANORMX) EIGEN055 C EIGEN056 14 ANORM=0.0 EIGEN057 DO 2 I=1,N EIGEN058 DO 2 J=I,N EIGEN059 IF(I-J) 15,2,15 EIGEN060 15 IA=I+(J*J-J)/2 EIGEN061 ANORM=ANORM+A(IA)*A(IA) EIGEN062 2 CONTINUE EIGEN063 IF(ANORM) 39, 39, 50 EIGEN064 50 ANORM=2.0*SQRT(ANORM) EIGEN065 ANRMX =ANORM*1.0E-8 EIGEN066 C EIGEN067 C INITIALIZE INDICATORS AND COMPUTE THRESHOLD, THR EIGEN068 C EIGEN069 IND=0 EIGEN070 THR=ANORM EIGEN071 3 THR=THR/FLOAT(N) EIGEN072 13 L=1 EIGEN073 4 M=L+1 EIGEN074 C EIGEN075 C COMPUTE SIN AND COS EIGEN076 C EIGEN077 5 MQ=(M*M-M)/2 EIGEN078 LQ=(L*L-L)/2 EIGEN079 LM=L+MQ EIGEN080 IF(ABS(A(LM))-THR) 7,32,32 EIGEN081 32 IND=1 EIGEN082 LL=L+LQ EIGEN083 MM=M+MQ EIGEN084 X=0.5*(A(LL)-A(MM)) EIGEN085 Y=-A(LM)/SQRT(A(LM)*A(LM)+X*X) EIGEN086 IF(X) 33,34,34 EIGEN087 33 Y=-Y EIGEN088 34 SINX=Y/SQRT(2.0*(1.0+(SQRT(1.0-Y*Y)))) EIGEN089 SINX2=SINX*SINX EIGEN090 COSX=SQRT(1.0-SINX2) EIGEN091 COSX2=COSX*COSX EIGEN092 SINCS =SINX*COSX EIGEN093 C EIGEN094 C ROTATE L AND M COLUMNS EIGEN095 C EIGEN096 ILQ=N*(L-1) EIGEN097 IMQ=N*(M-1) EIGEN098 DO 6 I=1,N EIGEN099 IQ=(I*I-I)/2 EIGEN100 IF(I-L) 16,11,16 EIGEN101 16 IF(I-M) 17,11,18 EIGEN102 17 IM=I+MQ EIGEN103 GO TO 19 EIGEN104 18 IM=M+IQ EIGEN105 19 IF(I-L) 37,23,23 EIGEN106 37 IL=I+LQ EIGEN107 GO TO 24 EIGEN108 23 IL=L+IQ EIGEN109 24 X=A(IL)*COSX-A(IM)*SINX EIGEN110 A(IM)=A(IL)*SINX+A(IM)*COSX EIGEN111 A(IL)=X EIGEN112 11 IF(MV-1) 26,6,26 EIGEN113 26 ILR=ILQ+I EIGEN114 IMR=IMQ+I EIGEN115 X=R(ILR)*COSX-R(IMR)*SINX EIGEN116 R(IMR)=R(ILR)*SINX+R(IMR)*COSX EIGEN117 R(ILR)=X EIGEN118 6 CONTINUE EIGEN119 X=2.0*A(LM)*SINCS EIGEN120 Y=A(LL)*COSX2+A(MM)*SINX2-X EIGEN121 X=A(LL)*SINX2+A(MM)*COSX2+X EIGEN122 A(LM)=(A(LL)-A(MM))*SINCS+A(LM)*(COSX2-SINX2) EIGEN123 A(LL)=Y EIGEN124 A(MM)=X EIGEN125 C EIGEN126 C TESTS FOR COMPLETION EIGEN127 C EIGEN128 C TEST FOR M = LAST COLUMN EIGEN129 7 IF(M-N) 35,8,35 EIGEN130 35 M=M+1 EIGEN131 GO TO 5 EIGEN132 C TEST FOR L = SECOND FROM LAST COLUMN EIGEN133 8 IF(L-(N-1)) 36,9,36 EIGEN134 36 L=L+1 EIGEN135 GO TO 4 EIGEN136 9 IF(IND-1) 10,38,10 EIGEN137 38 IND=0 EIGEN138 GO TO 13 EIGEN139 C COMPARE THRESHOLD WITH FINAL NORM EIGEN140 10 IF(THR-ANRMX ) 39,39,3 EIGEN141 C EIGEN142 C SORT EIGENVALUES AND EIGENVECTORS EIGEN143 C EIGEN144 39 IQ=-N EIGEN145 DO 30 I=1,N EIGEN146 IQ=IQ+N EIGEN147 LL=I+(I*I-I)/2 EIGEN148 JQ=N*(I-2) EIGEN149 DO 30 J=I,N EIGEN150 JQ=JQ+N EIGEN151 MM=J+(J*J-J)/2 EIGEN152 IF(A(LL)-A(MM)) 40,30,30 EIGEN153 40 X=A(LL) EIGEN154 A(LL)=A(MM) EIGEN155 A(MM)=X EIGEN156 IF(MV-1) 28,30,28 EIGEN157 28 DO 20 K=1,N EIGEN158 ILR=IQ+K EIGEN159 IMR=JQ+K EIGEN160 X=R(ILR) EIGEN161 R(ILR)=R(IMR) EIGEN162 20 R(IMR)=X EIGEN163 30 CONTINUE EIGEN164 RETURN EIGEN165 END EIGEN166 C ARRAY001 C ..................................................................ARRAY002 C ARRAY003 C SUBROUTINE ARRAY ARRAY004 C ARRAY005 C PURPOSE ARRAY006 C CONVERT DATA ARRAY FROM SINGLE TO DOUBLE DIMENSION OR VICE ARRAY007 C VERSA. THIS SUBROUTINE IS USED TO LINK THE USER PROGRAM ARRAY008 C WHICH HAS DOUBLE DIMENSION ARRAYS AND THE SSP SUBROUTINES ARRAY009 C WHICH OPERATE ON ARRAYS OF DATA IN A VECTOR FASHION. ARRAY010 C ARRAY011 C USAGE ARRAY012 C CALL ARRAY (MODE,I,J,N,M,S,D) ARRAY013 C ARRAY014 C DESCRIPTION OF PARAMETERS ARRAY015 C MODE - CODE INDICATING TYPE OF CONVERSION ARRAY016 C 1 - FROM SINGLE TO DOUBLE DIMENSION ARRAY017 C 2 - FROM DOUBLE TO SINGLE DIMENSION ARRAY018 C I - NUMBER OF ROWS IN ACTUAL DATA MATRIX ARRAY019 C J - NUMBER OF COLUMNS IN ACTUAL DATA MATRIX ARRAY020 C N - NUMBER OF ROWS SPECIFIED FOR THE MATRIX D IN ARRAY021 C DIMENSION STATEMENT ARRAY022 C M - NUMBER OF COLUMNS SPECIFIED FOR THE MATRIX D IN ARRAY023 C DIMENSION STATEMENT ARRAY024 C S - IF MODE=1, THIS VECTOR CONTAINS, AS INPUT, A DATA ARRAY025 C MATRIX OF SIZE I BY J IN CONSECUTIVE LOCATIONS ARRAY026 C COLUMN-WISE. IF MODE=2, IT CONTAINS A DATA MATRIX ARRAY027 C OF THE SAME SIZE AS OUTPUT. THE LENGTH OF VECTOR S ARRAY028 C IS IJ, WHERE IJ=I*J. ARRAY029 C D - IF MODE=1, THIS MATRIX (N BY M) CONTAINS, AS OUTPUT, ARRAY030 C A DATA MATRIX OF SIZE I BY J IN FIRST I ROWS AND ARRAY031 C J COLUMNS. IF MODE=2, IT CONTAINS A DATA MATRIX OF ARRAY032 C THE SAME SIZE AS INPUT. ARRAY033 C ARRAY034 C REMARKS ARRAY035 C VECTOR S CAN BE IN THE SAME LOCATION AS MATRIX D. VECTOR S ARRAY036 C IS REFERRED AS A MATRIX IN OTHER SSP ROUTINES, SINCE IT ARRAY037 C CONTAINS A DATA MATRIX. ARRAY038 C THIS SUBROUTINE CONVERTS ONLY GENERAL DATA MATRICES (STORAGEARRAY039 C MODE OF 0). ARRAY040 C ARRAY041 C SUBROUTINES AND FUNCTION SUBROUTINES REQUIRED ARRAY042 C NONE ARRAY043 C ARRAY044 C METHOD ARRAY045 C REFER TO THE DISCUSSION ON VARIABLE DATA SIZE IN THE SECTIONARRAY046 C DESCRIBING OVERALL RULES FOR USAGE IN THIS MANUAL. ARRAY047 C ARRAY048 C ..................................................................ARRAY049 C ARRAY050 SUBROUTINE ARRAY (MODE,I,J,N,M,S,D) DIMENSION S(1),D(1) C ARRAY053 NI=N-I ARRAY054 C ARRAY055 C TEST TYPE OF CONVERSION ARRAY056 C ARRAY057 IF(MODE-1) 100, 100, 120 ARRAY058 C ARRAY059 C CONVERT FROM SINGLE TO DOUBLE DIMENSION ARRAY060 C ARRAY061 100 IJ=I*J+1 ARRAY062 NM=N*J+1 ARRAY063 DO 110 K=1,J ARRAY064 NM=NM-NI ARRAY065 DO 110 L=1,I ARRAY066 IJ=IJ-1 ARRAY067 NM=NM-1 ARRAY068 110 D(NM)=S(IJ) ARRAY069 GO TO 140 ARRAY070 C ARRAY071 C CONVERT FROM DOUBLE TO SINGLE DIMENSION ARRAY072 C ARRAY073 120 IJ=0 ARRAY074 NM=0 ARRAY075 DO 130 K=1,J ARRAY076 DO 125 L=1,I ARRAY077 IJ=IJ+1 ARRAY078 NM=NM+1 ARRAY079 125 S(IJ)=D(NM) ARRAY080 130 NM=NM+NI ARRAY081 C ARRAY082 140 RETURN ARRAY083 END ARRAY084 SUBROUTINE SORT (A, N, B, C, D, K, SWITCH ) SORT 110 CSORT SORT 100 C SORT 120 C THIS ROUTINE SORTS INPUT ARRAY 'A' AND REARRANGES, OPTIONALLY, SORT 130 C ARRAYS 'B', 'C', AND 'D', IN ORDER CORRESPONDING TO 'A'. SORT 140 C N = NUMBER OF ITEMS IN 'A' (AND 'B', 'C', 'D', IF USED) SORT 150 C K = 0--SORT 'A' ONLY, 1--REARRANGE 'B', 2--REARRANGE 'B' AND 'C', SORT 160 C 3--REARRANGE 'B', 'C', AND 'D'. SORT 170 C IF 'SWITCH' IS POSITIVE, SORT WILL BE IN ASCENDING ORDER, SORT 180 C IF ZERO OR NEGATIVE, IN DESCENDING ORDER. SORT 190 C ALGORITHM FROM CACM, JULY 1959, PAGE 30 BY D. L. SHELL SORT 200 C SORT 210 DIMENSION A(N), B(N), C(N), D(N) INTEGER SWITCH SORT 230 105 KP1=K+1 SORT 240 IF(N.LE.1) GO TO 999 SORT 250 M = 1 SORT 260 106 M = M + M SORT 270 IF( M .LE. N ) GO TO 106 SORT 280 M = M - 1 SORT 290 994 M = M/2 SORT 300 IF(M.EQ.0) GO TO 999 SORT 310 KK = N-M SORT 320 J = 1 SORT 330 992 I = J SORT 340 996 IM = I + M SORT 350 IF(SWITCH) 810,810,800 SORT 360 800 IF(A(I).GT.A(IM)) GO TO 110 SORT 370 GO TO 995 SORT 380 810 IF(A(I).LT.A(IM)) GO TO 110 SORT 390 995 J = J+1 SORT 400 IF(J.GT.KK) GO TO 994 SORT 410 GO TO 992 SORT 420 110 TEMP=A(I) SORT 430 A(I) = A(IM) SORT 440 A(IM) = TEMP SORT 450 GO TO ( 140, 130, 120, 115), KP1 SORT 460 115 TEMP = D(I) SORT 470 D(I) = D(IM) SORT 480 D(IM) = TEMP SORT 490 120 TEMP=C(I) SORT 500 C(I) = C(IM) SORT 510 C(IM) = TEMP SORT 520 130 TEMP=B(I) SORT 530 B(I) = B(IM) SORT 540 B(IM) = TEMP SORT 550 140 I = I-M SORT 560 IF(I.LT.1) GO TO 995 SORT 570 GO TO 996 SORT 580 999 RETURN SORT 590 END SORT 600