C maxscal4.1.f C The author of this software is Yoshio Takane. 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 SECOND MDS Package of AT&T Bell Laboratories C For explanation of the method this software implements, see C "Nonmetric maximum likelihood multidimensional scaling from directional C rankings of similarities" by Y Takane and J D Carroll (1981) in C Psychometrika, vol. 46, pages 389-405. C Anyone who wants a copy of the program write-up may send a request to C takane@takane2.psych.mcgill.ca C----+----@----+----@----+----@----+----@----+----@----+----@----+----@ $ LOWLOAD 10 $ OPTION FORTRAN,RELMEM 20 $ FORTY FORM,NLSTIN 30 C MAXS4--MAXSCAL-4.1,JULY 1979.MAXIMUM LIKELIHOOD MDS,EUCLIDIAN MODEL 40 SUBROUTINE MAXS4 50 C MAXSCAL-4.1 (VERSION 1.01) 60 C MAXIMUM LIKELIHOOD MULTIDIMENSIONAL SCALING 70 C WITH THE SIMPLE EUCLIDIAN MODEL 80 C FROM VARIOUS TYPES OF DIRECTIONAL RANKINGS OF SIMILARITIES 90 C 100 C WRITTEN BY YOSHIO TAKANE 110 C JULY 1979 120 C 130 REAL*8 RX8(6500) 140 DIMENSION RX(13000), IQ(13000) 150 EQUIVALENCE (RX(1),RX8(1)) 160 COMMON /OPT/ IN, LOUT, NPH, TIT(20), FM(20), ME, INT, NCON, MXIT, 170 * CVCR, ITIE, INS, IDATA, INTPR, IDIST, IVEST, IPLTC, IPHC, IPHV, 180 * NBDMX1, NDMX1, NB2, NNB2, NBSQ2, NTIE, NBS, NMIS, NBDMX, ICODE, 190 * NCODE, NB1, FNB1, FNB, FNBI, FNB2, ILLT 200 C 210 NDI = 13000 220 NDR = 13000 230 IN = 5 240 LOUT = 6 250 NPH = 43 260 C 270 10 CONTINUE 280 C 290 C JOB PARAMETERS I/O 300 C 310 C JOB TITLE 320 READ(IN,9999,END=270) TIT 330 9999 FORMAT (20A4) 340 C 350 C PROBLEM SPECIFICATIONS 360 C 370 C NB: NUMBER OF STIMULI 380 C NS: NUMBER OF DATA SETS (DEFAULT=1) 390 C (NUMBER OF SUBJECTS IN CASE OF MATRIX INPUT) 400 C ME: ERROR MODEL 410 C 0= ADDITIVE (DEFAULT) 420 C 1= MULTIPLICATIVE 430 C NDMX: MAXIMUM NUMBER OF DIMENSIONS (DEFAULT=NDMN) 440 C NDMN: MINIMUM NUMBER OF DIMENSIONS (DEFAULT=1) 450 C INT: INITIAL CONFIGURTION 460 C 0= COMPUTE (DEFAULT) 470 C 1= INPUT FROM CARDS 480 C NCON: EXTERNAL CONSTRAINTS 490 C 0= NO (DEFAULT) 500 C 1= YES 510 C MXIT: MAXIMUM NUMBER OF ITERATIONS (DEFAULT=200) 520 C CVCR: CONVERGENCE CRITERION 530 C (DEFAULT= MAXIMUM RELATIVE CHANGE IN LIKELIHOOD 540 C LESS THAN 0.000002) 550 C 560 READ (IN,9998) NB, NS, ME, NDMX, NDMN, INT, NCON, MXIT, CVCR 570 9998 FORMAT (8I4, 8X, F10.0) 580 C 590 WRITE (LOUT,9997) 600 9997 FORMAT (54H1MAXSCAL-4 MAXIMUM LIKELIHOOD MULTIDIMENSIONAL SCALIN, 610 * 61HG FROM VARIOUS TYPES OF DIRECTIONAL RANKING PROCEDURES OF SIM, 620 * 15HILARITIES (SEM)) 630 C 640 WRITE (LOUT,9996) TIT 650 9996 FORMAT (/////13H JOB TITLE: , 20A4///) 660 C 670 IF (NS.EQ.0) NS = 1 680 IF (ME.NE.1) ME = 0 690 IF (NDMN.LT.1) NDMN = 1 700 IF (NDMX.LT.NDMN) NDMX = NDMN 710 IF (INT.NE.1) INT = 0 720 IF (NCON.NE.1) NCON = 0 730 IF (MXIT.LT.1) MXIT = 200 740 IF (CVCR.LE.0.0) CVCR = 2.0E-06 750 IF (NB.LT.NDMX+1) GO TO 210 760 IF (NS.LT.0) GO TO 220 770 C 780 NDMX1 = NDMX + 1 790 NBDMX1 = NB*NDMX1 800 NB2 = NB*NB 810 NBS = NB*NS 820 NNB2 = NB2*NS 830 NBSQ2 = (NB2-NB)/2 840 NBDMX = NB*NDMX 850 FNB = NB 860 FNBI = 1.0/FNB 870 FNB2 = NB2 880 NB1 = NB - 1 890 FNB1 = NB1 900 C 910 WRITE (LOUT,9995) NB, NS, NDMX, NDMN 920 9995 FORMAT (24H PROBLEM SPECIFICATIONS-//12H THERE ARE, I4, 2X, 930 * 11HSTIMULI AND, I4, 2X, 11HDATA SET(S)/22H THE SOLUTIONS ARE O, 940 * 32HBTAINED FOR ALL DIMENSIONALITIES, 8H BETWEEN, I3, 2X, 3HAND, 950 * I3) 960 IF (ME.EQ.0) WRITE (LOUT,9994) 970 9994 FORMAT (30H THE ERROR MODEL IS ADDITIVE) 980 IF (ME.EQ.1) WRITE (LOUT,9993) 990 9993 FORMAT (36H THE ERROR MODEL IS MULTIPLICATIVE) 1000 IF (INT.EQ.0) WRITE (LOUT,9992) 1010 9992 FORMAT (37H INITIAL ESTIMATES WILL BE COMPUTED) 1020 IF (INT.EQ.1) WRITE (LOUT,9991) 1030 9991 FORMAT (34H INITIAL ESTIMATES WILL BE INPUT) 1040 IF (NCON.EQ.0) WRITE (LOUT,9990) 1050 9990 FORMAT (36H THERE ARE NO EXTERNAL CONSTRAINTS) 1060 IF (NCON.EQ.1) WRITE (LOUT,9989) 1070 9989 FORMAT (38H THERE ARE SOME EXTERNAL CONSTRAINTS) 1080 WRITE (LOUT,9988) MXIT, CVCR 1090 9988 FORMAT (45H THE MAXIMUM NUMBER OF ITERATIONS IS SET TO, 1100 * I4/57H THE ITERATION WILL BE TERMINATED WHEN THE MAXIMUM RELA, 1110 * 40HTIVE CHANGE IN LIKELIHOOD GETS LESS THAN, F11.7) 1120 C 1130 C OUTPUT OPTIONS 1140 C 1150 C IDATA 1= PRINT INPUT DATA 1160 C INTPR 1= PRINT INITIAL ESTIMATES 1170 C IDIST 1= PRINT DERIVED DISTANCES 1180 C IVEST 1= PRINT VARIANCE-COVARIANCE ESTIMATES 1190 C ILLT 1= PRINT LOG-LIKELIHOOD OF EACH FIRST CHOICE 1200 C IPLTC 1= PLOT FINAL ESTIMATES 1210 C IPHC 1= PUNCH FINAL ESTIMATES 1220 C IPHV 1= PUNCH VARIANCE COVARIANCE ESTIMATES 1230 C 1240 READ (IN,9987) IDATA, INTPR, IDIST, IVEST, ILLT, IPLTC, IPHC, IPHV 1250 9987 FORMAT (6I4, 4X, 2I4) 1260 C 1270 IF (IDATA.NE.1) IDATA = 0 1280 IF (INTPR.NE.1) INTPR = 0 1290 IF (IDIST.NE.1) IDIST = 0 1300 IF (IVEST.NE.1) IVEST = 0 1310 IF (ILLT.NE.1) ILLT = 0 1320 IF (IPLTC.NE.1) IPLTC = 0 1330 IF (IPHC.NE.1) IPHC = 0 1340 IF (IPHV.NE.1) IPHV = 0 1350 IOUT = IDATA + INTPR + IDIST + IVEST + IPLTC + IPHC + IPHV 1360 IF (IOUT.NE.0) WRITE (LOUT,9986) 1370 9986 FORMAT (///16H OUTPUT OPTIONS-) 1380 WRITE (LOUT,9985) 1390 9985 FORMAT (1H ) 1400 IF (IDATA.EQ.1) WRITE (LOUT,9984) 1410 9984 FORMAT (19H PRINT INPUT DATA) 1420 IF (INTPR.EQ.1) WRITE (LOUT,9983) 1430 9983 FORMAT (26H PRINT INITIAL ESTIMATES) 1440 IF (IDIST.EQ.1) WRITE (LOUT,9982) 1450 9982 FORMAT (26H PRINT DERIVED DISTANCES) 1460 IF (IVEST.EQ.1) WRITE (LOUT,9981) 1470 9981 FORMAT (38H PRINT VARIANCE-COVARIANCE ESTIMATES) 1480 IF (ILLT.EQ.1) WRITE (LOUT,9980) 1490 9980 FORMAT (44H PRINT LOG-LIKELIHOOD OF EACH FIRST CHOICE) 1500 IF (IPLTC.EQ.1) WRITE (LOUT,9979) 1510 9979 FORMAT (24H PLOT FINAL ESTIMATES) 1520 IF (IPHC.EQ.1) WRITE (LOUT,9978) 1530 9978 FORMAT (24H PUNCH FINAL ESTIMATES) 1540 IF (IPHV.EQ.1) WRITE (LOUT,9977) 1550 9977 FORMAT (38H PUNCH VARIANCE-COVARINACE ESTIMATES) 1560 C 1570 C INPUT FORMAT 1580 C 1590 READ (IN,9999) FM 1600 WRITE (LOUT,9976) FM 1610 9976 FORMAT (///14H INPUT FORMAT-//6X, 20A4) 1620 C 1630 C DATA SPECIFICATIONS 1640 C 1650 C ICODE: DATA COLLECTION METHOD 1660 C 1=TETRADS 1670 C 2=TRIADS 1680 C 3=TRIADIC COMBINATIONS 1690 C 4=CONDITIONAL RANK ORDERS (DEFAULT) 1700 C 5=GENERAL DIRECTIONAL RANK ORDERS 1710 C NCODE: INPUT FORMAT CODE 1720 C WHEN ICODE=1, 1730 C 1=INPUT I,J,K,L,F(>); SET F(<)=NR-F(>),F(=)=0. 1740 C 2=INPUT I,J,K,L,F(>),F(<); SET F(=)=0. 1750 C 3=INPUT I,J,K,L,F(>),F(<); SET F(=)=NR-F(>)-F(<). 1760 C 4=INPUT I,J,K,L,F(>),F(<),F(=); (DEFAULT). 1770 C WHEN ICODE=2, 1780 C THE SAME THING APPLIES AS WHEN ICODE=1, EXECPT 1790 C THAT ONLY THREE INDICES (I,J & K) WILL BE INPUT 1800 C INSTEAD OF FOUR, WHERE I INDICATES THE COMMON ST. 1810 C WHEN ICODE=3, 1820 C NCODE IS IRRELEVANT; INPUT I,J,K,F(IJ)+F(<)+F(=). (RELEVANT ONLY WHEN ICODE=1 OR 2 AND THE 2450 C NUMBER OF REPLICATED OBSERVATIONS (NR) IS EQUAL FOR ALL 2460 C TETRADS OR TRIADS.) 2470 C 2480 READ (IN,9975) ICODE, NCODE, ITIE, INS, NMIS, NTIE, NR 2490 9975 FORMAT (20I4) 2500 C 2510 IF (ICODE.EQ.0) ICODE = 4 2520 IF (ICODE.EQ.1 .AND. NCODE.EQ.0) NCODE = 4 2530 IF (ICODE.EQ.2 .AND. NCODE.EQ.0) NCODE = 4 2540 IF (ICODE.EQ.4 .AND. NCODE.EQ.0) NCODE = 1 2550 IF (ICODE.EQ.5 .AND. NCODE.EQ.0) NCODE = 1 2560 IF (ITIE.NE.1) ITIE = 0 2570 IF (INS.NE.1) INS = 0 2580 IF (NS.EQ.1) INS = 0 2590 IF (ICODE.LT.0 .OR. ICODE.GT.5) GO TO 230 2600 IF (NCODE.LT.0) GO TO 240 2610 IF (ICODE.EQ.3) GO TO 30 2620 IF (ICODE.GT.3) GO TO 20 2630 IF (NCODE.GT.4) GO TO 240 2640 GO TO 30 2650 20 IF (ICODE+NCODE.GT.8) GO TO 240 2660 C 2670 30 WRITE (LOUT,9974) 2680 9974 FORMAT (///21H DATA SPECIFICATIONS-) 2690 WRITE (LOUT,9985) 2700 IF (ICODE.EQ.3) GO TO 40 2710 IF (ICODE.GT.3) GO TO 50 2720 IF (ICODE.EQ.1) WRITE (LOUT,9973) NCODE 2730 9973 FORMAT (53H THE DATA ARE COLLECTED BY THE METHOD OF TETRADS , 2740 * 19H(INPUT FORMAT CODE=, I2, 1H)) 2750 IF (ICODE.EQ.2) WRITE (LOUT,9972) NCODE 2760 9972 FORMAT (52H THE DATA ARE COLLECTED BY THE METHOD OF TRIADS , 2770 * 19H(INPUT FORMAT CODE=, I2, 1H)) 2780 GO TO 70 2790 40 WRITE (LOUT,9971) 2800 9971 FORMAT (51H THE DATA ARE COLLECTED BY THE METHOD OF TRIADIC , 2810 * 12HCOMBINATIONS) 2820 GO TO 70 2830 50 IF (ICODE.EQ.5) GO TO 60 2840 IF (NCODE.EQ.1) WRITE (LOUT,9970) 2850 9970 FORMAT (54H THE DATA ARE COLLECTED BY THE METHOD OF CONDITIONAL, 2860 * 1H , 25HRANK ORDERS: MATRIX INPUT) 2870 IF (NCODE.EQ.2) WRITE (LOUT,9969) 2880 9969 FORMAT (54H THE DATA ARE COLLECTED BY THE METHOD OF CONDITIONAL, 2890 * 1H , 28HRANK ORDERS: FREQUENCY INPUT) 2900 GO TO 70 2910 60 IF (NCODE.EQ.1) WRITE (LOUT,9968) 2920 9968 FORMAT (52H THE DATA ARE COLLECTED BY THE METHOD OF GENERAL , 2930 * 57HDIRECTIONAL RANK ORDERS: MATRIX INPUT, MATRIX CONDITIONAL) 2940 IF (NCODE.EQ.2) WRITE (LOUT,9967) 2950 9967 FORMAT (51H THE DATA ARE COLLECTED BY THE METHOD OF GENERAL , 2960 * 61HDIRECTIONAL RANK ORDERS: MATRIX INPUT, ARBITRARY PARTITIONING, 2970 * 1HS) 2980 IF (NCODE.EQ.3) WRITE (LOUT,9966) 2990 9966 FORMAT (51H THE DATA ARE COLLECTED BY THE METHOD OF GENERAL , 3000 * 40HDIRECTIONAL RANK ORDERS: FREQUENCY INPUT) 3010 70 IF (ITIE.EQ.0) WRITE (LOUT,9965) 3020 9965 FORMAT (52H TIED OBSERVATIONS, IF THEY EXIST, ARE DISREGARDED) 3030 IF (ITIE.EQ.1) WRITE (LOUT,9964) 3040 9964 FORMAT (50H TIED OBSERVATIONS, IF THEY EXIST, ARE RESPECTED) 3050 IF (NS.EQ.1) GO TO 80 3060 IF (INS.EQ.0) WRITE (LOUT,9963) 3070 9963 FORMAT (54H DISPERSIONS ARE ASSUMED CONSTANT ACROSS REPLICATION, 3080 * 1HS) 3090 IF (INS.EQ.1) WRITE (LOUT,9962) 3100 9962 FORMAT (52H DISPERSIONS ARE ASSUMED TO VARY OVER REPLICATIONS) 3110 80 CONTINUE 3120 IF (NMIS.EQ.1) WRITE (LOUT,9961) 3130 9961 FORMAT (54H IF MISSING DATA EXIST, RANKINGS ARE ASSUMED MADE WI, 3140 * 29HTH RESTRICTED SETS OF STIMULI) 3150 IF (ICODE.LT.4) NTIE = 0 3160 IF (ICODE.EQ.4 .AND. NCODE.EQ.1) NTIE = 0 3170 IF (ICODE.EQ.5 .AND. NCODE.LT.3) NTIE = 0 3180 IF (NTIE.NE.1) NTIE = 0 3190 IF (NTIE.EQ.1) WRITE (LOUT,9960) 3200 9960 FORMAT (54H INFORMATION REGARDING TIED RANKS (IN FREQUENCY INPU, 3210 * 2HT), 14H WILL BE INPUT) 3220 IF (ICODE.GE.3) GO TO 90 3230 IF (NCODE/2*2.EQ.NCODE) GO TO 90 3240 IF (NR.LE.0) GO TO 250 3250 WRITE (LOUT,9959) NR 3260 9959 FORMAT (54H THE NUMBER OF REPLICATIONS (ASSUMED EQUAL) PER TETR, 3270 * 2HAD, 3H IS, I4) 3280 C 3290 C IF NOT MATRIX INPUT, READ 3300 C 3310 90 IF (ICODE.EQ.4 .AND. NCODE.EQ.1) GO TO 120 3320 IF (ICODE.EQ.5 .AND. NCODE.LT.3) GO TO 120 3330 READ (IN,9958) (IQ(I),I=1,NS) 3340 9958 FORMAT (20I4) 3350 DO 100 I=1,NS 3360 IF (IQ(I).LE.0) GO TO 260 3370 100 CONTINUE 3380 WRITE (LOUT,9957) 3390 9957 FORMAT (///51H THE NUMBER OF OBSERVED CONDITIONAL RANK ORDERS (OR, 3400 * 22H TETRADS) PER DATA SET//9X, 3HSET, 6X, 3HNUM) 3410 WRITE (LOUT,9985) 3420 DO 110 I=1,NS 3430 WRITE (LOUT,9956) I, IQ(I) 3440 110 CONTINUE 3450 9956 FORMAT (6X, I5, I10) 3460 120 CONTINUE 3470 C 3480 C DATA INPUT 3490 C 3500 NTS = 0 3510 DO 130 K=1,NS 3520 NTS = NTS + IQ(K) 3530 130 CONTINUE 3540 NQ31 = NS + 1 3550 IF (ICODE.GT.3) GO TO 160 3560 IF (ICODE.EQ.3) GO TO 140 3570 NTS3 = NTS*3 3580 NQ32 = NQ31 + NTS3 3590 NQ33 = NQ32 + NTS3 3600 NQ34 = NQ33 + NTS3*4 3610 NQ35 = NQ34 + NTS3*2 3620 NQI = NQ35 + NTS*4 3630 IF (NQI.GT.NDI) GO TO 200 3640 NQ1 = NTS3 + 1 3650 NQ2 = NQ1 + NB2 3660 IF (INT.EQ.1) NQ2 = NQ1 + 1 3670 NQR = NQ2 + NTS3 3680 IF (NQR.GT.NDR) GO TO 190 3690 NQ3 = NQR 3700 C 3710 CALL TETRA(IQ(NQ35), RX(NQ2), IQ(1), IQ(NQ31), IQ(NQ32), 3720 * IQ(NQ33), IQ(NQ34), RX(1), NB, NS, NR, NTS) 3730 GO TO 150 3740 C 3750 140 NTS6 = NTS*6 3760 NQ32 = NQ31 + NTS6 3770 NQ33 = NQ32 + NTS6 3780 NQ34 = NQ33 + NTS6*6 3790 NQ35 = NQ34 + NTS6*3 3800 NQI = NQ35 + NTS6*2 3810 IF (INT.EQ.1) NQI = NQ35 + 1 3820 IF (NQI.GT.NDI) GO TO 200 3830 NQ1 = NTS6 + 1 3840 NQ2 = NQ1 + NB2 3850 IF (INT.EQ.1) NQ2 = NQ1 + 1 3860 NQ3 = NQ2 + NTS*9 3870 IF (INT.EQ.1) NQ3 = NQ2 + 1 3880 NQR = NQ3 + 6 3890 IF (NQR.GT.NDR) GO TO 190 3900 C 3910 CALL TRIADC(IQ(NQ35), RX(NQ2), IQ(1), IQ(NQ31), IQ(NQ32), 3920 * IQ(NQ33), IQ(NQ34), RX(1), RX(NQ3), NB, NS, NTS) 3930 C 3940 150 IF (INT.EQ.1) GO TO 180 3950 NQ4 = NQ3 + NBSQ2 3960 NQ5 = NQ4 + NBSQ2 3970 M = NBSQ2*(NBSQ2+1)/2 3980 NQR = NQ5 + M 3990 IF (NQR.GT.NDR) GO TO 190 4000 C 4010 CALL THURS(RX(NQ4), RX(NQ5), RX(NQ2), IQ(NQ35), RX(NQ1), RX(NQ3), 4020 * NB, NTS, M) 4030 GO TO 180 4040 C 4050 160 IF (ICODE.EQ.5) GO TO 170 4060 NBSX = NBS 4070 IF (NCODE.EQ.2) NBSX = NTS 4080 NQ32 = NQ31 + NBSX 4090 NQ33 = NQ32 + NBSX 4100 NBS1 = NBSX*NB1 4110 NQ34 = NQ33 + NBS1*2 4120 NQ35 = NQ34 + NBS1 4130 NQ36 = NQ35 + NB 4140 NQ37 = NQ36 + NB 4150 NQI = NQ37 + NB2 4160 IF (NQI.GT.NDI) GO TO 200 4170 NQ1 = NBSX + 1 4180 NQ2 = NQ1 + NB2 4190 IF (INT.EQ.1) NQ2 = NQ1 + 1 4200 NQ3 = NQ2 + NB 4210 NQR = NQ3 + NB 4220 IF (NQR.GT.NDR) GO TO 190 4230 C 4240 CALL CRANK(IQ(NQ37), IQ(1), IQ(NQ31), IQ(NQ32), RX(NQ2), RX(NQ3), 4250 * IQ(NQ35), IQ(NQ36), IQ(NQ34), IQ(NQ33), RX(NQ1), RX(1), NS, NB) 4260 GO TO 180 4270 C 4280 170 NBSQ4 = NBSQ2/2*NS 4290 IF (NCODE.EQ.3) NBSQ4 = NTS 4300 NQ32 = NQ31 + NBSQ4 4310 NQ33 = NQ32 + NBSQ4 4320 NBSQS = NBSQ2*NS 4330 IF (NCODE.EQ.3) NBSQS = NTS*NBSQ2 4340 NQ34 = NQ33 + NBSQS*2 4350 NQ35 = NQ34 + NBSQS 4360 NQ36 = NQ35 + NBSQ2 4370 NQ37 = NQ36 + NBSQ2 4380 NQ38 = NQ37 + NBSQ2*2 4390 NQI = NQ38 + NBSQ2 4400 IF (NQI.GT.NDI) GO TO 200 4410 NQ1 = NBSQ4 + 1 4420 NQ2 = NQ1 + NB2 4430 IF (INT.EQ.1) NQ2 = NQ1 + 1 4440 NQ3 = NQ2 + NBSQ2 4450 NQR = NQ3 + NBSQ2 4460 IF (NQR.GT.NDR) GO TO 190 4470 C 4480 CALL GRANK(IQ(NQ38), IQ(NQ35), IQ(NQ36), IQ(1), IQ(NQ31), 4490 * IQ(NQ32), RX(NQ1), IQ(NQ34), IQ(NQ33), RX(NQ2), RX(NQ3), 4500 * IQ(NQ37), RX(1), NB, NS) 4510 C 4520 180 NP0 = NB*NDMX + NS 4530 NQ3 = NQ2 + NB 4540 NQ5 = NQ3 + NB 4550 NP1 = MAX0(NBDMX1,NP0) 4560 NQ6 = NQ5 + NP1 4570 NQ7 = NQ6 + NP1 4580 NQ8 = NQ7 + NP1 4590 NP2 = NDMX1*NDMX1 4600 NQ9 = NQ8 + NP2 4610 NQ10 = NQ9 + NP2 4620 NQ11 = NQ10 + NDMX1 4630 NQ12 = NQ11 + NP0 4640 NP02 = NP0*(NP0+1)/2 4650 NQ13 = NQ12 + NP02 4660 NQ14 = NQ13 + NP0 4670 NQ15 = NQ14 + NP0 4680 NQ16 = NQ15 + NP02 4690 NQ17 = NQ16 + NBSQ2 4700 NQ18 = NQ17 + NBSQ2 4710 NQ19 = NQ18 + NBSQ2 4720 NQ20 = NDR/2 - NBDMX + 1 4730 NQ21 = NQ20 - NP0 4740 NQR = NQ19 + NP02 + (NBDMX+NP0)*2 4750 IF (NQR.GT.NDR) GO TO 190 4760 NQ36 = NQ35 + NP0 4770 NQI = NQ36 + NP0 4780 IF (NQI.GT.NDI) GO TO 200 4790 C 4800 CALL EXEC4(RX8(NQ21), RX(NQ1), RX(NQ2), RX(NQ3), RX(NQ5), 4810 * RX(NQ6), RX(NQ7), RX(NQ8), RX(NQ9), RX(NQ10), RX(NQ11), 4820 * RX(NQ12), RX(NQ13), RX(NQ14), RX(NQ15), RX(NQ16), RX(NQ17), 4830 * RX(NQ18), RX(NQ19), IQ(NQ34), IQ(NQ33), IQ(NQ35), IQ(NQ36), NB, 4840 * NS, NDMX, NDMN, RX8(NQ20), IQ(1), IQ(NQ31), IQ(NQ32), RX(1)) 4850 C 4860 GO TO 10 4870 C 4880 190 WRITE (LOUT,9955) NQR 4890 9955 FORMAT (///51H STORAGE REQUIREMENT FOR REAL ARRAYS EXCEEDS THE LI, 4900 * 3HMIT/37H INCREASE THE SIZE OF NDR TO AT LEAST, I8) 4910 270 STOP 4920 200 WRITE (LOUT,9954) NQI 4930 9954 FORMAT (///51H STORAGE REQUIREMENT FOR INTEGER ARRAYS EXCEEDS THE, 4940 * 6H LIMIT/37H INCREASE THE SIZE OF NDI TO AT LEAST, I8) 4950 STOP 4960 210 WRITE (LOUT,9953) NB, NDMX 4970 9953 FORMAT (///22H THE NUMBER OF STIMULI, I4, 2X, 16HIS TOO SMALL FOR, 4980 * 5H THE , 17HDIMENSIONALITY OF, I4) 4990 STOP 5000 220 WRITE (LOUT,9952) NS 5010 9952 FORMAT (///36H THE NUMBER OF DATA SETS IS NEGATIVE, I4) 5020 STOP 5030 230 WRITE (LOUT,9951) ICODE 5040 9951 FORMAT (///7H ICODE=, I3, 2X, 18HIS NOT PERMISSIBLE) 5050 STOP 5060 240 WRITE (LOUT,9950) NCODE 5070 9950 FORMAT (///51H THE DATA INPUT FORMAT CODE IS OUT OF RANGE (ICODE=, 5080 * I3, 1X, 6HNCODE=, I3, 1H)) 5090 STOP 5100 250 WRITE (LOUT,9949) NR 5110 9949 FORMAT (///51H THE NUMBER OF REPLICATED OBSERVATIONS PER TETRAD I, 5120 * 1HS, 14H INAPPROPRIATE, I4) 5130 STOP 5140 260 WRITE (LOUT,9948) (IQ(I),I=1,NS) 5150 9948 FORMAT (///49H SOME OF THE NUMBERS OF OBSERVED CONDITIONAL RANK, 5160 * 32H ORDERS (OR TETRADS) IS NEGATIVE, 10I4/(81X, 10I4/)) 5170 STOP 5180 END 5190 $ FORTY FORM 5200 C 5210 SUBROUTINE EXEC4(AG,P,DT,DU,U,X,Y,B,BK,ALAM,G,H,WK1,WK2,WK3,D,DD, 5220 1DS,WK4,IY,IZ,IPST,IAP,NB,NS,NDMX,NDMN,AH,NCR,NSC,NSR,F) 5230 C 5240 C EXECUTING ROUTINE 5250 C 5260 REAL*8 TM1,FLL4,FL 5270 REAL*8 AG(1),AH(1) 5280 DIMENSION P(NB,1),DT(1),DU(1),U(1),X(1),Y(1),B(1),BK(1), 5290 1ALAM(1),G(1),H(1),WK1(1),WK2(1),WK3(1),D(1),DD(1),DS(1), 5300 2WK4(1),F(1) 5310 DIMENSION IY(1),IZ(2,1),IPST(1),IAP(1),NCR(1),NSC(1),NSR(1) 5320 COMMON/OPT/IN,LOUT,NPH,TIT(20),FM(20),ME,INT,NCON,MXIT,CVCR,ITIE, 5330 1INS,IDATA,INTPR,IDIST,IVEST,IPLTC,IPHC,IPHV,NBDMX1,NDMX1,NB2,NNB2, 5340 2NBSQ2,NTIE,NBS,NMIS,NBDMX,ICODE,NCODE,NB1,FNB1,FNB,FNBI,FNB2,ILLT 5350 COMMON/NUM/NSQ2,STEPQ,IDDX 5360 C 5370 N=0 5380 NN=0 5390 DO 345 K=1,NS 5400 NCRK=NCR(K) 5410 DO 345 I=1,NCRK 5420 N=N+1 5430 NSCN=NSC(N) 5440 DO 345 J=1,NSCN 5450 NN=NN+1 5460 IF(IZ(1,NN).GT.IZ(2,NN)) GO TO 345 5470 II=IZ(1,NN) 5480 IZ(1,NN)=IZ(2,NN) 5490 IZ(2,NN)=II 5500 345 CONTINUE 5510 FNN=NN-N 5520 C 5530 READ(IN,8340) STEPQ 5540 8340 FORMAT(F10.5) 5550 IF(STEPQ.EQ.0.0) STEPQ=1.0 5560 WRITE(LOUT,8341) STEPQ 5570 8341 FORMAT(//28H INITIAL STEP SIZE IS SET TO,F10.4/) 5580 C 5590 C INITIAL ESTIMATES (COMPUTE) 5600 C 5610 IF(INT.EQ.1) GO TO 29 5620 C ADDITIVE CONSTANT 5630 CO=-1.0E10 5640 CX=1.0E10 5650 DO 26 I=2,NB 5660 I1=I-1 5670 DO 26 J=1,I1 5680 DO 26 K=1,NB 5690 IF(K.EQ.I) GO TO 26 5700 IF(K.EQ.J) GO TO 26 5710 C=P(I,J)-P(I,K)-P(J,K) 5720 IF(C.GT.CO) CO=C 5730 26 CONTINUE 5740 DO 27 I=2,NB 5750 I1=I-1 5760 DO 27 J=1,I1 5770 C=P(I,J) 5780 27 IF(C.LT.CX) CX=C 5790 CX=-CX 5800 IF(CX.GT.CO) CO=CX 5810 DO 28 I=2,NB 5820 I1=I-1 5830 DO 28 J=1,I1 5840 CX=P(I,J)+CO 5850 P(I,J)=CX*CX 5860 28 P(J,I)=P(I,J) 5870 C DOUBLE CENTERING 5880 T=0.0 5890 DO 21 I=1,NB 5900 DT(I)=0.0 5910 DO 22 J=1,NB 5920 22 DT(I)=DT(I)+P(J,I) 5930 T=T+DT(I) 5940 21 DT(I)=DT(I)*FNBI 5950 T=T/FNB2 5960 S=0.0 5970 DO 23 I=1,NB 5980 DO 23 J=1,I 5990 P(I,J)=DT(I)+DT(J)-P(I,J)-T 6000 IF(I.EQ.J) GO TO 24 6010 S=S+2.0*P(I,J)*P(I,J) 6020 GO TO 23 6030 24 S=S+P(I,I)*P(I,I) 6040 23 CONTINUE 6050 S=SQRT(FNB2/S) 6060 DO 25 I=1,NB 6070 DO 25 J=1,I 6080 P(I,J)=P(I,J)*S 6090 25 P(J,I)=P(I,J) 6100 C EIGEN DECOMPOSITION 6110 CALL CJEIG(P,U,X,NB,NDMX1,B,Y,ALAM,1,NB,NDMX1,BK) 6120 GO TO 30 6130 C 6140 C INITIAL ESTIMATES (INPUT) 6150 C 6160 29 DO 31 I=1,NB 6170 31 READ(IN,100) (U(J),J=I,NBDMX,NB) 6180 100 FORMAT(8F10.5) 6190 NN=0 6200 DO 32 J=1,NDMX 6210 C=0.0 6220 S=0.0 6230 N1=NN+1 6240 NN=NN+NB 6250 DO 33 N=N1,NN 6260 C=C+U(N) 6270 33 S=S+U(N)*U(N) 6280 C=C*FNBI 6290 ALAM(J)=S-FNB*C*C 6300 S=1.0/SQRT(ALAM(J)) 6310 DO 34 N=N1,NN 6320 34 U(N)=(U(N)-C)*S 6330 32 CONTINUE 6340 C 6350 C EXTERNAL CONSTRAINTS (DIMENSIONWISE EQUALITY CONSTRAINTS) 6360 C 6370 30 CONTINUE 6380 IF(NCON.EQ.0) GO TO 35 6390 DO 36 I=1,NB 6400 36 READ(IN,101) (IPST(J),J=I,NBDMX,NB) 6410 101 FORMAT(20I4) 6420 JJ=0 6430 IJA=0 6440 IJ=0 6450 DO 137 J=1,NDMX 6460 DO 37 I=1,NB 6470 IJ=IJ+1 6480 IF(IPST(IJ).GT.0) GO TO 38 6490 39 IJA=IJA+1 6500 IAP(IJ)=IJA 6510 GO TO 37 6520 38 IF(I.EQ.1) GO TO 39 6530 I1=JJ+1 6540 II=JJ+I-1 6550 DO 40 J1=I1,II 6560 IF(IPST(J1).EQ.IPST(IJ)) GO TO 41 6570 40 CONTINUE 6580 GO TO 39 6590 41 IAP(IJ)=IAP(J1) 6600 37 CONTINUE 6610 JJ=JJ+NB 6620 137 CONTINUE 6630 C 6640 WRITE(LOUT,355) 6650 355 FORMAT(1H1//21H STATUS OF PARAMETERS// 6660 110X,3HNO.,5X,8HSTIMULUS,3X,9HDIMENSION,5X,6HSTATUS,10X, 6670 217HACTIVE PARAMETERS) 6680 WRITE(LOUT,348) 6690 IJ=0 6700 DO 42 J=1,NDMX 6710 DO 42 I=1,NB 6720 IJ=IJ+1 6730 IF(IPST(IJ).GT.0) GO TO 43 6740 WRITE(LOUT,214) IJ,I,J,IAP(IJ) 6750 214 FORMAT(8X,I4,7X,I4,8X,I4,7X,4HFREE,20X,I3) 6760 GO TO 42 6770 43 WRITE(LOUT,215) IJ,I,J,IPST(IJ),IAP(IJ) 6780 215 FORMAT(8X,I4,7X,I4,8X,I4,7X,8HEQUALITY,I3,13X,I3) 6790 42 CONTINUE 6800 GO TO 57 6810 35 DO 58 I=1,NBDMX 6820 58 IAP(I)=I 6830 57 CONTINUE 6840 C 6850 C**********************************************************************C 6860 C 6870 C OBTAIN SOLUTIONS WITH VARYING DIMENSIONALITIES 6880 C 6890 NRP=NDMX-NDMN+1 6900 DO 987 IJKL=1,NRP 6910 C 6920 ND=NDMX-IJKL+1 6930 NBD=NB*ND 6940 C 6950 WRITE(LOUT,800) ND 6960 800 FORMAT(1H1//12H SOLUTION IN,I3,2X,12HDIMENSION(S)) 6970 C 6980 C INITIAL ESTIMATES FOR SPECIFIC DIMENSIONALITY 6990 C (AND UNDER SPECIFIC CONSTRAINTS) 7000 C 7010 S=0.0 7020 DO 44 J=1,ND 7030 44 S=S+ALAM(J) 7040 S=FNB/S 7050 IJ=0 7060 DO 45 J=1,ND 7070 C=SQRT(ALAM(J)*S) 7080 DO 45 I=1,NB 7090 IJ=IJ+1 7100 X(IJ)=U(IJ)*C 7110 IF(NCON.EQ.1) GO TO 45 7120 Y(IJ)=X(IJ) 7130 45 CONTINUE 7140 NT=NBD 7150 IF(NCON.EQ.0) GO TO 46 7160 NT=0 7170 JJ=0 7180 DO 47 J=1,ND 7190 JJ1=JJ+1 7200 NTC=IAP(JJ1) 7210 I1=JJ+2 7220 JJ=JJ+NB 7230 DO 48 II=I1,JJ 7240 48 IF(IAP(II).GT.NTC) NTC=IAP(II) 7250 NTC1=NTC-NT 7260 DO 49 I=1,NTC1 7270 DT(I)=0.0 7280 49 DU(I)=0.0 7290 DO 250 II=JJ1,JJ 7300 K=IAP(II)-NT 7310 DT(K)=DT(K)+X(II) 7320 250 DU(K)=DU(K)+1.0 7330 DO 251 I=1,NTC1 7340 251 DT(I)=DT(I)/DU(I) 7350 DO 252 II=JJ1,JJ 7360 KK=IAP(II) 7370 K=KK-NT 7380 X(II)=DT(K) 7390 252 Y(KK)=DT(K) 7400 NT=NT+NTC1 7410 47 CONTINUE 7420 C 7430 46 NT1=NT+1 7440 NT2=NT1 7450 IF(INS.EQ.1) NT2=NT+NS 7460 C=-5.0 7470 IF(ME.EQ.1) C=-8.0 7480 DO 56 I=NT1,NT2 7490 56 Y(I)=C 7500 NT22=NT2 7510 IF(ME.EQ.0) NT2=NT2-1 7520 WK1(NT22)=Y(NT22) 7530 I1=NBD+1 7540 I2=NBD+NS 7550 IF(INS.EQ.1) GO TO 156 7560 DO 157 II=I1,I2 7570 157 IAP(II)=NT1 7580 GO TO 158 7590 156 DO 159 II=I1,I2 7600 159 IAP(II)=NT+II-NBD 7610 158 CONTINUE 7620 C 7630 C PRINT INITIAL ESTIMATES 7640 C 7650 IF(INTPR.EQ.0) GO TO 53 7660 WRITE(LOUT,217) 7670 217 FORMAT(//6X,21HINITIAL CONFIGURATION) 7680 WRITE(LOUT,348) 7690 DO 54 I=1,NB 7700 54 WRITE(LOUT,415) I,(X(L),L=I,NBD,NB) 7710 C 7720 53 WRITE(LOUT,219) 7730 219 FORMAT(///6X,17HITERATION HISTORY//10X,9HITERATION,3X,8HLOG-LIKE, 7740 16HLIHOOD,3X,11HIMPROVEMENT,3X,22HDIRECTIONAL DERIVATIVE,3X, 7750 221HCONVERGENCE CRITERION) 7760 WRITE(LOUT,348) 7770 C 7780 C ITERATION STARTS HERE 7790 C 7800 NG=0 7810 EPS=1.0E-4 7820 KRANK=NB*ND-(ND+1)*ND/2 7830 IF(INS.EQ.1) KRANK=KRANK+NS-1 IDDX=0 7840 NSQ2=NT2*(NT2+1)/2 7850 C 7860 LLL=0 7870 1 LLL=LLL+1 7880 IF(LLL.GT.MXIT) GO TO 2 7890 C 7900 CALL GRHS4(D,DD,DS,IY,IZ,IAP,NCR,NSC,NSR,F,Y,AH,AG,G,H,TM1,NBD, 7910 1NB,ND,NS,NT,NT1,NT2,NG) 7920 C 7930 IF(NG.EQ.1) GO TO 980 7940 IF(LLL.GT.1) GO TO 700 7950 L1=LLL-1 7960 WRITE(LOUT,301) L1,TM1 7970 700 CONTINUE 7980 C 7990 C GAUSS-NEWTON ITERATION 8000 C 8010 DO 570 I=1,NT2 8020 570 WK2(I)=G(I) 8030 CALL MFSS(H,NT2,EPS,IRANK,KRANK,WK3) 8040 IF(IRANK.GT.0) GO TO 87 8050 C 8060 WRITE(LOUT,300) IRANK 8070 300 FORMAT(//42H *** INFORMATION MATRIX IS ILL CONDITIONED,5X, 8080 16HIRANK=,I3) 8090 EPS=EPS*10.0 8100 LLL=LLL-1 8110 GO TO 1 8120 C 8130 87 CONTINUE 8140 CALL MLSS(H,NT2,IRANK,WK3,0,WK2,IER) 8150 DDRV=0.0 8160 DO 8009 I=1,NT2 8170 8009 DDRV=DDRV+WK2(I)*G(I) 8180 9999 STEP=STEPQ 8190 LQ=0 8200 1000 LQ=LQ+1 8210 IF(LQ.GT.10) GO TO 3 8220 DO 89 I=1,NT2 8230 89 WK1(I)=Y(I)+STEP*WK2(I) 8240 FL=FLL4(D,DD,DS,WK1,IY,IZ,IAP,NBD,NB,ND,NS,NCR,NSC,NSR,F) 8250 IF(IDDX.NE.0) GO TO 9998 8260 IF(FL.GT.TM1) GO TO 90 8270 WRITE(LOUT,8020) FL 8280 8020 FORMAT(20X,F14.5) 8290 STEP=STEP*0.5 8300 GO TO 1000 8310 9998 IF(LLL.NE.1) GO TO 9997 8320 STEPQ=STEPQ*0.1 8330 WRITE(LOUT,9996) STEPQ 8340 9996 FORMAT(6X,48HSTEP SIZE HAS BEEN MULTIPLIED BY 0.1. IT IS NOW, 8350 1F10.4) 8360 GO TO 9999 8370 9997 WRITE(LOUT,9995) 8380 9995 FORMAT(///40H THE ERROR VARIANCE IS GETTING TOO SMALL/ 8390 151H CONSTRAINTS ARE NOT STRONG ENOUGH FOR THIS PROBLEM) 8400 GO TO 987 8410 C 8420 C CONVERGENCE CHECK 8430 C 8440 90 CONTINUE 8450 DIF=FL-TM1 8460 C=-DIF/FL 8470 TM1=FL 8480 WRITE(LOUT,301) LLL,FL,DIF,DDRV,C 8490 301 FORMAT(10X,I5,5X,F14.5,1X,F14.5,6X,F14.5,8X,F16.7) 8500 C 8510 DO 9000 I=1,NT2 8520 9000 Y(I)=WK1(I) 8530 IF(C.GT.CVCR) GO TO 1 8540 C 8550 WRITE(LOUT,450) LLL 8560 450 FORMAT(///24H CONVERGENCE ATTAINED IN,I4,2X,10HITERATIONS//) 8570 GO TO 4 8580 2 WRITE(LOUT,451) 8590 451 FORMAT(///26H MAXIMUM ITERATION REACHED/ 8600 125H CONVERGENCE NOT ATTAINED//) 8610 GO TO 4 8620 3 WRITE(LOUT,452) 8630 452 FORMAT(6X,47HNO IMPROVEMENT HAS BEEN OBTAINED BY THE ASCENT , 8640 19HDIRECTION) 8650 WRITE(LOUT,454) 8660 454 FORMAT(///21H ITERATION TERMINATES) 8670 C 8680 C PRINT RESULTS 8690 C 8700 4 CONTINUE 8710 DO 9001 I=1,NBD 8720 II=IAP(I) 8730 9001 X(I)=WK1(II) 8740 WRITE(LOUT,410) ND 8750 410 FORMAT(40H1FINAL ESTIMATES OF MODEL PARAMETERS (IN,I3,2X, 8760 113HDIMENSION(S))) 8770 WRITE(LOUT,411) 8780 411 FORMAT(//6X,30HDERIVED STIMULUS CONFIGURATION) 8790 WRITE(LOUT,348) 8800 DO 286 I=1,NB 8810 286 WRITE(LOUT,415) I,(X(J),J=I,NBD,NB) 8820 415 FORMAT(1H ,I10,10F12.5) 8830 C=0.0 8840 N=0 8850 DO 97 J=1,ND 8860 FMDM=0.0 8870 FMDS=0.0 8880 N1=N+1 8890 N=N+NB 8900 DO 98 I=N1,N 8910 FMDM=FMDM+X(I) 8920 98 FMDS=FMDS+X(I)*X(I) 8930 FMDM=FMDM*FNBI 8940 FMDS=FMDS*FNBI-FMDM*FMDM 8950 C=C+FMDS 8960 97 CONTINUE 8970 C=SQRT(C) 8980 WRITE(LOUT,981) C 8990 981 FORMAT(/6X,51HTHE OVERALL SIZE OF THE CONFIGURATION AS INDICATED , 9000 121HBY SQRT(TR(X'X)/N) IS,F12.5) 9010 NG=1 9020 LLL=0 9030 IF(ILLT.EQ.1) WRITE(LOUT,982) 9040 982 FORMAT(///52H LOG-LIKELIHOOD OF EACH SUCCESSIVE FIRST CHOICE (AND, 9050 114H A RANK ORDER)//6X,14HRANK ORDER NO.,3X,4HRANK,5X,7HSTIMULI,5X, 9060 218HLOG-LIKELIHOOD(LL),10X,6H-2(LL)) 9070 WRITE(LOUT,348) 9080 GO TO 1 9090 980 CONTINUE 9100 CALL GINV(H,NT2,KRANK,WK3,WK4,EPS,IPST,0,IER) 9110 WRITE(LOUT,400) KRANK 9120 400 FORMAT(38H1THE EFFECTIVE NUMBER OF PARAMETERS IS,I5) 9130 TM1=-2.0*TM1 9140 WRITE(LOUT,402) TM1 9150 402 FORMAT(/25H THE -2*LOG-LIKELIHOOD IS,F12.4) 9160 AIC=TM1+2.0*FLOAT(KRANK) 9170 WRITE(LOUT,403) AIC 9180 403 FORMAT(/34H THE VALUE OF THE AIC STATISTIC IS,F12.4) 9190 BIC=TM1+FLOAT(KRANK)*ALOG(FNN) 9200 WRITE(LOUT,404) BIC 9210 404 FORMAT(/34H THE VALUE OF THE BIC STATISTIC IS,F12.4) 9220 C 9230 WRITE(LOUT,499) 9240 499 FORMAT(///50H STANDARD ERRORS OF ESTIMATED STIMULUS COORDINATES) 9250 WRITE(LOUT,498) 9260 498 FORMAT(//6X,8HSTIMULUS,4X,9HDIMENSION,3X,12HACTIVE PARA.,5X, 9270 18HESTIMATE,4X,9HSTD ERROR) 9280 WRITE(LOUT,348) 9290 N=0 9300 DO 500 J=1,ND 9310 DO 500 I=1,NB 9320 N=N+1 9330 N1=IAP(N) 9340 N2=(N1+1)*N1/2 9350 STD=SQRT(H(N2)) 9360 500 WRITE(LOUT,501) I,J,N1,X(N),STD 9370 501 FORMAT(4X,I7,5X,I7,7X,I7,5X,2F13.5) 9380 C 9390 WRITE(LOUT,420) 9400 420 FORMAT(///37H ESTIMATES OF DISPERSION PARAMETER(S)) 9410 WRITE(LOUT,424) 9420 424 FORMAT(//6X,11HREPLICATION,13X,12HACTIVE PARA.,5X,8HESTIMATE,4X, 9430 19HSTD ERROR,4X,36HVARIANCE EQUIVALENT(IN NORMAL DIST.)) 9440 WRITE(LOUT,348) 9450 IF(NT2.LT.NT1) GO TO 985 9460 N=(NT+1)*NT/2 9470 DO 425 I=NT1,NT2 9480 N=N+I 9490 STD=SQRT(H(N)) 9500 C=-Y(I) 9510 CO=3.289868/C**2 9520 IF(INS.EQ.0) GO TO 423 9530 I1=I-NT 9540 WRITE(LOUT,426) I1,I,C,STD,CO 9550 426 FORMAT(4X,I7,19X,I7,5X,2F13.5,5X,F13.5) 9560 GO TO 425 9570 423 WRITE(LOUT,421) I,C,STD,CO 9580 421 FORMAT(9X,3HALL,18X,I7,5X,2F13.5,5X,F13.5) 9590 425 CONTINUE 9600 985 IF(ME.EQ.1) GO TO 989 9610 C=-Y(NT22) 9620 CO=3.289868/C**2 9630 IF(INS.EQ.0) GO TO 986 9640 WRITE(LOUT,970) NS,C,CO 9650 970 FORMAT(4X,I7,23X,5HFIXED,3X,F13.5,18X,F13.5) 9660 GO TO 989 9670 986 WRITE(LOUT,988) C,CO 9680 988 FORMAT(9X,3HALL,22X,5HFIXED,3X,F13.5,18X,F13.5) 9690 989 CONTINUE 9700 C 9710 IF(IDIST.EQ.0) GO TO 427 9720 WRITE(LOUT,428) 9730 428 FORMAT(18H1DERIVED DISTANCES) 9740 N=0 9750 DO 430 I=2,NB 9760 I1=I-1 9770 DO 430 J=1,I1 9780 N=N+1 9790 D(N)=0.0 9800 KK=0 9810 DO 431 K=1,ND 9820 K1=KK+I 9830 K2=KK+J 9840 KK=KK+NB 9850 431 D(N)=D(N)+(X(K1)-X(K2))*(X(K1)-X(K2)) 9860 430 D(N)=SQRT(D(N)) 9870 N1=NB1/10 9880 NRS=NB1-N1*10 9890 IF(NRS.GT.0) GO TO 350 9900 N1=N1-1 9910 350 N1=N1+1 9920 II=1 9930 DO 351 K=1,N1 9940 I2=II+9 9950 IF(I2.GT.NB1) I2=NB1 9960 WRITE(LOUT,339) (J,J=II,I2) 9970 WRITE(LOUT,348) 9980 I0=(II+1)*II/2-1 9990 II1=II+1 10000 DO 352 I=II1,NB 10010 I1=I0+1 10020 III=I-II1 10030 IF(III.GT.9) III=9 10040 I2=I1+III 10050 WRITE(LOUT,533) I,(D(J),J=I1,I2) 10060 352 I0=I0+I-1 10070 351 II=II+10 10080 C 10090 C PLOT RESULTS 10100 C 10110 427 IF(IPLTC.EQ.0) GO TO 189 10120 ABP=0.0 10130 DO 285 I=1,NBD 10140 WA=ABS(X(I)) 10150 285 IF(WA.GT.ABP) ABP=WA 10160 ABP=ABP*1.1 10170 ABM=-ABP 10180 IF(ND.EQ.1) GO TO 187 10190 DO 188 I=2,ND 10200 I1=I-1 10210 J2=(I-1)*NB+1 10220 DO 188 J=1,I1 10230 J1=(J-1)*NB+1 10240 WRITE(LOUT,260) TIT,J,I 10250 260 FORMAT(39H1PLOT OF DERIVED STIMULUS CONFIGURATION,5X,20A4/ 10260 120X,9HDIMENSION,I3,2X,21H(X-AXIS) VS DIMENSION,I3,2X,8H(Y-AXIS)) 10270 188 CALL PLOTR(X(J1),X(J2),ABP,ABP,ABM,ABM,NB,-2,NB) 10280 GO TO 189 10290 187 WRITE(LOUT,261) TIT 10300 261 FORMAT(39H1PLOT OF DERIVED STIMULUS CONFIGURATION,5X,20A4/ 10310 120X,20HONE DIMENSIONAL PLOT) 10320 CALL PLOTR(X(1),X(1),ABP,ABP,ABM,ABM,NB,-2,NB) 10330 C 10340 C PUNCH RESULTS 10350 C 10360 189 IF(IPHC.EQ.0) GO TO 174 10370 DO 175 I=1,NB 10380 175 WRITE(NPH,265) (X(J),J=I,NBD,NB) 10390 265 FORMAT(8E10.3) 10400 174 IF(IPHV.EQ.0) GO TO 313 10410 I2=0 10420 DO 277 I=1,NT2 10430 I1=I2+1 10440 I2=I2+I 10450 277 WRITE(NPH,265) (H(II),II=I1,I2) 10460 313 IF(IVEST.EQ.0) GO TO 987 10470 WRITE(LOUT,531) 10480 531 FORMAT(55H1ASYMPTOTIC VARIANCE/COVARIANCE ESTIMATES OF ESTIMATED , 10490 154HPARAMETERS (NUMBERS INDICATE ACTIVE PARAMETER NUMBERS)) 10500 N1=NT2/10 10510 NRS=NT2-N1*10 10520 IF(NRS.GT.0) GO TO 336 10530 N1=N1-1 10540 336 N1=N1+1 10550 II=1 10560 DO 337 K=1,N1 10570 I2=II+9 10580 IF(I2.GT.NT2) I2=NT2 10590 WRITE(LOUT,339) (J,J=II,I2) 10600 339 FORMAT(//9X,10I12) 10610 WRITE(LOUT,348) 10620 348 FORMAT(1H ) 10630 I0=(II+1)*II/2-1 10640 DO 338 I=II,NT2 10650 I1=I0+1 10660 III=I-II 10670 IF(III.GT.9) III=9 10680 I2=I1+III 10690 WRITE(LOUT,533) I,(H(J),J=I1,I2) 10700 533 FORMAT(1H ,I10,10F12.5) 10710 I0=I0+I 10720 338 CONTINUE 10730 II=II+10 10740 337 CONTINUE 10750 987 CONTINUE 10760 C 10770 RETURN 10780 END 10790 $ FORTY FORM 10800 C 10810 SUBROUTINE GRHS4(D, DD, DS, IY, IZ, IAP, NCR, NSC, NSR, F, Y, AH, 10820 * AG, G, H, TM1, NBD, NB, ND, NS, NT, NT1, NT2, NG) 10830 C 10840 C EVALUATE LOG-LIKELIHOOD, GRADIENT VECTOR AND/OR INFORMATION MATRIX 10850 C 10860 REAL*8 TM1,TM2,AL,SQ,SSQ,TSQ,TSSQ,AH(1),AG(1) 10870 DIMENSION D(1), DD(1), DS(1), IY(1), IZ(1), IAP(1), NCR(1), 10880 * NSC(1), NSR(1), F(1), Y(1), G(1), H(1) 10890 COMMON /OPT/ IN, LOUT, NPH, TIT(20), FM(20), ME, INT, NCON, MXIT, 10900 * CVCR, ITIE, INS, IDATA, INTPR, IDIST, IVEST, IPLTC, IPHC, IPHV, 10910 * NBDMX1, NDMX1, NB2, NNB2, NBSQ2, NTIE, NBS, NMIS, NBDMX, ICODE, 10920 * NCODE, NB1, FNB1, FNB, FNBI, FNB2, ILLT 10930 COMMON /NUM/ NSQ2 10940 C SUPPRESS UNDERFLOW MESSAGES 10950 CALL FXOPT(67, 1, 1, 0) 10960 C 10970 C CLEAR ACCUMULATORS 10980 TM1 = 0.0D0 10990 DO 10 I=1,NT2 11000 G(I) = 0.0 11010 10 CONTINUE 11020 DO 20 I=1,NSQ2 11030 H(I) = 0.0 11040 20 CONTINUE 11050 C CALCULATE DISTANCES 11060 N = 0 11070 DO 50 I=2,NB 11080 I1 = I - 1 11090 DO 40 J=1,I1 11100 N = N + 1 11110 TM2 = 0.0D0 11120 KK = 0 11130 DO 30 K=1,ND 11140 K1 = KK + I 11150 K2 = KK + J 11160 KK = KK + NB 11170 KK1 = IAP(K1) 11180 KK2 = IAP(K2) 11190 TM2 = TM2 + (Y(KK1)-Y(KK2))*(Y(KK1)-Y(KK2)) 11200 30 CONTINUE 11210 D(N) = TM2 11220 D(N) = SQRT(D(N)) 11230 40 CONTINUE 11240 50 CONTINUE 11250 IF (ME.EQ.0) GO TO 70 11260 DO 60 I=1,NBSQ2 11270 DS(I) = D(I) 11280 D(I) = ALOG(DS(I)) 11290 60 CONTINUE 11300 70 CONTINUE 11310 C 11320 JS = NBD 11330 N = 0 11340 IK = 0 11350 DO 350 K=1,NS 11360 IF (INS.EQ.0 .AND. K.GT.1) GO TO 110 11370 JS = JS + 1 11380 KK = IAP(JS) 11390 C = Y(KK) 11400 IF (ME.EQ.1) GO TO 90 11410 DO 80 I=1,NBSQ2 11420 DD(I) = EXP(C*D(I)) 11430 80 CONTINUE 11440 GO TO 110 11450 90 DO 100 I=1,NBSQ2 11460 DD(I) = DS(I)**C 11470 100 CONTINUE 11480 110 CONTINUE 11490 NCRK = NCR(K) 11500 DO 340 II=1,NCRK 11510 N = N + 1 11520 JJ = NSR(N) 11530 JJ1 = NSC(N) 11540 TSQ = 0.0D0 11550 TSSQ = 0.0D0 11560 DO 120 IP=1,NT 11570 AH(IP) = 0.0D0 11580 120 CONTINUE 11590 JJP1 = JJ + 1 11600 IF (JJP1.GT.JJ1) GO TO 150 11610 DO 140 J=JJP1,JJ1 11620 JY = IK + J 11630 JY = JY + JY 11640 J2 = IZ(JY) 11650 JY = JY - 1 11660 J1 = IZ(JY) 11670 JK = (J1-2)*(J1-1)/2 + J2 11680 TSQ = TSQ + DD(JK) 11690 TSSQ = TSSQ + D(JK)*DD(JK) 11700 DJK = D(JK) 11710 IF (ME.EQ.1) DJK = DS(JK)*DS(JK) 11720 YL = DD(JK)/DJK 11730 LL = 0 11740 DO 130 L=1,ND 11750 L1 = LL + J1 11760 L2 = LL + J2 11770 LL = LL + NB 11780 L11 = IAP(L1) 11790 L22 = IAP(L2) 11800 YLD = (Y(L11)-Y(L22))*YL 11810 AH(L11) = AH(L11) + YLD 11820 AH(L22) = AH(L22) - YLD 11830 130 CONTINUE 11840 140 CONTINUE 11850 C 11860 150 JX1 = 0 11870 ALT = 0.0 11880 DO 330 JA=1,JJ 11890 J = JJP1 - JA 11900 IF (NT2.LT.NT1) GO TO 170 11910 DO 160 L=NT1,NT2 11920 AG(L) = 0.0D0 11930 160 CONTINUE 11940 170 CONTINUE 11950 JY1 = IK + J 11960 JY = JY1 + JY1 11970 J2 = IZ(JY) 11980 JY = JY - 1 11990 J1 = IZ(JY) 12000 JK = (J1-2)*(J1-1)/2 + J2 12010 TSQ = TSQ + DD(JK) 12020 TSSQ = TSSQ + D(JK)*DD(JK) 12030 SQ = TSQ 12040 SSQ = TSSQ 12050 DJK = D(JK) 12060 IF (ME.EQ.1) DJK = DS(JK)*DS(JK) 12070 YL = DD(JK)/DJK 12080 LL = 0 12090 DO 180 L=1,ND 12100 L1 = LL + J1 12110 L2 = LL + J2 12120 LL = LL + NB 12130 L11 = IAP(L1) 12140 L22 = IAP(L2) 12150 YLD = (Y(L11)-Y(L22))*YL 12160 AH(L11) = AH(L11) + YLD 12170 AH(L22) = AH(L22) - YLD 12180 180 CONTINUE 12190 DO 190 IP=1,NT 12200 AG(IP) = AH(IP) 12210 190 CONTINUE 12220 C 12230 IF (NTIE.EQ.0) GO TO 250 12240 IF (IY(JY1).EQ.1) GO TO 250 12250 IF (ITIE.EQ.1) GO TO 220 12260 JX1 = JX1 + 1 12270 IF (JX1.EQ.1) GO TO 250 12280 I1 = J + 1 12290 I2 = J + JX1 - 1 12300 DO 210 M=I1,I2 12310 MY = IK + M 12320 MY = MY + MY 12330 M2 = IZ(MY) 12340 MY = MY - 1 12350 M1 = IZ(MY) 12360 IM = (M1-2)*(M1-1)/2 + M2 12370 SQ = SQ - DD(IM) 12380 SSQ = SSQ - D(IM)*DD(IM) 12390 DIM = D(IM) 12400 IF (ME.EQ.1) DIM = DS(IM)*DS(IM) 12410 YL = DD(IM)/DIM 12420 LL = 0 12430 DO 200 L=1,ND 12440 L1 = LL + M1 12450 L2 = LL + M2 12460 LL = LL + NB 12470 L11 = IAP(L1) 12480 L22 = IAP(L2) 12490 YLD = (Y(L11)-Y(L22))*YL 12500 AG(L11) = AG(L11) - YLD 12510 AG(L22) = AG(L22) + YLD 12520 200 CONTINUE 12530 210 CONTINUE 12540 IF (JX1.EQ.IY(JY1)) JX1 = 0 12550 GO TO 250 12560 220 JX1 = JX1 + 1 12570 IF (JX1.EQ.IY(JY1)) JX1 = 0 12580 IF (JX1.EQ.0) GO TO 250 12590 I1 = J - IY(JY1) + JX1 12600 I2 = J - 1 12610 DO 240 M=I1,I2 12620 MY = IK + M 12630 MY = MY + MY 12640 M2 = IZ(MY) 12650 MY = MY - 1 12660 M1 = IZ(MY) 12670 IM = (M1-2)*(M1-1)/2 + M2 12680 SQ = SQ + DD(IM) 12690 SSQ = SSQ + D(IM)*DD(IM) 12700 DIM = D(IM) 12710 IF (ME.EQ.1) DIM = DS(IM)*DS(IM) 12720 YL = DD(IM)/DIM 12730 LL = 0 12740 DO 230 L=1,ND 12750 L1 = LL + M1 12760 L2 = LL + M2 12770 LL = LL + NB 12780 L11 = IAP(L1) 12790 L22 = IAP(L2) 12800 YLD = (Y(L11)-Y(L22))*YL 12810 AG(L11) = AG(L11) + YLD 12820 AG(L22) = AG(L22) - YLD 12830 230 CONTINUE 12840 240 CONTINUE 12850 250 AL = C*D(JK) - DLOG(SQ) 12860 IF (NG.NE.1 .OR. ILLT.EQ.0) GO TO 260 12870 IF (J.EQ.JJ1 .AND. IY(JY1).EQ.1) GO TO 260 12880 AL2 = -2.0*AL 12890 ALT = ALT + AL2 12900 WRITE (LOUT,9999) N, J, J1, J2, AL, AL2 12910 9999 FORMAT (6X, I7, 8X, I4, 7X, 1H(, I2, 1H,, I2, 1H), 6X, F13.5, 7X, 12920 * F13.5) 12930 IF (AL2.GT.6.63) WRITE (LOUT,9998) 12940 9998 FORMAT (1H+, 77X, 3H***) 12950 IF (JA.NE.JJ) GO TO 260 12960 IF (ICODE.LT.3) GO TO 260 12970 WRITE (LOUT,9997) N, ALT 12980 9997 FORMAT (/6X, 20HTOTAL FOR RANK ORDER, I4, 35X, F13.5) 12990 WRITE (LOUT,9996) 13000 9996 FORMAT (1H ) 13010 260 CONTINUE 13020 AG(KK) = (D(JK)-SSQ/SQ)*F(N) 13030 DO 270 IP=1,NT 13040 AG(IP) = -AG(IP)/SQ 13050 270 CONTINUE 13060 C 13070 DJK = D(JK) 13080 IF (ME.EQ.1) DJK = DS(JK)*DS(JK) 13090 LL = 0 13100 DO 280 L=1,ND 13110 L1 = LL + J1 13120 L2 = LL + J2 13130 LL = LL + NB 13140 L11 = IAP(L1) 13150 L22 = IAP(L2) 13160 YLD = (Y(L11)-Y(L22))/DJK 13170 AG(L11) = AG(L11) + YLD 13180 AG(L22) = AG(L22) - YLD 13190 280 CONTINUE 13200 C 13210 DO 290 IP=1,NT 13220 AG(IP) = AG(IP)*C*F(N) 13230 290 CONTINUE 13240 C 13250 TM1 = TM1 + AL*F(N) 13260 AL = DEXP(AL) 13270 AL = AL/((1.0D0-AL)*F(N)) 13280 DO 300 L=1,NT2 13290 G(L) = G(L) + AG(L) 13300 300 CONTINUE 13310 L12 = 0 13320 DO 320 L1=1,NT2 13330 DO 310 L2=1,L1 13340 L12 = L12 + 1 13350 H(L12) = H(L12) + AG(L1)*AG(L2)*AL 13360 310 CONTINUE 13370 320 CONTINUE 13380 330 CONTINUE 13390 IK = IK + JJ1 13400 340 CONTINUE 13410 350 CONTINUE 13420 RETURN 13430 END 13440 $ FORTY FORM 13450 C 13460 DOUBLE PRECISION FUNCTION FLL4(D,DD,DS,Y,IY,IZ,IAP,NBD,NB,ND, 13470 1NS,NCR,NSC,NSR,F) 13480 C 13490 C EVALUATE THE LOG-LIKELIHOOD 13500 C 13510 REAL*8 FLL4,S,SS 13520 DIMENSION D(1),DD(1),DS(1),Y(1),IAP(1),IY(1),IZ(1) 13530 DIMENSION NCR(1),NSC(1),NSR(1),F(1) 13540 COMMON/OPT/IN,LOUT,NPH,TIT(20),FM(20),ME,INT,NCON,MXIT,CVCR,ITIE, 13550 1INS,IDATA,INTPR,IDIST,IVEST,IPLTC,IPHC,IPHV,NBDMX1,NDMX1,NB2,NNB2, 13560 2NBSQ2,NTIE,NBS,NMIS,NBDMX,ICODE,NCODE,NB1,FNB1,FNB,FNBI,FNB2,ILLT 13570 COMMON/NUM/NSQ2,STEPQ,IDDX 13580 FLL4=0.0D0 13590 C CALCULATE DISTANCES 13600 N=0 13610 DO 10 I=2,NB 13620 I1=I-1 13630 DO 10 J=1,I1 13640 N=N+1 13650 S=0.0D0 13660 KK=0 13670 DO 11 K=1,ND 13680 K1=KK+I 13690 K2=KK+J 13700 KK=KK+NB 13710 KK1=IAP(K1) 13720 KK2=IAP(K2) 13730 11 S=S+(Y(KK1)-Y(KK2))*(Y(KK1)-Y(KK2)) 13740 D(N)=S 13750 10 D(N)=SQRT(D(N)) 13760 IF(ME.EQ.0) GO TO 21 13770 DO 20 I=1,NBSQ2 13780 DS(I)=D(I) 13790 20 D(I)=ALOG(DS(I)) 13800 21 CONTINUE 13810 C 13820 N=0 13830 IK=0 13840 DO 12 K=1,NS 13850 IF(INS.EQ.0.AND.K.GT.1) GO TO 14 13860 JS=NBD+K 13870 KK=IAP(JS) 13880 C=Y(KK) 13890 IF(ME.EQ.1) GO TO 15 13900 DO 13 I=1,NBSQ2 13910 CD=C*D(I) 13920 IF(CD.LT.-88.0) GO TO 8000 13930 13 DD(I)=EXP(CD) 13940 IDDX=0 13950 GO TO 14 13960 8000 IDDX=IDDX+1 13970 RETURN 13980 15 DO 16 I=1,NBSQ2 13990 16 DD(I)=DS(I)**C 14000 14 CONTINUE 14010 NCRK=NCR(K) 14020 DO 17 II=1,NCRK 14030 N=N+1 14040 JJ=NSR(N) 14050 JJ1=NSC(N) 14060 S=0.0D0 14070 JJP1=JJ+1 14080 IF(JJP1.GT.JJ1) GO TO 201 14090 DO 200 J=JJP1,JJ1 14100 JY=IK+J 14110 JY=JY+JY 14120 J2=IZ(JY) 14130 JY=JY-1 14140 J1=IZ(JY) 14150 JK=(J1-2)*(J1-1)/2+J2 14160 200 S=S+DD(JK) 14170 201 JX=0 14180 DO 18 JA=1,JJ 14190 J=JJP1-JA 14200 JY1=IK+J 14210 JY=JY1+JY1 14220 J2=IZ(JY) 14230 JY=JY-1 14240 J1=IZ(JY) 14250 JK=(J1-2)*(J1-1)/2+J2 14260 FLL4=FLL4+C*D(JK)*F(N) 14270 S=S+DD(JK) 14280 SS=S 14290 IF(NTIE.EQ.0) GO TO 111 14300 IF(IY(JY1).EQ.1) GO TO 111 14310 IF(ITIE.EQ.1) GO TO 113 14320 JX=JX+1 14330 IF(JX.EQ.1) GO TO 111 14340 I1=J+1 14350 I2=J+JX-1 14360 DO 114 M=I1,I2 14370 MY=IK+M 14380 MY=MY+MY 14390 M2=IZ(MY) 14400 MY=MY-1 14410 M1=IZ(MY) 14420 IM=(M1-2)*(M1-1)/2+M2 14430 114 SS=SS-DD(IM) 14440 IF(JX.EQ.IY(JY1)) JX=0 14450 GO TO 111 14460 113 JX=JX+1 14470 IF(JX.EQ.IY(JY1)) JX=0 14480 IF(JX.EQ.0) GO TO 111 14490 I1=J-IY(JY1)+JX 14500 I2=J-1 14510 DO 115 M=I1,I2 14520 MY=IK+M 14530 MY=MY+MY 14540 M2=IZ(MY) 14550 MY=MY-1 14560 M1=IZ(MY) 14570 IM=(M1-2)*(M1-1)/2+M2 14580 115 SS=SS+DD(IM) 14590 111 FLL4=FLL4-DLOG(SS)*F(N) 14600 18 CONTINUE 14610 IK=IK+JJ1 14620 17 CONTINUE 14630 12 CONTINUE 14640 RETURN 14650 END 14660 $ FORTY FORM 14670 C 14680 SUBROUTINE TETRA(ID, FQ, NCR, NSC, NSR, IZ, IY, F, NB, NS, NR, 14690 * NTS) 14700 C 14710 C DATA INPUT FOR TETRAD/TRIAD JUDGMENTS 14720 C 14730 DIMENSION ID(1), FQ(1), NCR(1), NSC(1), NSR(1), IZ(2,1), IY(1), 14740 * F(1) 14750 COMMON /OPT/ IN, LOUT, NPH, TIT(20), FM(20), ME, INT, NCON, MXIT, 14760 * CVCR, ITIE, INS, IDATA, INTPR, IDIST, IVEST, IPLTC, IPHC, IPHV, 14770 * NBDMX1, NDMX1, NB2, NNB2, NBSQ2, NTIE, NBS, NMIS, NBDMX, ICODE, 14780 * NCODE, NB1, FNB1, FNB, FNBI, FNB2, ILLT 14790 C 14800 C DATA INPUT 14810 C 14820 FNR = NR 14830 I4 = 0 14840 J3 = 0 14850 DO 180 K=1,NS 14860 NTK = NCR(K) 14870 DO 170 I=1,NTK 14880 I1 = I4 + 1 14890 I2 = I4 + 2 14900 I3 = I4 + 3 14910 I4 = I4 + 4 14920 J1 = J3 + 1 14930 J2 = J3 + 2 14940 J3 = J3 + 3 14950 IF (ICODE.EQ.2) GO TO 10 14960 C TETRAD INPUT 14970 READ (IN,FM) (ID(J),J=I1,I4), (FQ(J),J=J1,J3) 14980 GO TO 20 14990 C TRIAD INPUT 15000 10 READ (IN,FM) ID(I1), ID(I2), ID(I4), (FQ(J),J=J1,J3) 15010 ID(I3) = ID(I1) 15020 20 IF (NCODE.EQ.4) GO TO 60 15030 IF (NCODE-2) 30, 40, 50 15040 30 FQ(J3) = 0.0 15050 FQ(J2) = FNR - FQ(J1) 15060 GO TO 60 15070 40 FQ(J3) = 0.0 15080 GO TO 60 15090 50 FQ(J3) = FNR - FQ(J1) - FQ(J2) 15100 C DATA CHECK 15110 60 DO 70 J=I1,I4 15120 IF (ID(J).LT.1 .OR. ID(J).GT.NB) GO TO 80 15130 70 CONTINUE 15140 GO TO 90 15150 80 WRITE (LOUT,9999) 15160 9999 FORMAT (///34H STIMULUS INDICES ARE OUT OF BOUND/13H EITHER NONPO, 15170 * 41HSITIVE INDICES OR INDICES GREATER THAN NB) 15180 GO TO 140 15190 90 IF (FQ(J1).LT.0.0) GO TO 100 15200 IF (FQ(J2).LT.0.0) GO TO 100 15210 IF (FQ(J3).LT.0.0) GO TO 100 15220 GO TO 110 15230 100 WRITE (LOUT,9998) 15240 9998 FORMAT (///21H NEGATIVE FREQUENCIES/26H EITHER NR), NR)+F, 15250 * 31H(<) OR NEGATIVE FREQUENCY INPUT) 15260 GO TO 140 15270 110 IF (ID(I1).EQ.ID(I2)) GO TO 130 15280 IF (ID(I3).EQ.ID(I4)) GO TO 130 15290 IF (ICODE.EQ.2) GO TO 120 15300 IF (ID(I1).EQ.ID(I3)) GO TO 120 15310 GO TO 160 15320 120 IF (ID(I2).NE.ID(I4)) GO TO 160 15330 130 WRITE (LOUT,9997) 15340 9997 FORMAT (///48H NONPERMISSIBLE COMBINATIONS OF STIMULUS INDICES/ 15350 * 31H EITHER I=J, K=L OR (I,J)=(K,L)) 15360 140 IF (NS.EQ.1) GO TO 150 15370 WRITE (LOUT,9995) K 15380 150 WRITE (LOUT,9994) 15390 WRITE (LOUT,9992) I, (ID(J),J=I1,I4), (FQ(J),J=J1,J3) 15400 STOP 15410 160 IF (FQ(J3).NE.0) NTIE = 1 15420 170 CONTINUE 15430 180 CONTINUE 15440 C 15450 C DATA OUTPUT 15460 C 15470 IF (IDATA.EQ.0) GO TO 220 15480 WRITE (LOUT,9996) 15490 9996 FORMAT (1H1, //36H INPUT DATA (TETRAD/TRIAD JUDGMENTS)) 15500 I4 = 0 15510 J3 = 0 15520 DO 210 K=1,NS 15530 IF (NS.EQ.1) GO TO 190 15540 WRITE (LOUT,9995) K 15550 9995 FORMAT (//6X, 10H DATA SET=, I3) 15560 190 WRITE (LOUT,9994) 15570 9994 FORMAT (/11X, 3HNO., 3X, 16HI J K L, 5X, 5HIJ>KL, 5X, 15580 * 5HIJ 0) OR 30990 C STORE IN WK (ISTORE =0) 31000 IP = 0 31010 M4 = 0 31020 IF (ISTORE.GT.0) REWIND ISTORE 31030 DO 460 I=1,K 31040 JP = 0 31050 M3 = NPKPK 31060 DO 390 J=1,K 31070 SUM = 0 31080 LP = 0 31090 M1 = IP 31100 DO 380 L=1,I 31110 M1 = M1 + 1 31120 IF (L-J) 360, 350, 350 31130 350 M2 = LP + J 31140 GO TO 370 31150 360 M2 = JP + L 31160 370 SUM = SUM + A(M1)*B(M2) 31170 LP = LP + L 31180 380 CONTINUE 31190 M3 = M3 + 1 31200 B(M3) = SUM 31210 JP = JP + J 31220 390 CONTINUE 31230 M3 = NPKPK 31240 M1 = IP 31250 DO 400 J=1,I 31260 M1 = M1 + 1 31270 M3 = M3 + 1 31280 A(M1) = B(M3) 31290 400 CONTINUE 31300 IF (K-I) 450, 450, 410 31310 410 IP1 = I + 1 31320 IF (ISTORE) 420, 420, 440 31330 420 DO 430 J=IP1,K 31340 M3 = M3 + 1 31350 M4 = M4 + 1 31360 WK(M4) = B(M3) 31370 430 CONTINUE 31380 GO TO 450 31390 440 M3 = M3 + 1 31400 WRITE (ISTORE) (B(M),M=M3,NPKP2K) 31410 450 IP = IP + I 31420 460 CONTINUE 31430 C STORE UPPER TRIANGLE OF R IN B 31440 IF (ISTORE) 470, 470, 490 31450 470 DO 480 M=1,NPK 31460 B(M) = WK(M) 31470 480 CONTINUE 31480 GO TO 510 31490 490 REWIND ISTORE 31500 M2 = 0 31510 DO 500 I=1,KM 31520 M1 = M2 + 1 31530 M2 = M2 + K - I 31540 READ (ISTORE) (B(M),M=M1,M2) 31550 500 CONTINUE 31560 C COMPUTE Y = R-TRANSPOSE R AND STORE IN A 31570 510 DO 610 J=1,K 31580 DO 590 I=J,K 31590 SUM = 0 31600 LPA = 0 31610 LPB = 0 31620 DO 580 L=1,K 31630 LMI = L - I 31640 IF (LMI) 520, 560, 560 31650 520 M1 = LPB - LMI 31660 LMJ = L - J 31670 IF (LMJ) 530, 540, 540 31680 530 M2 = LPB - LMJ 31690 SUM = SUM + B(M1)*B(M2) 31700 GO TO 550 31710 540 M2 = LPA + J 31720 SUM = SUM + B(M1)*A(M2) 31730 550 LPB = LPB + K - L 31740 GO TO 570 31750 560 M1 = LPA + I 31760 M2 = LPA + J 31770 SUM = SUM + A(M1)*A(M2) 31780 570 LPA = LPA + L 31790 580 CONTINUE 31800 WK(I) = SUM 31810 590 CONTINUE 31820 IPA = J*(J-1)/2 31830 DO 600 I=J,K 31840 M = IPA + J 31850 A(M) = WK(I) 31860 IPA = IPA + I 31870 600 CONTINUE 31880 610 CONTINUE 31890 IF (K.EQ.N) RETURN 31900 C COMPUTE U Y 31910 IP = NPKPK 31920 DO 670 I=KP,N 31930 M3 = IP 31940 JP = 0 31950 DO 660 J=1,K 31960 SUM = 0 31970 M1 = IP 31980 LP = 0 31990 DO 650 L=1,K 32000 M1 = M1 + 1 32010 IF (L-J) 630, 620, 620 32020 620 M2 = LP + J 32030 GO TO 640 32040 630 M2 = JP + L 32050 640 SUM = SUM + A(M1)*A(M2) 32060 LP = LP + L 32070 650 CONTINUE 32080 M3 = M3 + 1 32090 B(M3) = SUM 32100 JP = JP + J 32110 660 CONTINUE 32120 IP = IP + I 32130 670 CONTINUE 32140 C COMPUTE U Y U-TRANSPOSE 32150 IP = NPKPK 32160 DO 700 I=KP,N 32170 JP = NPKPK 32180 M3 = IP + K 32190 DO 690 J=KP,I 32200 SUM = 0 32210 M1 = IP 32220 M2 = JP 32230 DO 680 L=1,K 32240 M1 = M1 + 1 32250 M2 = M2 + 1 32260 SUM = SUM + B(M1)*A(M2) 32270 680 CONTINUE 32280 M3 = M3 + 1 32290 A(M3) = SUM 32300 JP = JP + J 32310 690 CONTINUE 32320 IP = IP + I 32330 700 CONTINUE 32340 IP = NPKPK 32350 DO 720 I=KP,N 32360 M = IP 32370 DO 710 J=1,K 32380 M = M + 1 32390 A(M) = B(M) 32400 710 CONTINUE 32410 IP = IP + I 32420 720 CONTINUE 32430 C RE-ORDER ROWS AND COLUMNS AND RETURN MOORE-PENROSE INVERSE IN A 32440 730 DO 750 IR=1,NM1 32450 I = N - IR 32460 ISUB = I*(I+1)/2 32470 KPIV = IPIV(I) 32480 IF (KPIV-I) 760, 750, 740 32490 740 KSUB = KPIV*(KPIV+1)/2 32500 CALL PERMUT(A, N, KPIV, KSUB, I, ISUB) 32510 750 CONTINUE 32520 RETURN 32530 760 WRITE (6,9998) KPIV, I 32540 9998 FORMAT (29H0TROUBLE WITH INTERCHANGE ..., 2I3) 32550 STOP 32560 END 32570 $ FORTY FORM 32580 SUBROUTINE PERMUT(A,N,KPIV,KSUB,I,ISUB) 32590 DIMENSION A(1) 32600 KMI=KPIV-I 32610 JI=KSUB-KMI 32620 IDC=JI-ISUB 32630 JJ=ISUB-I+1 32640 DO 1 K=JJ,ISUB 32650 KK=K+IDC 32660 HOLD=A(K) 32670 A(K)=A(KK) 32680 1 A(KK)=HOLD 32690 KK=KSUB 32700 DO 2 K=KPIV,N 32710 II=KK-KMI 32720 HOLD=A(KK) 32730 A(KK)=A(II) 32740 A(II)=HOLD 32750 2 KK=KK+K 32760 JJ=KPIV-1 32770 II=ISUB 32780 DO 3 K=I,JJ 32790 HOLD=A(II) 32800 A(II)=A(JI) 32810 A(JI)=HOLD 32820 II=II+K 32830 3 JI=JI+1 32840 RETURN 32850 END 32860 $ FORTY FORM 32870 SUBROUTINE TINV(A,N) 32880 C THIS SUBROUTINE INVERTS A TRIANGULAR MATRIX 32890 DIMENSION A(1) 32900 REAL*8 X,Y 32910 INTEGER P,R,S,T,U 32920 M=0 32930 DO 1 I=1,N 32940 M=M+I 32950 1 A(M)=1.0/A(M) 32960 P=0 32970 R=0 32980 T=0 32990 IF(N.EQ.1)RETURN 33000 DO 2 I=2,N 33010 P=P+1 33020 R=R+I 33030 Y=A(R+1) 33040 DO 2 J=2,I 33050 P=P+1 33060 T=T+1 33070 S=T 33080 U=I-2 33090 X=0 33100 DO 3 L=P,R 33110 K=R-L+P 33120 X=X-A(K)*A(S) 33130 S=S-U 33140 3 U=U-1 33150 2 A(P)=X*Y 33160 RETURN 33170 END 33180 $ FORTY FORM 33190 C 33200 SUBROUTINE SHEL9(A,C,NITEM) 33210 C 33220 C SORTING 33230 C 33240 C C(1) : C(I)=I WHEN ENTERING THE ROUTINE. 33250 C C(I) WILL HOLD THE INDEX NUMBER K OF THE I'TH 33260 C SMALLEST OBSERVATION A(K) WHEN EXIT. 33270 C A(1) : A VECTOR OF ITEMS ON WHICH SORTING IS 33280 C PERFORMED. 33290 C NITEM: THE NUMBER OF ITEMS TO BE SORTED. 33300 C 33310 C REMARKS: SORTING IS IN AN ASCENDING ORDER. 33320 C A WILL BE DESTROYED. 33330 C 33340 INTEGER C(1) 33350 REAL A(1) 33360 C 33370 M=NITEM 33380 20 M=M/2 33390 IF(M) 30,40,30 33400 30 K=NITEM-M 33410 J=1 33420 41 I=J 33430 49 L=I+M 33440 IF(A(L)-A(I)) 50,60,60 33450 50 B=A(I) 33460 A(I)=A(L) 33470 A(L)=B 33480 KK=C(I) 33490 C(I)=C(L) 33500 C(L)=KK 33510 I=I-M 33520 IF(I-1) 60,49,49 33530 60 J=J+1 33540 IF(J-K) 41,41,20 33550 40 CONTINUE 33560 RETURN 33570 END 33580 $ FORTY FORM 33590 C 33600 SUBROUTINE BLOC1(A,L,LL,N) 33610 C 33620 C IDENTIFY TIED BLOCKS 33630 C 33640 C A(1) : A VECTOR OF ITEMS ON WHICH BLOCKS OF TIED 33650 C OBSERVATIONS ARE IDENTIFIED. 33660 C L(1) : L(I) CONTAINS THE INDEX NUMBER OF THE I'TH 33670 C SMALLEST OBSERVATION A(K). 33680 C LL(1): LL(I) WILL BE FILLED WITH THE NUMBER OF 33690 C OBSERVATIONS EQUAL TO THE I'TH SMALLEST 33700 C OBSERVATION (INCLUDING ITSELF). 33710 C 33720 C REMARK: LL WILL BE CHANGED. 33730 C 33740 DIMENSION A(1),L(1),LL(1) 33750 C 33760 I=1 33770 3 II=I 33780 2 II=II+1 33790 IF(II.GT.N) GO TO 7 33800 K=L(I) 33810 KK=L(II) 33820 IF(A(K).EQ.A(KK)) GO TO 2 33830 5 ID=II-1 33840 DO 10 J=I,ID 33850 10 LL(J)=II-I 33860 IF(II.EQ.N) GO TO 6 33870 I=II 33880 GO TO 3 33890 6 LL(II)=1 33900 RETURN 33910 7 ID=II-1 33920 DO 11 J=I,ID 33930 11 LL(J)=II-I 33940 RETURN 33950 END 33960 $ FORTY FORM 33970 C 33980 SUBROUTINE CJEIG(A,U,V,N,ND,B,W,ALAM,NFT,NA,NB,BK) 33990 C 34000 C 34010 C EIGENVALUES AND VECTORS OF REAL SYMMETRIC MATRIX 34020 C BY CLINT AND JENNINGS 34030 C 34040 DIMENSION A(NA,NA),U(NA,NB),V(NA,NB),B(NB,NB),W(NA,NB),ALAM(NB) 34050 DIMENSION BK(NB,NB) 34060 C 34070 IF(NFT.NE.1) GO TO 11 34080 DO 35 J=1,ND 34090 DO 36 I=1,N 34100 36 U(I,J)=0.0 34110 35 U(J,J)=1.0 34120 11 CONTINUE 34130 C 34140 LL=0 34150 1 LL=LL+1 34160 IF(LL.GT.50) RETURN 34170 DO 10 I=1,N 34180 DO 10 J=1,ND 34190 V(I,J)=0.0 34200 DO 10 K=1,N 34210 10 V(I,J)=V(I,J)+A(I,K)*U(K,J) 34220 DO 13 I=1,ND 34230 DO 13 J=1,ND 34240 BK(I,J)=0.0 34250 DO 13 K=1,N 34260 13 BK(I,J)=BK(I,J)+U(K,I)*V(K,J) 34270 C 34280 CALL JAC (ND,BK,B,ALAM,NB) 34290 C 34300 DO 14 I=1,N 34310 DO 14 J=1,ND 34320 W(I,J)=0.0 34330 DO 14 K=1,ND 34340 14 W(I,J)=W(I,J)+V(I,K)*B(K,J) 34350 DO 15 I=1,ND 34360 DO 15 J=1,ND 34370 B(I,J)=0.0 34380 DO 15 K=1,N 34390 15 B(I,J)=B(I,J)+W(K,I)*W(K,J) 34400 C 34410 C CHOLESKY FACTORIZATION 34420 C 34430 DO 16 I=1,ND 34440 B(I,I)=SQRT(B(I,I)) 34450 IF(I-ND) 17,18,18 34460 17 J1=I+1 34470 DO 16 J=J1,ND 34480 B(I,J)=B(I,J)/B(I,I) 34490 DO 16 K=J1,J 34500 16 B(K,J)=B(K,J)-B(I,J)*B(I,K) 34510 18 CONTINUE 34520 C 34530 C INVERSE OF CHOLESKY FACTOR 34540 C 34550 DO 60 I=1,ND 34560 B(I,I)=1.0/B(I,I) 34570 IF(I-1) 60,60,61 34580 61 J2=I-1 34590 DO 62 J=1,J2 34600 B(I,J)=0.0 34610 DO 63 K=J,J2 34620 63 B(I,J)=B(I,J)+B(K,I)*B(K,J) 34630 B(I,J)=-B(I,J)*B(I,I) 34640 62 CONTINUE 34650 60 CONTINUE 34660 C 34670 DO 30 J=1,ND 34680 DO 30 I=1,N 34690 30 V(I,J)=W(I,J) 34700 DO 19 I=1,N 34710 DO 19 J=1,ND 34720 W(I,J)=0.0 34730 DO 19 K=1,J 34740 19 W(I,J)=W(I,J)+V(I,K)*B(J,K) 34750 C 34760 C TEST OF CONVERGENCE 34770 C 34780 DO 20 J=1,ND 34790 IF(ALAM(J).EQ.0.0) GO TO 20 34800 ABSA=ABS(ALAM(J)) 34810 SG=ABSA/ALAM(J) 34820 DO 29 I=1,N 34830 DIF=SG*W(I,J)-U(I,J) 34840 IF(ABS(DIF).GT.0.00005) GO TO 21 34850 29 CONTINUE 34860 20 CONTINUE 34870 DO 22 I=1,N 34880 DO 22 J=1,ND 34890 22 U(I,J)=W(I,J) 34900 RETURN 34910 21 DO 23 I=1,N 34920 DO 23 J=1,ND 34930 23 U(I,J)=W(I,J) 34940 GO TO 1 34950 END 34960 $ FORTY FORM 34970 C 34980 SUBROUTINE JAC(N,R,RD,D,NBMX) 34990 C SIMULTANEOUS METHOD 35000 C EIGENVALUES AND VECTORS BY JACOBI METHOD 35010 DIMENSION R(NBMX,NBMX),RD(NBMX,NBMX),D(NBMX) 35020 NL=100 35030 P=0.00001 35040 N11=N-1 35050 N2=N*2 35060 DO 10 I=1,N 35070 DO 10 J=1,N 35080 10 RD(I,J)=0.0 35090 DO 11 I=1,N 35100 11 RD(I,I)=1.0 35110 L=0 35120 14 DO 12 I=1,N 35130 12 D(I)=R(I,I) 35140 DO 13 I=1,N11 35150 I1=I+1 35160 DO 13 J=I1,N 35170 DR=R(I,I)-R(J,J) 35180 A=SQRT(DR**2+4.0*R(I,J)**2) 35190 IF((A+DR).LT.0.0) GO TO 33 35200 A=SQRT((A+DR)/(2.0*A)) 35210 GO TO 32 35220 33 A=0.0 35230 32 CONTINUE 35240 B=SQRT(1.0-A**2) 35250 C=SIGN(1.0,R(I,J)) 35260 DO 15 K=1,N2 35270 IF(K-N) 46,46,16 35280 46 CONTINUE 35290 U=R(K,I)*A*C+R(K,J)*B 35300 R(K,J)=-R(K,I)*B*C+R(K,J)*A 35310 R(K,I)=U 35320 GO TO 15 35330 16 K1=K-N 35340 U=RD(K1,I)*A*C+RD(K1,J)*B 35350 RD(K1,J)=-RD(K1,I)*B*C+RD(K1,J)*A 35360 RD(K1,I)=U 35370 15 CONTINUE 35380 DO 13 K=1,N 35390 U=R(I,K)*A*C+R(J,K)*B 35400 R(J,K)=-R(I,K)*B*C+R(J,K)*A 35410 13 R(I,K)=U 35420 DO 30 I=1,N 35430 30 D(I)=ABS((D(I)-R(I,I))/D(I)) 35440 S=0.0 35450 DO 17 I=1,N 35460 17 S=AMAX1(S,D(I)) 35470 IF(S-P) 38,38,35 35480 35 L=L+1 35490 IF(L-NL) 14,38,38 35500 38 CONTINUE 35510 DO 40 I=1,N 35520 40 D(I)=R(I,I) 35530 DO 513 I=1,N 35540 IF(D(I).LT.0.0) D(I)=0.0 35550 513 CONTINUE 35560 RETURN 35570 END 35580 $ FORTY FORM 35590 C 35600 SUBROUTINE PLOTR(X, Y, XA, YA, XI, YI, NPOI, ID, LONG) 35610 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 35620 C 35630 C ROUTINE TO GENERATE A ONE PAGE PLOT OF ARRAY -X- VS. -Y-. 35640 C 35650 C THE PARAMETERS -XMAX- -XMIN- -YMAX- -YMIN- INDICATE THE UPPER 35660 C AND LOWER BOUNDS FOR EACH AXIS. IF XMAX=XMIN THE ROUTINE GENERATES 35670 C ITS OWN BOUNDS FOR THE X AXIS, AND SIMILARLY IF YMAX=YMIN. 35680 C 35690 C IT IS ASSUMED THAT THE ARRAYS HAVE -NPOI- ENTRIES. 35700 C 35710 C THE PLOTTING IS DONE ON TAPE -LOUT-. 35720 C 35730 C IF -ID- IS POSITIVE, AXES WILL BE INCLUDED ON THE PLOT, 35740 C IF -ID- IS NEGATIVE, NO AXES WILL APPEAR. 35750 C IF ID IS PLUS OR MINUS 2, THE POINTS WILL BE IDENTIFIED, 35760 C OTHERWISE THEY WILL BE COUNTED. 35770 C 35780 C -LONG- IS THE LENGTH OF ARRAYS -X- AND -Y- IN THE CALLING PROGRAM. 35790 C 35800 C IF VECTOR X EQUALS VECTOR Y THEN POINTS WILL BE PLOTTED ALONG 35810 C THE HORIZONTAL AXIS (NO AXES WILL BE PLOTTED) 35820 C 35830 C THE PLOT IS CLEARED DURING COMPILATION 35840 C IF IN THE EXECUTION OF THE CALLING PROGRAM THE PLOT IS 35850 C LEFT UNCLEAR IT MAY BE CLEARED BY ... CALL CLEAR 35860 C 35870 C WRITTEN BY FORREST W. YOUNG, NOVEMBER, 1965 35880 C VERSION 7 (SEPT 1971) 35890 C 35900 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 35910 INTEGER*2 ITEM,HOLL,PTID,AIE,AMINUS,DD,PLUS,BLANK,A,B 35920 DIMENSION X(LONG), Y(LONG), ITEM(55,101), SMALL(21), HOLL(11), 35930 * PTID(35) 35940 COMMON /OPT/ IN, LOUT, NPH, TIT(20), FM(20), ME, INT, NCON, MXIT, 35950 * CVCR, ITIE, INS, IDATA, INTPR, IDIST, IVEST, IPLTC, IPHC, IPHV, 35960 * NBDMX1, NDMX1, NB2, NNB2, NBSQ2, NTIE, NBS, NMIS, NBDMX, ICODE, 35970 * NCODE, NB1, FNB1, FNB, FNBI, FNB2, ILLT 35980 DATA HOLL/1H ,1HX,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9,1HM/, 35990 XAIE,AMINUS,DD,PLUS,BLANK/1HI,1H-,1H.,1H*,1H / 36000 DATA PTID/"1","2","3","4","5","6","7","8","9","A","B","C","D", 36010 1"E","F","G","H","I","J","K","L","M","N","O","P","Q","R","S", 36020 2"T","U","V","W","X","Y","Z"/ 36030 C 36040 XMAX = XA 36050 YMAX = YA 36060 XMIN = XI 36070 YMIN = YI 36080 C 36090 IF (ITEM(1,1).EQ.BLANK) GO TO 30 36100 C 36110 C INITIAL CLEAR OF PLOT 36120 C 36130 DO 20 I=1,55 36140 DO 10 J=1,101 36150 ITEM(I,J) = BLANK 36160 10 CONTINUE 36170 20 CONTINUE 36180 C 36190 C CHECK TO SEE IF ONE DIMENSIONAL PLOT 36200 C 36210 30 IF1DIM = 0 36220 DO 40 I=1,NPOI 36230 IF (X(I)-Y(I)) 50, 40, 50 36240 40 CONTINUE 36250 IF1DIM = 1 36260 GO TO 80 36270 C 36280 C IF DESIRED, PREPARE AXES TO BE PLOTTED 36290 C AXES PLOTTED ON CENTER OF PAGE, NOT AT ORIGIN OF PLOT 36300 C 36310 50 IF (ID.LE.0) GO TO 80 36320 DO 60 I=1,55 36330 ITEM(I,51) = AIE 36340 60 CONTINUE 36350 DO 70 I=1,101 36360 ITEM(28,I) = AMINUS 36370 70 CONTINUE 36380 80 CONTINUE 36390 C 36400 C DETERMINE MINIMUM AND MAXIMUM OF X AXIS, IF NECESSARY 36410 C 36420 IF (XMAX.NE.XMIN) GO TO 100 36430 XMAX = -1.0E7 36440 XMIN = +1.0E7 36450 DO 90 I=1,NPOI 36460 IF (X(I).GT.XMAX) XMAX = X(I) 36470 IF (X(I).LT.XMIN) XMIN = X(I) 36480 90 CONTINUE 36490 C 36500 C DETERMINE MINIMUM AND MAXIMUM OF Y AXIS, IF NECESSARY 36510 C 36520 100 IF (YMAX.NE.YMIN) GO TO 120 36530 YMAX = -1.0E7 36540 YMIN = +1.0E7 36550 DO 110 I=1,NPOI 36560 IF (Y(I).GT.YMAX) YMAX = Y(I) 36570 IF (Y(I).LT.YMIN) YMIN = Y(I) 36580 110 CONTINUE 36590 120 CONTINUE 36600 C 36610 C DETERMINE RANGE AND INCREMENT OF BOTH AXES 36620 C 36630 SPANX = XMAX - XMIN 36640 SPANY = YMAX - YMIN 36650 DELX = SPANX/100.0 36660 DELY = SPANY/54.0 36670 VALUE = YMAX + DELY 36680 C 36690 C PREPARE LABELING OF AXES 36700 C 36710 SMALL(1) = XMIN 36720 DO 130 I=1,20 36730 SMALL(I+1) = SMALL(I) + 5.0*DELX 36740 130 CONTINUE 36750 C 36760 C PREPARE PLOT 36770 C 36780 DO 170 II=1,NPOI 36790 I = (YMAX-Y(II))/DELY + 1.5 36800 IF (IF1DIM.EQ.1) I = 28 36810 J = (X(II)-XMIN)/DELX + 1.5 36820 IF (I.GT.55 .OR. I.LT.1 .OR. J.GT.101 .OR. J.LT.1) GO TO 170 36830 IF (IABS(ID).NE.2) GO TO 140 36840 KK = MOD(II-1,35) + 1 36850 ITEM(I,J) = PTID(KK) 36860 GO TO 170 36870 140 DO 150 JJ=1,10 36880 IF (ITEM(I,J).EQ.HOLL(JJ)) GO TO 160 36890 150 CONTINUE 36900 IF (ITEM(I,J).EQ.AIE .OR. ITEM(I,J).EQ.AMINUS) ITEM(I,J) = HOLL(2) 36910 GO TO 170 36920 160 ITEM(I,J) = HOLL(JJ+1) 36930 170 CONTINUE 36940 C 36950 C LABEL PLOT 36960 C 36970 WRITE (LOUT,9999) DD, (SMALL(I),DD,I=2,20,2) 36980 9999 FORMAT (15X, A1, 10(F9.4, A1)) 36990 WRITE (LOUT,9998) (SMALL(I),I=1,21,2) 37000 9998 FORMAT (10X, 11F10.4) 37010 WRITE (LOUT,9997) 37020 9997 FORMAT (14X, 49H*.****.****.****.****.****.****.****.****.****.**, 37030 * 54H**.****.****.****.****.****.****.****.****.****.****.*) 37040 C 37050 C PRINT PLOT 37060 C 37070 DO 210 I=1,55 37080 VALUE = VALUE - DELY 37090 A = BLANK 37100 B = PLUS 37110 L = I + 2 37120 IF (L/3*3-L) 200, 180, 200 37130 180 B = DD 37140 IF (L/2*2-L) 200, 190, 200 37150 190 A = DD 37160 200 WRITE (LOUT,9996) VALUE, A, B, (ITEM(I,J),J=1,101), B, A, VALUE 37170 210 CONTINUE 37180 9996 FORMAT (1H , F11.2, 1X, 105A1, F11.2) 37190 C 37200 C LABEL PLOT 37210 C 37220 WRITE (LOUT,9997) 37230 C 37240 C CLEAR AXES, IF NECESSARY 37250 C 37260 IF (ID.LE.0) GO TO 240 37270 DO 220 I=1,55 37280 ITEM(I,51) = BLANK 37290 220 CONTINUE 37300 DO 230 I=1,101 37310 ITEM(28,I) = BLANK 37320 230 CONTINUE 37330 C 37340 C CLEAR PLOT 37350 C 37360 240 DO 250 II=1,NPOI 37370 I = (YMAX-Y(II))/DELY + 1.5 37380 IF (IF1DIM.EQ.1) I = 28 37390 J = (X(II)-XMIN)/DELX + 1.5 37400 IF (I.GT.55 .OR. I.LT.1 .OR. J.GT.101 .OR. J.LT.1) GO TO 250 37410 ITEM(I,J) = BLANK 37420 250 CONTINUE 37430 RETURN 37440 END 37450