C kyst.f C The authors of this software are Joseph B Kruskal, Forrest Young, C 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. There are few differences between kyst and kyst2a. 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----+----@----+----@----+----@----+----@----+----@----+----@----+----@ $ FORTREX STAB CMAIN MAIN 1 C ********************* KYST JANUARY,1973 *********MAIN 2 C MAIN 3 C KRUSKAL-YOUNG-SHEPARD-TORGERSON MULTIDIMENSIONAL SCALING PROGRAM MAIN 4 C MAIN 5 C MAIN 6 C WRITTEN BY MAIN 7 C MAIN 8 C DR. J. B. KRUSKAL MAIN 9 C BELL TELEPHONE LABORATORIES MAIN 10 C MURRAY HILL, N. J. MAIN 11 C MAIN 12 C DR. F. W. YOUNG MAIN 13 C PSYCHOMETRIC LABORATORY MAIN 14 C UNIVERSITY OF NORTH CAROLINA MAIN 15 C CHAPEL HILL, N. C. MAIN 16 C ASSISTED BY MAIN 17 C JUDITH B. SEERY MAIN 18 C BELL TELEPHONE LABORATORIES MAIN 19 C MURRAY HILL, N. J. MAIN 20 C MAIN 21 C MAIN 22 C MAIN 23 C **********************************************************************MAIN 24 C MAIN 25 C GENERAL REMARKS. MAIN 26 C MAIN 27 C MAIN 28 C KYST INCLUDES THE FOLLOWING SUBROUTINES MAIN 29 C MAIN 30 C BLDA EQSOLV PLOT SERCH MAIN 31 C BSEC1 FITM REGR SGEV MAIN 32 C CCACT FITMS RM1POW SORT MAIN 33 C CONFIG FITP RPOWER ST03 MAIN 34 C DATAPR INICON RROOT WTRAN MAIN 35 C DFLT NEWSTP RUNIFV XITEM MAIN 36 C DTRAN NRMLZ SBK NVIT1 MAIN 37 C MAIN 38 C ALL ROUTINES ARE WRITTEN IN FORTRAN IV, AND PUNCHED IN THE IBMEL 029 MAIN 39 C CHARACTER SET MAIN 40 C MAIN 41 C NO USE IS MADE OF SPECIAL OR NON-STANDARD SOFTWARE. MAIN 42 C MAIN 43 C ALL INPUT AND OUTPUT IS ONTO FORTRAN LOGICAL UNITS WITH NUMBERS MAIN 44 C CONTROLLED BY THESE VARIABLES LREAD,LPRINT,LPUNCH,LSCRAT. MAIN 45 C UNIT NUMBERS 5,6,7,8 HAVE BEEN USED RESPECTIVELY. TO CHANGE THESEMAIN 46 C ASSIGNMENTS, CHANGE THE VALUES FOR THE FOUR VARIABLES SET IN THE MAIN 47 C BLOCK DATA SUBPROGRAM. MAIN 48 C MAIN 49 C THE SCRATCH UNIT IS USED IN A MINOR WAY BY CCACT ONLY. MANY MAIN 50 C INSTALLATIONS WILL HAVE ALTERNATE METHODS OF DOING THE SAME THING.MAIN 51 C MAIN 52 C KYST USES THREE MACHINE DEPENDENT CONSTANTS PRECSN (RELATIVE MACHINE MAIN 53 C PRECISION), XMAG (LEAST POSITIVE MACHINE NUMBER), XMAXN (MINIMUM MAIN 54 C OF LARGEST POSITIVE MACHINE NUMBER AND MINUS LARGEST NEGATIVE MAIN 55 C MACHINE NUMBER). VALUES 1.5E-8, 1.0E-38, 1.0E38 HAVE BEEN USED MAIN 56 C RESPECTIVELY. TO CHANGE THESE ASSIGNMENTS, CHANGE THE VALUES MAIN 57 C FOR THESE VARIABLES SET IN THE BLOCK DATA SUBPROGRAM. MAIN 58 C MAIN 59 C **********************************************************************MAIN 60 C MAIN 61 C MAIN 62 C MAIN PROGRAM FOR KYST JANUARY,1973 MAIN 63 C WRITTEN BY J.KRUSKAL MAIN 64 C MODIFIED FOR KYST BY J.KRUSKAL AND J.SEERY JANUARY,1973 MAIN 65 DIMENSION DATA(1800), IJ(1800), DIST(1800), DHAT(1800) MAIN 66 DIMENSION WW(1800),LDIST(1800) MAIN 67 DIMENSION X(100,6), GR(100,6), GL(100,6) MAIN 68 DIMENSION CTITLE(80),FMAT(80),TITLE(80),LBLOCK(1800) MAIN 69 INTEGER GRNO(101), NOGRPS, LSPLIT, SPLITH MAIN 70 INTEGER GRSDIS(100), SDSWIT,SDSWT1 MAIN 71 REAL GRSTRS(100),PCOEFF(5) MAIN 72 REAL ITEM(101),PTID(100),PMAT(1800,2),RVEC(1800),XNUM(6) MAIN 73 DIMENSION STPL(10),DIMSV(10) MAIN 74 INTEGER DUMMY(25),RTYPE MAIN 75 C MAIN 76 C MAIN 77 COMMON /ACCUR/ PRECSN,XMAG,XMAXN MAIN 78 COMMON /KYST1/ DATA,WW,IJ,X MAIN 79 COMMON /KYST2/ GR,GL,RVEC,PCOEFF,DUMMY,PMAT MAIN 80 COMMON /MDSCL1/ LREAD, LPRINT, LPUNCH, LSCRAT MAIN 81 COMMON /MDSCL2/ MAIN 82 1 LDIMX, LDIMN, LDIMD, CUTOFF, STRMIN MAIN 83 2 , SFGRMN, COSAVW, ACSAVW, IFIRST, MATSWP MAIN 84 3 , SDSWIT, LCSWIT, LFITSW, R, NOIT MAIN 85 4 , SRATST, LSCH, LPUNSW, LSPL, LRANDC MAIN 86 5 , LDAPRT, LDIPRT, LREG, LHIPRT, NUDASW MAIN 87 6 , LPOLSP, LCONSW, LNFIXZ, DCON1, DCON2 MAIN 88 7 , DCON3, DCON4, DCON5, WCON1, WCON2 MAIN 89 8 , WCON3, WCON4, WCON5, NOITIN, LPLCON MAIN 90 9 , LPLSCT, LCOORS MAIN 91 COMMON /METRIC/ RTYPE,RM1,RECR,RR MAIN 92 COMMON /PLTCHR/ PTID,ITEM MAIN 93 C MAIN 94 EQUIVALENCE (DIST,LDIST) , (LBLOCK,RVEC) MAIN 95 EQUIVALENCE (PMAT(1,1),DIST(1)), (PMAT(1,2),DHAT(1)) MAIN 96 EQUIVALENCE (XNUM(1),ONE),(XNUM(2),TWO),(XNUM(3),THREE), MAIN 97 . (XNUM(4),FOUR),(XNUM(5),FIVE),(XNUM(6),SIX) MAIN 98 EXTERNAL WTRAN , DTRAN MAIN 99 C MAIN 100 DATA AAA,BEE,CEE,DEE,EEE,FFF /1HA,1HB,1HC,1HD,1HE,1HF/ MAIN 101 DATA ONE,TWO,THREE,FOUR,FIVE,SIX /1H1,1H2,1H3,1H4,1H5,1H6/ MAIN 102 DATA KPACK,IJPACK,NPREVZ /10000,100,0/ , AIE/1H-/ MAIN 103 C MAIN 104 C INITIALIZE PARAMETERS MAIN 105 C MAIN 106 100 CONTINUE MAIN 107 C MAIN 108 C ALPHABETICAL ORDER MAIN 109 ACSAVW=0.66 MAIN 110 COSAVW=0.66 MAIN 111 CUTOFF=-1.23E+20 MAIN 112 DCON1=0.0 MAIN 113 DCON2=1.0 MAIN 114 DCON3=1.0 MAIN 115 DCON4=0.0 MAIN 116 DCON5=1.0 MAIN 117 GRNO(1) = 1 MAIN 118 IFIRST=1 MAIN 119 ISTC=0 MAIN 120 LCONSW=4 MAIN 121 LCOORS=2 MAIN 122 LCSWIT=1 MAIN 123 LDAPRT = 1 MAIN 124 LDIMD=1 MAIN 125 LDIMN=2 MAIN 126 LDIMX=2 MAIN 127 LDIPRT = 1 MAIN 128 LFITSW=1 MAIN 129 LHIPRT = 2 MAIN 130 LNFIXZ = 0 MAIN 131 LPOLSP = 100 MAIN 132 LPLCON=2 MAIN 133 LPLSCT=1 MAIN 134 LPUNSW = 2 MAIN 135 LRANDC=-101 MAIN 136 LREG=0 MAIN 137 LSCH=1 MAIN 138 LSPL=302 MAIN 139 MATSWP=101 MAIN 140 NOIT=50 MAIN 141 NOITIN=1 MAIN 142 NUDASW=1 MAIN 143 R=2.0 MAIN 144 SDSIN=1.0 MAIN 145 SDSWIT = 10 MAIN 146 SFGRMN=0.0 MAIN 147 SRATST=0.999 MAIN 148 STRESS=1.0 MAIN 149 STRMIN=0.01 MAIN 150 WCON1=0.0 MAIN 151 WCON2=1.0 MAIN 152 WCON3=1.0 MAIN 153 WCON4=0.0 MAIN 154 WCON5=1.0 MAIN 155 WRITE(LPRINT, 17) MAIN 156 17 FORMAT(1H1) MAIN 157 C MAIN 158 C CONTROL CARD READ MAIN 159 C MAIN 160 1000 LCSWIT=1 MAIN 161 CALL CCACT MAIN 162 GO TO (1000, 1100, 1165, 1182, 1200, 2000, 9000, 1190 ), LCSWIT MAIN 163 C MAIN 164 C DATA READ MAIN 165 C MAIN 166 1100 CONTINUE MAIN 167 NREPL1=1 MAIN 168 NREPL3 = 1 MAIN 169 LSPLIT=LSPL/100 MAIN 170 SPLITH = MOD(LSPL,100) MAIN 171 MATSW = MOD(MATSWP,100) MAIN 172 LBLKDS = MATSWP/100 MAIN 173 LPOSW = MOD(LPOLSP,100) MAIN 174 LFITRM = LPOLSP/100 MAIN 175 IF(LREG.EQ.0) GO TO 1104 MAIN 176 SDSWIT=LREG MAIN 177 IF(LREG.LT.10) SDSWIT = SDSWIT+LPOSW MAIN 178 1104 CONTINUE MAIN 179 SDSWIT = SDSWIT+100*LFITRM MAIN 180 C MAIN 181 IF(NUDASW.NE.1) MLASTD=MM MAIN 182 IF(NUDASW.NE.1) GO TO 1106 MAIN 183 NUDASW=2 MAIN 184 MLASTD=0 MAIN 185 MMIN=0 MAIN 186 NOGRPS=0 MAIN 187 1106 M=MLASTD MAIN 188 MA = M+1 MAIN 189 C MAIN 190 READ(LREAD,10) TITLE MAIN 191 WRITE(LPRINT,11) TITLE MAIN 192 IF(MATSW.GE.4) GO TO 1102 MAIN 193 1101 READ (LREAD,12) NPART,NREPL2,NREPL4 MAIN 194 WRITE (LPRINT,13) NPART,NREPL2,NREPL4 MAIN 195 NROWS = NPART MAIN 196 NCOLS = NPART MAIN 197 GO TO 1103 MAIN 198 1102 READ (LREAD,12) NROWS,NCOLS,NREPL2,NREPL4 MAIN 199 WRITE (LPRINT,13) NROWS,NCOLS,NREPL2,NREPL4 MAIN 200 NPART = NROWS+NCOLS MAIN 201 C MAIN 202 1103 N=NPART MAIN 203 IF(NREPL4.NE.0) NREPL3 = NREPL4 MAIN 204 IF(NREPL2.NE.0) NREPL1 = NREPL2 MAIN 205 IF(LBLKDS.EQ.2) N = NPART*NREPL3 MAIN 206 IF(NREPL2.EQ.0 .OR. NREPL4.EQ.0 ) WRITE (LPRINT,16) MAIN 207 1105 READ(LREAD,10) FMAT MAIN 208 WRITE(LPRINT,11) FMAT MAIN 209 C MAIN 210 DO 1162 NREPL=1,NREPL3 MAIN 211 IJBLKD=0 MAIN 212 IF(LBLKDS.EQ.2) IJBLKD =(NREPL-1)*NPART MAIN 213 MLASTG=M MAIN 214 IA = 1 MAIN 215 IB = NROWS MAIN 216 IF(MATSW.EQ.2) IA = IFIRST MAIN 217 IF(MATSW.EQ.3) IB = NROWS-(IFIRST-1) MAIN 218 C MAIN 219 DO 1160 I = IA,IB MAIN 220 MLASTR=M MAIN 221 LTEMP = NCOLS MAIN 222 IF(MATSW.EQ.2) LTEMP = I-IFIRST+1 MAIN 223 IF(MATSW.EQ.3) LTEMP = NCOLS-I-IFIRST+2 MAIN 224 MB = MA + LTEMP * NREPL1 - 1 MAIN 225 ITRUE = I+IJBLKD MAIN 226 IF(MATSW.EQ.4) ITRUE = ITRUE+NCOLS MAIN 227 C MAIN 228 1130 READ (LREAD, FMAT) (DATA(MP),MP=MA,MB) MAIN 229 C MAIN 230 DO 1150 MP=MA,MB MAIN 231 IF( DATA(MP)-CUTOFF ) 1150, 1150, 1140 MAIN 232 1140 CONTINUE MAIN 233 M = M + 1 MAIN 234 DATA(M)=DATA(MP) MAIN 235 WW(M) = 1.0 MAIN 236 J=((MP-MA)/NREPL1)+1 MAIN 237 IF(MATSW.EQ.3) J = J + (I-1) + (IFIRST-1) MAIN 238 J = J+IJBLKD MAIN 239 IF(MATSW.EQ.5) J = J + NROWS MAIN 240 IJ(M)=IJPACK*(ITRUE-1)+J-1 MAIN 241 LDIST(M) = MP MAIN 242 1150 CONTINUE MAIN 243 C MAIN 244 IF(LSPLIT.EQ.1 .AND. M.GT.MLASTR) GO TO 1155 MAIN 245 GO TO 1160 MAIN 246 1155 NOGRPS=NOGRPS+1 MAIN 247 GRNO(NOGRPS+1) = M+1 MAIN 248 GRSDIS(NOGRPS)=SDSWIT MAIN 249 1160 MA = M + 1 MAIN 250 C MAIN 251 IF(LSPLIT.EQ.2 .AND. M.GT.MLASTG) GO TO 1161 MAIN 252 GO TO 1162 MAIN 253 1161 NOGRPS=NOGRPS+1 MAIN 254 GRNO(NOGRPS+1) = M+1 MAIN 255 GRSDIS(NOGRPS)=SDSWIT MAIN 256 1162 CONTINUE MAIN 257 C MAIN 258 IF(LSPLIT.EQ.3 .AND. M.GT.MLASTD) GO TO 1163 MAIN 259 IF(LSPLIT.EQ.4 .AND. SPLITH.EQ.2 .AND. M.GT.MLASTD) GO TO 1163 MAIN 260 GO TO 1164 MAIN 261 1163 NOGRPS=NOGRPS+1 MAIN 262 GRSDIS(NOGRPS)=SDSWIT MAIN 263 C MAIN 264 1164 CONTINUE MAIN 265 GRNO(NOGRPS+1) = M+1 MAIN 266 IF(LSPLIT.LE.3) SPLITH=2 MAIN 267 IF(LSPLIT.EQ.4) SPLITH=1 MAIN 268 IF((NOGRPS.EQ.1).OR.(MMIN.EQ.0)) MMIN=M MAIN 269 C MAIN 270 LSPL = 100 * LSPLIT + SPLITH MAIN 271 C MAIN 272 MM = M MAIN 273 GO TO 1000 MAIN 274 C MAIN 275 C WEIGHTS READ MAIN 276 C MAIN 277 1165 CONTINUE MAIN 278 M=MLASTD MAIN 279 MA=M+1 MAIN 280 DO 1181 NREPL=1,NREPL3 MAIN 281 IA = 1 MAIN 282 IB = NROWS MAIN 283 IF(MATSW.EQ.2) IA = IFIRST MAIN 284 IF(MATSW.EQ.3) IB = NROWS-(IFIRST-1) MAIN 285 DO 1180 I = IA,IB MAIN 286 LTEMP = NCOLS MAIN 287 IF(MATSW.EQ.2) LTEMP = I-IFIRST+1 MAIN 288 IF(MATSW.EQ.3) LTEMP = NCOLS-I-IFIRST+2 MAIN 289 MB = MA + (LTEMP-1)*NREPL1 MAIN 290 1172 READ (LREAD,FMAT)(WW(MP),MP=MA,MB) MAIN 291 DO 1177 MP=MA,MB MAIN 292 IF(LDIST(M+1).GT.MP) GO TO 1177 MAIN 293 M = M + 1 MAIN 294 WW(M)=WW(MP) MAIN 295 1177 CONTINUE MAIN 296 1180 MA = M + 1 MAIN 297 1181 CONTINUE MAIN 298 GO TO 1000 MAIN 299 C MAIN 300 C WEIGHT FORMATION BY WFUNCTION MAIN 301 C MAIN 302 1182 CONTINUE MAIN 303 WRITE(LPRINT,18) WCON1,WCON2,WCON3,WCON4,WCON5 MAIN 304 18 FORMAT(1H0,10X,69H** WEIGHTS FORMED ACCORDING TO RULE WW(I,J)=S+TMAIN 305 .*((A+B*DATA(I,J))**P)/,15X, 9H WHERE A=,F15.7,5X,2HB=,F15.7,5X, MAIN 306 .2HP=,F15.7,5X,2HS=,F15.7,5X,2HT=,F15.7) MAIN 307 MA=MLASTD+1 MAIN 308 DO 1185 M=MA,MM MAIN 309 TEMP1=DATA(M) MAIN 310 WW(M) = WTRAN(TEMP1) MAIN 311 1185 CONTINUE MAIN 312 GO TO 1000 MAIN 313 C MAIN 314 C DATA TRANSFORMATION BY DFUNCTION MAIN 315 C MAIN 316 1190 CONTINUE MAIN 317 WRITE(LPRINT,19) DCON1,DCON2,DCON3,DCON4,DCON5 MAIN 318 19 FORMAT(1H0,10X,73H** DATA TRANSFORMED ACCORDING TO RULE DATA(I,J)MAIN 319 .=S+T*((A=B*DATA(I,J))**P),/,15X, 9H WHERE A=,F15.7,5X,2HB=, MAIN 320 .F15.7,5X,2HP=,F15.7,5X,2HS=,F15.7,5X,2HT=,F15.7) MAIN 321 MA=MLASTD+1 MAIN 322 DO 1195 M=MA,MM MAIN 323 DATA(M)=DTRAN(DATA(M)) MAIN 324 1195 CONTINUE MAIN 325 GO TO 1000 MAIN 326 C MAIN 327 C CONFIGURATION READ MAIN 328 C MAIN 329 1200 LCONSW=2 MAIN 330 READ(LREAD,10) CTITLE MAIN 331 WRITE(LPRINT,11)CTITLE MAIN 332 READ (LREAD, 12) NCON, LDIMCO MAIN 333 WRITE (LPRINT, 13) NCON, LDIMCO MAIN 334 READ(LREAD,10) FMAT MAIN 335 WRITE(LPRINT,11) FMAT MAIN 336 DO 1210 I=1,NCON MAIN 337 READ (LREAD, FMAT) (X(I,L),L=1,LDIMCO) MAIN 338 WRITE(LPRINT,4008) I,(X(I,L),L=1,LDIMCO) MAIN 339 4008 FORMAT(1X,I2,10F7.3) MAIN 340 1210 CONTINUE MAIN 341 GO TO 1000 MAIN 342 C MAIN 343 C SOME INPUT FORMATS MAIN 344 C MAIN 345 10 FORMAT (80A1) MAIN 346 11 FORMAT (1H0,80A1) MAIN 347 12 FORMAT(24I3) MAIN 348 13 FORMAT(1X,24I3) MAIN 349 15 FORMAT(20F4.0) MAIN 350 16 FORMAT(118H THE NUMBER OF INTERNAL REPLICATES, OR THE NUMBER OF EMAIN 351 1XTERNAL REPLICATES HAS NOT BEEN SPECIFIED. IT IS TAKEN TO BE 1 . )MAIN 352 C MAIN 353 C COMPUTATION *************************************MAIN 354 2000 CONTINUE MAIN 355 C MAIN 356 FN=FLOAT (N) MAIN 357 SQRTN=SQRT (FN) MAIN 358 FNGRPS = FLOAT(NOGRPS) MAIN 359 LDIM=LDIMX MAIN 360 RR=R MAIN 361 RTYPE=3 MAIN 362 IF(R.EQ.1.0) RTYPE=1 MAIN 363 IF(R .EQ.2.0) RTYPE=2 MAIN 364 RM1=R-1.0 MAIN 365 RECR=1.0/R MAIN 366 IF(LDAPRT.EQ.2) MAIN 367 1 CALL DATAPR(GRNO,MM,NOGRPS,1) MAIN 368 WRITE (LPRINT,23) MAIN 369 C MAIN 370 C FINISH STARTING CONFIGURATION. LCONSW = 2,3,4. MAIN 371 C 4=USE INICON, 3=SAVED, 2=WAS READ IN. MAIN 372 C TO FILL IN ADDITIONAL POINTS, IF NEEDED MAIN 373 C LRANDC<0, USE ARBITRARY, ELSE USE RANDOM. MAIN 374 C MAIN 375 2100 CONTINUE MAIN 376 WRITE(LPRINT,11) TITLE MAIN 377 IF((LCONSW.EQ.4) .AND. (LRANDC.EQ.(-101)))GO TO 2140 MAIN 378 ISTCON = 1 MAIN 379 IF( LCONSW .EQ. 2 ) ISTCON = NCON + 1 MAIN 380 IF( LCONSW .EQ. 3 ) ISTCON = NPREVZ + 1 MAIN 381 IF( ISTCON .GT. N ) GO TO 2200 MAIN 382 IF(LRANDC.LE.0) GO TO 2110 MAIN 383 DO 2105 K=1,LRANDC MAIN 384 2105 TEMP=RUNIFV(1) MAIN 385 2110 TEMP1=0.01 MAIN 386 DO 2130 I=ISTCON,N MAIN 387 DO 2120 L=1,LDIMX MAIN 388 2120 X(I,L)=0.0 MAIN 389 K= MOD (I-1,LDIM) +1 MAIN 390 X(I,K)=TEMP1 MAIN 391 IF(LRANDC.LT.0) GO TO 2130 MAIN 392 DO 2125 L=1,LDIMX MAIN 393 TEMP = RUNIFV(1) MAIN 394 2125 X(I,L) = ALOG( TEMP/(1.0-TEMP)) MAIN 395 2130 TEMP1=TEMP1+0.01 MAIN 396 GO TO 2200 MAIN 397 2140 SDSWIT=MOD(GRSDIS(1),100) MAIN 398 IF(SDSWIT.EQ.11) SDSIN=-1.0 MAIN 399 IF(N.GT.60) NOITIN=0 MAIN 400 CALL INICON(N,LDIM,MMIN,NOITIN,LHIPRT,SDSIN,LSCH) MAIN 401 C MAIN 402 C SORT DATA AND IJ AND WW. MAIN 403 C ALSO RECORD BLOCKS OF EQUAL DATA VALUES. MAIN 404 C MAIN 405 2200 CONTINUE MAIN 406 NPREVZ = N MAIN 407 DO 2250 NG = 1,NOGRPS MAIN 408 MX = GRNO(NG) MAIN 409 MY = GRNO(NG+1) - 1 MAIN 410 MZ = MY - MX + 1 MAIN 411 SDSWIT = GRSDIS(NG) MAIN 412 SDSWIT = MOD(SDSWIT,100) MAIN 413 IF(SDSWIT.EQ.11)SDSWIT = -10 MAIN 414 CALL SORT( DATA(MX),MZ,IJ(MX),WW(MX),DUMMY,2,SDSWIT) MAIN 415 C MAIN 416 C SORT WILL SORT THE MM ELEMNTS OF DATA IN ALGEBRAIC MAIN 417 C ORDER, ASCENDING OR DESCENDING ACCORDING TO WHETHER MAIN 418 C SDSWIT IS + OR -. MAIN 419 C AT THE SAME TIME, THE ELEMENTS IN IJ AND IN WW WILL BMAIN 420 C REARRANGED IN EXACTLY THE SAME ORDER. THUS THE MAIN 421 C CORRESPONDENCE BETWEEN THE ELEMENTS OF DATA AND IJ MAIN 422 C AND WW IS PRESERVED. MAIN 423 C MAIN 424 DO 2240 MB = MX,MY MAIN 425 IF((DATA(MB+1).EQ.DATA(MB)).AND.(MB.NE.MY)) GO TO 2240 MAIN 426 IJ(MB)=MOD(IJ(MB),KPACK) MAIN 427 IJ(MB)=IJ(MB)+KPACK MAIN 428 2240 CONTINUE MAIN 429 2250 CONTINUE MAIN 430 C MAIN 431 C START COMPUTATION IN CURRENT DIMENSION MAIN 432 C MAIN 433 2300 FLDIM=FLOAT (LDIM) MAIN 434 ITNO=0 MAIN 435 COSAV=0.0 MAIN 436 SRATAV=0.8 MAIN 437 ACSAV=0.0 MAIN 438 STEP = 0.0 MAIN 439 NBAKUP = 0 MAIN 440 C MAIN 441 C INITIALIZE GRADIENT MAIN 442 C MAIN 443 2400 TEMP1=SQRT (1.0/FLDIM) MAIN 444 DO 2410 I=1,N MAIN 445 DO 2410 L=1,LDIM MAIN 446 2410 GR(I,L)=TEMP1 MAIN 447 SFGR=SQRTN MAIN 448 C MAIN 449 C PRINT HEADING FOR INFORMATION ABOUT SCALING IN CURRENT DIMENSION MAIN 450 C MAIN 451 2500 WRITE (LPRINT, 20) N, MM, NOGRPS, LDIM MAIN 452 IF(LHIPRT.EQ.2) WRITE (LPRINT,21) MAIN 453 WRITE (LPRINT, 22) MAIN 454 20 FORMAT(28H0HISTORY OF COMPUTATION. N= , I4, MAIN 455 1 15H. THERE ARE , I6, MAIN 456 2 26H DATA VALUES, SPLIT INTO , I4, MAIN 457 3 27H LISTS. DIMENSION = , I4 ) MAIN 458 21 FORMAT(52H0ITERATION STRESS SRAT SRATAV CAGRGL COSAV ACSAV, MAIN 459 1 16H SFGR STEP ) MAIN 460 22 FORMAT(1X) MAIN 461 23 FORMAT(1H1) MAIN 462 C MAIN 463 C START CURRENT ITERATION *************************************MAIN 464 C MAIN 465 C NORMALIZE CONFIGURATION. MOVE AND CLEAR GRADIENT. MAIN 466 C MAIN 467 3000 TEMP1=1.0 MAIN 468 IF(LCOORS.EQ.0) GO TO 3035 MAIN 469 TEMP1=0.0 MAIN 470 DO 3030 L=1,LDIM MAIN 471 TEMP2=0.0 MAIN 472 DO 3010 I=1,N MAIN 473 3010 TEMP2=TEMP2+X(I,L) MAIN 474 TEMP2=TEMP2/FN MAIN 475 DO 3020 I=1,N MAIN 476 X(I,L)=X(I,L)-TEMP2 MAIN 477 3020 TEMP1=TEMP1+X(I,L)**2 MAIN 478 3030 CONTINUE MAIN 479 TEMP1=SQRT (FN/TEMP1) MAIN 480 3035 DO 3040 L=1,LDIM MAIN 481 DO 3040 I=1,N MAIN 482 X(I,L)=TEMP1*X(I,L) MAIN 483 GL(I,L)=TEMP1*GR(I,L) MAIN 484 3040 GR(I,L)=0.0 MAIN 485 SFGL=TEMP1*SFGR MAIN 486 C MAIN 487 STBAMU = TEMP1 MAIN 488 C LOOP THROUGH THE DATA GROUPS MAIN 489 C MAIN 490 STRLST = STRESS MAIN 491 STRESS = 0.0 MAIN 492 DO 3340 NG = 1,NOGRPS MAIN 493 MX = GRNO(NG) MAIN 494 MY = GRNO(NG+1) - 1 MAIN 495 MZ = MY - MX + 1 MAIN 496 SDSWIT = GRSDIS(NG) MAIN 497 LFITRM = SDSWIT/100 MAIN 498 SDSWIT = MOD(SDSWIT,100) MAIN 499 IF(SDSWIT.EQ.11)SDSWIT = -10 MAIN 500 C MAIN 501 C COMPUTE DISTANCES AND FIND BEST-FITTING MONOTONE PSEUDO-DISTANCESMAIN 502 C MAIN 503 SUMW = 0.0 MAIN 504 DBAR = 0.0 MAIN 505 DO 3120 M=MX,MY MAIN 506 LTEMP1=MOD(IJ(M),KPACK) MAIN 507 I=LTEMP1/IJPACK+1 MAIN 508 J=MOD(LTEMP1,IJPACK)+1 MAIN 509 TEMP1=0.0 MAIN 510 DO 3110 L=1,LDIM MAIN 511 3110 TEMP1=TEMP1+RPOWER (X(I,L)-X(J,L)) MAIN 512 DIST(M)=RROOT (TEMP1) MAIN 513 DBAR=DBAR+DIST(M)*WW(M) MAIN 514 SUMW = SUMW + WW(M) MAIN 515 3120 CONTINUE MAIN 516 DBAR=DBAR/SUMW MAIN 517 C DBAS IS EITHER 0 OR DBAR ACCORDING TO WHETHER MAIN 518 C STRESS FORMULA 1 OR 2 IS BEING USED. MAIN 519 DBAS = 0.0 MAIN 520 IF(LSCH.EQ.2) DBAS = DBAR MAIN 521 IF(IABS(SDSWIT).GE.10) MAIN 522 1 CALL FITM(MX,MY,LFITSW) MAIN 523 IF(IABS(SDSWIT).LT.10) MAIN 524 1 CALL FITP(MX,MY,SDSWIT,LFITRM) MAIN 525 C MAIN 526 C CALCULATE U, T, AND STRESS MAIN 527 C MAIN 528 3200 U=0.0 MAIN 529 T=0.0 MAIN 530 DO 3210 M=MX,MY MAIN 531 U=U+(DIST(M)-DHAT(M))**2*WW(M) MAIN 532 3210 T=T+(DIST(M)-DBAS)**2*WW(M) MAIN 533 3215 U=SQRT (U) MAIN 534 TEMP1=T MAIN 535 T=SQRT (T) MAIN 536 GRSTRS(NG) = U/T MAIN 537 STRESS = STRESS + GRSTRS(NG)**2 MAIN 538 IF(U.EQ.0.0) GO TO 3340 MAIN 539 3220 RUT=1.0/(U*T) MAIN 540 UOT3=U/(TEMP1*T) MAIN 541 C MAIN 542 C CALCULATE THE (NEGATIVE) GRADIENT MAIN 543 C MAIN 544 3300 DO 3330 M = MX,MY MAIN 545 IF(DIST(M).EQ.0.0) GO TO 3330 MAIN 546 LTEMP1=MOD(IJ(M),KPACK) MAIN 547 I=LTEMP1/IJPACK+1 MAIN 548 J=MOD(LTEMP1,IJPACK)+1 MAIN 549 FACTOR=UOT3*(DIST(M)-DBAS)-RUT*(DIST(M)-DHAT(M)) MAIN 550 FACTOR = (FACTOR/RM1POW(DIST(M)) ) * WW(M) MAIN 551 FACTOR = FACTOR * GRSTRS(NG) MAIN 552 DO 3320 L=1,LDIM MAIN 553 TEMP1 = FACTOR * RM1POW(X(I,L)-X(J,L)) MAIN 554 GR(I,L)=GR(I,L)+TEMP1 MAIN 555 3320 GR(J,L)=GR(J,L)-TEMP1 MAIN 556 3330 CONTINUE MAIN 557 3340 CONTINUE MAIN 558 IF(STRESS .EQ. 0.0 ) GO TO 3700 MAIN 559 STRESS = SQRT( STRESS / FNGRPS ) MAIN 560 C MAIN 561 C AVOID MOVING FIXED POINTS MAIN 562 C MAIN 563 IF( LNFIXZ .EQ. 0) GO TO 3400 MAIN 564 DO 3360 I=1,LNFIXZ MAIN 565 DO 3360 L=1,LDIM MAIN 566 3360 GR(I,L) = 0.0 MAIN 567 C MAIN 568 C FIND SCALE FACTOR OF GRADIENT, ANGLE-COSINE BETWEEN GRADIENT MAIN 569 C AND PREVIOUS GRADIENT. MAIN 570 C MAIN 571 3400 SFGR=0.0 MAIN 572 CAGRGL=0.0 MAIN 573 DO 3410 I=1,N MAIN 574 DO 3410 L=1,LDIM MAIN 575 SFGR=SFGR+GR(I,L)**2 MAIN 576 3410 CAGRGL=CAGRGL+GR(I,L)*GL(I,L) MAIN 577 SFGR=SQRT (SFGR/FN) MAIN 578 C IF GRADIENT 0.0, SKIP AHEAD. MAIN 579 IF(SFGR) 3420,3700,3420 MAIN 580 3420 TEMP1=SFGR*SFGL*FN MAIN 581 CAGRGL=CAGRGL/TEMP1 MAIN 582 C MAIN 583 IF(ITNO.EQ.0 .OR. NBAKUP.GE.4) GO TO 3500 MAIN 584 IF(CAGRGL.GT.(-0.95) .AND. STRESS/STRLST.LT. 1.2001 ) GOTO 3500MAIN 585 C MAIN 586 C BACK UP ALONG LAST GRADIENT MAIN 587 C MAIN 588 NBAKUP = NBAKUP + 1 MAIN 589 IF(NBAKUP.EQ.1) STBASC=1.0 MAIN 590 STBASC=STBASC*STBAMU MAIN 591 TEMP1 = STEP MAIN 592 STEP = STEP / 10.0 MAIN 593 WRITE (LPRINT,38) STRESS, CAGRGL, STEP MAIN 594 38 FORMAT(10X, F7.3, 14X, F7.3, 22X, F8.4 ) MAIN 595 TEMP1 = (TEMP1 - STEP) / SFGL MAIN 596 TEMP1=STBASC*TEMP1 MAIN 597 DO 3430 I = 1,N MAIN 598 DO 3430 L = 1,LDIM MAIN 599 X(I,L) = X(I,L) - TEMP1*GL(I,L) MAIN 600 3430 GR(I,L) = GL(I,L) MAIN 601 SFGR = SFGL MAIN 602 STRESS = STRLST MAIN 603 GO TO 3000 MAIN 604 C MAIN 605 C STEP SIZE CALCULATIONS MAIN 606 C MAIN 607 3500 IF(ITNO) 9999, 3510, 3520 MAIN 608 3510 SRAT=0.8 MAIN 609 GO TO 3530 MAIN 610 3520 SRAT=STRESS/STRLST MAIN 611 3530 CALL NEWSTP( STEP, ITNO, SFGR, STRESS, MAIN 612 1 CAGRGL, COSAV, ACSAV, COSAVW, ACSAVW, SRAT, SRATAV ) MAIN 613 NBAKUP = 0 MAIN 614 C MAIN 615 C PRINT CURRENT STATUS OF COMPUTATION MAIN 616 C MAIN 617 3700 IF(LHIPRT.EQ.2) WRITE (LPRINT,30) ITNO,STRESS,SRAT,SRATAV, MAIN 618 1 CAGRGL,COSAV,ACSAV,SFGR,STEP MAIN 619 30 FORMAT(I10,F7.3,F7.3,F7.3,F7.3,F7.3,F7.3,F8.4, F8.4) MAIN 620 C MAIN 621 C DECIDE WHETHER TO CONTINUE ITERATING MAIN 622 C MAIN 623 3800 IF(STRESS) 9999, 3840, 3810 MAIN 624 3810 IF(SFGR-SFGRMN) 3850, 3850, 3815 MAIN 625 3815 TEMP1 = 0.5 * (1.0+SRATST) MAIN 626 TEMP2 = 1.0 - TEMP1 MAIN 627 IF( ABS (SRAT-TEMP1) - TEMP2 ) 3816, 3816, 3820 MAIN 628 3816 IF( ABS (SRATAV-TEMP1) - TEMP2 ) 3850, 3850, 3820 MAIN 629 3820 IF(STRESS-STRMIN) 3860, 3860, 3830 MAIN 630 3830 IF(ITNO-NOIT) 3900, 3870, 9999 MAIN 631 3840 CONTINUE MAIN 632 WRITE (LPRINT, 31) MAIN 633 31 FORMAT(24H0ZERO STRESS WAS REACHED ) MAIN 634 GO TO 4000 MAIN 635 3850 CONTINUE MAIN 636 WRITE (LPRINT, 32) MAIN 637 32 FORMAT(21H0MINIMUM WAS ACHIEVED ) MAIN 638 GO TO 4000 MAIN 639 3860 CONTINUE MAIN 640 WRITE (LPRINT, 33) MAIN 641 33 FORMAT(32H0SATISFACTORY STRESS WAS REACHED ) MAIN 642 GO TO 4000 MAIN 643 3870 CONTINUE MAIN 644 WRITE (LPRINT, 34) MAIN 645 34 FORMAT(39H0MAXIMUM NUMBER OF ITERATIONS WERE USED ) MAIN 646 GO TO 4000 MAIN 647 C MAIN 648 C CONTINUE ITERATING MAIN 649 C MAIN 650 3900 ITNO=ITNO+1 MAIN 651 TEMP1=STEP/SFGR MAIN 652 DO 3910 I=1,N MAIN 653 DO 3910 L=1,LDIM MAIN 654 3910 X(I,L)=X(I,L)+TEMP1*GR(I,L) MAIN 655 GO TO 3000 MAIN 656 C MAIN 657 C STOP ITERATING *************************************MAIN 658 C MAIN 659 C ROTATE FINAL CONFIGURATION TO PRINCIPAL COMPONENTS MAIN 660 C MAIN 661 4000 IF(LCOORS.LT.2) GO TO 4002 MAIN 662 C COMPUTE X TRANSPOSE TIMES X AND STORE UPPER HALF IN RVEC MAIN 663 KK=0 MAIN 664 DO 4410 K=1,LDIM MAIN 665 DO 4410 J=1,K MAIN 666 SUM=0.0 MAIN 667 DO 4405 I=1,N MAIN 668 4405 SUM=SUM+X(I,J)*X(I,K) MAIN 669 KK=KK+1 MAIN 670 4410 RVEC(KK)=SUM MAIN 671 C COMPUTE MATRIX OF EIGENVECTORS, GL MAIN 672 CALL SGEV(LDIM,LDIM,GL) MAIN 673 C COMPUTE X TIMES GL MAIN 674 DO 4420 K=1,LDIM MAIN 675 DO 4420 J=1,N MAIN 676 SUM=0.0 MAIN 677 DO 4415 I=1,LDIM MAIN 678 4415 SUM=SUM+X(J,I)*GL(I,K) MAIN 679 4420 GR(J,K)=SUM MAIN 680 DO 4425 I=1,N MAIN 681 DO 4425 J=1,LDIM MAIN 682 4425 X(I,J)=GR(I,J) MAIN 683 WRITE(LPRINT,152) MAIN 684 152 FORMAT(66H0THE FINAL CONFIGURATION HAS BEEN ROTATED TO PRINCIPAL CMAIN 685 1OMPONENTS.) MAIN 686 GO TO 4004 MAIN 687 C MAIN 688 4002 IF(LCOORS.EQ.1) WRITE(LPRINT,153) MAIN 689 153 FORMAT(110H0THE CONFIGURATION HAS BEEN NORMALIZED DURING THE ITERAMAIN 690 1TIONS BUT HAS NOT BEEN ROTATED TO PRINCIPAL COMPONENTS.) MAIN 691 IF(LCOORS.EQ.0) WRITE(LPRINT,154) MAIN 692 154 FORMAT(70H0NO NORMALIZATION WAS DONE TO THE CONFIGURATION DURING TMAIN 693 1HE ITERATIONS.) MAIN 694 4004 WRITE (LPRINT, 40)N,LDIM,STRESS,LSCH,(L,L=1,LDIM) MAIN 695 40 FORMAT(27H0THE FINAL CONFIGURATION OF,I4, MAIN 696 1 10H POINTS IN,I3, 22H DIMENSIONS HAS STRESS,F7.3,8H FORMULA ,I2 MAIN 697 2 /1X/29HLABEL FOR CONFIGURATION PLOTS,10X,19HFINAL CONFIGURATION, MAIN 698 3 /,39X,10I7) MAIN 699 IF(LPUNSW.EQ.2) GO TO 4005 MAIN 700 WRITE (LPUNCH,41) (TITLE(I),I=1,80),N,LDIM MAIN 701 41 FORMAT (14H CONFIGURATION/80A1/2I3) MAIN 702 WRITE (LPUNCH,52) MAIN 703 52 FORMAT (12H (2X,10F7.3)) MAIN 704 4005 DO 4010 I=1,N MAIN 705 WRITE(LPRINT,42) PTID(I),I,(X(I,L),L=1,LDIM) MAIN 706 IF(LPUNSW.EQ.2) GO TO 4010 MAIN 707 WRITE (LPUNCH, 43)I,(X(I,L),L=1,LDIM) MAIN 708 4010 CONTINUE MAIN 709 42 FORMAT(20X,A1,18X,I2,10F7.3) MAIN 710 43 FORMAT(I2,10F7.3) MAIN 711 WRITE (LPRINT,46) MAIN 712 DO 4020 NG=1,NOGRPS MAIN 713 MZ = GRNO(NG+1) - GRNO(NG) MAIN 714 SDSWT1=MOD(GRSDIS(NG),100) +1 MAIN 715 IF(SDSWT1-11) 4016,4013,4014 MAIN 716 4013 WRITE(LPRINT,60) NG,MZ,GRSTRS(NG) MAIN 717 60 FORMAT(1X,I5,2X,I5,F7.3,11H ASCENDING) MAIN 718 GO TO 4020 MAIN 719 4014 WRITE(LPRINT,62) NG,MZ,GRSTRS(NG) MAIN 720 62 FORMAT(1X,I5,2X,I5,F7.3,11H DESCENDING) MAIN 721 GO TO 4020 MAIN 722 4016 WRITE(LPRINT,61) NG,MZ,GRSTRS(NG),(PCOEFF(I),I=1,SDSWT1) MAIN 723 61 FORMAT(1X,I5,2X,I5,F7.3,12H POLYNOMIAL,5F15.7) MAIN 724 4020 CONTINUE MAIN 725 46 FORMAT(14H0DATA GROUP(S) ,/73H0SERIAL COUNT STRESS REGRESSION COEMAIN 726 .FFICIENTS (FROM DEGREE 0 TO MAX OF 4)) MAIN 727 4030 CONTINUE MAIN 728 IF(LDIPRT.EQ.2) MAIN 729 1 CALL DATAPR(GRNO,MM,NOGRPS,2) MAIN 730 C MAIN 731 C PLOTTING SECTION MAIN 732 C MAIN 733 C PLOT DIST AND DHAT VS. DATA MAIN 734 IF(LPLSCT.EQ.0) GO TO 5100 MAIN 735 WRITE(LPRINT,44) LDIM,LSCH,STRESS,TITLE MAIN 736 44 FORMAT(52H1DIST(D) AND DHAT(-) (Y-AXIS) VS. DATA (X-AXIS), FOR,I3,MAIN 737 . 27H DIMENSIONS. STRESS,FORMULA,I2,2H,=F8.4,/,30X,80A1) MAIN 738 PTID(1)=DEE MAIN 739 PTID(2)=AIE MAIN 740 CALL PLOT(DATA,PMAT,0.,0.,0.,0.,MM,2,-1,1800,1) MAIN 741 C PLOT ADDITIONAL SCATTER DIAGRAMS BY GROUPS--FIRST 5 GROUPS MAIN 742 C IF LPLSCT=1, ALL GROUPS IF LPLSCT=2 MAIN 743 IF(NOGRPS.EQ.1) GO TO 5090 MAIN 744 KPLT=NOGRPS MAIN 745 IF(LPLSCT.EQ.1) KPLT=MIN0(KPLT,5) MAIN 746 DO 5060 NG=1,KPLT MAIN 747 MX=GRNO(NG) MAIN 748 MY=GRNO(NG+1)-1 MAIN 749 MZ=MY-MX+1 MAIN 750 WRITE(LPRINT,4062) NG,LSCH,GRSTRS(NG),LDIM,TITLE MAIN 751 4062 FORMAT(43H1DIST(D) AND DHAT(-) VS. DATA FOR GROUP NO.,I3,17H. STRMAIN 752 .ESS,FORMULA,I2,2H,=,F8.4,12H DIMENSION=,I3,/,30X,80A1) MAIN 753 5060 CALL PLOT(DATA,PMAT,0.,0.,0.,0.,MZ,2,-1,1800,MX) MAIN 754 5090 CONTINUE MAIN 755 PTID(1)=AAA MAIN 756 PTID(2)=BEE MAIN 757 C MAIN 758 C PLOT CONFIGURATION MAIN 759 5100 IF(LPLCON.EQ.0) GO TO 5200 MAIN 760 IF(LDIM-1) 5190,5190,5125 MAIN 761 5125 IF(LPLCON.EQ.1) GO TO 5160 MAIN 762 C PLOT ALL PAIRS OF DIMENSIONS MAIN 763 DO 5145 J=2,LDIM MAIN 764 JMIN1=J-1 MAIN 765 DO 5145 I=1,JMIN1 MAIN 766 WRITE(LPRINT,5015) J,I,TITLE MAIN 767 5015 FORMAT(31H1CONFIGURATION PLOT DIMENSION,I2,23H (Y-AXIS) VS. DIMEMAIN 768 .NSION,I3,9H (X-AXIS),/,30X,80A1) MAIN 769 5145 CALL PLOT(X(1,I),X(1,J),0.,0.,0.,0.,N,1,2,100,1) MAIN 770 GO TO 5200 MAIN 771 C PLOT PAIRS OF DIMENSIONS SO THAT ALL DIMENSIONS ARE PLOTTED ONCMAIN 772 5160 DO 5180 J=2,LDIM,2 MAIN 773 IF((LDIM.EQ.5).AND.(J.EQ.4)) GO TO 5185 MAIN 774 JMIN1=J-1 MAIN 775 WRITE(LPRINT,5015) J,JMIN1,TITLE MAIN 776 CALL PLOT(X(1,JMIN1),X(1,J),0.,0.,0.,0.,N,1,2,100,1) MAIN 777 IF((J.NE.2) .OR. (MOD(LDIM,2).EQ.0)) GO TO 5180 MAIN 778 JP1=J+1 MAIN 779 WRITE(LPRINT,5015) JP1,JMIN1,TITLE MAIN 780 CALL PLOT(X(1,JMIN1),X(1,JP1),0.,0.,0.,0.,N,1,2,100,1) MAIN 781 5180 CONTINUE MAIN 782 GO TO 5200 MAIN 783 5185 WRITE(LPRINT,5015) LDIM,J,TITLE MAIN 784 CALL PLOT(X(1,J),X(1,LDIM),0.,0.,0.,0.,N,1,2,100,1) MAIN 785 GO TO 5200 MAIN 786 C ONE DIMENSIONAL CONFIGURATION PLOT MAIN 787 5190 WRITE(LPRINT,5195) TITLE MAIN 788 5195 FORMAT(1H1,20X,34HPLOT OF CONFIGURATION, DIMENSION 1,/,30X,80A1) MAIN 789 CALL PLOT(X(1,1),X(1,1),0.,0.,0.,0.,N,1,3,100,1) MAIN 790 C MAIN 791 C CHANGE DIMENSION MAIN 792 C MAIN 793 5200 ISTC=ISTC+1 MAIN 794 STPL(ISTC)= STRESS MAIN 795 DIMSV(ISTC)=LDIM MAIN 796 4100 LDIM=LDIM-LDIMD MAIN 797 IF(LDIM-LDIMN)4110,2300,2300 MAIN 798 C MAIN 799 C PLOT STRESS VERSUS DIMENSION IF MORE THAN TWO POINTS MAIN 800 C MAIN 801 4110 CONTINUE MAIN 802 IF(ISTC.LT.3) GO TO 100 MAIN 803 WRITE(LPRINT, 50) TITLE MAIN 804 50 FORMAT (1H1,10X,31HPLOT OF STRESS VERSUS DIMENSION, /,30X,80A1) MAIN 805 C INVERT TO IMPROVE READING OF PLOT MAIN 806 NVT=ISTC/2 MAIN 807 DO 8050 I=1,NVT MAIN 808 IC=ISTC+1-I MAIN 809 TEMP = DIMSV(I) MAIN 810 SEMP = STPL(I) MAIN 811 DIMSV(I) = DIMSV(IC) MAIN 812 STPL(I) =STPL(IC) MAIN 813 DIMSV(IC) = TEMP MAIN 814 STPL(IC) = SEMP MAIN 815 8050 CONTINUE MAIN 816 C FIND MAX STRESS FOR Y-AXIS , AND FILL IN PLOTTING CHARACTERS MAIN 817 YMA=0.0 MAIN 818 DO 8231 I=1,ISTC MAIN 819 IDIM=DIMSV(I) MAIN 820 PTID(I)=XNUM(IDIM) MAIN 821 8231 IF(YMA .LE. STPL(I)) YMA=STPL(I) MAIN 822 CALL PLOT(DIMSV,STPL,10.,YMA,0.,0.,ISTC,1,-2,10,1) MAIN 823 PTID(1)=AAA MAIN 824 PTID(2)=BEE MAIN 825 PTID(3)=CEE MAIN 826 PTID(4)=DEE MAIN 827 PTID(5)=EEE MAIN 828 PTID(6)=FFF MAIN 829 C MAIN 830 C REINITIALIZE, AND RETURN FOR MORE INPUT. MAIN 831 GO TO 100 MAIN 832 C MAIN 833 C NORMAL TERMINATION, AFTER READING STOP ONCONTROL CARD MAIN 834 9000 STOP MAIN 835 C MAIN 836 C TROUBLE EXIT MAIN 837 C MAIN 838 9999 WRITE (LPRINT, 99) MAIN 839 99 FORMAT(52H0KRUSKAL. IMPOSSIBLE BRANCH TAKEN FROM IF STATEMENT. ) MAIN 840 STOP MAIN 841 C MAIN 842 END MAIN 843 $ FORTRAN CBLOCK DATA BLOCK 1 BLOCK DATA BLOCK 2 C BLDA FOR KYST JANUARY,1973 BLOCK 3 C WRITTEN BY J.KRUSKAL BLOCK 4 C MODIFIED FOR KYST BY J.KRUSKAL AND J.SEERY JANUARY,1973 BLOCK 5 C BLOCK 6 C THIS PROGRAM HOLDS THE TABLE OF WORDS WHICH CCACT CONSULTS BLOCK 7 C LOGICALLY IT CONSISTS OF SEVERAL TABLES BLOCK 8 C THE VALUES IN MTAB INDICATE WHICH ROWS START NEW TABLES BLOCK 9 C EACH ENTRY IN THE TABLE CONTAINS 4 ITEMS BLOCK 10 C THE FIRST IS THE CODED ALPHABETIC WORD BLOCK 11 C (ONE WORD MAY BE ENTERED SEVERAL TIMES ) BLOCK 12 C THE SECOND INDICATES THE NATURE OF THIS ENTRY BLOCK 13 C 1 INDICATES INTEGER PARAMETERS WHICH MUST BE EXPLIC9TLY BLOCK 14 C READ IN BLOCK 15 C 2 INDICATES REAL PARAMETERS WHICH MUST BE EXPLICITLY BLOCK 16 C READ IN BLOCK 17 C 3 INDICATES AN INTEGER PARAMETER WHICH IS IMPLICITLY BLOCK 18 C DEFINED BLOCK 19 C 4 INDICATES A REAL PARAMETER WHICH IS IMPLICITLY BLOCK 20 C DEFINED BLOCK 21 C 5 INDICATES A PARAMETER WHICH BELONGS TO CCACT ONLY BLOCK 22 C 6 INDICATES AN IMPLICITLY DEFINED INTEGER PARAMETER BLOCK 23 C EFFECTING PROGRAM EXECUTION UPON RETURN FROM CCACT BLOCK 24 C THE THIRD INDICATES WHICH PARAMTER IS INVOLVED BLOCK 25 C SEPARATE NUMBERING FOR MAIN PROGRAM PARAMETERS AND CCACT BLOCK 26 C THE FOURTH GIVES THE VALUE FOR ANY INTERNALLY DEFINED PARAMETERS BLOCK 27 C IF THE PARAMETER IS OF TYPE 3,4 OR 6.IF OF TYPE 1 OR 2, THIS BLOCK 28 C ENTRY GIVES THE NUMBER OF ITEMS TO BE READ IN. BLOCK 29 C BLOCK 30 INTEGER LPAR( 42 ) BLOCK 31 REAL PAR( 42 ) BLOCK 32 EQUIVALENCE (LPAR,PAR) BLOCK 33 C BLOCK 34 INTEGER MTAB(13) BLOCK 35 INTEGER TAB1(4,16),TAB2(4,17),TAB3(4,16),TAB4(4,19),TAB5(4,18) BLOCK 36 REAL PTID(100),ITEM(101) BLOCK 37 C BLOCK 38 COMMON /ACCUR/ PRECSN,XMAG,XMAXN BLOCK 39 COMMON /MDSCL1/ LREAD, LPRINT, LPUNCH, LSCRAT BLOCK 40 COMMON /MDSCL2/ LPAR BLOCK 41 COMMON /MDSCL3/ MTAB,TAB1,TAB2,TAB3,TAB4,TAB5 BLOCK 42 COMMON /PLTCHR/ PTID,ITEM BLOCK 43 C BLOCK 44 DATA PRECSN,XMAG,XMAXN /1.5E-8,1.0E-38,1.0E38/ , BLOCK 45 . LREAD,LPRINT,LPUNCH,LSCRAT /5,6,7,8/ BLOCK 46 DATA MTAB /1,50,52,57,63,65,74,76,78,81,84,87,87/ BLOCK 47 C BLOCK 48 C TAB1 ENTRIES PERTAIN SUCCESSIVELY TO THESE CONTROL WORDS BLOCK 49 C DIMMAX, DIMMIN, DIMDIF, CUTOFF, STRMIN, BLOCK 50 C SFGRMN, COSAVW, ACSAVW, DIAGON, MATRIX, BLOCK 51 C HALFMA, LOWERH, UPPERH, LOWERC, ARBITR, BLOCK 52 C TORSCA BLOCK 53 DATA TAB1/ BLOCK 54 1 6989, 1, 1, 1, BLOCK 55 2 5375, 1, 2, 1, BLOCK 56 3 6923, 1, 3, 1, BLOCK 57 4 13520, 2, 4, 1, BLOCK 58 5 390, 2, 5, 1, BLOCK 59 6 2553, 2, 6, 1, BLOCK 60 7 7303, 2, 7, 1, BLOCK 61 8 11226, 2, 8, 1, BLOCK 62 9 11136, 5, 1, 2, BLOCK 63 A 9277, 3, 10, 1, BLOCK 64 B 3527, 3, 10, 2, BLOCK 65 C 6069, 3, 10, 2, BLOCK 66 D 8938, 3, 10, 3, BLOCK 67 E 1754, 3, 10, 4, BLOCK 68 F 10499, 3, 20, -99, BLOCK 69 G 2048, 3, 27, 4/ BLOCK 70 C BLOCK 71 C TAB2 ENTRIES PERTAIN SUCCESSIVELY TO THESE CONTROL WORDS BLOCK 72 C UPPERC, BLOCKD, DISSIM, DISSIM, SIMILA, BLOCK 73 C SIMILA, CONFIG, COMPUT, PRIMAR, SECOND, BLOCK 74 C R , ITERAT, FIX , SRATST, STOP , BLOCK 75 C PRE-IT, COORDI BLOCK 76 DATA TAB2/ BLOCK 77 1 4623, 3, 10, 5, BLOCK 78 2 8683, 5, 1, 5, BLOCK 79 3 15096, 6, 12, 2, BLOCK 80 4 15096, 3, 11, 10, BLOCK 81 5 6868, 6, 12, 2, BLOCK 82 6 6868, 3, 11, 11, BLOCK 83 7 14846, 6, 12, 5, BLOCK 84 8 11754, 6, 12, 6, BLOCK 85 9 765, 3, 13, 1, BLOCK 86 A 4119, 3, 13, 2, BLOCK 87 B 16326, 2, 14, 1, BLOCK 88 C 15170, 1, 15, 1, BLOCK 89 D 1855, 1, 28, 1, BLOCK 90 E 3720, 2, 16, 1, BLOCK 91 F 13171, 6, 12, 7, BLOCK 92 G 10738, 1, 39, 1, BLOCK 93 H 7261, 5, 1, 11/ BLOCK 94 C BLOCK 95 C TAB3 ENTRIES PERTAIN SUCCESSIVELY TO THESE CONTROL WORDS BLOCK 96 C WEIGHT, WFUNCT, SFORM1, SFORM2, CARDS, BLOCK 97 C NOCARD, SPLIT , DATA , SAVE , RANDOM, BLOCK 98 C PRINT , DFUNCT, REGRES, DCONST, WCONST, BLOCK 99 C PLOT BLOCK100 DATA TAB3/ BLOCK101 1 14543, 6, 12, 3, BLOCK102 2 11427, 6, 12, 4, BLOCK103 3 5318, 3, 17, 1, BLOCK104 4 6181, 3, 17, 2, BLOCK105 5 6927, 3, 18, 1, BLOCK106 6 16009, 3, 18, 2, BLOCK107 7 1966, 5, 1, 3, BLOCK108 8 6675, 6, 12, 2, BLOCK109 9 9189, 5, 1, 7, BLOCK110 A 8330, 1, 20, 1, BLOCK111 B 2775, 5, 1, 4, BLOCK112 C 10575, 6, 12, 8, BLOCK113 D 14439, 5, 1, 6, BLOCK114 E 267, 2, 29, 5, BLOCK115 F 1119, 2, 34, 5, BLOCK116 G 6878, 5, 1, 8/ BLOCK117 C BLOCK118 C TAB4 ENTRIES PERTAIN SUCCESSIVELY TO THESE CONTROL WORDS BLOCK119 C ABSENT, PRESEN, BYROWS, BYGROU, BYDECK, BLOCK120 C NOMORE, NOMORE, DATA , DISTAN, HISTOR, BLOCK121 C NODATA, NODIST, NOHIST, NO , YES , BLOCK122 C MONOTO, ASCEND, DESCEN, POLYNO BLOCK123 DATA TAB4/ BLOCK124 1 4258, 3, 9, 2, BLOCK125 2 2575, 3, 9, 1, BLOCK126 3 7761, 3, 19, 100, BLOCK127 4 11782, 3, 19, 200, BLOCK128 5 11288, 3, 19, 300, BLOCK129 6 5274, 3, 19, 400, BLOCK130 7 5274, 3, 19, 2, BLOCK131 8 6675, 3, 21, 2, BLOCK132 9 9824, 3, 22, 2, BLOCK133 A 12801, 3, 24, 2, BLOCK134 B 16057, 3, 21, 1, BLOCK135 C 5863, 3, 22, 1, BLOCK136 D 9395, 3, 24, 1, BLOCK137 E 9622, 3, 10, 100, BLOCK138 F 11125, 3, 10, 200, BLOCK139 G 15634, 3, 23, 0, BLOCK140 H 7782, 3, 23, 10, BLOCK141 I 11188, 3, 23, 11, BLOCK142 J 3756, 3, 26, 1/ BLOCK143 C BLOCK144 C TAB5 ENTRIES PERTAIN SUCCESSIVELY TO THESE CONTROL WORDS BLOCK145 C POLYNO, CONSTA, NOCONS, MULTIV, MULTIV, BLOCK146 C DATA , CONFIG, CONFIG, SCATTE, ALL , BLOCK147 C SOME , NONE , ALL , SOME , NONE , BLOCK148 C ROTATE, STANDA, AS-IS BLOCK149 DATA TAB5/ BLOCK150 1 3756, 1, 23, 1, BLOCK151 2 14387, 3, 26, 100, BLOCK152 3 5018, 3, 26, 200, BLOCK153 4 3608, 3, 26, 0, BLOCK154 5 3608, 1, 23, 1, BLOCK155 6 6675, 3, 25, 2, BLOCK156 7 14846, 3, 27, 3, BLOCK157 8 14846, 5, 1, 9, BLOCK158 9 11109, 5, 1, 10, BLOCK159 A 5766, 3, 40, 2, BLOCK160 B 13660, 3, 40, 1, BLOCK161 C 10008, 3, 40, 0, BLOCK162 D 5766, 3, 41, 2, BLOCK163 E 13660, 3, 41, 1, BLOCK164 F 10008, 3, 41, 0, BLOCK165 G 4503, 3, 42, 2, BLOCK166 H 3418, 3, 42, 1, BLOCK167 I 9499, 3, 42, 0/ BLOCK168 C BLOCK169 DATA PTID /1HA,1HB,1HC,1HD,1HE,1HF,1HG,1HH,1HI,1HJ,1HK,1HL,1HM, BLOCK170 1 1HN,1HO,1HP,1HQ,1HR,1HS,1HT,1HU,1HV,1HW,1HX,1HY, BLOCK171 2 1HZ,1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9,1HA,1HB,1HC, BLOCK172 3 1HD,1HE,1HF,1HG,1HH,1HI,1HJ,1HK,1HL,1HM,1HN,1HO, BLOCK173 4 1HP,1HQ,1HR,1HS,1HT,1HU,1HV,1HW,1HX,1HY,1HZ,1H1,1H2, BLOCK174 5 1H3,1H4,1H5,1H6,1H7,1H8,1H9,1HA,1HB,1HC,1HD,1HE, BLOCK175 6 1HF,1HG,1HH,1HI,1HJ,1HK,1HL,1HM,1HN,1HO,1HP,1HQ, BLOCK176 7 1HR,1HS,1HT,1HU,1HV,1HW,1HX,1HY,1HZ,1H1,1H2,1H3,1H4/ BLOCK177 DATA ITEM/101*1H / BLOCK178 END BLOCK179 $ FORTREX STAB CBSEC1 BSEC1 1 SUBROUTINE BSEC1(IB,N,ETA,A,B,W,E,S) BSEC1 2 C BSEC1 FOR KYST JANUARY,1973 BSEC1 3 C WRITTEN AND TECHNIQUES ADVOCATED IN BSEC1 4 C REFERENCE IMPLEMENTED FOR SGEV BSEC1 5 C BY P.BUSINGER JANUARY,1972 BSEC1 6 DIMENSION A(1),B(1),W(1) BSEC1 7 REAL ETA,E,S BSEC1 8 C BSEC1 9 C W.KAHAN, ACCURATE EIGENVALUES OF A SYM- BSEC1 10 C METRIC TRIDIAGONAL MATRIX. TECH.REPORT BSEC1 11 C NO.CS41, JULY 22,1966, COMP.SC.DEPT., BSEC1 12 C STANFORD UNIVERSITY. BSEC1 13 C BSEC1 14 C BSEC1 IS CALLED BY SGEV, AND BSEC1 15 C USES IB BISECTION STEPS TO FIND THE ALGEBRAICALLY BSEC1 16 C GREATEST EIGENVALUE E OF THE SYMMETRIC BSEC1 17 C TRIDIAGONAL MATRIX WHOSE DIAGONAL ELEMENTS ARE BSEC1 18 C A(1) THROUGH A(N) BSEC1 19 C AND WHOSE SUB-AND-SUPERDIAGONAL ELEMENTS ARE BSEC1 20 C B(1) THROUGH B(N-1) BSEC1 21 C FURTHERMORE, BSEC1 22 C W(1) THROUGH W(N) BSEC1 23 C ARE SET EQUAL TO THE DIAGONAL ELEMENTS OF THE BSEC1 24 C MATRIX U IN THE LU-DECOMPOSITION REQUIRED BSEC1 25 C FOR SUBSEQUENT INVERSE ITERATION. A AND B BSEC1 26 C ARE RESCALED, THE CHANGE IF SCALE BEING BSEC1 27 C REFLECTED IN S. ETA IS LEAST POSITIVE MA- BSEC1 28 C CHINE NUMBER. BSEC1 29 C BSEC1 30 REAL BB,D,T,X BSEC1 31 IF(N.GT.1)GOTO 5 BSEC1 32 E=A(1) BSEC1 33 W(1)=ETA BSEC1 34 GOTO 70 BSEC1 35 5 B(N)=0.E0 BSEC1 36 T=0.E0 BSEC1 37 DO 10 I=1,N BSEC1 38 T=AMAX1(T,ABS(A(I))) BSEC1 39 10 T=AMAX1(T,ABS(B(I))) BSEC1 40 IF(T.NE.0.E0)GOTO 20 BSEC1 41 E=0.E0 BSEC1 42 DO 19 I=1,N BSEC1 43 19 W(I)=ETA BSEC1 44 GOTO 70 BSEC1 45 20 S=S*T BSEC1 46 DO 30 I=1,N BSEC1 47 A(I)=A(I)/T BSEC1 48 30 B(I)=B(I)/T BSEC1 49 E=3.E0 BSEC1 50 D=-E BSEC1 51 DO 50 K=1,IB BSEC1 52 X=(D+E)/2.E0 BSEC1 53 U=1.E0 BSEC1 54 NU=0 BSEC1 55 BB=0.E0 BSEC1 56 DO 40 I=1,N BSEC1 57 U=(A(I)-BB/U)-X BSEC1 58 IF(U.GT.ETA)GOTO 35 BSEC1 59 U=AMIN1(U,-ETA) BSEC1 60 NU=NU+1 BSEC1 61 35 BB=B(I)**2 BSEC1 62 40 W(I)=U BSEC1 63 IF(NU.LT.N)D=X BSEC1 64 IF(NU.EQ.N)E=X BSEC1 65 50 CONTINUE BSEC1 66 IF(NU.EQ.N)GOTO 70 BSEC1 67 U=1.E0 BSEC1 68 BB=0.E0 BSEC1 69 DO 60 I=1,N BSEC1 70 U=(A(I)-BB/U)-E BSEC1 71 IF(U.GT.ETA)GOTO 55 BSEC1 72 U=AMIN1(U,-ETA) BSEC1 73 55 BB=B(I)**2 BSEC1 74 60 W(I)=U BSEC1 75 70 RETURN BSEC1 76 END BSEC1 77 $ FORTREX STAB CCACT CACT 1 SUBROUTINE CCACT CACT 2 C CCACT FOR KYST JANUARY,1973 CACT 3 C WRITTEN BY J.KRUSKAL CACT 4 C MODIFIED FOR KYST BY J.KRUSKAL AND J.SEERY JANUARY,1973 CACT 5 C CCACT--CONTROL CARD ACTIVITY. READS AND INTERPRETS CONTROL CARDS.CACT 6 C CACT 7 REAL MAP(40),DUMMY(27),ATAB(4,86) CACT 8 INTEGER TAB(4,86),IFAC(6) CACT 9 INTEGER MTAB(13) CACT 10 EQUIVALENCE (ATAB,TAB) CACT 11 C CACT 12 INTEGER LPAR( 42 ) CACT 13 REAL PAR( 42 ) CACT 14 EQUIVALENCE (MAP(1),DUMMY(1)), (MAP(28),CHTAB(1)) CACT 15 EQUIVALENCE (LPAR,PAR) CACT 16 C CACT 17 INTEGER IPARAM(5) CACT 18 REAL C(81) CACT 19 REAL WORD(18) CACT 20 REAL CHTAB(13) CACT 21 C CACT 22 INTEGER DFLTSW CACT 23 REAL BLANK, DOT CACT 24 INTEGER BLSW,NUMSW,DECSW,TYPE,XTYPE,T,TA, PARNO, TABNO CACT 25 C CACT 26 COMMON /MDSCL1/ LREAD, LPRINT, LPUNCH, LSCRAT CACT 27 COMMON /MDSCL2/ LPAR CACT 28 COMMON /MDSCL3/ MTAB, TAB CACT 29 C CACT 30 DATA BLANK, DOT, EQUALS, COMMA, C(81) /1H ,1H.,1H=,1H,,1H / CACT 31 DATA DOLLAR/1H$/ CACT 32 DATA MODP,IFAC/ 16381, 907, 887, 883, 881, 877, 863/ CACT 33 C CACT 34 DATA MAP/1H ,1HA,1HB,1HC,1HD,1HE,1HF,1HG,1HH,1HI,1HJ,1HK,1HL, CACT 35 . 1HM,1HN,1HO,1HP,1HQ,1HR,1HS,1HT,1HU,1HV,1HW,1HX,1HY, CACT 36 . 1HZ,1H0,1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9,1H-,1H+,1H./CACT 37 C CACT 38 C CACT 39 C CACT 40 1 FORMAT(80A1) CACT 41 11 FORMAT(1H0,80A1) CACT 42 2 FORMAT(1X,I1) CACT 43 3 FORMAT(1X,80A1) CACT 44 4 FORMAT(1X,I18) CACT 45 5 FORMAT(1X,F18.3) CACT 46 9 FORMAT(51H CCACT DIAGNOSTIC. CONTROL CARD ABOVE IS IMPROPER.) CACT 47 99 FORMAT( 12H ITEM NUMBER, I5, 20H HAS EXPECTED TYPE , I5, CACT 48 1 18H AND ACTUAL TYPE , I5, 16H . THIS ITEM IS , 18A1) CACT 49 98 FORMAT(12H ITEM NUMBER,I5,39H HAS ILLEGAL CHARACTER. THIS ITEM ISCACT 50 . ,18A1) CACT 51 C CACT 52 C READ AND PRINT CONTROL CARD CACT 53 C CACT 54 100 READ (LREAD,1) (C(I),I=1,80) CACT 55 WRITE(LPRINT,11) (C(I),I=1,80) CACT 56 C ANALYZE CONTROL CARD INTO TOKENS CACT 57 C CACT 58 C EACH TOKEN IS DELIMITED BY BLANKS CACT 59 C THERE ARE FOUR TYPES OF TOKENS. CACT 60 C ALPHABETIC, INTEGER, DECIMAL, AND DEFAULT CACT 61 C ALPHABETIC UNLESS FIRST CHARACTER IS DIGIT OR DEC POINT CACT 62 C OR PLUS OR MINUS OR DOLLAR SIGN CACT 63 C DECIMALS DISTINGUISHED BY DECIMAL POINT CACT 64 C DEFAULT = $ CACT 65 C CACT 66 REWIND LSCRAT CACT 67 BLSW=1 CACT 68 NUMSW=1 CACT 69 DECSW=1 CACT 70 DFLTSW=1 CACT 71 K=0 CACT 72 C CACT 73 DO 300 I=1,81 CACT 74 C CACT 75 X=C(I) CACT 76 GO TO (110,120), BLSW CACT 77 C CACT 78 110 IF(X.EQ.BLANK .OR. X.EQ.EQUALS .OR. X.EQ. COMMA ) GO TO 300CACT 79 BLSW=2 CACT 80 JA=I CACT 81 DO 115 KX=1,13 CACT 82 115 IF(X.EQ.CHTAB(KX))NUMSW=2 CACT 83 IF(X.EQ.DOLLAR) DFLTSW=2 CACT 84 IF(X.EQ.DOT) DECSW=2 CACT 85 GO TO 300 CACT 86 C CACT 87 120 IF(X.EQ.BLANK .OR. X.EQ.EQUALS .OR. X.EQ. COMMA ) GO TO 130CACT 88 IF(X.EQ.DOT) DECSW=2 CACT 89 GO TO 300 CACT 90 C CACT 91 130 K=K+1 CACT 92 JB = MIN0 (I-1,JA+16) CACT 93 JC=18-(JB-JA+1) CACT 94 TYPE=1 CACT 95 IF(NUMSW.EQ.2) TYPE=NUMSW+DECSW-1 CACT 96 IF(DFLTSW.EQ.2) TYPE=4 CACT 97 C CACT 98 WRITE (LSCRAT,2) TYPE CACT 99 IF(TYPE.EQ.4) GO TO 160 CACT 100 GO TO (140,150),NUMSW CACT 101 140 WRITE (LSCRAT,3) (C(J),J=JA,JB), (BLANK,J=1,JC) CACT 102 GO TO 160 CACT 103 150 WRITE (LSCRAT,3) (BLANK,J=1,JC), (C(J),J=JA,JB) CACT 104 160 BLSW=1 CACT 105 NUMSW=1 CACT 106 DECSW=1 CACT 107 DFLTSW=1 CACT 108 GO TO 300 CACT 109 C CACT 110 300 CONTINUE CACT 111 C CACT 112 C ANALYZE TOKENS AND SET PARAMETER VALUES ACCORDINGLY CACT 113 C CACT 114 KB=K CACT 115 IF(KB.EQ.0) RETURN CACT 116 REWIND LSCRAT CACT 117 XTYPE=1 CACT 118 IPARAM(1) = 1 CACT 119 NT=0 CACT 120 NTOKNS=0 CACT 121 NOPCC=0 CACT 122 C CACT 123 DO 1000 K=1,KB CACT 124 C CACT 125 READ (LSCRAT,2) TYPE CACT 126 IF( (XTYPE.EQ.1) .AND. (TYPE.NE.1) ) GO TO 995 CACT 127 IF( (XTYPE.NE.TYPE). AND. (TYPE.NE.4) ) GO TO 995 CACT 128 350 GO TO (400,410,420,430), TYPE CACT 129 C CACT 130 400 READ (LSCRAT,3) (WORD(L),L=1,18) CACT 131 C CONVERT FIRST SIX LETTERS OF ALPHABETIC WORD INTO CODE NUMBER. CACT 132 C LARGEST POSSIBLE CODE NUMBER IS LESS THAN 2**15. CACT 133 ICODE=0 CACT 134 DO 408 M1=1,6 CACT 135 X=WORD(M1) CACT 136 DO 402 M2=1,38 CACT 137 IF( X.EQ.MAP(M2) ) GO TO 404 CACT 138 402 CONTINUE CACT 139 C ILLEGAL CHARACTER CACT 140 WRITE (LPRINT,9) CACT 141 WRITE (LPRINT,98) K,(WORD(L),L=1,18) CACT 142 GO TO 999 CACT 143 404 M3=(M2-1)*IFAC(M1) CACT 144 IF( M3.GE.MODP ) M3=M3-MODP CACT 145 ICODE=ICODE+M3 CACT 146 IF( ICODE.GE.MODP ) ICODE=ICODE-MODP CACT 147 408 CONTINUE CACT 148 GO TO 510 CACT 149 C CACT 150 410 READ (LSCRAT,4) INTPAR CACT 151 ISUB=PARNO+NT CACT 152 LPAR(ISUB)=INTPAR CACT 153 GO TO 430 CACT 154 C CACT 155 420 READ (LSCRAT,5) DECPAR CACT 156 ISUB=PARNO+NT CACT 157 PAR(ISUB)=DECPAR CACT 158 C CACT 159 430 NT=NT+1 CACT 160 IF(NT.EQ.NTOKNS) XTYPE=1 CACT 161 GO TO 1000 CACT 162 C CACT 163 510 TABNO = IPARAM(1) CACT 164 MA = MTAB(TABNO) CACT 165 MB = MTAB(TABNO+1) - 1 CACT 166 C CACT 167 LMISSW = 0 CACT 168 DO 700 M=MA,MB CACT 169 IF( ICODE.NE.TAB(1,M) ) GO TO 700 CACT 170 LMISSW = 1 CACT 171 600 XTYPE=1 CACT 172 NT=0 CACT 173 IPARAM(1) = 1 CACT 174 PARNO=TAB(3,M) CACT 175 LTEMP=TAB(2,M) CACT 176 NTOKNS=TAB(4,M) CACT 177 IF(LTEMP.GT.2) NTOKNS=0 CACT 178 GO TO(610,620,630,640,650,660), LTEMP CACT 179 C CACT 180 C NAME OF INTEGER PARAMETER CACT 181 610 XTYPE=2 CACT 182 GO TO 700 CACT 183 C CACT 184 C NAME OF DECIMAL PARAMETER CACT 185 620 XTYPE=3 CACT 186 GO TO 700 CACT 187 C CACT 188 C IMPLICITLY SPECIFIED INTEGER PARAMETER CACT 189 C A SINGLE IMPLICITLY SPECIFIED PARAMTER CAN IN PRINCIPLE HOLD CACT 190 C AS MANY AS THREE PARAMTERS WHICH THE MAIN PROGRAM CONSIDERS CACT 191 C CONCEPTUALLY DISTINCT. ONE IS IN THE UNITS POSITION, ANOTHER IN CACT 192 C IN THE HUNDREDS POSITION, AND ANOTHER IN THE TEN-THOUSANDS CACT 193 C POSITION. CACT 194 630 LP = 1 CACT 195 LT=TAB(4,M) CACT 196 IF(LT.GE.100) LP = 100 CACT 197 IF(LT.GE.10000) LP = 10000 CACT 198 LQ = 100*LP CACT 199 LA = LPAR(PARNO) CACT 200 LPAR(PARNO)= LA-(MOD(LA,LQ)/LP)*LP +TAB(4,M) CACT 201 GO TO 700 CACT 202 C CACT 203 C IMPLICITLY SPECIFIED DECIMAL PARAMETER CACT 204 640 PAR(PARNO)=ATAB(4,M) CACT 205 GO TO 700 CACT 206 C CACT 207 C INTERNAL PARAMETER OF CCACT PROGRAM CACT 208 650 IPARAM(PARNO)=TAB(4,M) CACT 209 GO TO 700 CACT 210 C CACT 211 C IMPLICITLY SPECIFIED INTEGER PARAMETER CACT 212 C DETERMINES PROGRAM FLOW UPON RETURN FROM CCACT--SO ONLY ONE OF CACT 213 C THIS TYPE PER CONTROL CARD ALLOWED. CACT 214 660 IF(NOPCC)9999,665,996 CACT 215 665 NOPCC=1 CACT 216 LPAR(PARNO)=TAB(4,M) CACT 217 GO TO 700 CACT 218 C CACT 219 700 CONTINUE CACT 220 IF(LMISSW.EQ.0) GO TO 994 CACT 221 1000 CONTINUE CACT 222 IF(NT.LT.NTOKNS) GO TO 997 CACT 223 C CACT 224 1001 RETURN CACT 225 C CACT 226 994 WRITE(LPRINT,9) CACT 227 WRITE(LPRINT,998) (WORD(L),L=1,18) CACT 228 998 FORMAT(27H UNDEFINED CONTROL PHRASE ,18A1) CACT 229 GO TO 999 CACT 230 995 WRITE (LPRINT,9) CACT 231 WRITE (LPRINT,99) K,XTYPE,TYPE, (WORD(L),L=1,18) CACT 232 GO TO 999 CACT 233 996 WRITE(LPRINT,9) CACT 234 WRITE(LPRINT,6) CACT 235 6 FORMAT(69H0NO MORE THAN ONE CONTROL PHRASE OF TYPE 6 ALLOWED ON A CACT 236 .CONTROL CARD.) CACT 237 GO TO 999 CACT 238 997 WRITE(LPRINT,9) CACT 239 WRITE(LPRINT,97) CACT 240 97 FORMAT(47H CONTROL PHRASE NOT COMPLETED ON A SINGLE CARD.) CACT 241 999 CALL EXIT CACT 242 C CACT 243 9999 WRITE(LPRINT,7) CACT 244 7 FORMAT(23H0ERROR EXIT FROM CCACT.) CACT 245 GO TO 999 CACT 246 END CACT 247 $ FORTREX STAB CCONFIG CONFIG 1 SUBROUTINE CONFIG(N,ND) CONFIG 2 C CONFIG FOR KYST JANUARY,1973 CONFIG 3 C WRITTEN BY F.YOUNG CONFIG 4 C MODIFIED FOR KYST BY J.KRUSKAL AND J.SEERY JANUARY,1973 CONFIG 5 C ROUTINE TO OBTAIN A METRIC MULTIDIMENSIONAL SCALING SOLUTION CONFIG 6 C CONFIG 7 DIMENSION DATA(1800),WW(1800),IJ(1800),X(100,6) CONFIG 8 DIMENSION STORE(494),RWMEAN(100),ROOTS(6),XBEST(100,6), CONFIG 9 . DATAIN(5430) CONFIG10 REAL MEAN CONFIG11 COMMON /KYST1/ DATA,WW,IJ,X CONFIG12 COMMON /KYST2/ STORE,RWMEAN,ROOTS,XBEST,DATAIN CONFIG13 C CONFIG14 C PARAMETERS CONFIG15 C X - CONFIGURATION AT OUTPUT CONFIG16 C N - NUMBER OF POINTS IN SCALING CONFIG17 C ND - NUMBER OF DIMENSIONS OF SOLUTION CONFIG18 C DATAIN - ON INPUT, DATA FROM INICON. DESTROYED BY CONFIG CONFIG19 C RWMEAN - ROW AVERAGES CONFIG20 C ROOTS - ON RETURN FROM SGEV, THE FIRST ND EIGENVALUES CONFIG21 C CONFIG22 ISUB(I1,I2)=((I2-1)*(I2-2))/2+I1 CONFIG23 MEAN=0.0 CONFIG24 FN=N CONFIG25 C CONFIG26 C COMPUTE ROW,COLUMN, AND OVERALL AVERAGES OF DATIN VALUES SQUARED CONFIG27 C CONFIG28 DO 130 J=1,N CONFIG29 RWMEAN(J)=0.0 CONFIG30 DO 129 I=1,N CONFIG31 IF(I-J) 126,129,127 CONFIG32 126 IJC=ISUB(I,J) CONFIG33 GO TO 128 CONFIG34 127 IJC=ISUB(J,I) CONFIG35 128 RWMEAN(J)=DATAIN(IJC)**2+RWMEAN(J) CONFIG36 129 CONTINUE CONFIG37 130 MEAN=RWMEAN(J)+MEAN CONFIG38 MEAN=MEAN/FN**2 CONFIG39 C CONFIG40 C CONVERT DATAIN VALUES TO SCALAR-PRODUCT-LIKE VALUES BY THE CONFIG41 C HOUSEHOLDER-TORGERSON PROCEDURE. WORK FROM LAST ROW AND CONFIG42 C COLUMN BACKWARDS. CONFIG43 C CONFIG44 DO 135 JJ=1,N CONFIG45 J=N-JJ+1 CONFIG46 JM1=J-1 CONFIG47 J1=(JM1*(J-2))/2 CONFIG48 J2=(J*JM1)/2 CONFIG49 DO 135 II=1,J CONFIG50 I=J-II+1 CONFIG51 IJ1=J1+I CONFIG52 IJ2=J2+I CONFIG53 IF(I.EQ.J) GO TO 133 CONFIG54 C HOUSEHOLDER-TORGERSON FORMULAS CONFIG55 DATAIN(IJ2)=((RWMEAN(I)+RWMEAN(J))/FN-MEAN-DATAIN(IJ1)**2)/2.0 CONFIG56 GO TO 135 CONFIG57 133 DATAIN(IJ2)=((2.0*RWMEAN(I))/FN-MEAN)/2.0 CONFIG58 135 CONTINUE CONFIG59 C CONFIG60 C OBTAIN EIGENROOTS AND EIGENVECTORS CONFIG61 CALL SGEV(N,ND,X) CONFIG62 C CONFIG63 C COMPUTE CONFIGURATION CONFIG64 C CONFIG65 DO 136 I=1,ND CONFIG66 136 ROOTS(I)=SQRT(ABS(ROOTS(I))) CONFIG67 DO 137 I=1,N CONFIG68 DO 137 J=1,ND CONFIG69 137 X(I,J)=X(I,J)*ROOTS(J) CONFIG70 RETURN CONFIG71 END CONFIG72 $ FORTRAN CDATAPR DATAPR 1 SUBROUTINE DATAPR(GRNO,MM,NOGRPS,LCODE) DATAPR 2 C DATAPR FOR KYST JANUARY,1973 DATAPR 3 C WRITTEN BY J.KRUSKAL DATAPR 4 C DATAPR 5 C DATA PRINT DATAPR 6 C DATAPR 7 DIMENSION DATA(1800),WW(1800),IJ(1800),X(100,6) DATAPR 8 DIMENSION GR(100,6),GL(100,6),STORE(1830),DIST(1800),DHAT(1800), DATAPR 9 . P(50) DATAPR10 INTEGER Q(50),GRNO(1) DATAPR11 EQUIVALENCE(P(1),Q(1),STORE(1)) DATAPR12 COMMON /KYST1/ DATA,WW,IJ,X DATAPR13 COMMON /KYST2/ GR,GL,STORE,DIST,DHAT DATAPR14 COMMON /MDSCL1/ LREAD, LPRINT, LPUNCH, LSCRAT DATAPR15 DATA KPACK,IJPACK /10000,100/ DATAPR16 C DATAPR17 1 FORMAT( 1H1, 4(30H I J DATA WGHT GP NO ) ) DATAPR18 2 FORMAT(1H1, 3(40H I J DATA DIST DHAT WGHT GP NO ) ) DATAPR19 3 FORMAT( 1X, 5( 2I4, F9.3, F5.2, I3, I4, 1H, ) ) DATAPR20 4 FORMAT( 1X, 5( 2I4, F9.3, 2F5.2, F5.2, I3, I4, 1H, ) ) DATAPR21 C DATAPR22 100 NCOPA=3 DATAPR23 IF(LCODE.EQ.1) NCOPA=4 DATAPR24 C DATAPR25 MMLEFT=MM DATAPR26 200 MMNOW=MIN0( MMLEFT, 50*NCOPA ) DATAPR27 MMUSED=MM-MMLEFT DATAPR28 MMLEFT=MMLEFT-MMNOW DATAPR29 IF(LCODE.EQ.1) WRITE (LPRINT,1) DATAPR30 IF(LCODE.EQ.2) WRITE (LPRINT,2) DATAPR31 C DATAPR32 DO 300 LINENO=1,50 DATAPR33 NCONOW=(MMNOW-LINENO+50)/50 DATAPR34 IF(NCONOW.LE.0) GO TO 310 DATAPR35 C DATAPR36 L=0 DATAPR37 DO 270 LC=1,NCONOW DATAPR38 M = MMUSED + 50*(LC-1) + LINENO DATAPR39 C DATAPR40 IJM=MOD(IJ(M),KPACK) DATAPR41 I=IJM/IJPACK+1 DATAPR42 J=MOD(IJM,IJPACK)+1 DATAPR43 C DATAPR44 L=L+1 DATAPR45 Q(L)=I DATAPR46 L=L+1 DATAPR47 Q(L)=J DATAPR48 C DATAPR49 L=L+1 DATAPR50 P(L)=DATA(M) DATAPR51 C DATAPR52 IF(LCODE.EQ.1) GO TO 250 DATAPR53 C DATAPR54 L=L+1 DATAPR55 P(L)=DIST(M) DATAPR56 L=L+1 DATAPR57 P(L)=DHAT(M) DATAPR58 C DATAPR59 250 CONTINUE DATAPR60 C DATAPR61 L=L+1 DATAPR62 P(L)=WW(M) DATAPR63 C DATAPR64 DO 220 NG=1,NOGRPS DATAPR65 IF(M.GE.GRNO(NG)) NGM=NG DATAPR66 220 CONTINUE DATAPR67 L=L+1 DATAPR68 Q(L)=NGM DATAPR69 C DATAPR70 L=L+1 DATAPR71 Q(L)=M-GRNO(NGM)+1 DATAPR72 C DATAPR73 270 CONTINUE DATAPR74 C DATAPR75 LL=L DATAPR76 IF(LCODE.EQ.1) WRITE (LPRINT,3) (P(L),L=1,LL) DATAPR77 IF(LCODE.EQ.2) WRITE (LPRINT,4) (P(L),L=1,LL) DATAPR78 300 CONTINUE DATAPR79 C DATAPR80 310 IF(MMLEFT.GT.0) GO TO 200 DATAPR81 C DATAPR82 RETURN DATAPR83 C DATAPR84 END DATAPR85 $ FORTRAN CDFLT DFLT 1 SUBROUTINE DFLT(N,A,B,E,X) DFLT 2 C DFLT FOR KYST JANUARY,1973 DFLT 3 C WRITTEN AND DEFLATION TECHNIQUE DFLT 4 C ADVOCATED IN REFERENCE IMPLEMENTED DFLT 5 C FOR SGEV BY P.BUSINGER JANUARY,1972 DFLT 6 C DFLT 7 C G.GOLUB AND W.KAHAN, CALCULATING THE DFLT 8 C SINGULAR VALUES AND PSEUDO-INVERSE OF A DFLT 9 C MATRIX. J.SIAM NUMER. ANAL. 2,205-224(1965). DFLT 10 C DFLT IS CALLED BY SGEV, AND DFLT 11 REAL A(1),B(1),E,X(1) DFLT 12 C DFLT 13 C DEFLATES THE SYMMETRIC TRIDIAGONAL MATRIX DFLT 14 C WITH DIAGONAL ELEMENTS A(1) THROUGH A(N) DFLT 15 C AND SUB- AND SUPER-DIAGONAL ELEMENTS B(1) DFLT 16 C THROUGH B(N-1) INTO A SYMMETRIC TRIDIAGONAL DFLT 17 C MATRIX OF ORDER N-1. E IS THE EIGENVALUE DFLT 18 C BEING REMOVED AND X AN ASSOCIATED EIGEN- DFLT 19 C VECTOR. THE TANGENTS OF HALF THE ROTATION DFLT 20 C ANGLES OVERWRITE THE FIRST (N-1) COMPONENTS DFLT 21 C OF X. DFLT 22 C DFLT 23 REAL H,W,F,V,T,Q,S,C,AINEW DFLT 24 H=A(1)-E DFLT 25 W=B(1) DFLT 26 N1=N-1 DFLT 27 IF(N1.LT.1)GOTO 30 DFLT 28 DO 20 I=1,N1 DFLT 29 F=X(I) DFLT 30 V=X(I+1) DFLT 31 IF(AMAX1(ABS(H),ABS(W)).LE.AMAX1(ABS(F),ABS(V)))GOTO 10 DFLT 32 F=W DFLT 33 V=-H DFLT 34 10 T=0.E0 DFLT 35 IF(F*F.EQ.0.E0)GOTO 20 DFLT 36 T=F DFLT 37 IF(V.GE.0.E0)T=-T DFLT 38 T=T/(SQRT(F*F+V*V)+ABS(V)) DFLT 39 S=(T+T)/(1.E0+T*T) DFLT 40 C=(1.E0-T)*(1.E0+T)/(1.E0+T*T) DFLT 41 IF(I.GT.1)B(I-1)=C*H+S*W DFLT 42 F=C*A(I)+S*B(I) DFLT 43 Q=C*B(I)+S*A(I+1) DFLT 44 W=S*B(I+1) DFLT 45 B(I+1)=C*B(I+1) DFLT 46 AINEW=F*C+Q*S DFLT 47 H=Q*C-F*S DFLT 48 A(I+1)=A(I)+A(I+1)-AINEW DFLT 49 A(I)=AINEW DFLT 50 X(I+1)=C*X(I+1)-S*X(I) DFLT 51 20 X(I)=T DFLT 52 B(N-1)=0.E0 DFLT 53 30 RETURN DFLT 54 END DFLT 55 $ FORTRAN CDTRAN DTRAN 1 FUNCTION DTRAN(XX) DTRAN 2 C DTRAN FOR KYST JANUARY,1973 DTRAN 3 C WRITTEN FOR KYST BY J.SEERY JANUARY,1973 DTRAN 4 C DTRAN 5 C DTRAN--DATA TRANSFORMATION. DTRAN 6 DIMENSION DUMMY(28),DUM(4) DTRAN 7 COMMON /MDSCL2/ DUMMY,A,B,P,S,T,WCON1,WCON2,WCON3,WCON4,WCON5,DUM DTRAN 8 DTRAN=S+T*((A+B*XX)**P) DTRAN 9 RETURN DTRAN 10 END DTRAN 11 $ FORTRAN CEQSOLV EQSOLV 1 SUBROUTINE EQSOLV(A,B,C,NT,NF) EQSOLV 2 C EQSOLV FOR KYST JANUARY,1973 EQSOLV 3 C WRITTEN BY J.KRUSKAL EQSOLV 4 CEQSOLV SOLVES SIMULTANEOUS LINEAR EQUATIONS EQSOLV 5 REAL A(5,5), B(5), C(5) EQSOLV 6 C EQSOLV 7 NTM1=NT-1 EQSOLV 8 IF(NT.EQ.NF) GO TO 60 EQSOLV 9 DO 50 II=NF,NTM1 EQSOLV10 IIP1=II+1 EQSOLV11 DO 40 I=IIP1,NT EQSOLV12 E = A(I,II)/A(II,II) EQSOLV13 DO 30 J=NF,NT EQSOLV14 A(I,J)=A(I,J)-E*A(II,J) EQSOLV15 30 CONTINUE