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