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
