C mdscal.f
C The authors of this software are Joseph B Kruskal and Judith B Seery.

C Copyright (c) 1993 by AT&T.
C Permission to use, copy, modify, and distribute this software for any
C purpose without fee is hereby granted, provided that this entire
C notice is included in all copies of any software which is or includes
C a copy or modification of this software and in all copies of the
C supporting documentation for such software.
C THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMP-
C LIED WARRANTY.  IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY
C REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
C OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.

C This software comes from the FIRST MDS Package of AT&T Bell Laboratories

C The manual of how to use kyst2a is available in several different formats:
C see mds/kyst2a_manual.  It should be quite helpful in using mdscal.

C For explanation of the method this software implements, see
C "Multidimensional Scaling" (1978) by Joseph B. Kruskal and Myron Wish,
C Sage Publications: Beverly Hills, Calif.,

C "Multidimensional Scaling and Other Methods for Discovering Structure"
C (1977) by Joseph B. Kruskal, pp 296-339 in
C "Statistical Methods for Digital Computers" by Kurt Enslein,
C Anthony Ralston, and Herbert S. Wilf, Wiley: New York,

C "Multidimensional Scaling by Optimizing Goodness of Fit to a Nonmetric
C Hypothesis" (1964) by Joseph B. Kruskal in Psychometrika 29, pp 1-27,

C "Nonmetric Multidimensional Scaling: A Numerical Method" (1964) by
C Joseph B. Kruskal in Psychometrika 29(2) pp 115-129.
C----+----@----+----@----+----@----+----@----+----@----+----@----+----@

$      FORTRAN
C5MS                                                                    5MS    1
C                                                                       5MS    2
CMDSCAL            MDSCAL  MAIN PROGRAM                                 5MS    3
C MULTIDIMENSIONAL SCALING PROGRAM, VERSION 5MS, OCTOBER, 1971          5MS    4
C        DECEMBER, 1971, RELEASE                                        5MS    5
C        SUBROUTINES CCACT AND BLOCK DATA CORRECTED DECEMBER, 1971      5MS    6
C       ADAPTED BY J. KRUSKAL AND J. SEERY, NOVEMBER, 1971              5MS    7
C                                                                       5MS    8
C MULTIDIMENSIONAL SCALING PROGRAM. VERSION 5M                          5MS    9
C                                                                       5MS   10
C                                                                       5MS   11
C        THIS PROGRAM WAS WRITTEN BY                                    5MS   12
C                                                                       5MS   13
C             DR. J. B. KRUSKAL                                         5MS   14
C             BELL TELEPHONE LABORATORIES                               5MS   15
C             MURRAY HILL, N. J.                                        5MS   16
C                                                                       5MS   17
C                                                                       5MS   18
C        ADAPTED FOR IBM S/360 BY                                       5MS   19
C                                                                       5MS   20
C             F. J. CARMONE, JR.                                        5MS   21
C             DEPARTMENT OF MARKETING                                   5MS   22
C             DREXEL UNIVERSITY                                         5MS   23
C             PHILADELPHIA, PENNSYLVANIA    19104                       5MS   24
C             UNITED STATES OF AMERICA                                  5MS   25
C                                                                       5MS   26
C        DISTRIBUTED BY                                                 5MS   27
C                                                                       5MS   28
C             MARKETING SCIENCE INSTITUTE                               5MS   29
C             16 STORY STREET                                           5MS   30
C             CAMBRIDGE, MASS. 02138                                    5MS   31
C                                                                       5MS   32
C        OCTOBER 1, 1970                                                5MS   33
C                                                                       5MS   34
C                                                                       5MS   35
C             CHANGES MADE FOR "FIX" OPTION JAN 14, 1970                5MS   36
C             ROUTINES CHANGED  MAIN, CCACT, BLDA                       5MS   37
C                                                                       5MS   38
C                                                                       5MS   39
C                                                                       5MS   40
C             CHANGES MADE FOR SAVE CONFIGURATION AND RANDOM START      5MS   41
C             ON MAY 17, 1971. ROUTINES CHANGED  MAIN                   5MS   42
C                                                                       5MS   43
C                                                                       5MS   44
C GENERAL REMARKS.                                                      5MS   45
C                                                                       5MS   46
C THIS PROGRAM CONSISTS OF THE FOLLOWING ROUTINES.                      5MS   47
C      MDSCAL  MAIN ROUTINE                                             5MS   48
C      FIT     PERFORMS WEIGHTED-LEAST-SQUARES REGRESSION.              5MS   49
C      NEWSTP  FINDS NEW STEP SIZE FOR ITERATIVE MINIMIZATION.          5MS   50
C      CCACT   READS AND INTERPRETS CONTROL CARDS.                      5MS   51
C      SORT    SORTS ARRAYS. SIMPLE AND FAST,BEST AVAILABLE.            5MS   52
C      PLOTR   PLOTS DIST AND DHAT VERSUS DATA                          5MS   53
C      PLOT1   PLOTS THE TWO SPACE CONFIGURATION AND STRESS VS. DIMEN.  5MS   54
C NOTE   ALTHOUGH TWO PLOTTING ROUTINES ARE USED ARRAYS ARE EQUIVALENCED5MS   55
C        FOR MINIMUM SPACE.                                             5MS   56
C                                                                       5MS   57
C                                                                       5MS   58
C ALL ARE WRITTEN IN FORTRAN IV.                                        5MS   59
C                                                                       5MS   60
C NO USE IS MADE OF SPECIAL OR NON-STANDARD SOFTWARE.                   5MS   61
C                                                                       5MS   62
C ALL NORMAL INPUT-OUTPUT HANDLED BY  MDSCAL  AND  CCACT.               5MS   63
C                                                                       5MS   64
C FIT  HAS EMERGENCY DIAGNOSTIC OUTPUT.                                 5MS   65
C ALL INPUT AND OUTPUT IS ONTO FORTRAN LOGICAL UNITS WITH NUMBERS       5MS   66
C CONTROLLED BY THESE VARIABLES  LREAD,LPRINT,LPUNCH,LSCRAT.            5MS   67
C UNIT NUMBERS 5,6,43,2 HAVE BEEN USED RESPECTIVELY                     5MS   68
C TO CHANGE THESE ASSIGNMENTS, IT SUFFICES TO CHANGE THE VALUES         5MS   69
C FOR THE FOUR VARIABLES SET IN THE BLOCK DATA SUBPROGRAM.              5MS   70
C                                                                       5MS   71
C NOTE THAT THE SCRATCH UNIT IS USED IN A VERY MINOR WAY.               5MS   72
C IT IS USED ONLY BY  CCACT                                             5MS   73
C MANY INSTALLATIONS WILL HAVE ALTERNATE METHODS OF DOING THE SAME THING5MS   74
C                                                                       5MS   75
C      MDSCAL  THIS IS THE MAIN ROUTINE.                                5MS   76
C      IT DIRECTLY USES FIT, NEWSTP, CCACT, SORT.                       5MS   77
C                                                                       5MS   78
C                                                                       5MS   79
C      CARDS PUNCHED ACCORDING TO THE IBMEL 029 CHARACTER SET           5MS   80
C       NOTE   CONTROL PHRASE "TITLE" HAS BEEN DELETED                  5MS   81
C      NOTE   FUNCTIONS RPOWER, RM1POW, RROOT REPLACE RFUNCT            5MS   82
C                                                                       5MS   83
       DIMENSION   DATA(1800), IJ(1800), DIST(1800), DHAT(1800)         5MS   84
      DIMENSION                WW(1800),LDIST(1800)                     5MS   85
       DIMENSION    X(100,6), GR(100,6),  GL(100,6)                     5MS   86
      DIMENSION CTITLE(80),FMAT(80),TITLE(80)                           5MS   87
       INTEGER    GRNO(101), NOGRPS, LSPLIT,         SPLITH             5MS   88
       INTEGER GRSDIS(100), SDSWIT,SDSWT1                               5MS   89
       REAL    GRSTRS(100),PCOEFF(5)                                    5MS   90
      REAL SMALL(25),HOLL(11),CHAR(50),LABELS(55),ITEM(55,121)          5MS   91
      DIMENSION  STPL(10),DIMSV(10)                                     5MS   92
       INTEGER TWO9, TWO18 ,OUT                                         5MS   93
      INTEGER RTYPE                                                     5MS   94
C                                                                       5MS   95
C                                                                       5MS   96
       COMMON  /MDSCL1/  LREAD, LPRINT, LPUNCH, LSCRAT                  5MS   97
       COMMON  /MDSCL2/                                                 5MS   98
     1         LDIMX, LDIMN, LDIMD, CUTOFF, STRMIN                      5MS   99
     2 ,       SFGRMN, COSAVW, ACSAVW, IFIRST, MATSWP                   5MS  100
     3 ,       SDSWIT,  LCSWIT, LFITSW, R, NOIT                         5MS  101
     4 ,       SRATST, LSCH, LPUNSW, LSPL, LRANDC                       5MS  102
     5 ,       LDAPRT, LDIPRT, LREG,   LHIPRT, NUDASW                   5MS  103
     6 ,        LPOLSP, LCONSW, LNFIXZ, DCON1, DCON2                    5MS  104
     7 ,       DCON3, DCON4, DCON5, WCON1, WCON2                        5MS  105
     8 ,       WCON3, WCON4, WCON5                                      5MS  106
      COMMON  /MDSCL4/ PCOEFF                                           5MS  107
      COMMON /METRIC/ RTYPE,RM1,RECR,RR                                 5MS  108
      COMMON  /PLTSS/ SMALL,HOLL,CHAR,LABELS,ITEM,DD,BLANK,AIE,AMINUS   5MS  109
C                                                                       5MS  110
       EQUIVALENCE (DIST,LDIST)                                         5MS  111
       EXTERNAL                    WTRAN       , DTRAN                  5MS  112
C                                                                       5MS  113
       DATA    TWO9, TWO18, NPREVZ  /512,  262144,  0/                  5MS  114
C                                                                       5MS  115
C      INITIALIZE PARAMETERS                                            5MS  116
C                                                                       5MS  117
 100   CONTINUE                                                         5MS  118
C                                                                       5MS  119
C               ALPHABETICAL ORDER                                      5MS  120
       ACSAVW=0.66                                                      5MS  121
       COSAVW=0.66                                                      5MS  122
       CUTOFF=-1.23E+20                                                 5MS  123
      DCON1=1.0                                                         5MS  124
      DCON2=1.0                                                         5MS  125
      DCON3=0.0                                                         5MS  126
      DCON4=1.0                                                         5MS  127
      DCON5=0.0                                                         5MS  128
       GRNO(1) = 1                                                      5MS  129
       IFIRST=1                                                         5MS  130
      ISTC=0                                                            5MS  131
       LCONSW=1                                                         5MS  132
       LCSWIT=1                                                         5MS  133
       LDAPRT = 1                                                       5MS  134
       LDIMD=1                                                          5MS  135
       LDIMN=2                                                          5MS  136
       LDIMX=2                                                          5MS  137
       LDIPRT = 1                                                       5MS  138
       LFITSW=1                                                         5MS  139
       LHIPRT = 2                                                       5MS  140
       LNFIXZ = 0                                                       5MS  141
       LPOLSP = 100                                                     5MS  142
       LPUNSW = 2                                                       5MS  143
       LRANDC=-99                                                       5MS  144
       LREG=0                                                           5MS  145
       LSCH=2                                                           5MS  146
       LSPL=302                                                         5MS  147
       MATSWP=101                                                       5MS  148
       NOIT=50                                                          5MS  149
       NPREVZ = 0                                                       5MS  150
       NREPL1=1                                                         5MS  151
       NREPL3 = 1                                                       5MS  152
       NUDASW=1                                                         5MS  153
       R=2.0                                                            5MS  154
       SDSWIT = 10                                                      5MS  155
       SFGRMN=0.0                                                       5MS  156
       SRATST=0.999                                                     5MS  157
       STRESS=1.0                                                       5MS  158
       STRMIN=0.01                                                      5MS  159
      WCON1=1.0                                                         5MS  160
      WCON2=1.0                                                         5MS  161
      WCON3=0.0                                                         5MS  162
      WCON4=1.0                                                         5MS  163
      WCON5=0.0                                                         5MS  164
       WRITE(LPRINT, 17)                                                5MS  165
 17   FORMAT(1H1)                                                       5MS  166
C                                                                       5MS  167
C      CONTROL CARD READ                                                5MS  168
C                                                                       5MS  169
 1000  LCSWIT=1                                                         5MS  170
       CALL CCACT                                                       5MS  171
       GO TO (1000, 1100, 1165, 1182, 1200, 2000, 9000, 1190 ), LCSWIT  5MS  172
C                                                                       5MS  173
C      DATA READ                                                        5MS  174
C                                                                       5MS  175
 1100  CONTINUE                                                         5MS  176
       LSPLIT=LSPL/100                                                  5MS  177
       SPLITH = MOD(LSPL,100)                                           5MS  178
       MATSW = MOD(MATSWP,100)                                          5MS  179
       LBLKDS = MATSWP/100                                              5MS  180
      LPOSW = MOD(LPOLSP,100)                                           5MS  181
      LFITRM = LPOLSP/100                                               5MS  182
      IF(LREG.EQ.0) GO TO 1104                                          5MS  183
      SDSWIT=LREG                                                       5MS  184
      IF(LREG.LT.10) SDSWIT = SDSWIT+LPOSW                              5MS  185
1104  CONTINUE                                                          5MS  186
      SDSWIT = SDSWIT+100*LFITRM                                        5MS  187
C                                                                       5MS  188
       IF(NUDASW.NE.1) MLASTD=MM                                        5MS  189
       IF(NUDASW.NE.1) GO TO 1106                                       5MS  190
       NUDASW=2                                                         5MS  191
       MLASTD=0                                                         5MS  192
       NOGRPS=0                                                         5MS  193
1106   M=MLASTD                                                         5MS  194
       MA = M+1                                                         5MS  195
C                                                                       5MS  196
      READ(LREAD,10)   TITLE                                            5MS  197
      WRITE(LPRINT,11) TITLE                                            5MS  198
       IF(MATSW.GE.4) GO TO 1102                                        5MS  199
 1101  READ     (LREAD,12)     NPART,NREPL2,NREPL4                      5MS  200
       WRITE    (LPRINT,13)    NPART,NREPL2,NREPL4                      5MS  201
       NROWS = NPART                                                    5MS  202
       NCOLS = NPART                                                    5MS  203
       GO TO 1103                                                       5MS  204
 1102  READ     (LREAD,12)     NROWS,NCOLS,NREPL2,NREPL4                5MS  205
       WRITE    (LPRINT,13)    NROWS,NCOLS,NREPL2,NREPL4                5MS  206
       NPART = NROWS+NCOLS                                              5MS  207
C                                                                       5MS  208
 1103  N=NPART                                                          5MS  209
       IF(NREPL4.NE.0) NREPL3 = NREPL4                                  5MS  210
       IF(NREPL2.NE.0) NREPL1 = NREPL2                                  5MS  211
       IF(LBLKDS.EQ.2) N = NPART*NREPL3                                 5MS  212
       IF(NREPL2.EQ.0  .OR.  NREPL4.EQ.0  )  WRITE  (LPRINT,16)         5MS  213
 1105 READ(LREAD,10)     FMAT                                           5MS  214
      WRITE(LPRINT,11)   FMAT                                           5MS  215
C                                                                       5MS  216
       DO 1162 NREPL=1,NREPL3                                           5MS  217
       IJBLKD=0                                                         5MS  218
       IF(LBLKDS.EQ.2) IJBLKD =(NREPL-1)*NPART                          5MS  219
       MLASTG=M                                                         5MS  220
       IA = 1                                                           5MS  221
       IB = NROWS                                                       5MS  222
       IF(MATSW.EQ.2) IA = IFIRST                                       5MS  223
       IF(MATSW.EQ.3) IB = NROWS-(IFIRST-1)                             5MS  224
C                                                                       5MS  225
       DO 1160 I = IA,IB                                                5MS  226
       MLASTR=M                                                         5MS  227
       LTEMP = NCOLS                                                    5MS  228
       IF(MATSW.EQ.2) LTEMP = I-IFIRST+1                                5MS  229
       IF(MATSW.EQ.3) LTEMP = NCOLS-I-IFIRST+2                          5MS  230
       MB = MA + LTEMP * NREPL1 - 1                                     5MS  231
       ITRUE = I+IJBLKD                                                 5MS  232
       IF(MATSW.EQ.4) ITRUE = ITRUE+NCOLS                               5MS  233
C                                                                       5MS  234
 1130  READ           (LREAD, FMAT) (DATA(MP),MP=MA,MB)                 5MS  235
C                                                                       5MS  236
       DO 1150 MP=MA,MB                                                 5MS  237
       IF( DATA(MP)-CUTOFF )  1150, 1150, 1140                          5MS  238
 1140  CONTINUE                                                         5MS  239
       M = M + 1                                                        5MS  240
       DATA(M)=DATA(MP)                                                 5MS  241
       WW(M) = 1.0                                                      5MS  242
       J=((MP-MA)/NREPL1)+1                                             5MS  243
       IF(MATSW.EQ.3) J = J + (I-1) + (IFIRST-1)                        5MS  244
       J = J+IJBLKD                                                     5MS  245
       IF(MATSW.EQ.5) J = J + NROWS                                     5MS  246
       IJ(M) = TWO9 * ITRUE + J                                         5MS  247
       LDIST(M) = MP                                                    5MS  248
 1150  CONTINUE                                                         5MS  249
C                                                                       5MS  250
       IF(LSPLIT.EQ.1 .AND. M.GT.MLASTR) GO TO 1155                     5MS  251
       GO TO 1160                                                       5MS  252
 1155  NOGRPS=NOGRPS+1                                                  5MS  253
       GRNO(NOGRPS+1) = M+1                                             5MS  254
       GRSDIS(NOGRPS)=SDSWIT                                            5MS  255
 1160  MA = M + 1                                                       5MS  256
C                                                                       5MS  257
       IF(LSPLIT.EQ.2 .AND. M.GT.MLASTG) GO TO 1161                     5MS  258
       GO TO 1162                                                       5MS  259
 1161  NOGRPS=NOGRPS+1                                                  5MS  260
       GRNO(NOGRPS+1) = M+1                                             5MS  261
       GRSDIS(NOGRPS)=SDSWIT                                            5MS  262
 1162  CONTINUE                                                         5MS  263
C                                                                       5MS  264
       IF(LSPLIT.EQ.3 .AND. M.GT.MLASTD) GO TO 1163                     5MS  265
       IF(LSPLIT.EQ.4 .AND. SPLITH.EQ.2 .AND. M.GT.MLASTD) GO TO 1163   5MS  266
       GO TO 1164                                                       5MS  267
 1163  NOGRPS=NOGRPS+1                                                  5MS  268
       GRSDIS(NOGRPS)=SDSWIT                                            5MS  269
C                                                                       5MS  270
 1164  CONTINUE                                                         5MS  271
       GRNO(NOGRPS+1) = M+1                                             5MS  272
       IF(LSPLIT.LE.3) SPLITH=2                                         5MS  273
       IF(LSPLIT.EQ.4) SPLITH=1                                         5MS  274
C                                                                       5MS  275
       LSPL = 100 * LSPLIT + SPLITH                                     5MS  276
C                                                                       5MS  277
       MM = M                                                           5MS  278
       GO TO 1000                                                       5MS  279
C                                                                       5MS  280
C      WEIGHTS READ                                                     5MS  281
C                                                                       5MS  282
 1165  CONTINUE                                                         5MS  283
       M=MLASTD                                                         5MS  284
       MA=M+1                                                           5MS  285
       DO 1181 NREPL=1,NREPL3                                           5MS  286
       IA = 1                                                           5MS  287
       IB = NROWS                                                       5MS  288
       IF(MATSW.EQ.2) IA = IFIRST                                       5MS  289
       IF(MATSW.EQ.3) IB = NROWS-(IFIRST-1)                             5MS  290
       DO 1180 I = IA,IB                                                5MS  291
       LTEMP = NCOLS                                                    5MS  292
       IF(MATSW.EQ.2) LTEMP = I-IFIRST+1                                5MS  293
       IF(MATSW.EQ.3) LTEMP = NCOLS-I-IFIRST+2                          5MS  294
       MB = MA + (LTEMP-1)*NREPL1                                       5MS  295
 1172  READ           (LREAD,FMAT)(WW(MP),MP=MA,MB)                     5MS  296
       DO 1177 MP=MA,MB                                                 5MS  297
       IF(LDIST(M+1).GT.MP) GO TO 1177                                  5MS  298
       M = M + 1                                                        5MS  299
       WW(M)=WW(MP)                                                     5MS  300
 1177  CONTINUE                                                         5MS  301
 1180  MA = M + 1                                                       5MS  302
 1181  CONTINUE                                                         5MS  303
       GO TO 1000                                                       5MS  304
C                                                                       5MS  305
C      WEIGHT FORMATION BY WFUNCTION                                    5MS  306
C                                                                       5MS  307
 1182  CONTINUE                                                         5MS  308
      WRITE(LPRINT,18) WCON1,WCON2,WCON3,WCON4,WCON5                    5MS  309
 18   FORMAT(62H0WEIGHTS FORMED ACCORDING TO RULE  WW(I)=A*(B*DATA(I)+C)5MS  310
     .**D +E,/ 10H WHERE  A=,F15.7,5X,2HB=,F15.7,5X,2HC=,F15.7,5X,      5MS  311
     .2HD=,F15.7,5X,2HE=,F15.7)                                         5MS  312
       MA=MLASTD+1                                                      5MS  313
       DO 1185 M=MA,MM                                                  5MS  314
       TEMP1=DATA(M)                                                    5MS  315
       WW(M) = WTRAN(TEMP1)                                             5MS  316
 1185  CONTINUE                                                         5MS  317
       GO TO 1000                                                       5MS  318
C                                                                       5MS  319
C     DATA TRANSFORMATION BY DFUNCTION                                  5MS  320
C                                                                       5MS  321
 1190 CONTINUE                                                          5MS  322
      WRITE(LPRINT,19) DCON1,DCON2,DCON3,DCON4,DCON5                    5MS  323
 19   FORMAT(67H0DATA TRANSFORMED ACCORDING TO RULE   DATA(I)=A*(B*DATA(5MS  324
     .I)+C)**D +E,/ 10H WHERE  A=,F15.7,5X,2HB=,F15.7,5X,2HC=,F15.7,5X, 5MS  325
     .2HD=,F15.7,5X,2HE=,F15.7)                                         5MS  326
      MA=MLASTD+1                                                       5MS  327
      DO 1195 M=MA,MM                                                   5MS  328
      DATA(M)=DTRAN(DATA(M))                                            5MS  329
 1195 CONTINUE                                                          5MS  330
      GO TO 1000                                                        5MS  331
C                                                                       5MS  332
C      CONFIGURATION READ                                               5MS  333
C                                                                       5MS  334
 1200  LCONSW=2                                                         5MS  335
      READ(LREAD,10)  CTITLE                                            5MS  336
      WRITE(LPRINT,11)CTITLE                                            5MS  337
       READ           (LREAD, 12) NCON, LDIMCO                          5MS  338
       WRITE            (LPRINT, 13) NCON, LDIMCO                       5MS  339
      READ(LREAD,10)     FMAT                                           5MS  340
      WRITE(LPRINT,11)   FMAT                                           5MS  341
       DO 1210 I=1,NCON                                                 5MS  342
       READ           (LREAD, FMAT) (X(I,L),L=1,LDIMCO)                 5MS  343
       WRITE(LPRINT,42) I,(X(I,L),L=1,LDIMCO)                           5MS  344
 1210  CONTINUE                                                         5MS  345
       GO TO 1000                                                       5MS  346
C                                                                       5MS  347
C      SOME INPUT FORMATS                                               5MS  348
C                                                                       5MS  349
 10   FORMAT (80A1)                                                     5MS  350
 11   FORMAT (1H0,80A1)                                                 5MS  351
   12  FORMAT(24I3)                                                     5MS  352
   13  FORMAT(1X,24I3)                                                  5MS  353
  15  FORMAT(20F4.0)                                                    5MS  354
   16  FORMAT(118H THE NUMBER OF INTERNAL REPLICATES, OR THE NUMBER OF E5MS  355
     1XTERNAL REPLICATES HAS NOT BEEN SPECIFIED. IT IS TAKEN TO BE 1 . )5MS  356
C                                                                       5MS  357
C      COMPUTATION                 *************************************5MS  358
C                                                                       5MS  359
 2000 CONTINUE                                                          5MS  360
C   *****  PRINT GRAPH IDENTIFIERS     *****                            5MS  361
      OUT=LPRINT                                                        5MS  362
      WRITE(OUT,156)                                                    5MS  363
 156  FORMAT(1H0,10X,100HON THE SHEPARD DIAGRAM THE ORIGINAL DATA (DATA)5MS  364
     1 ARE PLOTTED ON THE Y AXIS AND DISTANCES (DIST,0) AND,/,10X       5MS  365
     2102HESTIMATED DISTANCES (DHAT,X) ON THE X AXIS. A + INDICATES TWO 5MS  366
     3VALUES ARE PLOTTED ON TOP OF EACH OTHER.,//////,10X,              5MS  367
     4 73HIDENTIFIERS FOR THE CONFIGURATION PLOT IN TWO DIMENSIONS ARE G5MS  368
     5IVEN BELOW ,/)                                                    5MS  369
      WRITE(OUT,151)  (J, J=1,25)                                       5MS  370
 151   FORMAT(1H0,40X,24HPOINT IDENTIFICATION KEY,//,5H PT #,25I5)      5MS  371
      WRITE(OUT,152)   (CHAR(I), I=1,25)                                5MS  372
 152  FORMAT(1X,4HCHAR,25(4X,A1)/)                                      5MS  373
      WRITE(OUT,153)   (J, J=26,50)                                     5MS  374
 153  FORMAT(5H0PT #,25I5)                                              5MS  375
      WRITE(OUT,154)  (CHAR(I), I=26,50)                                5MS  376
 154  FORMAT(1X,4HCHAR,25(4X,A1))                                       5MS  377
      WRITE(OUT,155)                                                    5MS  378
 155  FORMAT(60H0 ALL POINT NUMBERS ABOVE 50 ARE DESIGNATED BY THE SYMBO5MS  379
     2L  >,//,68HMULTIPLE POINTS AT THE SAME LOCATION ARE DESIGNATED BY 5MS  380
     3THE SYMBOL  ;)                                                    5MS  381
       FN=FLOAT (N)                                                     5MS  382
       SQRTN=SQRT (FN)                                                  5MS  383
       FNGRPS = FLOAT(NOGRPS)                                           5MS  384
       LDIM=LDIMX                                                       5MS  385
      RR=R                                                              5MS  386
      RTYPE=3                                                           5MS  387
      IF(R.EQ.1.0) RTYPE=1                                              5MS  388
      IF(R    .EQ.2.0) RTYPE=2                                          5MS  389
      RM1=R-1.0                                                         5MS  390
      RECR=1.0/R                                                        5MS  391
       IF(LDAPRT.EQ.2)                                                  5MS  392
     1 CALL DATAPR(DATA,IJ,WW,GRNO,MM,NOGRPS,DIST,DHAT,  1  )           5MS  393
C                                                                       5MS  394
       WRITE  (LPRINT,23)                                               5MS  395
C                                                                       5MS  396
C        FINISH STARTING CONFIGURATION. LCONSW = 1,2,3.                 5MS  397
C             2=WAS READ IN, 3=SAVED, 1=NEITHER. IF LRANDC=-99          5MS  398
C             USE ARBITRARY, ELSE USE RANDOM                            5MS  399
C                                                                       5MS  400
 2100  CONTINUE                                                         5MS  401
       ISTCON = 1                                                       5MS  402
       IF( LCONSW .EQ. 2 ) ISTCON = NCON + 1                            5MS  403
       IF( LCONSW .EQ. 3 ) ISTCON = NPREVZ + 1                          5MS  404
       IF( ISTCON .GT. N ) GO TO 2200                                   5MS  405
       IF(LRANDC.LE.0) GO TO 2110                                       5MS  406
       DO 2105  K=1,LRANDC                                              5MS  407
 2105  TEMP=RUNIFV(1)                                                   5MS  408
C2110  TEMP1=1.0                                                        5MS  409
 2110  TEMP1=0.01                                                       5MS  410
       DO 2130 I=ISTCON,N                                               5MS  411
       DO 2120 L=1,LDIMX                                                5MS  412
 2120  X(I,L)=0.0                                                       5MS  413
       K= MOD (I-1,LDIM) +1                                             5MS  414
       X(I,K)=TEMP1                                                     5MS  415
      IF(LRANDC.EQ.(-99)) GO TO 2130                                    5MS  416
       DO 2125 L=1,LDIMX                                                5MS  417
       TEMP = RUNIFV(1)                                                 5MS  418
 2125  X(I,L) = ALOG( TEMP/(1.0-TEMP))                                  5MS  419
C2130  TEMP1=TEMP1+1.0                                                  5MS  420
 2130  TEMP1=TEMP1+0.01                                                 5MS  421
C                                                                       5MS  422
C      SORT DATA AND IJ AND WW.                                         5MS  423
C      ALSO RECORD BLOCKS OF EQUAL DATA VALUES.                         5MS  424
C                                                                       5MS  425
 2200  CONTINUE                                                         5MS  426
       NPREVZ = N                                                       5MS  427
       DO 2250 NG = 1,NOGRPS                                            5MS  428
       MX = GRNO(NG)                                                    5MS  429
       MY = GRNO(NG+1) - 1                                              5MS  430
       MZ = MY - MX + 1                                                 5MS  431
       SDSWIT = GRSDIS(NG)                                              5MS  432
      SDSWIT = MOD(SDSWIT,100)                                          5MS  433
      IF(SDSWIT.EQ.11)SDSWIT = -10                                      5MS  434
       CALL SORT( DATA(MX),MZ,IJ(MX),WW(MX),DUMMY,2,SDSWIT)             5MS  435
C                                                                       5MS  436
C               SORT  WILL SORT THE MM ELEMNTS OF  DATA  IN ALGEBRAIC   5MS  437
C              ORDER, ASCENDING OR DESCENDING ACCORDING TO WHETHER      5MS  438
C              SDSWIT IS + OR -.                                        5MS  439
C              AT THE SAME TIME, THE ELEMENTS IN  IJ  AND IN  WW  WILL B5MS  440
C              REARRANGED IN EXACTLY THE SAME ORDER. THUS THE           5MS  441
C              CORRESPONDENCE BETWEEN THE ELEMENTS OF  DATA  AND  IJ    5MS  442
C              AND  WW  IS PRESERVED.                                   5MS  443
C                                                                       5MS  444
C                                                                       5MS  445
       K=1                                                              5MS  446
       MA = MX                                                          5MS  447
       DO 2240 MB = MX,MY                                               5MS  448
       IF(DATA(MB+1).NE.DATA(MB) .OR. MB .EQ. MY ) GO TO 2220           5MS  449
 2210  K=K+1                                                            5MS  450
       GO TO 2240                                                       5MS  451
 2220  DO 2230 M=MA,MB                                                  5MS  452
      IJ(M) = MOD(IJ(M),TWO 18)                                         5MS  453
 2230  IJ(M) = IJ(M) + TWO18 * K                                        5MS  454
       K=1                                                              5MS  455
       MA=MB+1                                                          5MS  456
 2240  CONTINUE                                                         5MS  457
 2250  CONTINUE                                                         5MS  458
C                                                                       5MS  459
C      START COMPUTATION IN CURRENT DIMENSION                           5MS  460
C                                                                       5MS  461
 2300  FLDIM=FLOAT (LDIM)                                               5MS  462
       ITNO=0                                                           5MS  463
       COSAV=0.0                                                        5MS  464
       SRATAV=0.8                                                       5MS  465
       ACSAV=0.0                                                        5MS  466
       STEP = 0.0                                                       5MS  467
       NBAKUP = 0                                                       5MS  468
C                                                                       5MS  469
C      INITIALIZE GRADIENT                                              5MS  470
C                                                                       5MS  471
 2400  TEMP1=SQRT (1.0/FLDIM)                                           5MS  472
       DO 2410 I=1,N                                                    5MS  473
       DO 2410 L=1,LDIM                                                 5MS  474
 2410  GR(I,L)=TEMP1                                                    5MS  475
       SFGR=SQRTN                                                       5MS  476
C                                                                       5MS  477
C      PRINT HEADING FOR INFORMATION ABOUT SCALING IN CURRENT DIMENSION 5MS  478
C                                                                       5MS  479
 2500 WRITE(LPRINT,11) TITLE                                            5MS  480
       WRITE  (LPRINT, 20) N, MM, NOGRPS, LDIM                          5MS  481
       IF(LHIPRT.EQ.2)  WRITE  (LPRINT,21)                              5MS  482
       WRITE            (LPRINT, 22)                                    5MS  483
 20    FORMAT(28H0HISTORY OF COMPUTATION. N= , I4,                      5MS  484
     1 15H.     THERE ARE  , I6,                                        5MS  485
     2 26H   DATA VALUES, SPLIT INTO  , I4,                             5MS  486
     3 27H    LISTS.     DIMENSION =  , I4     )                        5MS  487
 21    FORMAT(52H0ITERATION STRESS   SRAT SRATAV CAGRGL  COSAV  ACSAV,  5MS  488
     1 16H    SFGR    STEP )                                            5MS  489
   22  FORMAT(1X)                                                       5MS  490
 23    FORMAT(1H1)                                                      5MS  491
C                                                                       5MS  492
C      START CURRENT ITERATION     *************************************5MS  493
C                                                                       5MS  494
C      NORMALIZE CONFIGURATION. MOVE AND CLEAR GRADIENT.                5MS  495
C                                                                       5MS  496
 3000  TEMP1=0.0                                                        5MS  497
       DO 3030 L=1,LDIM                                                 5MS  498
       TEMP2=0.0                                                        5MS  499
       DO 3010 I=1,N                                                    5MS  500
 3010  TEMP2=TEMP2+X(I,L)                                               5MS  501
       TEMP2=TEMP2/FN                                                   5MS  502
       DO 3020 I=1,N                                                    5MS  503
       X(I,L)=X(I,L)-TEMP2                                              5MS  504
 3020  TEMP1=TEMP1+X(I,L)**2                                            5MS  505
 3030  CONTINUE                                                         5MS  506
       TEMP1=SQRT (FN/TEMP1)                                            5MS  507
       DO 3040 L=1,LDIM                                                 5MS  508
       DO 3040 I=1,N                                                    5MS  509
       X(I,L)=TEMP1*X(I,L)                                              5MS  510
       GL(I,L)=TEMP1*GR(I,L)                                            5MS  511
 3040  GR(I,L)=0.0                                                      5MS  512
       SFGL=TEMP1*SFGR                                                  5MS  513
C                                                                       5MS  514
       STBAMU = TEMP1                                                   5MS  515
C      LOOP THROUGH THE DATA GROUPS                                     5MS  516
C                                                                       5MS  517
       STRLST = STRESS                                                  5MS  518
       STRESS = 0.0                                                     5MS  519
       DO 3340 NG = 1,NOGRPS                                            5MS  520
       MX = GRNO(NG)                                                    5MS  521
       MY = GRNO(NG+1) - 1                                              5MS  522
       MZ = MY - MX + 1                                                 5MS  523
       SDSWIT = GRSDIS(NG)                                              5MS  524
      LFITRM = SDSWIT/100                                               5MS  525
      SDSWIT = MOD(SDSWIT,100)                                          5MS  526
      IF(SDSWIT.EQ.11)SDSWIT = -10                                      5MS  527
C                                                                       5MS  528
C      COMPUTE DISTANCES AND FIND BEST-FITTING MONOTONE PSEUDO-DISTANCES5MS  529
C                                                                       5MS  530
       SUMW = 0.0                                                       5MS  531
       DBAR = 0.0                                                       5MS  532
       DO 3120 M=MX,MY                                                  5MS  533
       LTEMP1 = MOD(IJ(M),TWO18)                                        5MS  534
       I = LTEMP1/TWO9                                                  5MS  535
       J = MOD(LTEMP1,TWO9)                                             5MS  536
       TEMP1=0.0                                                        5MS  537
       DO 3110 L=1,LDIM                                                 5MS  538
 3110  TEMP1=TEMP1+RPOWER (X(I,L)-X(J,L))                               5MS  539
       DIST(M)=RROOT (TEMP1)                                            5MS  540
      DBAR=DBAR+DIST(M)*WW(M)                                           5MS  541
       SUMW = SUMW + WW(M)                                              5MS  542
 3120  CONTINUE                                                         5MS  543
      DBAR=DBAR/SUMW                                                    5MS  544
C              DBAS IS EITHER  0  OR  DBAR  ACCORDING TO WHETHER        5MS  545
C              STRESS FORMULA 1 OR 2 IS BEING USED.                     5MS  546
       DBAS = 0.0                                                       5MS  547
       IF(LSCH.EQ.2) DBAS = DBAR                                        5MS  548
       IF(IABS(SDSWIT).GE.10)                                           5MS  549
     1 CALLFITM(DATA(MX),IJ(MX),DIST(MX),DHAT(MX),           WW(MX),    5MS  550
     2 MZ, SDSWIT, LFITSW  )                                            5MS  551
       IF(IABS(SDSWIT).LT.10)                                           5MS  552
     1 CALLFITP(DATA(MX),IJ(MX),DIST(MX),DHAT(MX),           WW(MX),    5MS  553
     2 MZ, SDSWIT, LFITRM  )                                            5MS  554
C                                                                       5MS  555
C      CALCULATE U, T, AND STRESS                                       5MS  556
C                                                                       5MS  557
 3200  U=0.0                                                            5MS  558
       T=0.0                                                            5MS  559
       DO 3210 M=MX,MY                                                  5MS  560
      U=U+(DIST(M)-DHAT(M))**2*WW(M)                                    5MS  561
 3210 T=T+(DIST(M)-DBAS)**2*WW(M)                                       5MS  562
 3215  U=SQRT (U)                                                       5MS  563
       TEMP1=T                                                          5MS  564
       T=SQRT (T)                                                       5MS  565
       GRSTRS(NG) = U/T                                                 5MS  566
       STRESS = STRESS + GRSTRS(NG)**2                                  5MS  567
       IF(U.EQ.0.0)   GO TO 3340                                        5MS  568
 3220  RUT=1.0/(U*T)                                                    5MS  569
       UOT3=U/(TEMP1*T)                                                 5MS  570
C                                                                       5MS  571
C      CALCULATE THE (NEGATIVE) GRADIENT                                5MS  572
C                                                                       5MS  573
 3300  DO 3330 M = MX,MY                                                5MS  574
       IF(DIST(M).EQ.0.0) GO TO 3330                                    5MS  575
       LTEMP1 = MOD(IJ(M),TWO18)                                        5MS  576
         I = LTEMP1/TWO9                                                5MS  577
       J = MOD(LTEMP1,TWO9)                                             5MS  578
       FACTOR=UOT3*(DIST(M)-DBAS)-RUT*(DIST(M)-DHAT(M))                 5MS  579
       FACTOR = (FACTOR/RM1POW(DIST(M)) ) * WW(M)                       5MS  580
       FACTOR = FACTOR * GRSTRS(NG)                                     5MS  581
       DO 3320 L=1,LDIM                                                 5MS  582
       TEMP1 = FACTOR * RM1POW(X(I,L)-X(J,L))                           5MS  583
       GR(I,L)=GR(I,L)+TEMP1                                            5MS  584
 3320  GR(J,L)=GR(J,L)-TEMP1                                            5MS  585
 3330  CONTINUE                                                         5MS  586
 3340  CONTINUE                                                         5MS  587
       IF(STRESS .EQ. 0.0 ) GO TO 3700                                  5MS  588
       STRESS = SQRT( STRESS / FNGRPS )                                 5MS  589
C                                                                       5MS  590
C             AVOID MOVING FIXED POINTS                                 5MS  591
C                                                                       5MS  592
       IF( LNFIXZ .EQ. 0) GO TO 3400                                    5MS  593
       DO 3360 I=1,LNFIXZ                                               5MS  594
       DO 3360 L=1,LDIM                                                 5MS  595
 3360  GR(I,L) = 0.0                                                    5MS  596
C                                                                       5MS  597
C      FIND SCALE FACTOR OF GRADIENT, ANGLE-COSINE BETWEEN GRADIENT     5MS  598
C      AND PREVIOUS GRADIENT.                                           5MS  599
C                                                                       5MS  600
 3400  SFGR=0.0                                                         5MS  601
       CAGRGL=0.0                                                       5MS  602
       DO 3410 I=1,N                                                    5MS  603
       DO 3410 L=1,LDIM                                                 5MS  604
       GR(I,L) = GR(I,L) / STRESS                                       5MS  605
       SFGR=SFGR+GR(I,L)**2                                             5MS  606
 3410  CAGRGL=CAGRGL+GR(I,L)*GL(I,L)                                    5MS  607
       SFGR=SQRT (SFGR/FN)                                              5MS  608
C              IF GRADIENT   0.0, SKIP AHEAD.                           5MS  609
       IF(SFGR) 3420,3700,3420                                          5MS  610
 3420  TEMP1=SFGR*SFGL*FN                                               5MS  611
       CAGRGL=CAGRGL/TEMP1                                              5MS  612
C                                                                       5MS  613
       IF(ITNO.EQ.0  .OR.  NBAKUP.GE.4) GO TO 3500                      5MS  614
       IF(CAGRGL.GT.(-0.95)  .AND.  STRESS/STRLST.LT. 1.2001 ) GOTO 35005MS  615
C                                                                       5MS  616
C      BACK UP ALONG LAST GRADIENT                                      5MS  617
C                                                                       5MS  618
       NBAKUP = NBAKUP + 1                                              5MS  619
       TEMP1 = STEP                                                     5MS  620
       STEP = STEP / 10.0                                               5MS  621
       WRITE  (LPRINT,38)   STRESS, CAGRGL,  STEP                       5MS  622
 38    FORMAT(10X, F7.3, 14X, F7.3, 22X, F8.4 )                         5MS  623
       TEMP1 = (TEMP1 - STEP) / SFGL                                    5MS  624
       TEMP1 = STBAMU * TEMP1                                           5MS  625
       DO 3430 I = 1,N                                                  5MS  626
       DO 3430 L = 1,LDIM                                               5MS  627
       X(I,L) = X(I,L) - TEMP1*GL(I,L)                                  5MS  628
 3430  GR(I,L) = GL(I,L)                                                5MS  629
       SFGR = SFGL                                                      5MS  630
       STRESS = STRLST                                                  5MS  631
       GO TO 3000                                                       5MS  632
C                                                                       5MS  633
C      STEP SIZE CALCULATIONS                                           5MS  634
C                                                                       5MS  635
3500   IF(ITNO) 9999, 3510, 3520                                        5MS  636
 3510  SRAT=0.8                                                         5MS  637
       GO TO 3530                                                       5MS  638
 3520  SRAT=STRESS/STRLST                                               5MS  639
 3530  CALL       NEWSTP( STEP, ITNO, SFGR, STRESS,                     5MS  640
     1 CAGRGL, COSAV, ACSAV, COSAVW, ACSAVW, SRAT, SRATAV )             5MS  641
       NBAKUP = 0                                                       5MS  642
C                                                                       5MS  643
C      PRINT CURRENT STATUS OF COMPUTATION                              5MS  644
C                                                                       5MS  645
 3700  IF(LHIPRT.EQ.2) WRITE (LPRINT,30) ITNO,STRESS,SRAT,SRATAV,       5MS  646
     1 CAGRGL,COSAV,ACSAV,SFGR,STEP                                     5MS  647
   30  FORMAT(I10,F7.3,F7.3,F7.3,F7.3,F7.3,F7.3,F8.4,   F8.4)           5MS  648
C                                                                       5MS  649
C      DECIDE WHETHER TO CONTINUE ITERATING                             5MS  650
C                                                                       5MS  651
 3800  IF(STRESS) 9999, 3840, 3810                                      5MS  652
 3810  IF(SFGR-SFGRMN) 3850, 3850, 3815                                 5MS  653
 3815  TEMP1 = 0.5 * (1.0+SRATST)                                       5MS  654
       TEMP2 = 1.0 - TEMP1                                              5MS  655
       IF( ABS (SRAT-TEMP1) - TEMP2 ) 3816, 3816, 3820                  5MS  656
 3816  IF( ABS (SRATAV-TEMP1) - TEMP2 ) 3850, 3850, 3820                5MS  657
 3820  IF(STRESS-STRMIN) 3860, 3860, 3830                               5MS  658
 3830  IF(ITNO-NOIT) 3900, 3870, 9999                                   5MS  659
 3840  CONTINUE                                                         5MS  660
       WRITE            (LPRINT, 31)                                    5MS  661
 31    FORMAT(24H0ZERO STRESS WAS REACHED )                             5MS  662
       GO TO 4000                                                       5MS  663
 3850  CONTINUE                                                         5MS  664
       WRITE            (LPRINT, 32)                                    5MS  665
 32    FORMAT(21H0MINIMUM WAS ACHIEVED )                                5MS  666
       GO TO 4000                                                       5MS  667
 3860  CONTINUE                                                         5MS  668
       WRITE            (LPRINT, 33)                                    5MS  669
 33    FORMAT(32H0SATISFACTORY STRESS WAS REACHED )                     5MS  670
       GO TO 4000                                                       5MS  671
 3870  CONTINUE                                                         5MS  672
       WRITE            (LPRINT, 34)                                    5MS  673
 34    FORMAT(39H0MAXIMUM NUMBER OF ITERATIONS WERE USED )              5MS  674
       GO TO 4000                                                       5MS  675
C                                                                       5MS  676
C      CONTINUE ITERATING                                               5MS  677
C                                                                       5MS  678
 3900  ITNO=ITNO+1                                                      5MS  679
       TEMP1=STEP/SFGR                                                  5MS  680
       DO 3910 I=1,N                                                    5MS  681
       DO 3910 L=1,LDIM                                                 5MS  682
 3910  X(I,L)=X(I,L)+TEMP1*GR(I,L)                                      5MS  683
       GO TO 3000                                                       5MS  684
C                                                                       5MS  685
C      STOP ITERATING              *************************************5MS  686
C                                                                       5MS  687
 4000  WRITE            (LPRINT, 40)N,LDIM,STRESS,LSCH,(L,L=1,LDIM)     5MS  688
 40    FORMAT(27H0THE FINAL CONFIGURATION OF,I4,                        5MS  689
     1 10H POINTS IN,I3, 22H DIMENSIONS HAS STRESS,F7.3,8H FORMULA ,I2  5MS  690
     2 /1X/20H FINAL CONFIGURATION/10I7)                                5MS  691
       IF(LPUNSW.EQ.2) GO TO 4005                                       5MS  692
      WRITE    (LPUNCH,41) (TITLE(I),I=1,80),N,LDIM                     5MS  693
 41   FORMAT (14H CONFIGURATION/80A1/2I3)                               5MS  694
      WRITE     (LPUNCH,52)                                             5MS  695
 52   FORMAT (12H (2X,10F7.3))                                          5MS  696
 4005  DO 4010 I=1,N                                                    5MS  697
       WRITE            (LPRINT, 42)I,(X(I,L),L=1,LDIM)                 5MS  698
       IF(LPUNSW.EQ.2) GO TO 4010                                       5MS  699
       WRITE            (LPUNCH, 43)I,(X(I,L),L=1,LDIM)                 5MS  700
 4010  CONTINUE                                                         5MS  701
   42  FORMAT(1X,I2,10F7.3)                                             5MS  702
   43  FORMAT(I2,10F7.3)                                                5MS  703
       WRITE  (LPRINT,46)                                               5MS  704
       DO 4020 NG=1,NOGRPS                                              5MS  705
       MZ = GRNO(NG+1) - GRNO(NG)                                       5MS  706
      SDSWT1=MOD(GRSDIS(NG),100) +1                                     5MS  707
      IF(SDSWT1-11) 4016,4013,4014                                      5MS  708
 4013 WRITE(LPRINT,60) NG,MZ,GRSTRS(NG)                                 5MS  709
 60   FORMAT(1X,I5,2X,I5,F7.3,11H  ASCENDING)                           5MS  710
      GO TO 4020                                                        5MS  711
 4014 WRITE(LPRINT,62) NG,MZ,GRSTRS(NG)                                 5MS  712
 62   FORMAT(1X,I5,2X,I5,F7.3,11H DESCENDING)                           5MS  713
      GO TO 4020                                                        5MS  714
 4016 WRITE(LPRINT,61) NG,MZ,GRSTRS(NG),(PCOEFF(I),I=1,SDSWT1)          5MS  715
 61   FORMAT(1X,I5,2X,I5,F7.3,12H  POLYNOMIAL,5F15.7)                   5MS  716
 4020  CONTINUE                                                         5MS  717
 46   FORMAT(14H0DATA GROUP(S) ,/73H0SERIAL  COUNT STRESS REGRESSION COE5MS  718
     .FFICIENTS (FROM DEGREE 0 TO MAX OF 4))                            5MS  719
 4030  CONTINUE                                                         5MS  720
       IF(LDIPRT.EQ.2)                                                  5MS  721
     1 CALL DATAPR(DATA,IJ,WW,GRNO,MM,NOGRPS,DIST,DHAT,  2  )           5MS  722
C                                                                       5MS  723
      WRITE  (LPRINT,51)                                                5MS  724
 51   FORMAT(10H0*********)                                             5MS  725
       WRITE            (LPRINT, 44)   LDIM  ,STRESS,TITLE              5MS  726
 44   FORMAT(1H1,10X,29HDIST AND DHAT VERSUS DATA FOR,I4,               5MS  727
     2     22H  DIMENSIONS. STRESS =,F8.4/11X,80A1)                     5MS  728
C                                                                       5MS  729
C      PLOT SHEPARD DIAGRAM AND TWO SPACE CONFIGURATION                 5MS  730
C                                                                       5MS  731
C      CHANGE DIMENSION                                                 5MS  732
C                                                                       5MS  733
       CALL PLOTR(DIST,DATA,DHAT,MM,LPRINT,+1)                          5MS  734
       IF( NOGRPS .EQ. 1 ) GO TO 4065                                   5MS  735
C                                                                       5MS  736
C      PLOT ADDITIONAL SHEPARD DIAGRAMS BY GROUPS -- LIMIT OF 5         5MS  737
C                                                                       5MS  738
       KPLT = NOGRPS                                                    5MS  739
       IF( NOGRPS .GT. 5 ) KPLT = 5                                     5MS  740
       DO 4060 NG=1,KPLT                                                5MS  741
       MX = GRNO(NG)                                                    5MS  742
       MY = GRNO(NG+ 1) - 1                                             5MS  743
       MZ = MY - MX + 1                                                 5MS  744
       WRITE(LPRINT, 4062) NG, GRSTRS(NG), LDIM, TITLE                  5MS  745
 4062 FORMAT (42H1 DIST AND DHAT VERSUS DATA FOR GROUP NO. , I3,3X,     5MS  746
     1 8HSTRESS =,F6.4,3X, 11HDIMENSION =,I3,/,10X,80A1)                5MS  747
       CALL PLOTR( DIST(MX), DATA(MX), DHAT(MX), MZ, LPRINT, +1 )       5MS  748
 4060  CONTINUE                                                         5MS  749
 4065  CONTINUE                                                         5MS  750
      IF(LDIM.EQ.2)  WRITE(LPRINT,48)     TITLE,STRESS                  5MS  751
 48   FORMAT (30H1 TWO SPACE CONFIGURATION FOR ,80A1, 9H STRESS =,      5MS  752
     1      F7.4///)                                                    5MS  753
      IF(LDIM.EQ.2) CALL PLOT1(X(1,1),X(1,2),N,+1,+1,+4.,-4.,3.,-3.)    5MS  754
      ISTC=ISTC+1                                                       5MS  755
      STPL(ISTC)= STRESS                                                5MS  756
      DIMSV(ISTC)=LDIM                                                  5MS  757
 4100  LDIM=LDIM-LDIMD                                                  5MS  758
       IF(LDIM-LDIMN)4110,2300,2300                                     5MS  759
 4110 CONTINUE                                                          5MS  760
C                                                                       5MS  761
C     PLOT STRESS VERSUS DIMENSION IF MORE THAN TWO POINTS              5MS  762
C                                                                       5MS  763
      IF(ISTC.LT.3) GO TO 100                                           5MS  764
      WRITE(LPRINT, 50)  TITLE                                          5MS  765
 50   FORMAT (1H1,10X,31HPLOT OF STRESS VERSUS DIMENSION,80A1/)         5MS  766
C  *****  INVERT TO IMPROVE READING OF PLOT  *****                      5MS  767
       NVT=ISTC/2                                                       5MS  768
       DO 8050 I=1,NVT                                                  5MS  769
       IC=ISTC+1-I                                                      5MS  770
       TEMP = DIMSV(I)                                                  5MS  771
       SEMP = STPL(I)                                                   5MS  772
       DIMSV(I) = DIMSV(IC)                                             5MS  773
       STPL(I) = STPL(IC)                                               5MS  774
       DIMSV(IC) = TEMP                                                 5MS  775
       STPL(IC) = SEMP                                                  5MS  776
 8050  CONTINUE                                                         5MS  777
      CALL PLOT1(DIMSV,STPL,ISTC,-1,+1,12.,0.,1.08,0.)                  5MS  778
       WRITE            (LPRINT, 45)                                    5MS  779
 45    FORMAT(50H0************************************************* )   5MS  780
       GO TO 100                                                        5MS  781
C              REINITIALIZE, AND RETURN FOR MORE INPUT.                 5MS  782
C                                                                       5MS  783
C      NORMAL TERMINATION, AFTER READING  STOP  ONCONTROL CARD          5MS  784
C                                                                       5MS  785
 9000  STOP                                                             5MS  786
C                                                                       5MS  787
C      TROUBLE EXIT                                                     5MS  788
C                                                                       5MS  789
 9999  WRITE            (LPRINT, 99)                                    5MS  790
 99    FORMAT(52H0KRUSKAL. IMPOSSIBLE BRANCH TAKEN FROM IF STATEMENT. ) 5MS  791
       STOP                                                             5MS  792
C                                                                       5MS  793
      END                                                               5MS  794
$      FORTRAN
CRROO                                                                   RROO   1
      FUNCTION RROOT(ZZ)                                                RROO   2
C      MDSCAL, VERSION 5MS, OCTOBER, 1971                               RROO   3
C      WRITTEN BY J. KRUSKAL AND J. SEERY, NOVEMBER, 1971               RROO   4
C                                                                       RROO   5
      INTEGER RTYPE                                                     RROO   6
      COMMON /METRIC/ RTYPE,RM1,RECR,R                                  RROO   7
      Z=ZZ                                                              RROO   8
      GO TO (310,320,330), RTYPE                                        RROO   9
 310  RROOT=Z                                                           RROO  10
      RETURN                                                            RROO  11
 320  RROOT=SQRT(Z)                                                     RROO  12
      RETURN                                                            RROO  13
 330  RROOT=Z**RECR                                                     RROO  14
      RETURN                                                            RROO  15
      END                                                               RROO  16
$      FORTRAN
CRM1P                                                                   RM1P   1
      FUNCTION RM1POW(YY)                                               RM1P   2
C      MDSCAL, VERSION 5MS, OCTOBER, 1971                               RM1P   3
C      WRITTEN BY J. KRUSKAL AND J. SEERY, NOVEMBER, 1971               RM1P   4
C                                                                       RM1P   5
      INTEGER RTYPE                                                     RM1P   6
      COMMON /METRIC/ RTYPE,RM1,RECR,R                                  RM1P   7
      Y=YY                                                              RM1P   8
      GO TO (210,220,230), RTYPE                                        RM1P   9
 210  RM1POW=0.0                                                        RM1P  10
      IF(Y.NE.0.0) RM1POW=SIGN(1.0,Y)                                   RM1P  11
      RETURN                                                            RM1P  12
 220  RM1POW=Y                                                          RM1P  13
      RETURN                                                            RM1P  14
 230  RM1POW=SIGN(ABS(Y)**RM1,Y)                                        RM1P  15
      RETURN                                                            RM1P  16
      END                                                               RM1P  17
$      FORTRAN
CRPOW                                                                   RPOW   1
      FUNCTION RPOWER(XX)                                               RPOW   2
C      MDSCAL, VERSION 5MS, OCTOBER, 1971                               RPOW   3
C      WRITTEN BY J. KRUSKAL AND J. SEERY, NOVEMBER, 1971               RPOW   4
C                                                                       RPOW   5
      INTEGER RTYPE                                                     RPOW   6
      COMMON /METRIC/ RTYPE,RM1,RECR,R                                  RPOW   7
      X=XX                                                              RPOW   8
      GO TO (110,120,130), RTYPE                                        RPOW   9
 110  RPOWER=ABS(X)                                                     RPOW  10
      RETURN                                                            RPOW  11
 120  RPOWER=X*X                                                        RPOW  12
      RETURN                                                            RPOW  13
 130  RPOWER=ABS(X)**R                                                  RPOW  14
      RETURN                                                            RPOW  15
      END                                                               RPOW  16
$      FORTRAN
CFITM                                                                   FITM   1
      SUBROUTINEFITM(DATA,IJ,DIST,DHAT,       WW,  MM,  SDSWIT, LFITSW )FITM   2
C      MDSCAL, VERSION 5MS, OCTOBER, 1971                               FITM   3
C      UNCHANGED FROM VERSION 4, JANUARY,1968                           FITM   4
C                                                                       FITM   5
C      FITM    PERFORMS WEIGHTED-LEAST-SQUARES MONOTONE REGRESSION      FITM   6
C      THIS ROUTINE FINDS THE VALUES FOR  DHAT  WHICH ARE MONOTONIC     FITM   7
C      AND WHICH BEST APPROXIMATE  THE VALUES OF  DIST,                 FITM   8
C      IN THE SENSE THAT THE SUM OF THE SQUARED DEVIATIONS,             FITM   9
C      EACH ONE WEIGHTED BY THE CORRESPONDING VALUE IN  WW ,            FITM  10
C      IS A MINIMUM.                                                    FITM  11
C      IT DIRECTLY USES  SORT.                                          FITM  12
C                                                                       FITM  13
       DIMENSION   DATA(1800), IJ(1800), DIST(1800), DHAT(1800)         FITM  14
      DIMENSION LBLOCK(1800),WW(1800)                                   FITM  15
       INTEGER TWO9, TWO18                                              FITM  16
       INTEGER SDSWIT                                                   FITM  17
C                                                                       FITM  18
       COMMON  /MDSCL1/  LREAD, LPRINT, LPUNCH, LSCRAT                  FITM  19
C                                                                       FITM  20
       DATA    TWO9, TWO18  /512,  262144/                              FITM  21
C                                                                       FITM  22
C      FORM FIRST APPROXIMATION TO CORRECT PARTITON                     FITM  23
C                                                                       FITM  24
C              IF LFITSW 1,  USE PRIMARY APPROACH.  THUS WE SORT EACH   FITM  25
C              BLOCK OF EQUAL VALUES OF DATA ACCORDING TO DIST VALUES.  FITM  26
C              THEN WE CREATE PARTITON INTO BLOCKS OF SIZE 1 TO START.  FITM  27
C                                                                       FITM  28
C              IF LFITSW 2,   USE SECONDARY APPROACH. THUS WE START WITHFITM  29
C              PARTITION INTO BLOCKS OF EQUAL DATA VALUES.              FITM  30
C                                                                       FITM  31
C              WITHIN EACH BLOCK  OF WHATEVER SIZE, THE FIRST DHAT VALUEFITM  32
C      GIVES THE WEIGHTED TOTAL OF THE DIST VALUES IN THAT BLOCK, WHILE FITM  33
C              THE LAST DHAT VALUE GIVES THE SUM OF THE WEIGHTS FOR THATFITM  34
C              BLOCK.  SIMILARLY, WITHIN EACH BLOCK, THE FIRST LBLOCK   FITM  35
C              VALUE AND THE LAST LBLOCK VALUE BOTH GIVE THE SIZE OF THEFITM  36
C              BLOCK.                                                   FITM  37
C      BLOCKS OF SIZE 1 FORM A PARTIAL EXCEPTION. EVERYTHIN IS THE SAME FITM  38
C              EXCEPT THAT THE SUM OF THE W VALUES IS NOT FOUND IN THE  FITM  39
C              LAST DHAT VALUE IN THE BLOCK.                            FITM  40
C                                                                       FITM  41
  100  MA=1                                                             FITM  42
 110   K = IJ(MA) / TWO18                                               FITM  43
       MB=MA+K-1                                                        FITM  44
       GO TO ( 200, 300 ), LFITSW                                       FITM  45
C                                                                       FITM  46
C              PRIMARY APPROACH                                         FITM  47
C                                                                       FITM  48
  200  IF(K-1) 9999, 220, 210                                           FITM  49
C              IF K 1, SAVE SORTING TIME                                FITM  50
  210  CALL SORT( DIST(MA), K, IJ(MA), DATA(MA),WW(MA),3,+1)            FITM  51
C                                                                       FITM  52
C               SORT  WILL SORT THE K ELEMENTS OF  DIST  IN ALGEBRAIC   FITM  53
C              ORDER.                                                   FITM  54
C              BECAUSE THE FINAL ARGUMENT IS ZERO, SORTING WILL BE      FITM  55
C              ASCENDING OR DESCENDING ACCORDING TO THE VALUE OF SDSWIT FITM  56
C              SUPPLIED DURING THE CALL FOR SORT IN MDSCAL.             FITM  57
C              THE ELEMENTS OF IJ  AND  DATA  WILL BE REARRANGED IN     FITM  58
C              EXACTLY THE SAME ORDER AS THOSE OF  DIST .               FITM  59
C      ALSO THE ELEMENTS OF  WW .                                       FITM  60
C                                                                       FITM  61
  220  DO 230 M=MA,MB                                                   FITM  62
       DHAT(M) = DIST(M) * WW(M)                                        FITM  63
  230  LBLOCK(M)=1                                                      FITM  64
       GO TO 400                                                        FITM  65
C                                                                       FITM  66
C              SECONDARY APPROACH                                       FITM  67
C                                                                       FITM  68
C                                                                       FITM  69
  300  LBLOCK(MA)=K                                                     FITM  70
       LBLOCK(MB)=K                                                     FITM  71
       TEMP1=0.0                                                        FITM  72
       TEMP2 = 0.0                                                      FITM  73
       DO 310 M=MA,MB                                                   FITM  74
       TEMP1 = TEMP1 + DIST(M) * WW(M)                                  FITM  75
       TEMP2 = TEMP2 + WW(M)                                            FITM  76
 310   CONTINUE                                                         FITM  77
       DHAT(MA)=TEMP1                                                   FITM  78
       IF(K-1) 9999, 400, 320                                           FITM  79
 320   DHAT(MB) = TEMP2                                                 FITM  80
C              PROCEED TO NEXT BOCK. QUERY. IS IT THE LAST              FITM  81
C                                                                       FITM  82
  400  MA=MA+K                                                          FITM  83
       IF(MA-MM-1) 110, 500, 9999                                       FITM  84
C                                                                       FITM  85
C      START MAIN COMPUTATION      *************************************FITM  86
C                                                                       FITM  87
  500  MA=1                                                             FITM  88
  510  LUD=2                                                            FITM  89
       NSATIS=0                                                         FITM  90
  520  K=LBLOCK(MA)                                                     FITM  91
       MB=MA+K-1                                                        FITM  92
       IF(K-1) 9999, 530, 540                                           FITM  93
 530   WT =WW(MB)                                                       FITM  94
       GO TO 550                                                        FITM  95
 540   WT = DHAT(MB)                                                    FITM  96
 550   DAV = DHAT(MA) / WT                                              FITM  97
       GO TO ( 600, 700 ), LUD                                          FITM  98
C                                                                       FITM  99
C              IS BLOCK DOWN-SATISFIED. IF NOT, MERGE                   FITM 100
C                                                                       FITM 101
  600  IF(MA-1) 9999, 630, 610                                          FITM 102
  610  MBD=MA-1                                                         FITM 103
       KD=LBLOCK(MBD)                                                   FITM 104
       MAD=MBD-KD+1                                                     FITM 105
       IF(KD-1) 9999, 613, 615                                          FITM 106
 613   WTD =WW(MBD)                                                     FITM 107
       GO TO 617                                                        FITM 108
 615   WTD = DHAT(MBD)                                                  FITM 109
 617   DAVD = DHAT(MAD) / WTD                                           FITM 110
       IF( DAVD-DAV ) 630, 620, 620                                     FITM 111
  620  KNEW=K+KD                                                        FITM 112
       LBLOCK(MAD)=KNEW                                                 FITM 113
       LBLOCK(MB)=KNEW                                                  FITM 114
       DTONEW = DHAT(MAD) + DHAT(MA)                                    FITM 115
       DHAT(MAD) = DTONEW                                               FITM 116
       DHAT(MB) = WT + WTD                                              FITM 117
       NSATIS=0                                                         FITM 118
       MA=MAD                                                           FITM 119
       GO TO 800                                                        FITM 120
  630  NSATIS=NSATIS+1                                                  FITM 121
       GO TO 800                                                        FITM 122
C                                                                       FITM 123
C              IS BLOCK UP-SATISFIED. IF NOT, MERGE                     FITM 124
C                                                                       FITM 125
  700  IF(MB-MM) 710, 730, 9999                                         FITM 126
  710  MAU=MB+1                                                         FITM 127
       KU=LBLOCK(MAU)                                                   FITM 128
       MBU=MAU+KU-1                                                     FITM 129
       IF(KU-1) 9999, 713, 715                                          FITM 130
 713   WTU =WW(MBU)                                                     FITM 131
       GO TO 717                                                        FITM 132
 715   WTU = DHAT(MBU)                                                  FITM 133
 717   DAVU = DHAT(MAU) / WTU                                           FITM 134
       IF( DAV-DAVU ) 730, 720, 720                                     FITM 135
  720  KNEW=K+KU                                                        FITM 136
       LBLOCK(MA)=KNEW                                                  FITM 137
       LBLOCK(MBU)=KNEW                                                 FITM 138
       DTONEW = DHAT(MA) + DHAT(MAU)                                    FITM 139
       DHAT(MA) = DTONEW                                                FITM 140
       DHAT(MBU) = WT + WTU                                             FITM 141
       NSATIS=0                                                         FITM 142
       GO TO 800                                                        FITM 143
  730  NSATIS=NSATIS+1                                                  FITM 144
       GO TO 800                                                        FITM 145
C                                                                       FITM 146
C              PROCEED TO NEXT BLOCK IF READY.                          FITM 147
C                                                                       FITM 148
  800  LUD = 3-LUD                                                      FITM 149
C              QUERY. IS BLOCK BOTH UP AND DOWN SATISFIED. IF NOT, RETURFITM 150
       IF(NSATIS-1) 520, 520, 810                                       FITM 151
C              QUERY. IS THIS LAST BLOCK. IF NOT, GO ON TO NEXT BLOCK.  FITM 152
  810  IF(MB-MM) 820, 900, 9999                                         FITM 153
  820  MA=MB+1                                                          FITM 154
       GO TO 510                                                        FITM 155
C                                                                       FITM 156
C      MAIN COMPUTATION COMPLETE. PLACE ANSWERS IN DHAT.                FITM 157
C                                                                       FITM 158
  900  MA=1                                                             FITM 159
  910  K=LBLOCK(MA)                                                     FITM 160
       MB=MA+K-1                                                        FITM 161
       IF(K-1) 9999, 940, 920                                           FITM 162
 920   TEMP1 = DHAT(MA) / DHAT(MB)                                      FITM 163
       DO 930 M=MA,MB                                                   FITM 164
  930  DHAT(M)=TEMP1                                                    FITM 165
       GO TO 945                                                        FITM 166
 940   DHAT(MA) = DIST(MA)                                              FITM 167
 945   MA = MB + 1                                                      FITM 168
       IF(MA-MM-1) 910, 950, 9999                                       FITM 169
  950  RETURN                                                           FITM 170
C                                                                       FITM 171
C      TROUBLE EXIT                                                     FITM 172
C                                                                       FITM 173
 9999  WRITE            (LPRINT, 99)                                    FITM 174
   99  FORMAT(50H0KRUSKAL. IMPOSSIBLE BRANCH TAKEN ON IF STATEMENT. )   FITM 175
       STOP                                                             FITM 176
C                                                                       FITM 177
      END                                                               FITM 178
$      FORTRAN
CEQSV                                                                   EQSV   1
       SUBROUTINE EQSOLV(A,B,C,NT,NF)                                   EQSV   2
C      MDSCAL, VERSION 5MS, OCTOBER, 1971                               EQSV   3
C      UNCHANGED FROM VERSION 4, JANUARY,1968                           EQSV   4
C                                                                       EQSV   5
CEQSOLV        SOLVES SIMULTANEOUS LINEAR EQUATIONS                     EQSV   6
       REAL    A(5,5), B(5), C(5)                                       EQSV   7
C                                                                       EQSV   8
       NTM1=NT-1                                                        EQSV   9
       IF(NT.EQ.NF) GO TO 60                                            EQSV  10
       DO 50 II=NF,NTM1                                                 EQSV  11
       IIP1=II+1                                                        EQSV  12
       DO 40 I=IIP1,NT                                                  EQSV  13
       E = A(I,II)/A(II,II)                                             EQSV  14
       DO 30 J=NF,NT                                                    EQSV  15
       A(I,J)=A(I,J)-E*A(II,J)                                          EQSV  16
 30    CONTINUE                                                         EQSV  17
       B(I)=B(I)-E*B(II)                                                EQSV  18
 40    CONTINUE                                                         EQSV  19
 50    CONTINUE                                                         EQSV  20
C                                                                       EQSV  21
 60    DO 90 IC=NF,NT                                                   EQSV  22
       I=NT+NF-IC                                                       EQSV  23
       E=B(I)                                                           EQSV  24
       IF(I.EQ.NT) GO TO 80                                             EQSV  25
       IP1=I+1                                                          EQSV  26
       DO 70 J=IP1,NT                                                   EQSV  27
       E=E-C(J)*A(I,J)                                                  EQSV  28
 70    CONTINUE                                                         EQSV  29
 80    C(I) = E / A(I,I)                                                EQSV  30
 90    CONTINUE                                                         EQSV  31
       RETURN                                                           EQSV  32
       END                                                              EQSV  33
$      FORTRAN
CRUFV                                                                   RUFV   1
       FUNCTION RUNIFV(A)                                               RUFV   2
C      MDSCAL, VERSION 5MS, OCTOBER, 1971                               RUFV   3
C      UNCHANGED FROM VERSION 4, JANUARY,1968                           RUFV   4
CRUNIFV        RANDOM UNIFROM VARIABLE FOR MDSCAL                       RUFV   5
C                                                                       RUFV   6
C      THIS IS A SIMPLE ROUTINE FOR GENERATING RANDOM NUMBERS.          RUFV   7
C      THEY ARE UNIFORMLY DISTRIBUTED BETWEEN 0 AND 1.                  RUFV   8
C      IT IS NOT FAST, NOR DOES IT PRODUCE VERY HIGH QUALITY NUMBERS.   RUFV   9
C      ITS MAIN VIRTUE IS THAT IT IS IN FORTRAN AND WILL WORK ON        RUFV  10
C      ALMOST ANY MACHINE ON WHICH FORTRAN IS AVAILABLE.                RUFV  11
C                                                                       RUFV  12
       DATA MODULO,FLMOD,K/8192,8192.0,1/                               RUFV  13
C      MODULO   2**13                                                   RUFV  14
C                                                                       RUFV  15
       DO 10 I=1,15                                                     RUFV  16
 10    K = MOD(5*K, MODULO)                                             RUFV  17
       Z = FLOAT(K)/FLMOD                                               RUFV  18
       RUNIFV = Z                                                       RUFV  19
       RETURN                                                           RUFV  20
C                                                                       RUFV  21
       END                                                              RUFV  22
$      FORTRAN
CREGR                                                                   REGR   1
       FUNCTION REGR(DA,I)                                              REGR   2
C      MDSCAL, VERSION 5MS, OCTOBER, 1971                               REGR   3
C      UNCHANGED FROM VERSION 4, JANUARY,1968                           REGR   4
CREGR          CALCULATES VALUES OF FUNCTIONS REGRESSED OVER            REGR   5
C                                                                       REGR   6
       J=MIN0(I,4)                                                      REGR   7
       GO TO (10,20,30,40), J                                           REGR   8
C                                                                       REGR   9
 10    REGR=1.0                                                         REGR  10
       RETURN                                                           REGR  11
C                                                                       REGR  12
 20    REGR=DA                                                          REGR  13
       RETURN                                                           REGR  14
C                                                                       REGR  15
 30    REGR=DA*DA                                                       REGR  16
       RETURN                                                           REGR  17
C                                                                       REGR  18
 40    REGR=DA**(I-1)                                                   REGR  19
       RETURN                                                           REGR  20
C                                                                       REGR  21
       END                                                              REGR  22
$      FORTRAN
CNEWS                                                                   NEWS   1
       SUBROUTINE NEWSTP( STEP, ITNO, SFGR, STRESS,                     NEWS   2
     1 CAGRGL, COSAV, ACSAV, COSAVW, ACSAVW, SRAT, SRATAV )             NEWS   3
C      MDSCAL, VERSION 5MS, OCTOBER, 1971                               NEWS   4
C      UNCHANGED FROM VERSION 4, JANUARY,1968                           NEWS   5
CNEWSTP            NEWSTP  FOR MDSCAL                                   NEWS   6
C                                                                       NEWS   7
C      NEWSTP  THIS SUBROUTINE COMPUTES THE STEP SIZE.                  NEWS   8
C                                                                       NEWS   9
C              THE MAIN PURPOSE OF THIS ROUTINE IS TO COMPUTE THE NEW   NEWS  10
C              VALUE OF  STEP .                                         NEWS  11
C              INCIDENTALLY, IT UPDATES  COSAV ,  ACSAV , AND  SRATAV . NEWS  12
C                                                                       NEWS  13
C      UPDATE THREE AVERAGE QUANTITIES                                  NEWS  14
C                                                                       NEWS  15
       COSAV = CAGRGL*COSAVW  +  COSAV*(1.0-COSAVW)                     NEWS  16
       ACSAV = ABS (CAGRGL)*ACSAVW  +  ACSAV*(1.0-ACSAVW)               NEWS  17
       SRATAV = (SRAT**0.33334)  *  (SRATAV**0.66666)                   NEWS  18
       IF(ITNO) 100, 100, 200                                           NEWS  19
C                                                                       NEWS  20
C              GUESS INITIAL STEP SIZE                                  NEWS  21
C                                                                       NEWS  22
 100   STEP= (25.0*STRESS) * SFGR                                       NEWS  23
       RETURN                                                           NEWS  24
C                                                                       NEWS  25
C              FIND NEW STEP SIZE                                       NEWS  26
C                                                                       NEWS  27
 200   ANG=4.0**COSAV                                                   NEWS  28
       TEMP1 = 1.0 + (AMIN1 (1.0,SRATAV) ) ** 5                         NEWS  29
       TEMP2 = 1.0 + ( ACSAV - ABS (COSAV) )                            NEWS  30
       BIAS = 1.6 / (TEMP1*TEMP2)                                       NEWS  31
       GOODLK = SQRT (AMIN1 (1.0,SRAT) )                                NEWS  32
       STEP = STEP * ANG * BIAS * GOODLK                                NEWS  33
       RETURN                                                           NEWS  34
       END                                                              NEWS  35
$      FORTRAN
CSORT                                                                   SORT   1
       SUBROUTINE SORT (A, N, B, C, D, K, SWITCH )                      SORT   2
C      MDSCAL, VERSION 5MS, OCTOBER, 1971                               SORT   3
C      UNCHANGED FROM VERSION 4, JANUARY,1968                           SORT   4
C                                                                       SORT   5
C     THIS ROUTINE SORTS INPUT ARRAY 'A' AND REARRANGES, OPTIONALLY,    SORT   6
C     ARRAYS 'B', 'C', AND 'D', IN ORDER CORRESPONDING TO 'A'.          SORT   7
C     N = NUMBER OF ITEMS IN 'A' (AND 'B', 'C', 'D', IF USED)           SORT   8
C     K = 0--SORT 'A' ONLY, 1--REARRANGE 'B', 2--REARRANGE 'B' AND 'C', SORT   9
C         3--REARRANGE 'B', 'C', AND 'D'.                               SORT  10
C     IF 'SWITCH' IS POSITIVE, SORT WILL BE IN ASCENDING ORDER,         SORT  11
C                 IF ZERO OR NEGATIVE, IN DESCENDING ORDER.             SORT  12
C     ALGORITHM FROM CACM, JULY 1959, PAGE 30 BY D. L. SHELL            SORT  13
C                                                                       SORT  14
       DIMENSION A(1800), B(1800), C(1800), D(1800)                     SORT  15
       INTEGER  SWITCH                                                  SORT  16
 105   KP1=K+1                                                          SORT  17
       IF(N.LE.1) GO TO 999                                             SORT  18
      M = 1                                                             SORT  19
  106  M = M + M                                                        SORT  20
       IF( M .LE. N ) GO TO 106                                         SORT  21
       M = M - 1                                                        SORT  22
 994       M = M/2                                                      SORT  23
       IF(M.EQ.0) GO TO 999                                             SORT  24
       KK = N-M                                                         SORT  25
       J = 1                                                            SORT  26
992     I = J                                                           SORT  27
996     IM = I + M                                                      SORT  28
       IF(SWITCH)  810,810,800                                          SORT  29
800     IF(A(I).GT.A(IM)) GO TO 110                                     SORT  30
       GO TO 995                                                        SORT  31
810    IF(A(I).LT.A(IM))  GO TO 110                                     SORT  32
995     J = J+1                                                         SORT  33
       IF(J.GT.KK)  GO TO 994                                           SORT  34
       GO TO 992                                                        SORT  35
 110   TEMP=A(I)                                                        SORT  36
       A(I) = A(IM)                                                     SORT  37
       A(IM) = TEMP                                                     SORT  38
       GO TO ( 140, 130, 120, 115), KP1                                 SORT  39
 115   TEMP = D(I)                                                      SORT  40
       D(I) = D(IM)                                                     SORT  41
        D(IM) = TEMP                                                    SORT  42
 120   TEMP=C(I)                                                        SORT  43
       C(I) = C(IM)                                                     SORT  44
       C(IM) = TEMP                                                     SORT  45
 130   TEMP=B(I)                                                        SORT  46
       B(I) = B(IM)                                                     SORT  47
       B(IM) = TEMP                                                     SORT  48
140       I = I-M                                                       SORT  49
       IF(I.LT.1)  GO TO 995                                            SORT  50
        GO TO 996                                                       SORT  51
999      RETURN                                                         SORT  52
      END                                                               SORT  53
$      FORTRAN
CDAPR                                                                   DAPR   1
C      MDSCAL, VERSION 4, JANUARY 1968, BELL TEL LABS                   DAPR   2
C      MDSCAL, VERSION 5MS, OCTOBER, 1971                               DAPR   3
CDATAPR ADAPTED BY J.KRUSKAL AND J.SEERY, JUNE, 1971                    DAPR   4
CDATAPR        DATA PRINT FOR MDSCAL                                    DAPR   5
C                                                                       DAPR   6
       SUBROUTINE DATAPR( DATA, IJ, WW,GRNO, MM, NOGRPS, DIST, DHAT, LCODAPR   7
     1 DE)                                                              DAPR   8
C                                                                       DAPR   9
       REAL    DATA(1), WW(1), DIST(1), DHAT(1),P(50)                   DAPR  10
C                                                                       DAPR  11
       INTEGER  IJ(1),GRNO(1),Q(50)                                     DAPR  12
       INTEGER TWO9, TWO18                                              DAPR  13
C                                                                       DAPR  14
       EQUIVALENCE (P,Q)                                                DAPR  15
C                                                                       DAPR  16
       COMMON   /MDSCL1/ LREAD, LPRINT, LPUNCH, LSCRAT                  DAPR  17
C                                                                       DAPR  18
C                                                                       DAPR  19
 1     FORMAT( 1H1, 4(30H   I   J   DATA   WGHT GP  NO  ) )             DAPR  20
 2     FORMAT(1H1, 3(40H   I   J   DATA   DIST DHAT WGHT GP  NO ) )     DAPR  21
 3     FORMAT( 1X, 5( 2I4, F9.3,        F5.2, I3, I4, 1H, ) )           DAPR  22
 4     FORMAT( 1X, 5( 2I4, F9.3, 2F5.2, F5.2, I3, I4, 1H, ) )           DAPR  23
C                                                                       DAPR  24
       DATA    TWO9, TWO18  /512, 262144 /                              DAPR  25
C                                                                       DAPR  26
 100   NCOPA=3                                                          DAPR  27
       IF(LCODE.EQ.1) NCOPA=4                                           DAPR  28
C                                                                       DAPR  29
       MMLEFT=MM                                                        DAPR  30
 200   MMNOW=MIN0( MMLEFT, 50*NCOPA )                                   DAPR  31
       MMUSED=MM-MMLEFT                                                 DAPR  32
       MMLEFT=MMLEFT-MMNOW                                              DAPR  33
       IF(LCODE.EQ.1) WRITE (LPRINT,1)                                  DAPR  34
       IF(LCODE.EQ.2) WRITE (LPRINT,2)                                  DAPR  35
C                                                                       DAPR  36
       DO 300 LINENO=1,50                                               DAPR  37
       NCONOW=(MMNOW-LINENO+50)/50                                      DAPR  38
       IF(NCONOW.LE.0) GO TO 310                                        DAPR  39
C                                                                       DAPR  40
       L=0                                                              DAPR  41
       DO 270 LC=1,NCONOW                                               DAPR  42
       M = MMUSED + 50*(LC-1) + LINENO                                  DAPR  43
C                                                                       DAPR  44
       IJM = MOD( IJ(M), TWO18 )                                        DAPR  45
       I=IJM/TWO9                                                       DAPR  46
       J = MOD( IJM, TWO9 )                                             DAPR  47
C                                                                       DAPR  48
       L=L+1                                                            DAPR  49
       Q(L)=I                                                           DAPR  50
       L=L+1                                                            DAPR  51
       Q(L)=J                                                           DAPR  52
C                                                                       DAPR  53
       L=L+1                                                            DAPR  54
       P(L)=DATA(M)                                                     DAPR  55
C                                                                       DAPR  56
       IF(LCODE.EQ.1) GO TO 250                                         DAPR  57
C                                                                       DAPR  58
       L=L+1                                                            DAPR  59
       P(L)=DIST(M)                                                     DAPR  60
       L=L+1                                                            DAPR  61
       P(L)=DHAT(M)                                                     DAPR  62
C                                                                       DAPR  63
 250   CONTINUE                                                         DAPR  64
C                                                                       DAPR  65
       L=L+1                                                            DAPR  66
       P(L)=WW(M)                                                       DAPR  67
C                                                                       DAPR  68
       DO 220 NG=1,NOGRPS                                               DAPR  69
       IF(M.GE.GRNO(NG)) NGM=NG                                         DAPR  70
 220   CONTINUE                                                         DAPR  71
       L=L+1                                                            DAPR  72
       Q(L)=NGM                                                         DAPR  73
C                                                                       DAPR  74
       L=L+1                                                            DAPR  75
       Q(L)=M-GRNO(NGM)+1                                               DAPR  76
C                                                                       DAPR  77
 270   CONTINUE                                                         DAPR  78
C                                                                       DAPR  79
       LL=L                                                             DAPR  80
       IF(LCODE.EQ.1) WRITE (LPRINT,3) (P(L),L=1,LL)                    DAPR  81
       IF(LCODE.EQ.2) WRITE (LPRINT,4) (P(L),L=1,LL)                    DAPR  82
 300   CONTINUE                                                         DAPR  83
C                                                                       DAPR  84
 310   IF(MMLEFT.GT.0) GO TO 200                                        DAPR  85
C                                                                       DAPR  86
       RETURN                                                           DAPR  87
C                                                                       DAPR  88
       END                                                              DAPR  89
$      FORTRAN
CPLTR                                                                   PLTR   1
      SUBROUTINE PLOTR(X,Y,YP,NPOI,OUT,ID)                              PLTR   2
C      MDSCAL, VERSION 5MS, OCTOBER, 1971                               PLTR   3
CPLOTR  ADAPTED BY J. KRUSKAL AND J. SEERY     JUNE, 1971               PLTR   4
C                                                                       PLTR   5
C ROUTINE TO GENERATE A ONE PAGE PLOT OF ARRAY -X- VS. -Y-.             PLTR   6
C                                                                       PLTR   7
C THE PARAMETERS -XMAX- -XMIN- YMAX- - YMIN- INDICATE THE UPPER         PLTR   8
C AND LOWER BOUNDS FOR EACH AXIS. IF XMAX=XMIN THE ROUTINE GENERATES    PLTR   9
C ITS OWN BOUNDS FOR THE X AXIS, AND SIMILARLY IF YMAX=YMIN.            PLTR  10
C                                                                       PLTR  11
C IT IS ASSUMED THAT THE ENTRIES HAVE -NPOI- ENTRIES                    PLTR  12
C                                                                       PLTR  13
C THE PLOTTING IS DONE ON TAPE -OUT-.                                   PLTR  14
C                                                                       PLTR  15
C IF -ID- IS NEGATIVE, AXES WILL BE INCLUDED ON THE PLOT.               PLTR  16
C IF -ID- IS POSITIVE, NO AXES WILL APPEAR.                             PLTR  17
C                                                                       PLTR  18
C -LONG- IS THE LENGTH OF THE ARRAYS -X- AND -Y- IN THE CALLING PROGRAM.PLTR  19
C                                                                       PLTR  20
C   X    PLOTTED ON X AXIS                                              PLTR  21
C   Y    PLOTTED ON Y AXIS                                              PLTR  22
C  YP    PLOTTED ON X AXIS                                              PLTR  23
C                                                                       PLTR  24
C WRITTEN BY FORREST W. YOUNG, NOVEMBER, 1965                           PLTR  25
C      VERSION 3, APRIL 1967                                            PLTR  26
C                                                                       PLTR  27
C--------ADAPTED BY F. J. CARMONE,JR.  5/1/67                           PLTR  28
C                                                                       PLTR  29
C                                                                       PLTR  30
      DIMENSION X(NPOI),Y(NPOI),ITEM(55,121),SMALL(25),HOLL(11),YP(NPOI)PLTR  31
      REAL LABELS(55),ITEM,CHAR(50)                                     PLTR  32
      COMMON  /PLTSS/ SMALL,HOLL,CHAR,LABELS,ITEM,DD,BLANK,AIE,AMINUS   PLTR  33
      INTEGER OUT                                                       PLTR  34
      DATA AAB/1HX/, DOO/1H0/, PL/1H+/, PLUS/1H*/                       PLTR  35
 3001 FORMAT(14X,103H*.****.****.****.****.****.****.****.****.****.****PLTR  36
     1.****.****.****.****.****.****.****.****.****.****.*)             PLTR  37
 3300 FORMAT(2H  ,A1,F9.2,1X,105A1,F11.2)                               PLTR  38
 3301 FORMAT(15X,A1,10(F9.4,A1))                                        PLTR  39
 3302 FORMAT(10X,11F10.4)                                               PLTR  40
      DO 115 I=1,55                                                     PLTR  41
      DO 115 J=1,101                                                    PLTR  42
  115 ITEM(I,J)=BLANK                                                   PLTR  43
  116 CONTINUE                                                          PLTR  44
      IF(ID.GT.0)GO TO 119                                              PLTR  45
      DO 117 I=1,55                                                     PLTR  46
  117 ITEM(I,51)=AIE                                                    PLTR  47
      DO 118 I=1,101                                                    PLTR  48
  118 ITEM(28,I)=AMINUS                                                 PLTR  49
  119 CONTINUE                                                          PLTR  50
      XMAX=X(1)                                                         PLTR  51
      XMIN =X(2)                                                        PLTR  52
      DO 121 I=1,NPOI                                                   PLTR  53
      IF(X(I).GT.XMAX)XMAX=X(I)                                         PLTR  54
      IF(YP(I).GT.XMAX)XMAX=YP(I)                                       PLTR  55
      IF(X(I).LT.XMIN)XMIN=X(I)                                         PLTR  56
      IF(YP(I).LT.XMIN)XMIN=YP(I)                                       PLTR  57
  121 CONTINUE                                                          PLTR  58
      YMAX=Y(1)                                                         PLTR  59
      YMIN=Y(2)                                                         PLTR  60
      DO 123 I=1,NPOI                                                   PLTR  61
      IF(Y(I).GT.YMAX)YMAX=Y(I)                                         PLTR  62
      IF(Y(I).LT.YMIN)YMIN=Y(I)                                         PLTR  63
  123 CONTINUE                                                          PLTR  64
      SPANX=XMAX-XMIN                                                   PLTR  65
      SPANY=YMAX-YMIN                                                   PLTR  66
      XMAX=XMAX+.05*SPANX                                               PLTR  67
      YMAX=YMAX+.05*SPANY                                               PLTR  68
      XMIN=XMIN-.05*SPANX                                               PLTR  69
      YMIN=YMIN-.05*SPANY                                               PLTR  70
      SPANX=XMAX-XMIN                                                   PLTR  71
      SPANY=YMAX-YMIN                                                   PLTR  72
      DELX=SPANX/100.0                                                  PLTR  73
      DELY=SPANY/54.0                                                   PLTR  74
      VALUE=YMAX+DELY                                                   PLTR  75
      SMALL(1)=XMIN                                                     PLTR  76
      DO 120 I=1,20                                                     PLTR  77
  120 SMALL(I+1)=SMALL(I)+5.0*DELX                                      PLTR  78
      DO 135 II=1,NPOI                                                  PLTR  79
      I=(YMAX-Y(II))/DELY+1.5                                           PLTR  80
      J=(X(II)-XMIN)/DELX+1.5                                           PLTR  81
      JP=(YP(II)-XMIN)/DELX+1.5                                         PLTR  82
      IF(I.GT.56.OR.I.LT.1.OR.J.LT.1.OR.J.GT.101)GO TO 135              PLTR  83
      IF(JP.LT.1.OR.JP.GT.101)GO TO 135                                 PLTR  84
      IF(ITEM(I,J).NE.BLANK)GO TO 200                                   PLTR  85
      ITEM( I,J)=DOO                                                    PLTR  86
      GO TO 201                                                         PLTR  87
  200 ITEM(I,J)=PL                                                      PLTR  88
  201 CONTINUE                                                          PLTR  89
      ITEM(I,JP)=AAB                                                    PLTR  90
      GO TO 135                                                         PLTR  91
C 202 CONTINUE                        (OMITTED FOR GE VERSION)          PLTR  92
  135 CONTINUE                                                          PLTR  93
  140  FORMAT(/)                                                        PLTR  94
      WRITE(OUT,3301)DD,(SMALL(I),DD,I=2,20,2)                          PLTR  95
      WRITE(OUT,3302)(SMALL(I),I=1,21,2)                                PLTR  96
      WRITE(OUT,3001)                                                   PLTR  97
      DO 330 I=1,55                                                     PLTR  98
      VALUE=VALUE-DELY                                                  PLTR  99
      A=BLANK                                                           PLTR 100
      B=PLUS                                                            PLTR 101
      L=I+2                                                             PLTR 102
      IF(L/3*3-L)330,310,330                                            PLTR 103
  310 B=DD                                                              PLTR 104
      IF(L/2*2-L)330,320,330                                            PLTR 105
  320 A=DD                                                              PLTR 106
  330 WRITE(OUT,3300)LABELS(I),VALUE,A,B,(ITEM(I,J),J=1,101),B,A,VALUE  PLTR 107
      WRITE(OUT,3001)                                                   PLTR 108
      WRITE(OUT,3301)DD,(SMALL(I),DD,I=2,20,2)                          PLTR 109
      WRITE(OUT,3302)(SMALL(I),I=1,21,2)                                PLTR 110
      CONTINUE                                                          PLTR 111
      RETURN                                                            PLTR 112
      END                                                               PLTR 113
$      FORTRAN
CDTRN                                                                   DTRN   1
      FUNCTION DTRAN(XX)                                                DTRN   2
C      MDSCAL, VERSION 5MS, OCTOBER, 1971                               DTRN   3
CDTRAN  WRITTEN BY J.KRUSKAL AND J. SEERY, AUGUST, 1971                 DTRN   4
      DIMENSION DUMMY(28)                                               DTRN   5
      COMMON /MDSCL2/ DUMMY,A,B,C,D,E,WCON1,WCON2,WCON3,WCON4,WCON5     DTRN   6
      DTRAN= A*(B*XX+C)**D +E                                           DTRN   7
      RETURN                                                            DTRN   8
      END                                                               DTRN   9
$      FORTRAN
CWTRN                                                                   WTRN   1
      FUNCTION WTRAN(XX)                                                WTRN   2
C      MDSCAL, VERSION 5MS, OCTOBER, 1971                               WTRN   3
CWTRAN WRITTEN BY J.KRUSKAL AND J.SEERY, AUGUST, 1971                   WTRN   4
C        (REPLACES OLD DUMMY VERSION)                                   WTRN   5
      DIMENSION DUMMY(28)                                               WTRN   6
      COMMON /MDSCL2/ DUMMY,DCON1,DCON2,DCON3,DCON4,DCON5,A,B,C,D,E     WTRN   7
      WTRAN= A*(B*XX+C)**D +E                                           WTRN   8
      RETURN                                                            WTRN   9
      END                                                               WTRN  10
$      FORTRAN
CFITP                                                                   FITP   1
       SUBROUTINEFITP(DATA,IJ,DIST,DHAT,WW,MM,NTI,K)                    FITP   2
C      MDSCAL, VERSION 5MS, OCTOBER, 1971                               FITP   3
CFITP   ADAPTED BY J.KRUSKAL AND J.SEERY, AUGUST, 1971                  FITP   4
CFITP          LINEAR REGRESSION FIT ROUTINE FOR MDSCAL                 FITP   5
C                                                                       FITP   6
       REAL DATA(1),DIST(1),DHAT(1),WW(1)                               FITP   7
       INTEGER IJ(1)                                                    FITP   8
       REAL A(5,5),B(5),V(5),C(5)                                       FITP   9
      COMMON  /MDSCL4/ C                                                FITP  10
C                                                                       FITP  11
       NT = MAX0( 1, MIN0( 5,NTI ) )                                    FITP  12
C                                                                       FITP  13
       SUMW=0.0                                                         FITP  14
       DO 1105 I=K,NT                                                   FITP  15
       B(I) = 0.0                                                       FITP  16
       DO 1105 J=K,I                                                    FITP  17
       A(I,J) = 0.0                                                     FITP  18
 1105  CONTINUE                                                         FITP  19
C                                                                       FITP  20
       DO 1150 M=1,MM                                                   FITP  21
       DA=DATA(M)                                                       FITP  22
       SUMW=SUMW+WW(M)                                                  FITP  23
C                                                                       FITP  24
       DO 1110 I=K,NT                                                   FITP  25
       V(I)=REGR(DA,I)                                                  FITP  26
 1110  CONTINUE                                                         FITP  27
C                                                                       FITP  28
       DO 1140 I=K,NT                                                   FITP  29
       DO 1130 J=K,I                                                    FITP  30
       A(I,J)=A(I,J)+V(I)*V(J)*WW(M)                                    FITP  31
 1130  CONTINUE                                                         FITP  32
       B(I)=B(I)+V(I)*DIST(M)*WW(M)                                     FITP  33
 1140  CONTINUE                                                         FITP  34
C                                                                       FITP  35
 1150  CONTINUE                                                         FITP  36
C                                                                       FITP  37
C                                                                       FITP  38
       DO 1155 I = K,NT                                                 FITP  39
       DO 1155 J= K,I                                                   FITP  40
       A(J,I) = A(I,J)                                                  FITP  41
 1155  CONTINUE                                                         FITP  42
C                                                                       FITP  43
       CALL EQSOLV(A,B,C,NT,K)                                          FITP  44
C                                                                       FITP  45
       DO 1170 M=1,MM                                                   FITP  46
       DA = DATA(M)                                                     FITP  47
       TEMP=0.0                                                         FITP  48
       DO 1160 I=K,NT                                                   FITP  49
       TEMP=TEMP+C(I)*REGR(DA,I)                                        FITP  50
 1160  CONTINUE                                                         FITP  51
       DHAT(M)=TEMP                                                     FITP  52
 1170  CONTINUE                                                         FITP  53
C                                                                       FITP  54
       RETURN                                                           FITP  55
       END                                                              FITP  56
$      FORTRAN
CPLT1                                                                   PLT1   1
      SUBROUTINE PLOT1 (Z,Y,NPOI,AXES,IDENT,XMAX,XMIN,YMAX,YMIN)        PLT1   2
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCPLT1   3
C      MDSCAL, VERSION 5MS, OCTOBER, 1971                               PLT1   4
CPLOT1  ADAPTED BY J. KRUSKAL AND J. SEERY   JUNE,1971                  PLT1   5
C                                                                       PLT1   6
C   ROUTINE TO GENERATE A SQUARE PLOT OF ARRAY X VS. Y.                 PLT1   7
C                                                                       PLT1   8
C   THE INPUT VECTORS HAVE NPOI ENTRIES (A NUMBER LESS THAN OR EQUAL TO PLT1   9
C  THE LENGTH OF THE ARGUMENT ARRAYS IN THE CALLING PROGRAM.            PLT1  10
C                                                                       PLT1  11
C   THE PLOTTING IS DONE ON TAPE -OUT-.                                 PLT1  12
C                                                                       PLT1  13
C   IF -AXES- IS POSITIVE, AXES WILL APPEAR ON THE PLOT.                PLT1  14
C   IF -AXES- IS NEGATIVE, NO AXES WILL APPEAR.                         PLT1  15
C                                                                       PLT1  16
C                                                                       PLT1  17
C   IF -XMAX- = -XMIN- THE PROGRAM WILL DETERMINE ITS OWN SCALE FOR     PLT1  18
C   THE X AXIS, AND SIMILARLY FOR THE Y AXIS IF -YMAX- = -YMIN-.        PLT1  19
C  IF -XMAX- IS NOT EQUAL TO -XMIN-, THE PROGRAM WILL USE THESE VALUES  PLT1  20
C   FOR THE UPPER AND LOWER BOUNDS OF THE X AXIS, AND SIMILARLY FOR THE PLT1  21
C   Y AXIS.                                                             PLT1  22
C                                                                       PLT1  23
C   IF -IDENT- IS POSITIVE, THE PLOTTED POINTS WILL BE IDENTIFIED BY    PLT1  24
C   NUMBER.    THERE ARE FIFTY (50) IDENTIFICATION CHARACTERS AVAILABLE.PLT1  25
C   IF -IDENT- IS NEGATIVE, THE PLOTTED POINTS WILL APPEAR AS "X", AND  PLT1  26
C   MULTIPLE POINTS WILL BE IDENTIFIED BY THE NUMBER OF SUPERIMPOSED    PLT1  27
C  POINTS AT THAT LOCATION.  ***NOTE*** THIS OPTION SHOULD BE USED WHEN PLT1  28
C   NPOI SIGNIFICANTLY EXCEEDS 50.                                      PLT1  29
C                                                                       PLT1  30
C   WRITTEN BY FORREST W. YOUNG, NOVEMBER 1965.                         PLT1  31
C         VERSION 3, APRIL 1967                                         PLT1  32
C                                                                       PLT1  33
C----------ADAPTED BY F.J. CARMONE, JR., 5/1/67.                        PLT1  34
C----------MODIFIED BY R.F. MCCRACKEN  11/15/67.                        PLT1  35
C                                                                       PLT1  36
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCPLT1  37
C                                                                       PLT1  38
      DIMENSION Z(1),Y(1),SMALL(25),CHAR(50),HOLL(11)                   PLT1  39
      REAL LABELS(55),ITEM(55,121)                                      PLT1  40
      COMMON  /PLTSS/ SMALL,HOLL,CHAR,LABELS,ITEM,DD,BLANK,AIE,AMINUS   PLT1  41
      INTEGER OUT,G,AXES                                                PLT1  42
       COMMON  /MDSCL1/  LREAD, LPRINT, LPUNCH, LSCRAT                  PLT1  43
      DATA ORIG/1H0/, AMULT/1H;/, GT50/1H>/, PLUS/1H*/                  PLT1  44
       OUT = LPRINT                                                     PLT1  45
      DO 115 I = 1,55                                                   PLT1  46
      DO 115 J=1,121                                                    PLT1  47
  115 ITEM(I,J) = BLANK                                                 PLT1  48
      IF(XMAX .NE. XMIN) GO TO 222                                      PLT1  49
      XMAX = -1.0E7                                                     PLT1  50
      XMIN = +1.0E7                                                     PLT1  51
      DO 221 I = 1,NPOI                                                 PLT1  52
      IF(Z(I) .GT. XMAX) XMAX = Z(I)                                    PLT1  53
      IF(Z(I) .LT. XMIN) XMIN = Z(I)                                    PLT1  54
  221 CONTINUE                                                          PLT1  55
  222 IF(YMAX .NE. YMIN) GO TO 224                                      PLT1  56
      YMAX = -1.0E7                                                     PLT1  57
      YMIN = +1.0E7                                                     PLT1  58
      DO 223 I = 1,NPOI                                                 PLT1  59
      IF(Y(I) .GT. YMAX) YMAX = Y(I)                                    PLT1  60
      IF(Y(I) .LT. YMIN) YMIN = Y(I)                                    PLT1  61
  223 CONTINUE                                                          PLT1  62
  224 CONTINUE                                                          PLT1  63
      SPANX = XMAX - XMIN                                               PLT1  64
      SPANY = YMAX - YMIN                                               PLT1  65
      DELX = SPANX/120.0                                                PLT1  66
      DELY = SPANY/54.0                                                 PLT1  67
      IF(AXES .LE. 0) GO TO 119                                         PLT1  68
      K = 0                                                             PLT1  69
      M = 0                                                             PLT1  70
      IF(XMIN .LT..0 .AND. XMAX .GT..0) GO TO 200                       PLT1  71
      GO TO 202                                                         PLT1  72
  200 K = (-XMIN/DELX) + 1.5                                            PLT1  73
      DO 201 I = 1,55                                                   PLT1  74
  201 ITEM(I,K) = AIE                                                   PLT1  75
  202 CONTINUE                                                          PLT1  76
      IF(YMIN .LT..0 .AND. YMAX .GT..0) GO TO 203                       PLT1  77
      GO TO 119                                                         PLT1  78
  203 M = (-YMIN/DELY) + 1.5                                            PLT1  79
      M = 56 - M                                                        PLT1  80
      DO 204 J = 1,121                                                  PLT1  81
  204 ITEM(M,J) = AMINUS                                                PLT1  82
      IF(K .EQ. 0 .OR. M .EQ. 0) GO TO 119                              PLT1  83
      ITEM(M,K) = ORIG                                                  PLT1  84
  119 CONTINUE                                                          PLT1  85
      VALUE = YMAX + DELY                                               PLT1  86
      SMALL(1) = XMIN                                                   PLT1  87
      DO 120 I=1,24                                                     PLT1  88
  120 SMALL(I+1) = SMALL(I) + 5.0*DELX                                  PLT1  89
      IF(IDENT .GT. 0) GO TO 140                                        PLT1  90
      DO 135 II=1,NPOI                                                  PLT1  91
       I = (YMAX - Y(II))/DELY + 1.5                                    PLT1  92
       J = (Z(II) - XMIN)/DELX + 1.5                                    PLT1  93
      IF(I.GT.55.OR.I.LT.1.OR.J.GT.121.OR.J.LT.1) GO TO 135             PLT1  94
      DO 134 JJ=1,10                                                    PLT1  95
      IF(ITEM(I,J) .EQ. HOLL(JJ)) GO TO 136                             PLT1  96
  134 CONTINUE                                                          PLT1  97
      IF(ITEM(I,J) .EQ. AIE .OR. ITEM(I,J) .EQ. AMINUS .OR. ITEM(I,J)   PLT1  98
     1 .EQ. ORIG)   ITEM(I,J) = HOLL(2)                                 PLT1  99
      GO TO 135                                                         PLT1 100
  136 ITEM(I,J) = HOLL(JJ+1)                                            PLT1 101
  135 CONTINUE                                                          PLT1 102
      GO TO 141                                                         PLT1 103
  140 CONTINUE                                                          PLT1 104
      DO 137 II = 1,NPOI                                                PLT1 105
      I = (YMAX - Y(II))/DELY + 1.5                                     PLT1 106
      J = (Z(II) - XMIN)/DELX + 1.5                                     PLT1 107
      IF(I.GT.55.OR.I.LT.1.OR.J.GT.121.OR.J.LT.1) GO TO 137             PLT1 108
      IF(ITEM(I,J) .EQ. BLANK .OR. ITEM(I,J) .EQ. AIE .OR. ITEM(I,J)    PLT1 109
     1 .EQ. AMINUS .OR. ITEM(I,J) .EQ. ORIG) GO TO  138                 PLT1 110
      ITEM(I,J) = AMULT                                                 PLT1 111
      GO TO 137                                                         PLT1 112
  138 CONTINUE                                                          PLT1 113
      IF(II .GT. 50) GO TO 139                                          PLT1 114
      ITEM(I,J) = CHAR(II)                                              PLT1 115
      GO TO 137                                                         PLT1 116
  139 ITEM(I,J) = GT50                                                  PLT1 117
  137 CONTINUE                                                          PLT1 118
  141 CONTINUE                                                          PLT1 119
      WRITE(OUT,3001)                                                   PLT1 120
 3001 FORMAT( 8X,123H.*....*....*....*....*....*....*....*....*....*....PLT1 121
     1*....*....*....*....*....*....*....*....*....*....*....*....*....*PLT1 122
     2....*.)                                                           PLT1 123
      DO 330 I = 1,55                                                   PLT1 124
      VALUE = VALUE - DELY                                              PLT1 125
      A = BLANK                                                         PLT1 126
      B = DD                                                            PLT1 127
      L = I + 2                                                         PLT1 128
      IF(L/5*5-L) 330,310,330                                           PLT1 129
  310 B = PLUS                                                          PLT1 130
      IF(L/2*2-L) 330,320,330                                           PLT1 131
  320 A = PLUS                                                          PLT1 132
  330 WRITE(OUT,3300) VALUE,A,B,(ITEM(I,J), J=1,121)  ,B,A              PLT1 133
 3300 FORMAT(1X,F6.3,125A1)                                             PLT1 134
      WRITE(OUT,3001)                                                   PLT1 135
      WRITE(OUT,3301) DD,(SMALL(I),DD, I=2,24,2)                        PLT1 136
 3301 FORMAT( 9X,A1,12(F9.4,A1))                                        PLT1 137
      WRITE(OUT,3302) (SMALL(I), I=1,23,2)                              PLT1 138
 3302 FORMAT(4X,12F10.4)                                                PLT1 139
      IF(IDENT .LE. 0) GO TO 150                                        PLT1 140
  150 CONTINUE                                                          PLT1 141
       WRITE(6,160)                                                     PLT1 142
  160 FORMAT(1H0)                                                       PLT1 143
      RETURN                                                            PLT1 144
      END                                                               PLT1 145
$      FORTRAN
CCACT                                                                   CACT   1
       SUBROUTINE CCACT                                                 CACT   2
C      MDSCAL, VERSION 5MS, OCTOBER, 1971                               CACT   3
CCCACT  ADAPTED BY J. KRUSKAL AND J. SEERY    AUGUST, 1971              CACT   4
C      CORRECTED DECEMBER, 1971                                         CACT   5
C                                                                       CACT   6
C             CHANGES MADE FOR "FIX" OPTION JAN 14, 1970                CACT   7
C             ROUTINES CHANGED  MAIN, CCACT, BLDA                       CACT   8
C                                                                       CACT   9
      REAL MAP(40),DUMMY(27),ATAB(4,70)                                 CACT  10
      INTEGER TAB(4,70),IFAC(6)                                         CACT  11
      INTEGER    MTAB(10)                                               CACT  12
       EQUIVALENCE                         (ATAB,TAB)                   CACT  13
C                                                                       CACT  14
       INTEGER                         LPAR(   38 )                     CACT  15
       REAL                             PAR(   38 )                     CACT  16
      EQUIVALENCE (MAP(1),DUMMY(1)), (MAP(28),CHTAB(1))                 CACT  17
       EQUIVALENCE           (LPAR,PAR)                                 CACT  18
C                                                                       CACT  19
       INTEGER  IPARAM(5)                                               CACT  20
       REAL    C(81)                                                    CACT  21
      REAL WORD(18)                                                     CACT  22
       REAL CHTAB(13)                                                   CACT  23
C                                                                       CACT  24
      INTEGER DFLTSW                                                    CACT  25
       REAL BLANK, DOT                                                  CACT  26
       INTEGER BLSW,NUMSW,DECSW,TYPE,XTYPE,T,TA, PARNO, TABNO           CACT  27
C                                                                       CACT  28
       COMMON  /MDSCL1/  LREAD, LPRINT, LPUNCH, LSCRAT                  CACT  29
       COMMON  /MDSCL2/  LPAR                                           CACT  30
       COMMON  /MDSCL3/  MTAB,  TAB                                     CACT  31
C                                                                       CACT  32
        DATA    BLANK, DOT, EQUALS, COMMA, C(81) /1H ,1H.,1H=,1H,,1H  / CACT  33
      DATA DOLLAR/1H$/                                                  CACT  34
      DATA MODP,IFAC/ 16381, 907, 887, 883, 881, 877, 863/              CACT  35
C                                                                       CACT  36
      DATA  MAP/1H ,1HA,1HB,1HC,1HD,1HE,1HF,1HG,1HH,1HI,1HJ,1HK,1HL,    CACT  37
     .          1HM,1HN,1HO,1HP,1HQ,1HR,1HS,1HT,1HU,1HV,1HW,1HX,1HY,    CACT  38
     .          1HZ,1H0,1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9,1H.,1H+,1H-/CACT  39
C                                                                       CACT  40
C                                                                       CACT  41
C                                                                       CACT  42
 1     FORMAT(80A1)                                                     CACT  43
 11    FORMAT(1H0,80A1)                                                 CACT  44
 2     FORMAT(1X,I1)                                                    CACT  45
 3    FORMAT(1X,80A1)                                                   CACT  46
 4     FORMAT(1X,I18)                                                   CACT  47
 5     FORMAT(1X,F18.3)                                                 CACT  48
 9     FORMAT(51H MDSCAL DIAGNOSTIC. CONTROL CARD ABOVE IS IMPROPER.)   CACT  49
  99    FORMAT( 12H ITEM NUMBER, I5, 20H   HAS EXPECTED TYPE ,  I5,     CACT  50
     1 18H  AND ACTUAL TYPE , I5, 16H . THIS ITEM IS , 18A1)            CACT  51
 98   FORMAT(12H ITEM NUMBER,I5,39H  HAS ILLEGAL CHARACTER. THIS ITEM ISCACT  52
     .  ,18A1)                                                          CACT  53
C                                                                       CACT  54
C      READ AND PRINT CONTROL CARD                                      CACT  55
C                                                                       CACT  56
 100   READ  (LREAD,1)  (C(I),I=1,80)                                   CACT  57
       WRITE(LPRINT,11) (C(I),I=1,80)                                   CACT  58
C      ANALYZE CONTROL CARD INTO TOKENS                                 CACT  59
C                                                                       CACT  60
C              EACH  TOKEN  IS DELIMITED BY BLANKS                      CACT  61
C              THERE ARE  FOUR TYPES OF TOKENS.                         CACT  62
C              ALPHABETIC, INTEGER, DECIMAL, AND DEFAULT                CACT  63
C              ALPHABETIC UNLESS FIRST CHARACTER IS DIGIT OR DEC POINT  CACT  64
C                      OR PLUS OR MINUS  OR DOLLAR SIGN                 CACT  65
C              DECIMALS       DISTINGUISHED BY DECIMAL POINT            CACT  66
C              DEFAULT = $                                              CACT  67
C                                                                       CACT  68
       REWIND LSCRAT                                                    CACT  69
       BLSW=1                                                           CACT  70
       NUMSW=1                                                          CACT  71
       DECSW=1                                                          CACT  72
      DFLTSW=1                                                          CACT  73
       K=0                                                              CACT  74
C                                                                       CACT  75
       DO 300 I=1,81                                                    CACT  76
C                                                                       CACT  77
       X=C(I)                                                           CACT  78
       GO TO (110,120), BLSW                                            CACT  79
C                                                                       CACT  80
 110   IF(X.EQ.BLANK  .OR.   X.EQ.EQUALS  .OR.   X.EQ. COMMA ) GO TO 300CACT  81
       BLSW=2                                                           CACT  82
       JA=I                                                             CACT  83
       DO 115 KX=1,13                                                   CACT  84
 115   IF(X.EQ.CHTAB(KX))NUMSW=2                                        CACT  85
      IF(X.EQ.DOLLAR) DFLTSW=2                                          CACT  86
       IF(X.EQ.DOT) DECSW=2                                             CACT  87
       GO TO 300                                                        CACT  88
C                                                                       CACT  89
 120   IF(X.EQ.BLANK  .OR.   X.EQ.EQUALS  .OR.   X.EQ. COMMA ) GO TO 130CACT  90
       IF(X.EQ.DOT) DECSW=2                                             CACT  91
       GO TO 300                                                        CACT  92
C                                                                       CACT  93
 130   K=K+1                                                            CACT  94
       JB = MIN0 (I-1,JA+16)                                            CACT  95
       JC=18-(JB-JA+1)                                                  CACT  96
       TYPE=1                                                           CACT  97
       IF(NUMSW.EQ.2) TYPE=NUMSW+DECSW-1                                CACT  98
      IF(DFLTSW.EQ.2) TYPE=4                                            CACT  99
C                                                                       CACT 100
       WRITE (LSCRAT,2) TYPE                                            CACT 101
      IF(TYPE.EQ.4) GO TO 160                                           CACT 102
       GO TO (140,150),NUMSW                                            CACT 103
 140  WRITE (LSCRAT,3) (C(J),J=JA,JB), (BLANK,J=1,JC)                   CACT 104
       GO TO 160                                                        CACT 105
 150  WRITE (LSCRAT,3) (BLANK,J=1,JC), (C(J),J=JA,JB)                   CACT 106
 160   BLSW=1                                                           CACT 107
       NUMSW=1                                                          CACT 108
       DECSW=1                                                          CACT 109
      DFLTSW=1                                                          CACT 110
       GO TO 300                                                        CACT 111
C                                                                       CACT 112
 300   CONTINUE                                                         CACT 113
C                                                                       CACT 114
C      ANALYZE TOKENS AND SET PARAMETER VALUES ACCORDINGLY              CACT 115
C                                                                       CACT 116
       KB=K                                                             CACT 117
       IF(KB.EQ.0) RETURN                                               CACT 118
       REWIND LSCRAT                                                    CACT 119
       XTYPE=1                                                          CACT 120
       IPARAM(1) = 1                                                    CACT 121
      NT=0                                                              CACT 122
      NTOKNS=0                                                          CACT 123
C                                                                       CACT 124
       DO 1000 K=1,KB                                                   CACT 125
C                                                                       CACT 126
       READ  (LSCRAT,2) TYPE                                            CACT 127
      IF( (XTYPE.EQ.1) .AND. (TYPE.NE.1) ) GO TO 995                    CACT 128
      IF( (XTYPE.NE.TYPE). AND. (TYPE.NE.4) ) GO TO  995                CACT 129
 350  GO TO (400,410,420,430), TYPE                                     CACT 130
C                                                                       CACT 131
 400  READ (LSCRAT,3) (WORD(L),L=1,18)                                  CACT 132
C     CONVERT FIRST SIX LETTERS OF ALPHABETIC WORD INTO CODE NUMBER.    CACT 133
C     LARGEST POSSIBLE CODE NUMBER IS LESS THAN 2**15.                  CACT 134
      ICODE=0                                                           CACT 135
      DO 408 M1=1,6                                                     CACT 136
      X=WORD(M1)                                                        CACT 137
      DO 402 M2=1,37                                                    CACT 138
      IF( X.EQ.MAP(M2) ) GO TO 404                                      CACT 139
 402  CONTINUE                                                          CACT 140
C     ILLEGAL CHARACTER                                                 CACT 141
      WRITE (LPRINT,9)                                                  CACT 142
      WRITE (LPRINT,98) K,(WORD(L),L=1,18)                              CACT 143
      GO TO 999                                                         CACT 144
 404  M3=(M2-1)*IFAC(M1)                                                CACT 145
      IF( M3.GE.MODP ) M3=M3-MODP                                       CACT 146
      ICODE=ICODE+M3                                                    CACT 147
      IF( ICODE.GE.MODP ) ICODE=ICODE-MODP                              CACT 148
 408  CONTINUE                                                          CACT 149
       GO TO 510                                                        CACT 150
C                                                                       CACT 151
 410   READ  (LSCRAT,4) INTPAR                                          CACT 152
      ISUB=PARNO+NT                                                     CACT 153
      LPAR(ISUB)=INTPAR                                                 CACT 154
       GO TO 430                                                        CACT 155
C                                                                       CACT 156
 420   READ  (LSCRAT,5) DECPAR                                          CACT 157
      ISUB=PARNO+NT                                                     CACT 158
      PAR(ISUB)=DECPAR                                                  CACT 159
C                                                                       CACT 160
 430  NT=NT+1                                                           CACT 161
      IF(NT.EQ.NTOKNS) XTYPE=1                                          CACT 162
       GO TO 1000                                                       CACT 163
C                                                                       CACT 164
 510   TABNO = IPARAM(1)                                                CACT 165
       MA = MTAB(TABNO)                                                 CACT 166
       MB = MTAB(TABNO+1) - 1                                           CACT 167
C                                                                       CACT 168
       LMISSW = 0                                                       CACT 169
       DO 700 M=MA,MB                                                   CACT 170
      IF( ICODE.NE.TAB(1,M) ) GO TO 700                                 CACT 171
       LMISSW = 1                                                       CACT 172
 600   XTYPE=1                                                          CACT 173
      NT=0                                                              CACT 174
       IPARAM(1) = 1                                                    CACT 175
      PARNO=TAB(3,M)                                                    CACT 176
      LTEMP=TAB(2,M)                                                    CACT 177
      NTOKNS=TAB(4,M)                                                   CACT 178
      IF(LTEMP.GT.2) NTOKNS=0                                           CACT 179
       GO TO (610,620,630,640,650),  LTEMP                              CACT 180
C                                                                       CACT 181
C      NAME OF INTEGER PARAMETER                                        CACT 182
 610   XTYPE=2                                                          CACT 183
       GO TO 700                                                        CACT 184
C                                                                       CACT 185
C      NAME OF DECIMAL PARAMETER                                        CACT 186
 620   XTYPE=3                                                          CACT 187
       GO TO 700                                                        CACT 188
C                                                                       CACT 189
C      IMPLICITLY SPECIFIED INTEGER PARAMETER                           CACT 190
C      A SINGLE IMPLICITLY SPECIFIED PARAMTER CAN IN PRINCIPLE HOLD     CACT 191
C      AS MANY AS THREE PARAMTERS WHICH THE MAIN PROGRAM CONSIDERS      CACT 192
C      CONCEPTUALLY DISTINCT.  ONE IS IN THE UNITS POSITION, ANOTHER IN CACT 193
C      IN THE HUNDREDS POSITION, AND ANOTHER IN THE TEN-THOUSANDS       CACT 194
C      POSITION.                                                        CACT 195
 630   LP = 1                                                           CACT 196
      LT=TAB(4,M)                                                       CACT 197
       IF(LT.GE.100) LP = 100                                           CACT 198
       IF(LT.GE.1) LP = 10000                                       CACT 1990000
       LQ = 100*LP                                                      CACT 200
       LA = LPAR(PARNO)                                                 CACT 201
      LPAR(PARNO)= LA-(MOD(LA,LQ)/LP)*LP +TAB(4,M)                      CACT 202
       GO TO 700                                                        CACT 203
C                                                                       CACT 204
C      IMPLICITLY SPECIFIED DECIMAL PARAMETER                           CACT 205
 640  PAR(PARNO)=ATAB(4,M)                                              CACT 206
       GO TO 700                                                        CACT 207
C                                                                       CACT 208
C      INTERNAL PARAMETER OF CCACT PROGRAM                              CACT 209
 650  IPARAM(PARNO)=TAB(4,M)                                            CACT 210
       GO TO 700                                                        CACT 211
C                                                                       CACT 212
 700   CONTINUE                                                         CACT 213
      IF(LMISSW.EQ.0) GO TO 994                                         CACT 214
 1000  CONTINUE                                                         CACT 215
      IF(NT.LT.NTOKNS) GO TO 997                                        CACT 216
C                                                                       CACT 217
 1001  RETURN                                                           CACT 218
C                                                                       CACT 219
 994  WRITE(LPRINT,9)                                                   CACT 220
      WRITE(LPRINT,998) (WORD(L),L=1,18)                                CACT 221
 998  FORMAT(27H UNDEFINED CONTROL PHRASE  ,18A1)                       CACT 222
      GO TO 999                                                         CACT 223
 995  WRITE (LPRINT,9)                                                  CACT 224
      WRITE (LPRINT,99) K,XTYPE,TYPE, (WORD(L),L=1,18)                  CACT 225
      GO TO 999                                                         CACT 226
 997  WRITE(LPRINT,9)                                                   CACT 227
      WRITE(LPRINT,97)                                                  CACT 228
 97   FORMAT(47H CONTROL PHRASE NOT COMPLETED ON A SINGLE CARD.)        CACT 229
 999  CALL EXIT                                                         CACT 230
C                                                                       CACT 231
       END                                                              CACT 232
$      FORTRAN
CBLDA                                                                   BLDA   1
       BLOCK   DATA                                                     BLDA   2
C      MDSCAL, VERSION 5MS, OCTOBER, 1971                               BLDA   3
CBLOCK  ADAPTED BY J.KRUSKAL AND J.SEERY, AUGUST, 1971                  BLDA   4
C      CORRECTED DECEMBER, 1971                                         BLDA   5
C                                                                       BLDA   6
C      THIS PROGRAM HOLDS THE  TABLE OF WORDS WHICH CCACT CONSULTS      BLDA   7
C      LOGICALLY IT CONSISTS OF SEVERAL TABLES                          BLDA   8
C      THE VALUES IN MTAB INDICATE WHICH ROWS START NEW TABLES          BLDA   9
C      EACH ENTRY IN THE TABLE CONTAINS 4 ITEMS                         BLDA  10
C      THE FIRST IS THE CODED ALPHABETIC WORD                           BLDA  11
C             (ONE WORD MAY BE ENTERED SEVERAL TIMES )                  BLDA  12
C      THE SECOND INDICATES THE NATURE OF THIS ENTRY                    BLDA  13
C              1 INDICATES INTEGER PARAMETERS WHICH MUST BE  EXPLIC9TLY BLDA  14
C                  READ IN                                              BLDA  15
C              2 INDICATES REAL PARAMETERS WHICH MUST BE EXPLICITLY     BLDA  16
C                  READ IN                                              BLDA  17
C              3 INDICATES AN INTEGER PARAMETER WHICH IS IMPLICITLY     BLDA  18
C                  DEFINED                                              BLDA  19
C              4 INDICATES A REAL     PARAMETER WHICH IS IMPLICITLY     BLDA  20
C                  DEFINED                                              BLDA  21
C              5 INDICATES A PARAMETER WHICH BELONGS TO CCACT ONLY      BLDA  22
C      THE THIRD INDICATES WHICH PARAMTER IS INVOLVED                   BLDA  23
C                                                                       BLDA  24
C              SEPARATE NUMBERING FOR MAIN PROGRAM PARAMETERS AND CCACT BLDA  25
C      THE FOURTH GIVES THE VALUE FOR ANY INTERNALLY DEFINED PARAMETERS BLDA  26
C                                                                       BLDA  27
C                                                                       BLDA  28
C             CHANGES MADE FOR "FIX" OPTION AND BLK. DATA PROB. APRL,70 BLDA  29
C             ROUTINES CHANGED  MAIN, CCACT, BLDA                       BLDA  30
C       NOTE   CONTROL PHRASE "TITLE" HAS BEEN DELETED                  BLDA  31
C                                                                       BLDA  32
       INTEGER                         LPAR(   38 )                     BLDA  33
       REAL                             PAR(   38 )                     BLDA  34
       EQUIVALENCE           (LPAR,PAR)                                 BLDA  35
C                                                                       BLDA  36
      INTEGER    MTAB(10)                                               BLDA  37
      INTEGER TAB1(4,14),TAB2(4,15),TAB3(4,15),TAB4(4,13),TAB5(4,13)    BLDA  38
C                                                                       BLDA  39
      REAL SMALL(25),HOLL(11),CHAR(50),LABELS(55),ITEM(55,121)          BLDA  40
      COMMON  /PLTSS/ SMALL,HOLL,CHAR,LABELS,ITEM,DD,BLANK,AIE,AMINUS   BLDA  41
C                                                                       BLDA  42
       COMMON  /MDSCL1/  LREAD, LPRINT, LPUNCH, LSCRAT                  BLDA  43
       COMMON  /MDSCL2/  LPAR                                           BLDA  44
       COMMON  /MDSCL3/  MTAB,TAB1,TAB2,TAB3,TAB4,TAB5                  BLDA  45
C                                                                       BLDA  46
      DATA  LREAD,LPRINT,LPUNCH,LSCRAT /5,6,43,2/                       BLDA  47
C                                                                       BLDA  48
      DATA   MTAB /1,45,47,52,58,60,69,71,71,71/                        BLDA  49
C                                                                       BLDA  50
C     TAB1 ENTRIES PERTAIN SUCCESSIVELY TO THESE CONTROL WORDS          BLDA  51
C          DIMMAX,    DIMMIN,    DIMDIF,    CUTOFF,    STRMIN,          BLDA  52
C          SFGRMN,    COSAVW,    ACSAVW,    DIAGON,    MATRIX,          BLDA  53
C          HALFMA,    LOWERH,    UPPERH,    LOWERC                      BLDA  54
       DATA  TAB1/                                                      BLDA  55
     1  6989,    1,    1,    1,                                         BLDA  56
     2  5375,    1,    2,    1,                                         BLDA  57
     3  6923,    1,    3,    1,                                         BLDA  58
     4 13520,    2,    4,    1,                                         BLDA  59
     5   390,    2,    5,    1,                                         BLDA  60
     6  2553,    2,    6,    1,                                         BLDA  61
     7  7303,    2,    7,    1,                                         BLDA  62
     8 11226,    2,    8,    1,                                         BLDA  63
     9 11136,    5,    1,    2,                                         BLDA  64
     A  9277,    3,   10,    1,                                         BLDA  65
     B  3527,    3,   10,    2,                                         BLDA  66
     C  6069,    3,   10,    2,                                         BLDA  67
     D  8938,    3,   10,    3,                                         BLDA  68
     E  1754,    3,   10,    4/                                         BLDA  69
C                                                                       BLDA  70
C     TAB2 ENTRIES PERTAIN SUCCESSIVELY TO THESE CONTROL WORDS          BLDA  71
C          UPPERC,    BLOCKD,    DISSIM,    DISSIM,    SIMILA,          BLDA  72
C          SIMILA,    CONFIG,    COMPUT,    PRIMAR,    SECOND,          BLDA  73
C          R     ,    ITERAT,    FIX   ,    SRATST,    STOP             BLDA  74
      DATA  TAB2/                                                       BLDA  75
     1  4623,    3,   10,    5,                                         BLDA  76
     2  8683,    5,    1,    5,                                         BLDA  77
     3 15096,    3,   12,    2,                                         BLDA  78
     4 15096,    3,   11,   10,                                         BLDA  79
     5  6868,    3,   12,    2,                                         BLDA  80
     6  6868,    3,   11,   11,                                         BLDA  81
     7 14846,    3,   12,    5,                                         BLDA  82
     8 11754,    3,   12,    6,                                         BLDA  83
     9   765,    3,   13,    1,                                         BLDA  84
     A  4119,    3,   13,    2,                                         BLDA  85
     B 16326,    2,   14,    1,                                         BLDA  86
     C 15170,    1,   15,    1,                                         BLDA  87
     D  1855,    1,   28,    1,                                         BLDA  88
     E  3720,    2,   16,    1,                                         BLDA  89
     F 13171,    3,   12,    7/                                         BLDA  90
C                                                                       BLDA  91
C     TAB3 ENTRIES PERTAIN SUCCESSIVELY TO THESE CONTROL WORDS          BLDA  92
C          WEIGHT,     WFUNCT,   SFORM1,    SFORM2,    CARDS,           BLDA  93
C          NOCARD,    SPLIT ,    DATA  ,    SAVE  ,    RANDOM,          BLDA  94
C          PRINT ,    DFUNCT,    REGRES,     DCONST,   WCONST           BLDA  95
      DATA  TAB3/                                                       BLDA  96
     1 14543,    3,   12,    3,                                         BLDA  97
     2 11427,    3,   12,    4,                                         BLDA  98
     3  5318,    3,   17,    1,                                         BLDA  99
     4  6181,    3,   17,    2,                                         BLDA 100
     5  6927,    3,   18,    1,                                         BLDA 101
     6 16009,    3,   18,    2,                                         BLDA 102
     7  1966,    5,    1,    3,                                         BLDA 103
     8  6675,    3,   12,    2,                                         BLDA 104
     9  9189,    5,    1,    7,                                         BLDA 105
     A  8330,    1,   20,    1,                                         BLDA 106
     B  2775,    5,    1,    4,                                         BLDA 107
     C 10575,    3,   12,    8,                                         BLDA 108
     D 14439,    5,    1,    6,                                         BLDA 109
     E   267,    2,   29,    5,                                         BLDA 110
     F  2026,    2,   34,    5/                                         BLDA 111
C                                                                       BLDA 112
C     TAB4 ENTRIES PERTAIN SUCCESSIVELY TO THESE CONTROL WORDS          BLDA 113
C          ABSENT,    PRESEN,    BYROWS,    BYGROU,    BYDECK,          BLDA 114
C          NOMORE,    NOMORE,    DATA  ,    DISTAN,    HISTOR,          BLDA 115
C          NODATA,    NODIST,    NOHIST                                 BLDA 116
      DATA  TAB4/                                                       BLDA 117
     1  4258,    3,    9,    2,                                         BLDA 118
     2  2575,    3,    9,    1,                                         BLDA 119
     3  7761,    3,   19,  100,                                         BLDA 120
     4 11782,    3,   19,  200,                                         BLDA 121
     5 11288,    3,   19,  300,                                         BLDA 122
     6  5274,    3,   19,  400,                                         BLDA 123
     7  5274,    3,   19,    2,                                         BLDA 124
     8  6675,    3,   21,    2,                                         BLDA 125
     9  9824,    3,   22,    2,                                         BLDA 126
     A 12801,    3,   24,    2,                                         BLDA 127
     B 16057,    3,   21,    1,                                         BLDA 128
     C  5863,    3,   22,    1,                                         BLDA 129
     D  9395,    3,   24,    1/                                         BLDA 130
C                                                                       BLDA 131
C     TAB5 ENTRIES PERTAIN SUCCESSIVELY TO THESE CONTROL WORDS          BLDA 132
C          NO    ,    YES   ,    MONOTO,    ASCEND,    DESCEN,          BLDA 133
C          POLYNO,    POLYNO,    CONSTA,    NOCONS,    MULTIV,          BLDA 134
C          MULTIV,    DATA  ,    CONFIG                                 BLDA 135
      DATA  TAB5/                                                       BLDA 136
     1  9622,    3,   10,  100,                                         BLDA 137
     2 11125,    3,   10,  200,                                         BLDA 138
     3 15634,    3,   23,    0,                                         BLDA 139
     4  7782,    3,   23,   10,                                         BLDA 140
     5 11188,    3,   23,   11,                                         BLDA 141
     6  3756,    3,   26,    1,                                         BLDA 142
     7  3756,    1,   23,    1,                                         BLDA 143
     8 14387,    3,   26,  100,                                         BLDA 144
     9  5018,    3,   26,  200,                                         BLDA 145
     A  3608,    3,   26,    0,                                         BLDA 146
     B  3608,    1,   23,    1,                                         BLDA 147
     C  6675,    3,   25,    2,                                         BLDA 148
     D 14846,    3,   27,    3/                                         BLDA 149
C                                                                       BLDA 150
      DATA BLANK/1H /, DD/1H./,   AIE/1H /, AMINUS/1H-/                 BLDA 151
      DATA LABELS /12*1H ,1HS,1H ,1HH,1H ,1HE,1H ,1HP,1H ,1HA,1H ,1HR,  BLDA 152
     1                1H ,1HD,5*1H ,  1HD,1H ,1HI,1H ,1HA,1H ,1HG,1H ,  BLDA 153
     2                1HR,1H ,1HA,1H ,1HM,12*1H /                       BLDA 154
      DATA CHAR /1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9,1HA,1HB,1HC,1HD,   BLDA 155
     1   1HE,1HF,1HG,1HH,1HI,1HJ,1HK,1HL,1HM,1HN,1HO,1HP,1HQ,1HR,1HS,   BLDA 156
     2   1HT,1HU,1HV,1HW,1HX,1HY,1HZ,1H+,    1H/,1H=,1H*,1H&,1H$,1H@,   BLDA 157
     3   1H%,1H-  ,1H<,1H(,1H),1H",1H#,1H'/                             BLDA 158
      DATA HOLL /1H ,1HX,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9,1HM/           BLDA 159
       END                                                              BLDA 160