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<JK<IK), 1830
C F(IJ<IK<JK),F(JK<IJ<IK),F(JK<IK<IJ),F(IK<IJ<JK), 1840
C F(IK<JK<IJ). 1850
C WHEN ICODE=4, 1860
C 1=MATRIX INPUT; INPUT NS NBXNB DATA MATRICES; 1870
C (DEFAULT). 1880
C 2=FREQUENCY INPUT; INPUT F(R),N,M,I,I(1),...,I(NB-1), 1890
C WHERE, 1900
C F(R): FREQUENCY OF A CONDITIONAL RANK ORDER; 1910
C (DEFAULT=1.0). 1920
C N : NUMBER OF COMPARISON STIMULI; 1930
C (DEFAULT=NB-1). 1940
C M : NUMBER OF STIMULI ACTUALLY RANKED; 1950
C (DEFAULT=N). 1960
C I : INDEX OF A STANDARD STIMULUS. 1970
C I(J): INDEX OF A STIMULUS WHICH IS THE J'TH 1980
C MOST SIMILAR STIMULUS TO THE STANDARD. 1990
C IF NTIE=1, INPUT THE NUMBER OF TIES FOR EACH I(J) 2000
C ON THE FOLLOWING CARD IN FORMAT (20I4). 2010
C WHEN ICODE=5, 2020
C 1=MATRIX INPUT; MATRIX CONDITIONAL; INPUT LOWER TRI- 2030
C ANGULAR PORTIONS (W/O DIAGONALS) OF NS NBXNB DATA 2040
C MATRICES; (DEFAULT). 2050
C 2=MATRIX INPUT; ARBITRARY PARTITIONINGS WITHIN 2060
C MATRICES; INPUT LOWER TRIANGULAR PORTIONS (W/O 2070
C DIAG.) OF NS PAIRS OF NBXNB DATA MATRICES AND PARTI- 2080
C TIONING MATRICES; (DISSIMILARITIES BELONGING TO THE 2090
C SAME SET ARE INDICATED BY THE SAME NUMBER; INPUT 2100
C FORMAT OF PARTITIONING INFORMATION IS 40I2). 2110
C 3=FREQUENCY INPUT; INPUT TWO CARDS FOR EACH RANK 2120
C ORDER; 2130
C CARD 1: F(R),N,M, IN FORMAT(F10.0,2I5), WHERE 2140
C F(R): FREQUENCY OF A RANK ORDER; (DEFAULT=1.0). 2150
C N : NUMBER OF DISSIMILARITIES INVOLVED IN A 2160
C RANK ORDER; (DEFAULT=NB*(NB-1)/2). 2170
C M : NUMBER OF DISSIMILARITIES ACTUALLY RANKED; 2180
C (DEFAULT=N). 2190
C CARD 2: (I(1),J(1)),...,(I(N),J(N)), WHERE (I(K), 2200
C J(K)) IS THE INDEX OF THE K'TH SMALLEST DIS- 2210
C SIMILARITY. 2220
C IF NTIE=1, INPUT THE NUMBER OF TIES FOR EACH 2230
C (I(K),J(K)) ON THE FOLLOWING CARD IN FORMAT (20I4). 2240
C 2250
C ALL INFORMATION (EXECPT F(.)) SHOULD BE INPUT IN INTEGER 2260
C FORM; THE FORMAT SHOUL BE SPECIFIED IN FM EXCEPT FOR CARD 1 2270
C IN ICODE=5 & NCODE=3, WHICH IS FIXED TO (F10.0,2I5). 2280
C 2290
C ITIE: TIED RANKS (RELEVANT ONLY IN CASE THERE ARE TIED RANKS) 2300
C 0= WEAK TIES (DEFAULT) 2310
C 1= STRONG TIES 2320
C INS: DISPERSION PARAMETERS 2330
C 0= NO DIFFERENCE OVER DATA SETS (DEFAULT;NS=1) 2340
C 1= DIFFERENT DISPERSIONS OVER DATA SETS 2350
C NMIS: MISSING DATA OPTIONS (RELEVANT ONLY IN CASE OF 2360
C MATRIX INPUT) 2370
C 0= ALL COMPARISON STIMULI ARE PRESENT (DEFAULT) 2380
C 1= COMPARISON STIMULI ARE RESTRICTED 2390
C NTIE: INFORMATION REGARDING TIED RANKS (RELEVANT ONLY IN CASE 2400
C OF FREQUENCY INPUT; I.E., NCODE=3 IN ICODE=4 OR NCODE=2 2410
C IN ICODE=5) 2420
C 0= ASSUMED NO TIES (DEFAULT) 2430
C 1= READ FROM CARDS 2440
C NR: F(>)+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<F(>), NR<F(>)+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<KL, 5X, 5HIJ=KL) 15590
WRITE (LOUT,9993) 15600
9993 FORMAT (1H ) 15610
NTK = NCR(K) 15620
DO 200 I=1,NTK 15630
I1 = I4 + 1 15640
I4 = I4 + 4 15650
J1 = J3 + 1 15660
J3 = J3 + 3 15670
WRITE (LOUT,9992) I, (ID(J),J=I1,I4), (FQ(J),J=J1,J3) 15680
9992 FORMAT (10X, I3, 4I5, 3F10.0) 15690
200 CONTINUE 15700
210 CONTINUE 15710
C 15720
C CONVERT THE DATA FORMAT 15730
C 15740
220 N = 0 15750
I4 = 0 15760
J3 = 0 15770
DO 260 K=1,NS 15780
NCRK = 0 15790
NTK = NCR(K) 15800
DO 250 I=1,NTK 15810
I1 = I4 + 1 15820
I2 = I4 + 2 15830
I3 = I4 + 3 15840
I4 = I4 + 4 15850
J1 = J3 + 1 15860
IF (FQ(J1).EQ.0.0) GO TO 230 15870
NCRK = NCRK + 1 15880
N = N + 1 15890
F(N) = FQ(J1) 15900
NSC(N) = 2 15910
NSR(N) = 2 15920
NN = N + N 15930
IZ(1,NN) = ID(I1) 15940
IZ(2,NN) = ID(I2) 15950
IY(NN) = 1 15960
NN = NN - 1 15970
IZ(1,NN) = ID(I3) 15980
IZ(2,NN) = ID(I4) 15990
IY(NN) = 1 16000
230 J2 = J3 + 2 16010
IF (FQ(J2).EQ.0.0) GO TO 240 16020
NCRK = NCRK + 1 16030
N = N + 1 16040
F(N) = FQ(J2) 16050
NSC(N) = 2 16060
NSR(N) = 2 16070
NN = N + N 16080
IZ(1,NN) = ID(I3) 16090
IZ(2,NN) = ID(I4) 16100
IY(NN) = 1 16110
NN = NN - 1 16120
IZ(1,NN) = ID(I1) 16130
IZ(2,NN) = ID(I2) 16140
IY(NN) = 1 16150
240 J3 = J3 + 3 16160
IF (FQ(J3).EQ.0.0) GO TO 250 16170
NCRK = NCRK + 1 16180
N = N + 1 16190
F(N) = FQ(J3) 16200
NSC(N) = 2 16210
NSR(N) = 2 16220
NN = N + N 16230
IZ(1,NN) = ID(I3) 16240
IZ(2,NN) = ID(I4) 16250
IY(NN) = 2 16260
NN = NN - 1 16270
IZ(1,NN) = ID(I1) 16280
IZ(2,NN) = ID(I2) 16290
IY(NN) = 2 16300
250 CONTINUE 16310
NCR(K) = NCRK 16320
260 CONTINUE 16330
C 16340
IF (NTIE.EQ.1) GO TO 280 16350
AL = 0.0 16360
J3 = 0 16370
NN = 0 16380
FO = 0.0 16390
DO 270 K=1,NTS 16400
J1 = J3 + 1 16410
J2 = J3 + 2 16420
J3 = J3 + 3 16430
FQJ1 = FQ(J1) 16440
FQJ2 = FQ(J2) 16450
IF (FQJ1.EQ.0.0 .OR. FQJ2.EQ.0.0) GO TO 270 16460
NN = NN + 1 16470
FN = FQJ1 + FQJ2 16480
FO = FO + FN 16490
AL = AL + FQJ1*ALOG(FQJ1/FN) + FQJ2*ALOG(FQJ2/FN) 16500
270 CONTINUE 16510
AL = -2.0*AL 16520
NN = NTS - NN 16530
WRITE (LOUT,9991) AL, NTS, NN 16540
9991 FORMAT (///43H THE -2*LOG-LIKELIHOOD OF THE NULL MODEL IS, F12.4, 16550
* 3X, 4HWITH, I5, 2X, 20HPARAMETERS (OF WHICH, I4, 2X, 9HARE P=0 O, 16560
* 4HR 1)) 16570
AIC = AL + 2.0*NTS 16580
WRITE (LOUT,9990) AIC 16590
9990 FORMAT (/52H THE VALUE OF THE AIC STATISTIC OF THE NULL MODEL IS, 16600
* F12.4) 16610
BIC = AL + NTS*ALOG(FO) 16620
WRITE (LOUT,9989) BIC 16630
9989 FORMAT (/52H THE VALUE OF THE BIC STATISTIC OF THE NULL MODEL IS, 16640
* F12.4) 16650
C 16660
280 IF (INS.EQ.1) RETURN 16670
NCRK = 0 16680
DO 290 K=1,NS 16690
NCRK = NCRK + NCR(K) 16700
290 CONTINUE 16710
NCR(1) = NCRK 16720
NS = 1 16730
RETURN 16740
END 16750
$ FORTY FORM 16760
C 16770
SUBROUTINE TRIADC(ID, FQ, NCR, NSC, NSR, IZ, IY, F, WK, NB, NS, 16780
* NTS) 16790
C 16800
C DATA INPUT FOR TRIADIC COMBINATIONS 16810
C 16820
DIMENSION ID(1), FQ(1), NCR(1), NSC(1), NSR(1), IZ(2,1), IY(1), 16830
* F(1), WK(1) 16840
COMMON /OPT/ IN, LOUT, NPH, TIT(20), FM(20), ME, INT, NCON, MXIT, 16850
* CVCR, ITIE, INS, IDATA, INTPR, IDIST, IVEST, IPLTC, IPHC, IPHV, 16860
* NBDMX1, NDMX1, NB2, NNB2, NBSQ2, NTIE, NBS, NMIS, NBDMX, ICODE, 16870
* NCODE, NB1, FNB1, FNB, FNBI, FNB2, ILLT 16880
C 16890
IF (IDATA.EQ.0) GO TO 10 16900
WRITE (LOUT,9999) 16910
9999 FORMAT (1H1, //34H INPUT DATA (TRIADIC COMBINATIONS)) 16920
10 N = 0 16930
FO = 0.0 16940
AL = 0.0 16950
J3 = 0 16960
I4 = 0 16970
DO 210 K=1,NS 16980
IF (IDATA.EQ.0) GO TO 30 16990
IF (NS.EQ.1) GO TO 20 17000
WRITE (LOUT,9998) K 17010
9998 FORMAT (//6X, 9HDATA SET=, I3) 17020
20 WRITE (LOUT,9997) 17030
9997 FORMAT (/11X, 3HNO., 3X, 11HI J K, 4X, 17HIJ<JK<IK IJ<IK<J, 17040
* 3HK , 38HJK<IJ<IK JK<IK<IJ IK<IJ<JK IK<JK<IJ) 17050
WRITE (LOUT,9996) 17060
9996 FORMAT (1H ) 17070
30 CONTINUE 17080
NCRK = 0 17090
NTK = NCR(K) 17100
DO 200 I=1,NTK 17110
C INPUT 17120
READ (IN,FM) K1, K2, K3, (WK(J),J=1,6) 17130
C DATA CHECK 17140
IF (K1.LT.1 .OR. K1.GT.NB) GO TO 50 17150
IF (K2.LT.1 .OR. K2.GT.NB) GO TO 50 17160
IF (K3.LT.1 .OR. K3.GT.NB) GO TO 50 17170
DO 40 J=1,6 17180
IF (WK(J).LT.0.0) GO TO 60 17190
40 CONTINUE 17200
IF (K1.EQ.K2) GO TO 70 17210
IF (K1.EQ.K3) GO TO 70 17220
IF (K2.EQ.K3) GO TO 70 17230
GO TO 100 17240
50 WRITE (LOUT,9995) 17250
9995 FORMAT (///34H STIMULUS INDICES ARE OUT OF BOUND/13H EITHER NONPO, 17260
* 41HSITIVE INDICES OR INDICES GREATER THAN NB) 17270
GO TO 80 17280
60 WRITE (LOUT,9994) 17290
9994 FORMAT (///21H NEGATIVE FREQUENCIES) 17300
GO TO 80 17310
70 WRITE (LOUT,9993) 17320
9993 FORMAT (///47H NONPERMISSIBLE COMBINATIONS OF STIMULS INDICES/ 17330
* 23H EITHER I=J, I=K OR J=K) 17340
80 IF (NS.EQ.1) GO TO 90 17350
WRITE (LOUT,9998) K 17360
90 WRITE (LOUT,9997) 17370
WRITE (LOUT,9992) I, K1, K2, K3, (WK(J),J=1,6) 17380
STOP 17390
C OUTPUT 17400
100 IF (IDATA.EQ.1) WRITE (LOUT,9992) I, K1, K2, K3, (WK(J),J=1,6) 17410
9992 FORMAT (10X, I3, 3I5, 6F10.0) 17420
FN = 0.0 17430
DO 110 J=1,6 17440
FN = FN + WK(J) 17450
110 CONTINUE 17460
IF (FN.EQ.0.0) GO TO 130 17470
DO 120 J=1,6 17480
IF (WK(J).EQ.0.0) GO TO 120 17490
AL = AL + WK(J)*ALOG(WK(J)/FN) 17500
120 CONTINUE 17510
FO = FO + FN 17520
130 CONTINUE 17530
C CONVERT THE DATA FORMAT 17540
IF (WK(1).EQ.0.0) GO TO 140 17550
NCRK = NCRK + 1 17560
N = N + 1 17570
F(N) = WK(1) 17580
NSC(N) = 3 17590
NSR(N) = 3 17600
NN = N*3 17610
IZ(1,NN) = K1 17620
IZ(2,NN) = K3 17630
IY(NN) = 1 17640
NN = NN - 1 17650
IZ(1,NN) = K2 17660
IZ(2,NN) = K3 17670
IY(NN) = 1 17680
NN = NN - 1 17690
IZ(1,NN) = K1 17700
IZ(2,NN) = K2 17710
IY(NN) = 1 17720
140 IF (WK(2).EQ.0.0) GO TO 150 17730
NCRK = NCRK + 1 17740
N = N + 1 17750
F(N) = WK(2) 17760
NSC(N) = 3 17770
NSR(N) = 3 17780
NN = N*3 17790
IZ(1,NN) = K2 17800
IZ(2,NN) = K3 17810
IY(NN) = 1 17820
NN = NN - 1 17830
IZ(1,NN) = K1 17840
IZ(2,NN) = K3 17850
IY(NN) = 1 17860
NN = NN - 1 17870
IZ(1,NN) = K1 17880
IZ(2,NN) = K2 17890
IY(NN) = 1 17900
150 IF (WK(3).EQ.0.0) GO TO 160 17910
NCRK = NCRK + 1 17920
N = N + 1 17930
F(N) = WK(3) 17940
NSC(N) = 3 17950
NSR(N) = 3 17960
NN = N*3 17970
IZ(1,NN) = K1 17980
IZ(2,NN) = K3 17990
IY(NN) = 1 18000
NN = NN - 1 18010
IZ(1,NN) = K1 18020
IZ(2,NN) = K2 18030
IY(NN) = 1 18040
NN = NN - 1 18050
IZ(1,NN) = K2 18060
IZ(2,NN) = K3 18070
IY(NN) = 1 18080
160 IF (WK(4).EQ.0.0) GO TO 170 18090
NCRK = NCRK + 1 18100
N = N + 1 18110
F(N) = WK(4) 18120
NSC(N) = 3 18130
NSR(N) = 3 18140
NN = N*3 18150
IZ(1,NN) = K1 18160
IZ(2,NN) = K2 18170
IY(NN) = 1 18180
NN = NN - 1 18190
IZ(1,NN) = K1 18200
IZ(2,NN) = K3 18210
IY(NN) = 1 18220
NN = NN - 1 18230
IZ(1,NN) = K2 18240
IZ(2,NN) = K3 18250
IY(NN) = 1 18260
170 IF (WK(5).EQ.0.0) GO TO 180 18270
NCRK = NCRK + 1 18280
N = N + 1 18290
F(N) = WK(5) 18300
NSC(N) = 3 18310
NSR(N) = 3 18320
NN = N*3 18330
IZ(1,NN) = K2 18340
IZ(2,NN) = K3 18350
IY(NN) = 1 18360
NN = NN - 1 18370
IZ(1,NN) = K1 18380
IZ(2,NN) = K2 18390
IY(NN) = 1 18400
NN = NN - 1 18410
IZ(1,NN) = K1 18420
IZ(2,NN) = K3 18430
IY(NN) = 1 18440
180 IF (WK(6).EQ.0.0) GO TO 190 18450
NCRK = NCRK + 1 18460
N = N + 1 18470
F(N) = WK(6) 18480
NSC(N) = 3 18490
NSR(N) = 3 18500
NN = N*3 18510
IZ(1,NN) = K1 18520
IZ(2,NN) = K2 18530
IY(NN) = 1 18540
NN = NN - 1 18550
IZ(1,NN) = K2 18560
IZ(2,NN) = K3 18570
IY(NN) = 1 18580
NN = NN - 1 18590
IZ(1,NN) = K1 18600
IZ(2,NN) = K3 18610
IY(NN) = 1 18620
C DATA FOR INITIALIZATION 18630
190 IF (INT.EQ.1) GO TO 200 18640
J1 = J3 + 1 18650
J2 = J3 + 2 18660
J3 = J3 + 3 18670
FQ(J1) = WK(3) + WK(4) + WK(6) 18680
FQ(J2) = WK(1) + WK(2) + WK(5) 18690
FQ(J3) = 0.0 18700
I1 = I4 + 1 18710
I2 = I4 + 2 18720
I3 = I4 + 3 18730
I4 = I4 + 4 18740
ID(I1) = K1 18750
ID(I2) = K2 18760
ID(I3) = K2 18770
ID(I4) = K3 18780
J1 = J3 + 1 18790
J2 = J3 + 2 18800
J3 = J3 + 3 18810
FQ(J1) = WK(2) + WK(5) + WK(6) 18820
FQ(J2) = WK(1) + WK(3) + WK(4) 18830
FQ(J3) = 0.0 18840
I1 = I4 + 1 18850
I2 = I4 + 2 18860
I3 = I4 + 3 18870
I4 = I4 + 4 18880
ID(I1) = K2 18890
ID(I2) = K3 18900
ID(I3) = K1 18910
ID(I4) = K3 18920
J1 = J3 + 1 18930
J2 = J3 + 2 18940
J3 = J3 + 3 18950
FQ(J1) = WK(4) + WK(5) + WK(6) 18960
FQ(J2) = WK(1) + WK(2) + WK(3) 18970
FQ(J3) = 0.0 18980
I1 = I4 + 1 18990
I2 = I4 + 2 19000
I3 = I4 + 3 19010
I4 = I4 + 4 19020
ID(I1) = K1 19030
ID(I2) = K2 19040
ID(I3) = K1 19050
ID(I4) = K3 19060
200 CONTINUE 19070
NCR(K) = NCRK 19080
210 CONTINUE 19090
AL = -2.0*AL 19100
NTS5 = NTS*5 19110
NN = NTS*6 - N 19120
WRITE (LOUT,9991) AL, NTS5, NN 19130
9991 FORMAT (///43H THE -2*LOG-LIKELIHOOD OF THE NULL MODEL IS, F12.4, 19140
* 3X, 4HWITH, I5, 2X, 20HPARAMETERS (OF WHICH, I4, 2X, 9HARE P=0 O, 19150
* 4HR 1)) 19160
AIC = AL + 2.0*NTS5 19170
WRITE (LOUT,9990) AIC 19180
9990 FORMAT (/52H THE VALUE OF THE AIC STATISTIC OF THE NULL MODEL IS, 19190
* F12.4) 19200
BIC = AL + NTS5*ALOG(FO) 19210
WRITE (LOUT,9989) BIC 19220
9989 FORMAT (/52H THE VALUE OF THE BIC STATISTIC OF THE NULL MODEL IS, 19230
* F12.4) 19240
IF (INT.EQ.0) NTS = I4/4 19250
IF (INS.EQ.1) RETURN 19260
NCRK = 0 19270
DO 220 K=1,NS 19280
NCRK = NCRK + NCR(K) 19290
220 CONTINUE 19300
NCR(1) = NCRK 19310
NS = 1 19320
RETURN 19330
END 19340
$ FORTY FORM 19350
C 19360
SUBROUTINE CRANK(O, NCR, NSC, NSR, DT, DU, IX1, IX2, IY, IZ, P, 19370
* F, NS, NB) 19380
C 19390
C DATA INPUT FOR CONDITIONAL RANK ORDERS 19400
C 19410
INTEGER O(1) 19420
DIMENSION NCR(1), NSC(1), NSR(1), DT(1), DU(1), IX1(1), IX2(1), 19430
* IY(1), IZ(2,1), P(NB,1), F(1) 19440
COMMON /OPT/ IN, LOUT, NPH, TIT(20), FM(20), ME, INT, NCON, MXIT, 19450
* CVCR, ITIE, INS, IDATA, INTPR, IDIST, IVEST, IPLTC, IPHC, IPHV, 19460
* NBDMX1, NDMX1, NB2, NNB2, NBSQ2, NTIE, NBS, NMIS, NBDMX, ICODE, 19470
* NCODE, NB1, FNB1, FNB, FNBI, FNB2, ILLT 19480
C 19490
IF (INT.EQ.1) GO TO 30 19500
DO 20 I=1,NB 19510
DO 10 J=1,NB 19520
P(I,J) = 0.0 19530
10 CONTINUE 19540
20 CONTINUE 19550
C 19560
30 IF (NCODE.EQ.2) GO TO 240 19570
C 19580
C MATRIX INPUT 19590
IF (IDATA.EQ.0) GO TO 40 19600
WRITE (LOUT,9999) 19610
9999 FORMAT (1H1//49H INPUT DATA (CONDITIONAL RANK ORDERS): MATRIX INP, 19620
* 2HUT) 19630
40 N = 0 19640
I0 = 0 19650
DO 210 K=1,NS 19660
C INPUT 19670
IT = 0 19680
DO 50 I=1,NB 19690
II = IT + 1 19700
IT = IT + NB 19710
READ (IN,FM) (O(J),J=II,IT) 19720
50 CONTINUE 19730
C OUTPUT 19740
IF (IDATA.EQ.0) GO TO 110 19750
IF (NS.GT.1) WRITE (LOUT,9998) K 19760
9998 FORMAT (//6X, 8HSUBJECT=, I3) 19770
N1 = NB/30 19780
NRS = NB - N1*30 19790
IF (NRS.NE.0) GO TO 60 19800
NRS = 30 19810
N1 = N1 - 1 19820
60 N1 = N1 + 1 19830
I2 = 0 19840
DO 100 I=1,N1 19850
IF (I.LT.N1) GO TO 70 19860
J = NRS 19870
GO TO 80 19880
70 J = 30 19890
80 I1 = I2 + 1 19900
I2 = I2 + J 19910
WRITE(LOUT,9997) (KK,KK=I1,I2) 19920
9997 FORMAT (/11X, 30I4) 19930
WRITE (LOUT,9996) 19940
9996 FORMAT (1H ) 19950
II = (I-1)*30 19960
DO 90 L=1,NB 19970
I1 = II + 1 19980
IT = II + J 19990
WRITE (LOUT,9995) L, (O(LL),LL=I1,IT) 20000
9995 FORMAT (6X, I4, 1X, 30I4) 20010
II = II + NB 20020
90 CONTINUE 20030
100 CONTINUE 20040
110 CONTINUE 20050
C CONVERT THE DATA FORMAT 20060
NCR(K) = NB 20070
NN = 0 20080
DO 170 I=1,NB 20090
N = N + 1 20100
NSC(N) = 0 20110
DO 140 J=1,NB 20120
IX1(J) = J 20130
NN = NN + 1 20140
IF (J.EQ.I) GO TO 120 20150
IF (O(NN).GT.0) GO TO 130 20160
DT(J) = 1.0E10 20170
DU(J) = DT(J) 20180
NSC(N) = NSC(N) + 1 20190
GO TO 140 20200
120 DT(J) = 1.0E20 20210
O(NN) = 0 20220
DU(J) = DT(J) 20230
GO TO 140 20240
130 DT(J) = O(NN) 20250
DU(J) = DT(J) 20260
140 CONTINUE 20270
NSR(N) = NB1 - NSC(N) 20280
NSC(N) = NB1 20290
IF (NMIS.EQ.1) NSC(N) = NSR(N) 20300
C 20310
CALL SHEL9(DT, IX1, NB) 20320
CALL BLOC1(DU, IX1, IX2, NB) 20330
C 20340
F(N) = 1.0 20350
JJ = NSC(N) 20360
DO 150 J=1,JJ 20370
I0 = I0 + 1 20380
IZ(1,I0) = I 20390
IZ(2,I0) = IX1(J) 20400
IY(I0) = IX2(J) 20410
150 CONTINUE 20420
JJ = NSR(N) 20430
DO 160 J=1,JJ 20440
IF (IX2(J).GT.1) NTIE = 1 20450
160 CONTINUE 20460
170 CONTINUE 20470
C INITIAL DISTANCES 20480
NN = 0 20490
N = N - NB 20500
DO 200 I=1,NB 20510
N = N + 1 20520
C = FNB1/NSC(N) 20530
CX = C*0.5 20540
CO = (NSC(N)+NSR(N)-1)*CX 20550
IF (NMIS.EQ.1) CO = (NSR(N)+1)*CX 20560
DO 190 J=1,NB 20570
NN = NN + 1 20580
IF (I.EQ.J) GO TO 190 20590
IF (O(NN).EQ.0) GO TO 180 20600
P(I,J) = P(I,J) + O(NN)*C 20610
GO TO 190 20620
180 P(I,J) = P(I,J) + CO 20630
190 CONTINUE 20640
200 CONTINUE 20650
210 CONTINUE 20660
C SYMMETRIZE 20670
DO 230 I=1,NB 20680
P(I,I) = 0.0 20690
IF (I.EQ.1) GO TO 230 20700
I1 = I - 1 20710
DO 220 J=1,I1 20720
P(I,J) = P(I,J) + P(J,I) 20730
P(J,I) = P(I,J) 20740
220 CONTINUE 20750
230 CONTINUE 20760
GO TO 490 20770
C 20780
C FREQUENCY INPUT 20790
240 IF (IDATA.EQ.0) GO TO 250 20800
WRITE (LOUT,9994) 20810
9994 FORMAT (1H1, //47H INPUT DATA (CONDITIONAL RANK ORDERS): FREQUENC, 20820
* 7HY INPUT) 20830
250 N = 0 20840
NY = 0 20850
DO 430 K=1,NS 20860
IF (IDATA.EQ.0) GO TO 260 20870
IF (NS.GT.1) WRITE (LOUT,9993) K 20880
9993 FORMAT (//6X, 9HDATA SET=, I3) 20890
NM = 20 20900
IF (NB1.LT.NM) NM = NB1 20910
WRITE(LOUT,9992) (I,I=1,NM) 20920
9992 FORMAT (/6X, 27HNO. FREQ. N M I, 4X, 20I4) 20930
WRITE (LOUT,9996) 20940
260 NCRK = NCR(K) 20950
DO 420 JJ=1,NCRK 20960
N = N + 1 20970
C INPUT 20980
READ (IN,FM) F(N), NN, NSR(N), I, (O(J),J=1,NB1) 20990
IF (F(N).EQ.0.0) F(N) = 1.0 21000
IF (NN.EQ.0) NN = NB1 21010
IF (NSR(N).EQ.0) NSR(N) = NN 21020
NSC(N) = NN 21030
C DATA CHECK 21040
IF (F(N).LT.0.0) GO TO 290 21050
IF (NN.LT.2 .OR. NN.GT.NB1) GO TO 290 21060
IF (NSR(N).LT.1 .OR. NSR(N).GT.NN) GO TO 290 21070
IF (I.LT.1 .OR. I.GT.NB) GO TO 290 21080
DO 280 J=1,NN 21090
IF (O(J).LT.1 .OR. O(J).GT.NB) GO TO 290 21100
IF (I.EQ.O(J)) GO TO 300 21110
DO 270 J1=J,NN 21120
IF (J1.EQ.J) GO TO 270 21130
IF (O(J).EQ.O(J1)) GO TO 300 21140
270 CONTINUE 21150
280 CONTINUE 21160
GO TO 320 21170
290 WRITE (LOUT,9991) 21180
9991 FORMAT (///34H THE DATA OUT OF PERMISSIBLE RANGE) 21190
GO TO 310 21200
300 WRITE (LOUT,9990) 21210
9990 FORMAT (///48H NONPERMISSIBLE COMBINATIONS OF STIMULUS INDICES/ 21220
* 56H EITHER I=I(J) FOR SOME J, OR I(J)=I(K) FOR SOME J AND K) 21230
310 IF (NS.GT.1) WRITE (LOUT,9993) K 21240
WRITE (LOUT,9992) 21250
WRITE (LOUT,9989) JJ, F(N), NN, NSR(N), I, (O(J),J=1,NN) 21260
STOP 21270
C OUTPUT 21280
320 IF (IDATA.EQ.1) WRITE (LOUT,9989) JJ, F(N), NN, NSR(N), I, 21290
* (O(J),J=1,NN) 21300
9989 FORMAT (4X, I4, F10.0, 3I5, 4X, 20I4/(37X, 20I4)) 21310
C CONVERT THE DATA FORMAT 21320
DO 330 J=1,NN 21330
NY = NY + 1 21340
IZ(1,NY) = I 21350
IZ(2,NY) = O(J) 21360
IY(NY) = 1 21370
330 CONTINUE 21380
IF (NTIE.EQ.0) GO TO 370 21390
J1 = NY - NN + 1 21400
READ (IN,9988) (IY(J),J=J1,NY) 21410
9988 FORMAT (20I4) 21420
J1 = J1 - 1 21430
340 J1 = J1 + 1 21440
IF (J1.GT.NY) GO TO 370 21450
IF (IY(J1).LE.0) IY(J1) = 1 21460
IF (IY(J1).EQ.1) GO TO 340 21470
J2 = IY(J1) 21480
J11 = J1 + 1 21490
J22 = J1 + J2 - 1 21500
IF (J22.LE.NY) GO TO 350 21510
J22 = NY 21520
J2 = NY - J1 + 1 21530
350 DO 360 J=J11,J22 21540
IY(J2) = J2 21550
360 CONTINUE 21560
J1 = J22 21570
GO TO 340 21580
C INITIAL DISTANCES 21590
370 IF (INT.EQ.1) GO TO 420 21600
C = FNB1/NSC(N) 21610
CX = C*0.5 21620
CO = (NSC(N)+NSR(N)-1)*CX 21630
IF (NMIS.EQ.1) CO = (NSR(N)+1)*CX 21640
DO 410 J=1,NN 21650
I1 = I 21660
I2 = O(J) 21670
IF (I1.GT.I2) GO TO 380 21680
I1 = I2 21690
I2 = I 21700
380 IF (J.GT.NSR(N)) GO TO 390 21710
P(I1,I2) = P(I1,I2) + J*C*F(N) 21720
GO TO 400 21730
390 P(I1,I2) = P(I1,I2) + CO*F(N) 21740
400 P(I2,I1) = P(I2,I1) + F(N) 21750
410 CONTINUE 21760
420 CONTINUE 21770
430 CONTINUE 21780
C SYMMETRIZE 21790
IF (INT.EQ.1) GO TO 490 21800
N = 0 21810
FN = 0.0 21820
DO 450 I=2,NB 21830
I1 = I - 1 21840
DO 440 J=1,I1 21850
IF (P(J,I).EQ.0.0) GO TO 440 21860
P(I,J) = P(I,J)/P(J,I) 21870
N = N + 1 21880
FN = FN + P(I,J) 21890
440 CONTINUE 21900
450 CONTINUE 21910
FN = FN/N 21920
DO 480 I=2,NB 21930
I1 = I - 1 21940
DO 470 J=1,I1 21950
IF (P(J,I).NE.0.0) GO TO 460 21960
P(I,J) = FN 21970
460 P(J,I) = P(I,J) 21980
470 CONTINUE 21990
480 CONTINUE 22000
490 IF (INS.EQ.1) RETURN 22010
NCRK = 0 22020
DO 500 K=1,NS 22030
NCRK = NCRK + NCR(K) 22040
500 CONTINUE 22050
NCR(1) = NCRK 22060
NS = 1 22070
RETURN 22080
END 22090
$ FORTY FORM 22100
C 22110
SUBROUTINE GRANK(O, IX1, IX2, NCR, NSC, NSR, P, IY, IZ, WK1, WK2, 22120
* ID, F, NB, NS) 22130
C 22140
C DATA INPUT FOR GENERAL DIRECTIONAL RANK ORDERS 22150
C 22160
INTEGER O(1) 22170
DIMENSION IX1(1), IX2(1), NCR(1), NSC(1), NSR(1), IY(1), IZ(2,1), 22180
* P(NB,1), WK1(1), WK2(1), ID(2,1), F(1) 22190
COMMON /OPT/ IN, LOUT, NPH, TIT(20), FM(20), ME, INT, NCON, MXIT, 22200
* CVCR, ITIE, INS, IDATA, INTPR, IDIST, IVEST, IPLTC, IPHC, IPHV, 22210
* NBDMX1, NDMX1, NB2, NNB2, NBSQ2, NTIE, NBS, NMIS, NBDMX, ICODE, 22220
* NCODE, NB1, FNB1, FNB, FNBI, FNB2, ILLT 22230
C 22240
FNBSQ2 = NBSQ2 22250
IF (INT.EQ.1) GO TO 30 22260
DO 20 I=1,NB 22270
DO 10 J=1,NB 22280
P(I,J) = 0.0 22290
10 CONTINUE 22300
20 CONTINUE 22310
30 N = 0 22320
DO 50 I=2,NB 22330
I1 = I - 1 22340
DO 40 J=1,I1 22350
N = N + 1 22360
ID(1,N) = I 22370
ID(2,N) = J 22380
40 CONTINUE 22390
50 CONTINUE 22400
C 22410
IF (NCODE.EQ.3) GO TO 480 22420
C 22430
C MATRIX INPUT 22440
IF (IDATA.EQ.0) GO TO 60 22450
WRITE (LOUT,9999) 22460
9999 FORMAT (1H1, //46H INPUT DATA (GENERAL DIRECTIONAL RANK ORDERS):, 22470
* 13H MATRIX INPUT) 22480
IF (NCODE.EQ.2) WRITE (LOUT,9998) 22490
9998 FORMAT (/10X, 39HPARTITIONING INFORMATION IN PARENTHESES) 22500
60 CONTINUE 22510
N = 0 22520
NY = 0 22530
DO 390 K=1,NS 22540
C INPUT 22550
I2 = 0 22560
DO 70 I=2,NB 22570
I1 = I2 + 1 22580
I2 = I2 + I - 1 22590
READ (IN,FM) (O(J),J=I1,I2) 22600
70 CONTINUE 22610
IF (NCODE.EQ.1) GO TO 90 22620
I2 = 0 22630
DO 80 I=2,NB 22640
I1 = I2 + 1 22650
I2 = I2 + I - 1 22660
READ (IN,9997) (IX1(J),J=I1,I2) 22670
80 CONTINUE 22680
9997 FORMAT (40I2) 22690
90 CONTINUE 22700
C OUTPUT 22710
IF (IDATA.EQ.0) GO TO 150 22720
IF (NS.GT.1) WRITE (LOUT,9996) K 22730
9996 FORMAT (//6X, 8HSUBJECT=, I3) 22740
N1 = NB1/10 22750
NRS = NB1 - N1*10 22760
IF (NRS.GT.0) GO TO 100 22770
N1 = N1 - 1 22780
100 N1 = N1 + 1 22790
II = 1 22800
DO 140 L=1,N1 22810
I2 = II + 9 22820
IF (I2.GT.NB1) I2 = NB1 22830
WRITE(LOUT,9995) (J,J=II,I2) 22840
9995 FORMAT (//9X, 10I12) 22850
WRITE (LOUT,9994) 22860
9994 FORMAT (1H ) 22870
I0 = (II+1)*II/2 - 1 22880
II1 = II + 1 22890
DO 130 I=II1,NB 22900
I1 = I0 + 1 22910
III = I - II1 22920
IF (III.GT.9) III = 9 22930
I2 = I1 + III 22940
IF (NCODE.EQ.2) GO TO 110 22950
WRITE (LOUT,9993) I, (O(J),J=I1,I2) 22960
9993 FORMAT (1H , 2I10, 9I12) 22970
GO TO 120 22980
110 WRITE (LOUT,9992) I, (O(J),IX1(J),J=I1,I2) 22990
9992 FORMAT (1H , I10, 10(I6, 2X, 1H(, I2, 1H))) 23000
120 I0 = I0 + I - 1 23010
130 CONTINUE 23020
II = II + 10 23030
140 CONTINUE 23040
150 CONTINUE 23050
C 23060
C MATRIX CONDITIONAL 23070
C CONVERT THE DATA FORMAT 23080
IF (NCODE.EQ.2) GO TO 230 23090
NCR(K) = 1 23100
N = N + 1 23110
NSR(N) = 0 23120
DO 170 I=1,NBSQ2 23130
IX2(I) = I 23140
IF (O(I).EQ.0) GO TO 160 23150
WK1(I) = O(I) 23160
WK2(I) = WK1(I) 23170
GO TO 170 23180
160 WK1(I) = 1.0E10 23190
WK2(I) = WK1(I) 23200
NSR(N) = NSR(N) + 1 23210
170 CONTINUE 23220
C 23230
NY1 = NY + 1 23240
CALL SHEL9(WK1, IX2, NBSQ2) 23250
CALL BLOC1(WK2, IX2, IY(NY1), NBSQ2) 23260
C 23270
F(N) = 1.0 23280
NSR(N) = NBSQ2 - NSR(N) 23290
NSC(N) = NBSQ2 23300
IF (NMIS.EQ.1) NSC(N) = NSR(N) 23310
NK = NSR(N) 23320
DO 180 I=1,NK 23330
NY = NY + 1 23340
IF (IY(NY).GT.1) NTIE = 1 23350
180 CONTINUE 23360
NY = NY - NK 23370
NK = NSC(N) 23380
DO 190 I=1,NK 23390
NY = NY + 1 23400
II = IX2(I) 23410
IZ(1,NY) = ID(1,II) 23420
IZ(2,NY) = ID(2,II) 23430
190 CONTINUE 23440
C INITIAL DISTANCES 23450
C = FNBSQ2/NK 23460
CX = C*0.5 23470
CO = (NK+NSR(N)-1)*CX 23480
IF (NMIS.EQ.1) CO = (NSR(N)+1)*CX 23490
NN = 0 23500
DO 220 I=2,NB 23510
I1 = I - 1 23520
DO 210 J=1,I1 23530
NN = NN + 1 23540
IF (O(NN).EQ.0) GO TO 200 23550
P(I,J) = P(I,J) + O(NN)*C 23560
GO TO 210 23570
200 P(I,J) = P(I,J) + CO 23580
210 CONTINUE 23590
220 CONTINUE 23600
GO TO 390 23610
C 23620
C GENERAL PARTITIONINGS 23630
C CONVERT THE DATA FORMAT 23640
230 CONTINUE 23650
NW = 0 23660
NN = 0 23670
DO 250 I1=2,NB 23680
I11 = I1 - 1 23690
DO 240 J=1,I11 23700
NN = NN + 1 23710
IF (IX1(NN).EQ.0) IX1(NN) = 10000 23720
240 CONTINUE 23730
250 CONTINUE 23740
260 CONTINUE 23750
MN = 10000 23760
NN = 0 23770
DO 280 I1=2,NB 23780
I11 = I1 - 1 23790
DO 270 I2=1,I11 23800
NN = NN + 1 23810
IF (IX1(NN).LT.MN) MN = IX1(NN) 23820
270 CONTINUE 23830
280 CONTINUE 23840
IF (MN.EQ.10000) GO TO 380 23850
NN = 0 23860
NSCN = 0 23870
NSRN = 0 23880
DO 300 I1=2,NB 23890
I11 = I1 - 1 23900
DO 290 I2=1,I11 23910
NN = NN + 1 23920
IF (IX1(NN).NE.MN) GO TO 290 23930
NSCN = NSCN + 1 23940
WK1(NSCN) = O(NN) 23950
ID(1,NSCN) = I1 23960
ID(2,NSCN) = I2 23970
IF (O(NN).NE.0) NSRN = NSRN + 1 23980
290 CONTINUE 23990
300 CONTINUE 24000
IF (NSCN.LT.2) GO TO 260 24010
IF (NSRN.LT.1) GO TO 260 24020
NW = NW + 1 24030
N = N + 1 24040
NSC(N) = NSCN 24050
NSR(N) = NSRN 24060
F(N) = 1.0 24070
DO 310 I=1,NSCN 24080
WK2(I) = WK1(I) 24090
IX2(I) = I 24100
310 CONTINUE 24110
C 24120
NY1 = NY + 1 24130
CALL SHEL9(WK1, IX2, NSCN) 24140
CALL BLOC1(WK2, IX2, IY(NY1), NSCN) 24150
C 24160
DO 320 I=1,NSRN 24170
NY = NY + 1 24180
IF (IY(NY).GT.1) NTIE = 1 24190
320 CONTINUE 24200
NY = NY - NSRN 24210
DO 330 I=1,NSCN 24220
NY = NY + 1 24230
II = IX2(I) 24240
IZ(1,NY) = ID(1,II) 24250
IZ(2,NY) = ID(2,II) 24260
330 CONTINUE 24270
C INITIAL DISTANCES 24280
C = FNBSQ2/NSCN 24290
CX = C*0.5 24300
CO = (NSCN+NSRN-1)*CX 24310
IF (NMIS.EQ.1) CO = (NSRN+1)*CX 24320
NN = 0 24330
DO 370 I1=2,NB 24340
I11 = I1 - 1 24350
DO 360 J1=1,I11 24360
NN = NN + 1 24370
IF (IX1(NN).NE.MN) GO TO 360 24380
IX1(NN) = 10000 24390
IF (O(NN).EQ.0) GO TO 340 24400
P(I1,J1) = P(I1,J1) + O(NN)*C 24410
GO TO 350 24420
340 P(I1,J1) = P(I1,J1) + CO 24430
350 P(J1,I1) = P(J1,I1) + 1.0 24440
360 CONTINUE 24450
370 CONTINUE 24460
GO TO 260 24470
380 NCR(K) = NW 24480
390 CONTINUE 24490
C SYMMETRIZE 24500
IF (NCODE.EQ.2) GO TO 420 24510
DO 410 I=2,NB 24520
I1 = I - 1 24530
DO 400 J=1,I1 24540
P(J,I) = P(I,J) 24550
400 CONTINUE 24560
410 CONTINUE 24570
GO TO 770 24580
420 N = 0 24590
FN = 0.0 24600
DO 440 I=2,NB 24610
I1 = I - 1 24620
DO 430 J=1,I1 24630
IF (P(J,I).EQ.0.0) GO TO 430 24640
P(I,J) = P(I,J)/P(J,I) 24650
N = N + 1 24660
FN = FN + P(I,J) 24670
430 CONTINUE 24680
440 CONTINUE 24690
FN = FN/N 24700
DO 470 I=2,NB 24710
I1 = I - 1 24720
DO 460 J=1,I1 24730
IF (P(J,I).NE.0.0) GO TO 450 24740
P(I,J) = FN 24750
450 P(J,I) = P(I,J) 24760
460 CONTINUE 24770
470 CONTINUE 24780
GO TO 770 24790
C 24800
C FREQUENCY INPUT 24810
480 IF (IDATA.EQ.0) GO TO 490 24820
WRITE (LOUT,9991) 24830
9991 FORMAT (1H1, //41H INPUT (GENERAL DIRECTIONAL RANK ORDERS):, 24840
* 16H FREQUENCY INPUT) 24850
WRITE (LOUT,9990) 24860
9990 FORMAT (1H ) 24870
490 N = 0 24880
NY = 0 24890
DO 710 K=1,NS 24900
IF (IDATA.EQ.0) GO TO 500 24910
IF (NS.GT.1) WRITE (LOUT,9989) K 24920
9989 FORMAT (//6X, 9HDATA SET=, I3) 24930
WRITE(LOUT,9988) (I,I=1,10) 24940
9988 FORMAT (/6X, 22HNO. FREQ. N M, 4X, I4, 9I7) 24950
WRITE (LOUT,9990) 24960
500 NCRK = NCR(K) 24970
DO 700 JJ=1,NCRK 24980
N = N + 1 24990
C INPUT 25000
READ (IN,9987) F(N), NN, NSR(N) 25010
9987 FORMAT (F10.0, 2I5) 25020
IF (F(N).EQ.0.0) F(N) = 1.0 25030
IF (NN.EQ.0) NN = NBSQ2 25040
IF (NSR(N).EQ.0) NSR(N) = NN 25050
C DATA CHECK 25060
IF (F(N).LT.0.0) GO TO 510 25070
IF (NN.LT.2 .OR. NN.GT.NBSQ2) GO TO 510 25080
IF (NSR(N).LT.1 .OR. NSR(N).GT.NN) GO TO 510 25090
GO TO 520 25100
510 WRITE (LOUT,9986) 25110
9986 FORMAT (///34H THE DATA OUT OF PERMISSIBLE RANGE) 25120
GO TO 580 25130
520 NSC(N) = NN 25140
I = NY + 1 25150
NY = NY + NN 25160
C INPUT 25170
READ (IN,FM) (IZ(1,J),IZ(2,J),J=I,NY) 25180
C DATA CHECK 25190
DO 540 J=I,NY 25200
DO 530 I1=1,2 25210
IF (IZ(I1,J).LT.1 .OR. IZ(I1,J).GT.NB) GO TO 510 25220
530 CONTINUE 25230
540 CONTINUE 25240
NYM1 = NY - 1 25250
DO 560 J=I,NYM1 25260
JP1 = J + 1 25270
DO 550 J1=JP1,NY 25280
IF (IZ(1,J).EQ.IZ(1,J1) .AND. IZ(2,J).EQ.IZ(2,J1)) GO TO 570 25290
550 CONTINUE 25300
560 CONTINUE 25310
GO TO 590 25320
570 WRITE (LOUT,9985) 25330
9985 FORMAT (///51H SOME INDEX PAIRS APPEAR MORE THAN ONCE IN ONE RANK) 25340
580 IF (NS.GT.1) WRITE (LOUT,9989) K 25350
WRITE (LOUT,9988) 25360
WRITE (LOUT,9983) JJ, F(N), NN, NSR(N), (IZ(1,J),IZ(2,J),J=I,NY) 25370
STOP 25380
590 IF (NTIE.EQ.1) GO TO 610 25390
DO 600 J=I,NY 25400
IY(J) = 1 25410
600 CONTINUE 25420
GO TO 650 25430
610 READ (IN,9984) (IY(J),J=I,NY) 25440
9984 FORMAT (20I4) 25450
J1 = I - 1 25460
620 J1 = J1 + 1 25470
IF (J1.GT.NY) GO TO 650 25480
IF (IY(J1).LE.0) IY(J1) = 1 25490
IF (IY(J1).EQ.1) GO TO 620 25500
J2 = IY(J1) 25510
J11 = J1 + 1 25520
J22 = J1 + J2 - 1 25530
IF (J22.LE.NY) GO TO 630 25540
J22 = NY 25550
J2 = NY - J1 + 1 25560
630 DO 640 J=J11,J22 25570
IY(J) = J2 25580
640 CONTINUE 25590
J1 = J22 25600
GO TO 620 25610
650 IF (IDATA.EQ.1) WRITE (LOUT,9983) JJ, F(N), NN, NSR(N), 25620
* (IZ(1,J),IZ(2,J),J=I,NY) 25630
9983 FORMAT (4X, I4, F10.0, 2I5, 4X, 10(I3, 1H,, I2, 1H;)/(32X, 10(I3, 25640
* 1H,, I2, 1H;))) 25650
C INITIAL DISTANCES 25660
IF (INT.EQ.1) GO TO 700 25670
C = FNBSQ2/NN 25680
CX = C*0.5 25690
CO = (NN+NSR(N)-1)*CX 25700
IF (NMIS.EQ.1) CO = (NSR(N)+1)*CX 25710
DO 690 J=I,NY 25720
I1 = IZ(1,J) 25730
I2 = IZ(2,J) 25740
IF (I1.GT.I2) GO TO 660 25750
I0 = I1 25760
I1 = I2 25770
I2 = I0 25780
660 JX = J - I + 1 25790
IF (JX.GT.NSR(N)) GO TO 670 25800
P(I1,I2) = P(I1,I2) + JX*C*F(N) 25810
GO TO 680 25820
670 P(I1,I2) = P(I1,I2) + CO*F(N) 25830
680 P(I2,I1) = P(I2,I1) + F(N) 25840
690 CONTINUE 25850
700 CONTINUE 25860
710 CONTINUE 25870
C SYMMETRIZE 25880
IF (INT.EQ.1) GO TO 770 25890
N = 0 25900
FN = 0.0 25910
DO 730 I=2,NB 25920
I1 = I - 1 25930
DO 720 J=1,I1 25940
IF (P(J,I).EQ.0.0) GO TO 720 25950
P(I,J) = P(I,J)/P(J,I) 25960
N = N + 1 25970
FN = FN + P(I,J) 25980
720 CONTINUE 25990
730 CONTINUE 26000
FN = FN/N 26010
DO 760 I=2,NB 26020
I1 = I - 1 26030
DO 750 J=1,I1 26040
IF (P(J,I).NE.0.0) GO TO 740 26050
P(I,J) = FN 26060
740 P(J,I) = P(I,J) 26070
750 CONTINUE 26080
760 CONTINUE 26090
770 IF (INS.EQ.1) RETURN 26100
NCRK = 0 26110
DO 780 K=1,NS 26120
NCRK = NCRK + NCR(K) 26130
780 CONTINUE 26140
NCR(1) = NCRK 26150
NS = 1 26160
RETURN 26170
END 26180
$ FORTY FORM 26190
C 26200
SUBROUTINE THURS(S,G,FQ,ID,P,WK3,NB,NTS,M) 26210
C 26220
DIMENSION S(1),G(1),FQ(1),ID(1),P(NB,1),WK3(1) 26230
COMMON/OPT/IN,LOUT,NPH,TIT(20),FM(20),ME,INT,NCON,MXIT,CVCR,ITIE, 26240
1INS,IDATA,INTPR,IDIST,IVEST,IPLTC,IPHC,IPHV,NBDMX1,NDMX1,NB2,NNB2, 26250
2NBSQ2,NTIE,NBS,NMIS,NBDMX,ICODE,NCODE,NB1,FNB1,FNB,FNBI,FNB2 26260
C 26270
DO 9 I=1,NBSQ2 26280
9 S(I)=0.0 26290
DO 10 I=1,M 26300
10 G(I)=0.0 26310
C 26320
K3=0 26330
I4=0 26340
DO 11 J=1,NTS 26350
K1=K3+1 26360
K2=K3+2 26370
K3=K3+3 26380
A=FQ(K3)*0.5 26390
FN=FQ(K1)+A 26400
IF(FN.EQ.0.0) GO TO 12 26410
FD=FQ(K2)+A 26420
IF(FD.EQ.0.0) GO TO 12 26430
I1=I4+1 26440
I2=I4+2 26450
I3=I4+3 26460
I4=I4+4 26470
J1=ID(I1) 26480
J2=ID(I2) 26490
J3=(J1-2)*(J1-1)/2+J2 26500
IF(J2.GT.J1) J3=(J2-2)*(J2-1)/2+J1 26510
J1=ID(I3) 26520
J2=ID(I4) 26530
J4=(J1-2)*(J1-1)/2+J2 26540
IF(J2.GT.J1) J4=(J2-2)*(J2-1)/2+J1 26550
A=ALOG(FN/FD) 26560
S(J3)=S(J3)+A 26570
S(J4)=S(J4)-A 26580
I1=(J3+1)*J3/2 26590
I2=(J4+1)*J4/2 26600
I3=(J3-1)*J3/2+J4 26610
IF(J4.GT.J3) I3=(J4-1)*J4/2+J3 26620
G(I1)=G(I1)+1.0 26630
G(I2)=G(I2)+1.0 26640
G(I3)=G(I3)-1.0 26650
GO TO 11 26660
12 I4=I4+4 26670
11 CONTINUE 26680
C 26690
KRANK=0 26700
CALL MFSS(G,NBSQ2,0.00001,IRANK,KRANK,WK3) 26710
CALL MLSS(G,NBSQ2,IRANK,WK3,0,S,IER) 26720
C 26730
I4=0 26740
DO 13 I=2,NB 26750
I1=I-1 26760
DO 13 J=1,I1 26770
I4=I4+1 26780
P(I,J)=S(I4) 26790
13 P(J,I)=S(I4) 26800
DO 14 I=1,NB 26810
14 P(I,I)=0.0 26820
C 26830
RETURN 26840
END 26850
$ FORTY FORM 26860
SUBROUTINE MLSS(A,N,IRANK,TRAC,INC,RHS,IER) MLS26870
C MLS26880
C MLS26890
C DIMENSIONED DUMMY VARIABLES MLS26900
DIMENSION A(1),TRAC(1),RHS(1) MLS26910
DOUBLE PRECISION SUM MLS26920
C MLS26930
C TEST OF SPECIFIED DIMENSIONS MLS26940
IDEF=N-IRANK MLS26950
IF(N)33,33,1 MLS26960
1 IF(IRANK)33,33,2 MLS26970
2 IF(IDEF)33,3,3 MLS26980
C MLS26990
C CALCULATE AUXILIARY VALUES MLS27000
3 ITE=IRANK*(IRANK+1)/2 MLS27010
IX2=IRANK+1 MLS27020
NP1=N+1 MLS27030
IER=0 MLS27040
C MLS27050
C INTERCHANGE RIGHT HAND SIDE MLS27060
JJ=1 MLS27070
II=1 MLS27080
4 DO 6 I=1,N MLS27090
J=TRAC(II) MLS27100
IF(J)31,31,5 MLS27110
5 HOLD=RHS(II) MLS27120
RHS(II)=RHS(J) MLS27130
RHS(J)=HOLD MLS27140
6 II=II+JJ MLS27150
IF(JJ)32,7,7 MLS27160
C MLS27170
C PERFORM STEP 2 IF NECESSARY MLS27180
7 ISW=1 MLS27190
IF(INC*IDEF)8,28,8 MLS27200
C MLS27210
C CALCULATE X1 = X1 + U * X2 MLS27220
8 ISTA=ITE MLS27230
DO 10 I=1,IRANK MLS27240
ISTA=ISTA+1 MLS27250
JJ=ISTA MLS27260
SUM=0.D0 MLS27270
DO 9 J=IX2,N MLS27280
SUM=SUM+A(JJ)*RHS(J) MLS27290
9 JJ=JJ+J MLS27300
10 RHS(I)=RHS(I)+SUM MLS27310
GOTO(11,28,11),ISW MLS27320
C MLS27330
C CALCULATE X2 = TRANSPOSE(U) * X1 MLS27340
11 ISTA=ITE MLS27350
DO 15 I=IX2,N MLS27360
JJ=ISTA MLS27370
SUM=0.D0 MLS27380
DO 12 J=1,IRANK MLS27390
JJ=JJ+1 MLS27400
12 SUM=SUM+A(JJ)*RHS(J) MLS27410
GOTO(13,13,14),ISW MLS27420
13 SUM=-SUM MLS27430
14 RHS(I)=SUM MLS27440
15 ISTA=ISTA+I MLS27450
GOTO(16,29,30),ISW MLS27460
C MLS27470
C INITIALIZE STEP (4) OR STEP (8) MLS27480
16 ISTA=IX2 MLS27490
IEND=N MLS27500
JJ=ITE+ISTA MLS27510
C MLS27520
C DIVISION OF X1 BY TRANSPOSE OF TRIANGULAR MATRIX MLS27530
17 SUM=0.D0 MLS27540
DO 20 I=ISTA,IEND MLS27550
IF(A(JJ))18,31,18 MLS27560
18 RHS(I)=(RHS(I)-SUM)/A(JJ) MLS27570
IF(I-IEND)19,21,21 MLS27580
19 JJ=JJ+ISTA MLS27590
SUM=0.D0 MLS27600
DO 20 J=ISTA,I MLS27610
SUM=SUM+A(JJ)*RHS(J) MLS27620
20 JJ=JJ+1 MLS27630
C MLS27640
C DIVISION OF X1 BY TRIANGULAR MATRIX MLS27650
21 SUM=0.D0 MLS27660
II=IEND MLS27670
DO 24 I=ISTA,IEND MLS27680
RHS(II)=(RHS(II)-SUM)/A(JJ) MLS27690
IF(II-ISTA)25,25,22 MLS27700
22 KK=JJ-1 MLS27710
SUM=0.D0 MLS27720
DO 23 J=II,IEND MLS27730
SUM=SUM+A(KK)*RHS(J) MLS27740
23 KK=KK+J MLS27750
JJ=JJ-II MLS27760
24 II=II-1 MLS27770
25 IF(IDEF)26,30,26 MLS27780
26 GOTO(27,11,8),ISW MLS27790
C MLS27800
C PERFORM STEP (5) MLS27810
27 ISW=2 MLS27820
GOTO 8 MLS27830
C MLS27840
C PERFORM STEP (6) MLS27850
28 ISTA=1 MLS27860
IEND=IRANK MLS27870
JJ=1 MLS27880
ISW=2 MLS27890
GOTO 17 MLS27900
C MLS27910
C PERFORM STEP (8) MLS27920
29 ISW=3 MLS27930
GOTO 16 MLS27940
C MLS27950
C REINTERCHANGE CALCULATED SOLUTION MLS27960
30 II=N MLS27970
JJ=-1 MLS27980
GOTO 4 MLS27990
C MLS28000
C ERROR RETURN IN CASE OF ZERO DIVISOR MLS28010
31 IER=1 MLS28020
32 RETURN MLS28030
C MLS28040
C ERROR RETURN IN CASE OF ILLEGAL DIMENSION MLS28050
33 IER=-1 MLS28060
RETURN MLS28070
END MLS28080
$ FORTY FORM 28090
SUBROUTINE MFSS(A,N,EPS,IRANK,KRANK,TRAC) 28100
DIMENSION A(1),TRAC(1) 28110
DOUBLE PRECISION SUM 28120
IF(N)36,36,1 28130
1 IRANK=0 28140
IF(KRANK.LE.0) KRANK=N 28150
ISUB=0 28160
KPIV=0 28170
J=0 28180
PIV=0. 28190
DO 3 K=1,N 28200
J=J+K 28210
TRAC(K)=A(J) 28220
IF(A(J)-PIV)3,3,2 28230
2 PIV=A(J) 28240
KSUB=J 28250
KPIV=K 28260
3 CONTINUE 28270
DO 32 I=1,N 28280
ISUB=ISUB+I 28290
TEST=ABS(A(ISUB)*EPS) 28300
TNEG=-AMAX1(TEST,1E-6) 28310
IM1=I-1 28320
4 KMI=KPIV-I 28330
IF(KMI)35,9,5 28340
5 CALL PERMUT(A,N,KPIV,KSUB,I,ISUB) 28350
9 IF(IRANK)22,10,10 28360
10 TRAC(KPIV)=TRAC(I) 28370
TRAC(I)=KPIV 28380
KK=IM1-IRANK 28390
KMI=ISUB-KK 28400
PIV=0. 28410
IDC=IRANK+1 28420
JI=ISUB-1 28430
JK=KMI 28440
JJ=ISUB-I 28450
DO 19 K=I,N 28460
SUM=0.D0 28470
IF(KK)13,13,11 28480
11 DO 12 J=KMI,JI 28490
SUM=SUM-A(J)*A(JK) 28500
12 JK=JK+1 28510
13 JJ=JJ+K 28520
IF(K-I)14,14,16 28530
14 SUM=A(ISUB)+SUM 28540
IF(TNEG-SUM) 40,40,41 28550
41 IRANK=-3 28560
RETURN 28570
40 IF(SUM.LE.TEST.OR.(K.GT.KRANK.AND.IRANK.EQ.0)) GO TO 20 28580
15 A(ISUB)=DSQRT(SUM) 28590
KPIV=I+1 28600
GOTO 19 28610
16 SUM=(A(JK)+SUM)/A(ISUB) 28620
A(JK)=SUM 28630
IF(A(JJ))19,19,17 28640
17 TRAC(K)=TRAC(K)-SUM*SUM 28650
HOLD=TRAC(K)/A(JJ) 28660
IF(PIV-HOLD)18,19,19 28670
18 PIV=HOLD 28680
KPIV=K 28690
KSUB=JJ 28700
19 JK=JJ+IDC 28710
GOTO 32 28720
20 IF(IRANK)21,21,37 28730
21 IRANK=-1 28740
GOTO 4 28750
22 IRANK=IM1 28760
II=ISUB-IRANK 28770
JI=II 28780
DO 26 K=1,IRANK 28790
JI=JI-1 28800
JK=ISUB-1 28810
JJ=K-1 28820
DO 26 J=I,N 28830
IDC=IRANK 28840
SUM=0.D0 28850
KMI=JI 28860
KK=JK 28870
IF(JJ)25,25,23 28880
23 DO 24 L=1,JJ 28890
IDC=IDC-1 28900
SUM=SUM-A(KMI)*A(KK) 28910
KMI=KMI-IDC 28920
24 KK=KK-1 28930
25 A(KK)=(SUM+A(KK))/A(KMI) 28940
26 JK=JK+J 28950
JJ=ISUB-I 28960
PIV=0. 28970
KK=ISUB-1 28980
DO 31 K=I,N 28990
JJ=JJ+K 29000
IDC=0 29010
DO 28 J=K,N 29020
SUM=0.D0 29030
KMI=JJ+IDC 29040
DO 27 L=II,KK 29050
JK=L+IDC 29060
27 SUM=SUM+A(L)*A(JK) 29070
A(KMI)=SUM 29080
28 IDC=IDC+J 29090
A(JJ)=A(JJ)+1.D0 29100
TRAC(K)=A(JJ) 29110
IF(PIV-A(JJ))29,30,30 29120
29 KPIV=K 29130
KSUB=JJ 29140
PIV=A(JJ) 29150
30 II=II+K 29160
KK=KK+K 29170
31 CONTINUE 29180
GOTO 4 29190
32 CONTINUE 29200
33 IF(IRANK)35,34,35 29210
34 IRANK=N 29220
35 RETURN 29230
36 IRANK=-1 29240
RETURN 29250
37 IRANK=-2 29260
RETURN 29270
END 29280
$ FORTY FORM 29290
C 29300
SUBROUTINE GINV(A, N, K, B, WK, EPS, IPIV, ISTORE, IER) 29310
C SUBROUTINE FOR COMPUTING THE MOORE-PENROSE INVERSE OF A SYMMETRIC 29320
C POSITIVE SEMIDEFINITE MATRIX. 29330
DIMENSION A(1), WK(1), B(1), IPIV(1) 29340
REAL*8 SUM,WSUM 29350
IER = 0 29360
IF (N-1) 10, 20, 30 29370
10 IER = -2 29380
RETURN 29390
20 IF (A(1).GT.0.0) A(1) = 1.0/A(1) 29400
RETURN 29410
C MFSS REPLACES A BY T , U , AND X , WHERE T IS LOWER TRIANGULAR OF 29420
C ORDER K , U IS N-K BY K , AND X IS LOWER TRIANGULAR OF ORDER 29430
30 CALL MFSS(A, N, EPS, IRANK, K, WK) 29440
IF (IRANK) 40, 40, 70 29450
40 IF (3+IRANK) 50, 50, 60 29460
50 WRITE (6,9999) 29470
9999 FORMAT (52H0INPUT MATRIX FOR GINV IS NOT POSITIVE SEMIDEFINITE./ 29480
* 24H COMPUTATION TERMINATED.) 29490
IER = -3 29500
RETURN 29510
60 IER = -1 29520
RETURN 29530
70 IF (K.NE.IRANK .AND. K.GT.0) IER = 1 29540
DO 80 I=1,N 29550
IPIV(I) = WK(I) 29560
80 CONTINUE 29570
K = IRANK 29580
KM = K - 1 29590
KP = K + 1 29600
NMK = N - K 29610
NM1 = N - 1 29620
NPK = K*(K-1)/2 29630
NPKPK = NPK + K 29640
NPKP2K = NPKPK + 1 29650
NP = N*(N+1)/2 29660
C COMPUTE T-INVERSE 29670
CALL TINV(A, K) 29680
IF (N-K) 90, 90, 140 29690
C IF N = K , COMPUTE T-INVERSE T-INVERSE-TRANSPOSE 29700
90 JP = 0 29710
DO 130 J=1,K 29720
IP = JP 29730
DO 110 I=J,K 29740
SUM = 0 29750
LP = IP 29760
DO 100 L=I,K 29770
M1 = LP + I 29780
M2 = LP + J 29790
SUM = SUM + A(M1)*A(M2) 29800
LP = LP + L 29810
100 CONTINUE 29820
WK(I) = SUM 29830
IP = IP + I 29840
110 CONTINUE 29850
IP = JP 29860
DO 120 I=J,K 29870
M = IP + J 29880
A(M) = WK(I) 29890
IP = IP + I 29900
120 CONTINUE 29910
JP = JP + J 29920
130 CONTINUE 29930
GO TO 730 29940
C STORE X IN B 29950
140 M = 0 29960
IP = NPKPK 29970
DO 160 I=KP,N 29980
M1 = IP + K 29990
DO 150 J=KP,I 30000
M = M + 1 30010
M1 = M1 + 1 30020
B(M) = A(M1) 30030
150 CONTINUE 30040
IP = IP + I 30050
160 CONTINUE 30060
C COMPUTE X-INVERSE 30070
CALL TINV(B, NMK) 30080
C REPLACE X IN A BY X-INVERSE 30090
M = 0 30100
IP = NPKPK 30110
DO 180 I=KP,N 30120
M1 = IP + KP 30130
DO 170 J=KP,I 30140
M = M + 1 30150
A(M1) = B(M) 30160
M1 = M1 + 1 30170
170 CONTINUE 30180
IP = IP + I 30190
180 CONTINUE 30200
C COMPUTE W = X-INVERSE U AND STORE IN B 30210
IP = NPKPK 30220
DO 210 I=KP,N 30230
M3 = IP 30240
DO 200 J=1,K 30250
SUM = 0 30260
LP = NPKPK 30270
M1 = IP + KP 30280
DO 190 L=KP,I 30290
M2 = LP + J 30300
SUM = SUM + A(M1)*A(M2) 30310
M1 = M1 + 1 30320
LP = LP + L 30330
190 CONTINUE 30340
M3 = M3 + 1 30350
B(M3) = SUM 30360
200 CONTINUE 30370
IP = IP + I 30380
210 CONTINUE 30390
C COMPUTE R = T-INVERSE - T-INVERSE W-TRANSPOSE W AND STORE LOWER 30400
C TRIANGLE IN A AND UPPER TRIANGLE IN B 30410
IF (ISTORE) 220, 310, 310 30420
C IN-CORE STOREAGE ECONOMIZING SOLUTION 30430
220 MB = 0 30440
IP = 0 30450
DO 300 I=1,K 30460
DO 250 J=1,K 30470
SUM = 0 30480
LP = 0 30490
M1 = IP 30500
DO 240 L=1,I 30510
WSUM = 0 30520
MP = NPKPK 30530
DO 230 M=KP,N 30540
M3 = MP + L 30550
M2 = MP + J 30560
WSUM = WSUM + B(M3)*B(M2) 30570
MP = MP + M 30580
230 CONTINUE 30590
M1 = M1 + 1 30600
SUM = SUM + A(M1)*WSUM 30610
240 CONTINUE 30620
WK(J) = SUM 30630
250 CONTINUE 30640
M1 = IP 30650
DO 260 J=1,I 30660
M1 = M1 + 1 30670
A(M1) = A(M1) - WK(J) 30680
260 CONTINUE 30690
IF (K-I) 290, 290, 270 30700
270 IP1 = I + 1 30710
DO 280 J=IP1,K 30720
MB = MB + 1 30730
B(MB) = -WK(J) 30740
280 CONTINUE 30750
290 IP = IP + I 30760
300 CONTINUE 30770
GO TO 510 30780
C COMPUTE S = I - W-TRANSPOSE W AND STORE IN B 30790
310 IP = 0 30800
DO 340 I=1,K 30810
M = IP 30820
DO 330 J=1,I 30830
SUM = 0 30840
LP = NPKPK 30850
DO 320 L=KP,N 30860
M1 = LP + I 30870
M2 = LP + J 30880
SUM = SUM - B(M1)*B(M2) 30890
LP = LP + L 30900
320 CONTINUE 30910
M = M + 1 30920
B(M) = SUM 30930
330 CONTINUE 30940
B(M) = B(M) + 1.0 30950
IP = IP + I 30960
340 CONTINUE 30970
C COMPUTE R - T-INVERSE S AND STORE LOWER TRIANGLE IN A AND 30980
C WRITE UPPER TRIANGLE ONTO SCRATCH FILE ISTORE (ISTORE > 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
.