C kyst2a.f C The authors of this software are Joseph B Kruskal, Forrest W Young, and C 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 SECOND MDS Package of AT&T Bell Laboratories C The manual of how to use this program is available in several different C formats: see mds/kyst2a_manual. C For explanation of the method this software implements, see C "Multidimensional Scaling" by Joseph B. Kruskal and Myron Wish, C a monograph published by Sage Publications, Beverly Hills, Calif. C in 1978 as series 07, item 011, in the Sage University Papers, and C "Multidimensional Scaling by Optimizing Goodness of Fit to a Nonmetric C Hypothesis" by Joseph B. Kruskal in Psychometrika in 1964, vol 29, C nos. 1 and 2, pages 1-27 and 115-129. C----+----@----+----@----+----@----+----@----+----@----+----@----+----@ CMAIN MAIN 1 C KYST-2A AUGUST, 1977 MAIN 2 C MAIN 3 C KRUSKAL-YOUNG-SHEPARD-TORGERSON MULTIDIMENSIONAL SCALING PROGRAM MAIN 4 C MAIN 5 C VERSION 1 JANUARY, 1973 MAIN 6 C VERSION 2 OCTOBER, 1976 MAIN 7 C VERSION 2A AUGUST, 1977 MAIN 8 C (THIS VERSION INCLUDES REVISED MAIN PROGRAM MENTIONED MAIN 9 C IN THE ERRATA SHEET OF JULY 1, 1977, AND INCLUDES THE MAIN 10 C CORRECTIONS FROM THE ERRATA SHEET OF JULY 31, 1977.) MAIN 11 C MAIN 12 C WRITTEN BY: MAIN 13 C MAIN 14 C DR. J. B. KRUSKAL MAIN 15 C BELL TELEPHONE LABORATORIES MAIN 16 C MURRAY HILL, N. J. MAIN 17 C MAIN 18 C DR. F. W. YOUNG MAIN 19 C PSYCHOMETRIC LABORATORY MAIN 20 C UNIVERSITY OF NORTH CAROLINA MAIN 21 C CHAPEL HILL, N. C. MAIN 22 C ASSISTED BY: MAIN 23 C JUDITH B. SEERY MAIN 24 C BELL TELEPHONE LABORATORIES MAIN 25 C MURRAY HILL, N. J. MAIN 26 C MAIN 27 C MAIN 28 C MAIN 29 C **********************************************************************MAIN 30 C MAIN 31 C MAIN 31 C GENERAL REMARKS. MAIN 32 C MAIN 33 C MAIN 34 C KYST INCLUDES THE FOLLOWING SUBROUTINES: MAIN 35 C MAIN 36 C BLDA EQSOLV PLOT SERCH MAIN 37 C BSEC1 FITM REGR SGEV MAIN 38 C CCACT FITMS RM1POW SORT MAIN 39 C CONFIG FITP RPOWER ST03 MAIN 40 C DATAPR INICON RROOT WTRAN MAIN 41 C DFLT NEWSTP RUNIFV XITEM MAIN 42 C DTRAN NRMLZ SBK NVIT1 MAIN 43 C MAIN 44 C ALL ROUTINES ARE WRITTEN IN FORTRAN IV, AND PUNCHED IN THE IBMEL 029 MAIN 45 C CHARACTER SET MAIN 46 C MAIN 47 C NO USE IS MADE OF SPECIAL OR NON-STANDARD SOFTWARE. MAIN 48 C MAIN 49 C ALL INPUT AND OUTPUT IS ONTO FORTRAN LOGICAL UNITS WITH NUMBERS MAIN 50 C CONTROLLED BY THESE VARIABLES: LREAD,LPRINT,LPUNCH,LSCRAT. MAIN 51 C UNIT NUMBERS 5,6,7,8 HAVE BEEN USED RESPECTIVELY. TO CHANGE THESEMAIN 52 C ASSIGNMENTS, CHANGE THE VALUES FOR THE FOUR VARIABLES SET IN THE MAIN 53 C BLOCK DATA SUBPROGRAM. MAIN 54 C MAIN 55 C THE SCRATCH UNIT IS USED IN A MINOR WAY BY CCACT ONLY. MANY MAIN 56 C INSTALLATIONS WILL HAVE ALTERNATE METHODS OF DOING THE SAME THING.MAIN 57 C MAIN 58 C KYST USES THREE MACHINE DEPENDENT CONSTANTS: PRECSN (RELATIVE MACHINE MAIN 59 C PRECISION), XMAG (LEAST POSITIVE MACHINE NUMBER), XMAXN (MINIMUM MAIN 60 C OF LARGEST POSITIVE MACHINE NUMBER AND MINUS LARGEST NEGATIVE MAIN 61 C MACHINE NUMBER), FLPRIN (FLOATING LINES PRINTED PER INCH OF PAPER)MAIN 62 C , CHPRIN (CHARACTERS PRINTED PER INCH ON A LINE). VALUES ARE MAIN 63 C 1.5E-8, 1.0E-38, 1.0E38, 6.0, AND 10.0 MAIN 64 C RESPECTIVELY. TO CHANGE THESE ASSIGNMENTS, CHANGE THE VALUES MAIN 65 C FOR THESE VARIABLES SET IN THE BLOCK DATA SUBPROGRAM. MAIN 66 C MAIN 67 C **********************************************************************MAIN 68 C MAIN 69 C MAIN 70 C MAIN PROGRAM FOR KYST JANUARY,1973 MAIN 71 C WRITTEN BY J.KRUSKAL MAIN 72 C MODIFIED FOR KYST BY J.KRUSKAL AND J.SEERY JANUARY,1973 MAIN 73 DIMENSION DATA(1800), IJ(1800), DIST(1800), DHAT(1800) MAIN 74 DIMENSION WW(1800),LDIST(1800) MAIN 75 DIMENSION X(100,6), GR(100,6), GL(100,6) MAIN 76 REAL PMAT(1800,2),RVEC(1800) MAIN 77 INTEGER GRNO(101), NOGRPS, LSPLIT, SPLITH MAIN 78 INTEGER GRSDIS(100), SDSWIT,SDSWT1 MAIN 79 INTEGER FMAT(80),PTID(100),XNUM(6),ONE,TWO,THREE,FOUR,FIVE,SIX MAIN 80 INTEGER AAA,BEE,CEE,DEE,FFF,AIE,EEE MAIN 81 INTEGER TITLE(80),CTITLE(80),LBLOCK(1800),ITEM(101) MAIN 82 REAL GRSTRS(100),PCOEFF(5),DUMMY(25) MAIN 83 DIMENSION STPL(10),DIMSV(10) MAIN 84 DIMENSION X2(100),Y2(100),PTID2(100),DUMEST(100) MAIN 85 INTEGER RTYPE MAIN 86 C MAIN 87 EXTERNAL WTRAN , DTRAN MAIN 88 C MAIN 89 COMMON /ACCUR/ PRECSN,XMAG,XMAXN,FLPRIN,CHPRIN MAIN 90 COMMON /KYST1/ DATA,WW,X,IJ MAIN 91 COMMON /KYST2/ GR,GL,RVEC,PCOEFF,DUMMY,PMAT MAIN 92 COMMON /MDSCL1/ LREAD, LPRINT, LPUNCH, LSCRAT MAIN 93 COMMON /MDSCL2/ MAIN 94 1 LDIMX, LDIMN, LDIMD, CUTOFF, STRMIN MAIN 95 2 , SFGRMN, COSAVW, ACSAVW, IFIRST, MATSWP MAIN 96 3 , SDSWIT, LCSWIT, LFITSW, RCOM, NOIT MAIN 97 4 , SRATST, LSCH, LPUNSW, LSPL, LRANDC MAIN 98 5 , LDAPRT, LDIPRT, LREG, LHIPRT, NUDASW MAIN 99 6 , LPOLSP, LCONSW, LNFIXZ, DCON1, DCON2 MAIN 100 7 , DCON3, DCON4, DCON5, WCON1, WCON2 MAIN 101 8 , WCON3, WCON4, WCON5, NOITIN, LPLCON MAIN 102 9 , LPLSCT, LCOORS, LDATAF,SPREAD MAIN 103 COMMON /METRIC/ RM1,RECR,R ,RTYPE MAIN 104 COMMON /PLTCHR/ PTID,ITEM MAIN 105 C MAIN 106 EQUIVALENCE (DIST(1),LDIST(1)) , (LBLOCK(1),RVEC(1)) MAIN 107 EQUIVALENCE (PMAT(1,1),DIST(1)), (PMAT(1,2),DHAT(1)) MAIN 108 EQUIVALENCE (XNUM(1),ONE),(XNUM(2),TWO),(XNUM(3),THREE), MAIN 109 . (XNUM(4),FOUR),(XNUM(5),FIVE),(XNUM(6),SIX), MAIN 110 1 (GR(1,1),X2(1)),(GR(1,2),Y2(1)),(GR(1,3),PTID2(1)) MAIN 111 C MAIN 112 DATA AAA,BEE,CEE,DEE,EEE,FFF /1HA,1HB,1HC,1HD,1HE,1HF/ MAIN 113 DATA ONE,TWO,THREE,FOUR,FIVE,SIX /1H1,1H2,1H3,1H4,1H5,1H6/ MAIN 114 DATA KPACK,IJPACK,NPREVZ /10000,100,0/ , AIE/1H-/ MAIN 115 DATA IVRSN /0/ MAIN 116 C MAIN 117 C INITIALIZE PARAMETERS MAIN 118 C MAIN 119 100 CONTINUE MAIN 120 C MAIN 121 C ALPHABETICAL ORDER MAIN 122 ACSAVW=0.66 MAIN 123 COSAVW=0.66 MAIN 124 CUTOFF=-1.23E+20 MAIN 125 DCON1=0.0 MAIN 126 DCON2=1.0 MAIN 127 DCON3=1.0 MAIN 128 DCON4=0.0 MAIN 129 DCON5=1.0 MAIN 130 GRNO(1) = 1 MAIN 131 IFIRST=1 MAIN 132 ISTC=0 MAIN 133 LCONSW=4 MAIN 134 LCOORS=2 MAIN 135 LCSWIT=1 MAIN 136 LDAPRT = 1 MAIN 137 LDATAF=LREAD MAIN 138 LDIMD=1 MAIN 139 LDIMN=2 MAIN 140 LDIMX=2 MAIN 141 LDIPRT = 1 MAIN 142 LFITSW=1 MAIN 143 LHIPRT = 2 MAIN 144 LNFIXZ = 0 MAIN 145 LPOLSP = 100 MAIN 146 LPLCON=2 MAIN 147 LPLSCT=3 MAIN 148 LPUNSW = 2 MAIN 149 LRANDC=-101 MAIN 150 LREG=0 MAIN 151 LSCH=1 MAIN 152 LSPL=302 MAIN 153 MATSWP=101 MAIN 154 NOIT=50 MAIN 155 NOITIN=1 MAIN 156 NUDASW=1 MAIN 157 R=2.0 MAIN 158 RCOM=-10.0 MAIN 159 SDSIN=1.0 MAIN 160 SDSWIT = 10 MAIN 161 SFGRMN=0.0 MAIN 162 SRATST=0.999 MAIN 163 SPREAD=10.0 MAIN 164 STRESS=1.0 MAIN 165 STRMIN=0.01 MAIN 166 WCON1=0.0 MAIN 167 WCON2=1.0 MAIN 168 WCON3=1.0 MAIN 169 WCON4=0.0 MAIN 170 WCON5=1.0 MAIN 171 WRITE(LPRINT, 17) MAIN 172 17 FORMAT(1H1) MAIN 173 IF(IVRSN.EQ.0) WRITE(LPRINT,9) MAIN 174 9 FORMAT(1H0,25X,43HOUTPUT FROM KYST2A--VERSION OF AUGUST, 1977) MAIN 175 IVRSN=1 MAIN 176 C MAIN 177 C CONTROL CARD READ MAIN 178 C MAIN 179 1000 LCSWIT=1 MAIN 180 CALL CCACT MAIN 181 GO TO (1000, 1100, 1165, 1182, 1200, 2000, 9000, 1190 ), LCSWIT MAIN 182 C MAIN 183 C DATA READ MAIN 184 C MAIN 185 1100 CONTINUE MAIN 186 NREPL1=1 MAIN 187 NREPL3 = 1 MAIN 188 LSPLIT=LSPL/100 MAIN 189 SPLITH = MOD(LSPL,100) MAIN 190 MATSW = MOD(MATSWP,100) MAIN 191 LBLKDS = MATSWP/100 MAIN 192 LPOSW = MOD(LPOLSP,100) MAIN 193 LFITRM = LPOLSP/100 MAIN 194 IF(LREG.EQ.0) GO TO 1104 MAIN 195 SDSWIT=LREG MAIN 196 IF(LREG.LT.10) SDSWIT = SDSWIT+LPOSW MAIN 197 1104 CONTINUE MAIN 198 SDSWIT = SDSWIT+100*LFITRM MAIN 199 C MAIN 200 IF(NUDASW.NE.1) MLASTD=MM MAIN 201 IF(NUDASW.NE.1) GO TO 1106 MAIN 202 NUDASW=2 MAIN 203 MLASTD=0 MAIN 204 MMIN=0 MAIN 205 NOGRPS=0 MAIN 206 1106 M=MLASTD MAIN 207 MA = M+1 MAIN 208 MSERA=MA MAIN 209 C MAIN 210 READ(LREAD,10) TITLE MAIN 211 WRITE(LPRINT,11) TITLE MAIN 212 IF(MATSW.GE.4) GO TO 1102 MAIN 213 1101 READ (LREAD,12) NPART,NREPL2,NREPL4 MAIN 214 WRITE (LPRINT,13) NPART,NREPL2,NREPL4 MAIN 215 NROWS = NPART MAIN 216 NCOLS = NPART MAIN 217 GO TO 1103 MAIN 218 1102 READ (LREAD,12) NROWS,NCOLS,NREPL2,NREPL4 MAIN 219 WRITE (LPRINT,13) NROWS,NCOLS,NREPL2,NREPL4 MAIN 220 NPART = NROWS+NCOLS MAIN 221 C MAIN 222 1103 N=NPART MAIN 223 IF(NREPL4.NE.0) NREPL3 = NREPL4 MAIN 224 IF(NREPL2.NE.0) NREPL1 = NREPL2 MAIN 225 IF(LBLKDS.EQ.2) N = NPART*NREPL3 MAIN 226 IF(NREPL2.EQ.0 .OR. NREPL4.EQ.0 ) WRITE (LPRINT,16) MAIN 227 1105 READ(LREAD,10) FMAT MAIN 228 WRITE(LPRINT,11) FMAT MAIN 229 C MAIN 230 DO 1162 NREPL=1,NREPL3 MAIN 231 IJBLKD=0 MAIN 232 IF(LBLKDS.EQ.2) IJBLKD =(NREPL-1)*NPART MAIN 233 MLASTG=M MAIN 234 IA = 1 MAIN 235 IB = NROWS MAIN 236 IF(MATSW.EQ.2) IA = IFIRST MAIN 237 IF(MATSW.EQ.3) IB = NROWS-(IFIRST-1) MAIN 238 C MAIN 239 DO 1160 I = IA,IB MAIN 240 MLASTR=M MAIN 241 LTEMP = NCOLS MAIN 242 IF(MATSW.EQ.2) LTEMP = I-IFIRST+1 MAIN 243 IF(MATSW.EQ.3) LTEMP = NCOLS-I-IFIRST+2 MAIN 244 MB = MA + LTEMP * NREPL1 - 1 MAIN 245 MSERB=MSERA+LTEMP*NREPL1-1 MAIN 246 ITRUE = I+IJBLKD MAIN 247 IF(MATSW.EQ.4) ITRUE = ITRUE+NCOLS MAIN 248 C MAIN 249 1130 READ (LDATAF,FMAT) (DATA(MP),MP=MA,MB) MAIN 250 C MAIN 251 DO 1150 MP=MA,MB MAIN 252 IF( DATA(MP)-CUTOFF ) 1150, 1150, 1140 MAIN 253 1140 CONTINUE MAIN 254 M = M + 1 MAIN 255 DATA(M)=DATA(MP) MAIN 256 WW(M) = 1.0 MAIN 257 J=((MP-MA)/NREPL1)+1 MAIN 258 IF(MATSW.EQ.3) J = J + (I-1) + (IFIRST-1) MAIN 259 J = J+IJBLKD MAIN 260 IF(MATSW.EQ.5) J = J + NROWS MAIN 261 IJ(M)=IJPACK*(ITRUE-1)+J-1 MAIN 262 LDIST(M) = MSERA+MP-MA MAIN 263 1150 CONTINUE MAIN 264 C MAIN 265 IF(LSPLIT.EQ.1 .AND. M.GT.MLASTR) GO TO 1155 MAIN 266 GO TO 1158 MAIN 267 1155 NOGRPS=NOGRPS+1 MAIN 268 GRNO(NOGRPS+1) = M+1 MAIN 269 GRSDIS(NOGRPS)=SDSWIT MAIN 270 1158 MA = M + 1 MAIN 271 1160 MSERA=MSERB+1 MAIN 272 C MAIN 273 IF(LSPLIT.EQ.2 .AND. M.GT.MLASTG) GO TO 1161 MAIN 274 GO TO 1162 MAIN 275 1161 NOGRPS=NOGRPS+1 MAIN 276 GRNO(NOGRPS+1) = M+1 MAIN 277 GRSDIS(NOGRPS)=SDSWIT MAIN 278 1162 CONTINUE MAIN 279 C MAIN 280 IF(LSPLIT.EQ.3 .AND. M.GT.MLASTD) GO TO 1163 MAIN 281 IF(LSPLIT.EQ.4 .AND. SPLITH.EQ.2 .AND. M.GT.MLASTD) GO TO 1163 MAIN 282 GO TO 1164 MAIN 283 1163 NOGRPS=NOGRPS+1 MAIN 284 GRSDIS(NOGRPS)=SDSWIT MAIN 285 C MAIN 286 1164 CONTINUE MAIN 287 GRNO(NOGRPS+1) = M+1 MAIN 288 IF(LSPLIT.LE.3) SPLITH=2 MAIN 289 IF(LSPLIT.EQ.4) SPLITH=1 MAIN 290 IF((NOGRPS.EQ.1).OR.(MMIN.EQ.0)) MMIN=M MAIN 291 C MAIN 292 LSPL = 100 * LSPLIT + SPLITH MAIN 293 C MAIN 294 MM = M MAIN 295 LDATAF=LREAD MAIN 296 GO TO 1000 MAIN 297 C MAIN 298 C WEIGHTS READ MAIN 299 C MAIN 300 1165 CONTINUE MAIN 301 M=MLASTD MAIN 302 MA=M+1 MAIN 303 MSERA=MA MAIN 304 DO 1181 NREPL=1,NREPL3 MAIN 305 IA = 1 MAIN 306 IB = NROWS MAIN 307 IF(MATSW.EQ.2) IA = IFIRST MAIN 308 IF(MATSW.EQ.3) IB = NROWS-(IFIRST-1) MAIN 309 DO 1180 I = IA,IB MAIN 310 LTEMP = NCOLS MAIN 311 IF(MATSW.EQ.2) LTEMP = I-IFIRST+1 MAIN 312 IF(MATSW.EQ.3) LTEMP = NCOLS-I-IFIRST+2 MAIN 313 MB = MA + LTEMP*NREPL1-1 MAIN 314 MSERB=MSERA+LTEMP*NREPL1-1 MAIN 315 1172 READ (LDATAF,FMAT) (WW(MP),MP=MA,MB) MAIN 316 DO 1177 MP=MA,MB MAIN 317 IF(LDIST(M+1).GT.(MSERA+MP-MA)) GO TO 1177 MAIN 318 M = M + 1 MAIN 319 WW(M)=WW(MP) MAIN 320 1177 CONTINUE MAIN 321 MA = M + 1 MAIN 322 1180 MSERA=MSERB+1 MAIN 323 1181 CONTINUE MAIN 324 LDATAF = LREAD MAIN 325 GO TO 1000 MAIN 326 C MAIN 327 C WEIGHT FORMATION BY WFUNCTION MAIN 328 C MAIN 329 1182 CONTINUE MAIN 330 WRITE(LPRINT,18) WCON1,WCON2,WCON3,WCON4,WCON5 MAIN 331 18 FORMAT(1H0,10X,69H** WEIGHTS FORMED ACCORDING TO RULE, WW(I,J)=S+TMAIN 332 .*((A+B*DATA(I,J))**P)/ 15X, 9H WHERE A=,F15.7,5X,2HB=,F15.7,5X, MAIN 333 .2HP=,F15.7,5X,2HS=,F15.7,5X,2HT=,F15.7) MAIN 334 MA=MLASTD+1 MAIN 335 DO 1185 M=MA,MM MAIN 336 TEMP1=DATA(M) MAIN 337 WW(M) = WTRAN(TEMP1) MAIN 338 1185 CONTINUE MAIN 339 GO TO 1000 MAIN 340 C MAIN 341 C DATA TRANSFORMATION BY DFUNCTION MAIN 342 C MAIN 343 1190 CONTINUE MAIN 344 WRITE(LPRINT,19) DCON1,DCON2,DCON3,DCON4,DCON5 MAIN 345 19 FORMAT(1H0,10X,73H** DATA TRANSFORMED ACCORDING TO RULE, DATA(I,J)MAIN 346 .=S+T*((A=B*DATA(I,J))**P) / 15X, 9H WHERE A=,F15.7,5X,2HB=, MAIN 347 .F15.7,5X,2HP=,F15.7,5X,2HS=,F15.7,5X,2HT=,F15.7) MAIN 348 MA=MLASTD+1 MAIN 349 DO 1195 M=MA,MM MAIN 350 DATA(M)=DTRAN(DATA(M)) MAIN 351 1195 CONTINUE MAIN 352 GO TO 1000 MAIN 353 C MAIN 354 C CONFIGURATION READ MAIN 355 C MAIN 356 1200 LCONSW=2 MAIN 357 READ(LREAD,10) CTITLE MAIN 358 WRITE(LPRINT,11)CTITLE MAIN 359 READ (LREAD, 12) NCON, LDIMCO MAIN 360 WRITE (LPRINT, 13) NCON, LDIMCO MAIN 361 READ(LREAD,10) FMAT MAIN 362 WRITE(LPRINT,11) FMAT MAIN 363 DO 1210 I=1,NCON MAIN 364 READ (LDATAF,FMAT) (X(I,L),L=1,LDIMCO) MAIN 365 WRITE(LPRINT,4008) I,(X(I,L),L=1,LDIMCO) MAIN 366 4008 FORMAT(1X,I2,10F7.3) MAIN 367 1210 CONTINUE MAIN 368 LDATAF = LREAD MAIN 369 GO TO 1000 MAIN 370 C MAIN 371 C SOME INPUT FORMATS MAIN 372 C MAIN 373 10 FORMAT (80A1) MAIN 374 11 FORMAT (1H0,80A1) MAIN 375 12 FORMAT(24I3) MAIN 376 13 FORMAT(1X,24I3) MAIN 377 15 FORMAT(20F4.0) MAIN 378 16 FORMAT(118H THE NUMBER OF INTERNAL REPLICATES, OR THE NUMBER OF EMAIN 379 1XTERNAL REPLICATES HAS NOT BEEN SPECIFIED. IT IS TAKEN TO BE 1 . )MAIN 380 C MAIN 381 C COMPUTATION *************************************MAIN 382 2000 CONTINUE MAIN 383 C MAIN 384 FN=FLOAT (N) MAIN 385 SQRTN=SQRT (FN) MAIN 386 FNGRPS = FLOAT(NOGRPS) MAIN 387 LDIM=LDIMX MAIN 388 IF(RCOM.NE.(-10.0)) R=RCOM MAIN 389 RR=R MAIN 390 RTYPE=3 MAIN 391 IF(R.EQ.1.0) RTYPE=1 MAIN 392 IF(R .EQ.2.0) RTYPE=2 MAIN 393 RM1=R-1.0 MAIN 394 RECR=1.0/R MAIN 395 IF(LDAPRT.EQ.2) MAIN 396 1 CALL DATAPR(GRNO,MM,NOGRPS,1) MAIN 397 WRITE (LPRINT,23) MAIN 398 C MAIN 399 C FINISH STARTING CONFIGURATION. LCONSW = 2,3,4. MAIN 400 C 4=USE INICON, 3=SAVED, 2=WAS READ IN. MAIN 401 C TO FILL IN ADDITIONAL POINTS, IF NEEDED: MAIN 402 C LRANDC<0, USE ARBITRARY, ELSE USE RANDOM. MAIN 403 C MAIN 404 2100 CONTINUE MAIN 405 WRITE(LPRINT,11) TITLE MAIN 406 IF((LCONSW.EQ.4) .AND. (LRANDC.EQ.(-101)))GO TO 2140 MAIN 407 ISTCON = 1 MAIN 408 IF( LCONSW .EQ. 2 ) ISTCON = NCON + 1 MAIN 409 IF( LCONSW .EQ. 3 ) ISTCON = NPREVZ + 1 MAIN 410 IF( ISTCON .GT. N ) GO TO 2200 MAIN 411 IF(LRANDC.LE.0) GO TO 2110 MAIN 412 DO 2105 K=1,LRANDC MAIN 413 2105 TEMP=RUNIFV(1.0) MAIN 414 2110 TEMP1=0.01 MAIN 415 DO 2130 I=ISTCON,N MAIN 416 DO 2120 L=1,LDIMX MAIN 417 2120 X(I,L)=0.0 MAIN 418 K= MOD (I-1,LDIM) +1 MAIN 419 X(I,K)=TEMP1 MAIN 420 IF(LRANDC.LT.0) GO TO 2130 MAIN 421 DO 2125 L=1,LDIMX MAIN 422 TEMP = RUNIFV(1.0) MAIN 423 2125 X(I,L) = ALOG( TEMP/(1.0-TEMP)) MAIN 424 2130 TEMP1=TEMP1+0.01 MAIN 425 GO TO 2200 MAIN 426 2140 SDSWIT=MOD(GRSDIS(1),100) MAIN 427 IF(SDSWIT.EQ.11) SDSIN=-1.0 MAIN 428 IF(N.GT.60) NOITIN=0 MAIN 429 CALL INICON(N,LDIM,MMIN,NOITIN,LHIPRT,SDSIN,LSCH,SPREAD) MAIN 430 C MAIN 431 C SORT DATA AND IJ AND WW. MAIN 432 C ALSO RECORD BLOCKS OF EQUAL DATA VALUES. MAIN 433 C MAIN 434 2200 CONTINUE MAIN 435 NPREVZ = N MAIN 436 DO 2250 NG = 1,NOGRPS MAIN 437 MX = GRNO(NG) MAIN 438 MY = GRNO(NG+1) - 1 MAIN 439 MZ = MY - MX + 1 MAIN 440 SDSWIT = GRSDIS(NG) MAIN 441 SDSWIT = MOD(SDSWIT,100) MAIN 442 IF(SDSWIT.EQ.11)SDSWIT = -10 MAIN 443 CALL SORT( DATA(MX),MZ,IJ(MX),WW(MX),DUMMY,2,SDSWIT) MAIN 444 C MAIN 445 C 'SORT' WILL SORT THE MM ELEMNTS OF 'DATA' IN ALGEBRAIC MAIN 446 C ORDER, ASCENDING OR DESCENDING ACCORDING TO WHETHER MAIN 447 C SDSWIT IS + OR -. MAIN 448 C AT THE SAME TIME, THE ELEMENTS IN 'IJ' AND IN 'WW' WILL BMAIN 449 C REARRANGED IN EXACTLY THE SAME ORDER. THUS THE MAIN 450 C CORRESPONDENCE BETWEEN THE ELEMENTS OF 'DATA' AND 'IJ' MAIN 451 C AND 'WW' IS PRESERVED. MAIN 452 C MAIN 453 DO 2240 MB = MX,MY MAIN 454 IF((DATA(MB+1).EQ.DATA(MB)).AND.(MB.NE.MY)) GO TO 2240 MAIN 455 IJ(MB)=MOD(IJ(MB),KPACK) MAIN 456 IJ(MB)=IJ(MB)+KPACK MAIN 457 2240 CONTINUE MAIN 458 2250 CONTINUE MAIN 459 C MAIN 460 C START COMPUTATION IN CURRENT DIMENSION MAIN 461 C MAIN 462 2300 FLDIM=FLOAT (LDIM) MAIN 463 ITNO=0 MAIN 464 COSAV=0.0 MAIN 465 SRATAV=0.8 MAIN 466 ACSAV=0.0 MAIN 467 STEP = 0.0 MAIN 468 NBAKUP = 0 MAIN 469 C MAIN 470 C INITIALIZE GRADIENT MAIN 471 C MAIN 472 2400 TEMP1=SQRT (1.0/FLDIM) MAIN 473 DO 2410 I=1,N MAIN 474 DO 2410 L=1,LDIM MAIN 475 2410 GR(I,L)=TEMP1 MAIN 476 SFGR=SQRTN MAIN 477 C MAIN 478 C PRINT HEADING FOR INFORMATION ABOUT SCALING IN CURRENT DIMENSION MAIN 479 C MAIN 480 2500 WRITE (LPRINT, 20) N, MM, NOGRPS, LDIM MAIN 481 IF(LHIPRT.EQ.2) WRITE (LPRINT,21) MAIN 482 WRITE (LPRINT, 22) MAIN 483 20 FORMAT(28H0HISTORY OF COMPUTATION. N= , I4, MAIN 484 1 15H. THERE ARE , I6, MAIN 485 2 26H DATA VALUES, SPLIT INTO , I4, MAIN 486 3 27H LISTS. DIMENSION = , I4 ) MAIN 487 21 FORMAT(52H0ITERATION STRESS SRAT SRATAV CAGRGL COSAV ACSAV, MAIN 488 1 16H SFGR STEP ) MAIN 489 22 FORMAT(1X) MAIN 490 23 FORMAT(1H1) MAIN 491 C MAIN 492 C START CURRENT ITERATION *************************************MAIN 493 C MAIN 494 C NORMALIZE CONFIGURATION. MOVE AND CLEAR GRADIENT. MAIN 495 C MAIN 496 3000 TEMP1=1.0 MAIN 497 IF(LCOORS.EQ.0) GO TO 3035 MAIN 498 TEMP1=0.0 MAIN 499 DO 3030 L=1,LDIM MAIN 500 TEMP2=0.0 MAIN 501 DO 3010 I=1,N MAIN 502 3010 TEMP2=TEMP2+X(I,L) MAIN 503 TEMP2=TEMP2/FN MAIN 504 DO 3020 I=1,N MAIN 505 X(I,L)=X(I,L)-TEMP2 MAIN 506 3020 TEMP1=TEMP1+X(I,L)**2 MAIN 507 3030 CONTINUE MAIN 508 TEMP1=SQRT (FN/TEMP1) MAIN 509 3035 DO 3040 L=1,LDIM MAIN 510 DO 3040 I=1,N MAIN 511 X(I,L)=TEMP1*X(I,L) MAIN 512 GL(I,L)=TEMP1*GR(I,L) MAIN 513 3040 GR(I,L)=0.0 MAIN 514 SFGL=TEMP1*SFGR MAIN 515 C MAIN 516 STBAMU = TEMP1 MAIN 517 C LOOP THROUGH THE DATA GROUPS MAIN 518 C MAIN 519 STRLST = STRESS MAIN 520 STRESS = 0.0 MAIN 521 DO 3340 NG = 1,NOGRPS MAIN 522 MX = GRNO(NG) MAIN 523 MY = GRNO(NG+1) - 1 MAIN 524 MZ = MY - MX + 1 MAIN 525 SDSWIT = GRSDIS(NG) MAIN 526 LFITRM = SDSWIT/100 MAIN 527 SDSWIT = MOD(SDSWIT,100) MAIN 528 IF(SDSWIT.EQ.11)SDSWIT = -10 MAIN 529 C MAIN 530 C COMPUTE DISTANCES AND FIND BEST-FITTING MONOTONE PSEUDO-DISTANCESMAIN 531 C MAIN 532 SUMW = 0.0 MAIN 533 DBAR = 0.0 MAIN 534 DO 3120 M=MX,MY MAIN 535 LTEMP1=MOD(IJ(M),KPACK) MAIN 536 I=LTEMP1/IJPACK+1 MAIN 537 J=MOD(LTEMP1,IJPACK)+1 MAIN 538 TEMP1=0.0 MAIN 539 DO 3110 L=1,LDIM MAIN 540 3110 TEMP1=TEMP1+RPOWER (X(I,L)-X(J,L)) MAIN 541 DIST(M)=RROOT (TEMP1) MAIN 542 DBAR=DBAR+DIST(M)*WW(M) MAIN 543 SUMW = SUMW + WW(M) MAIN 544 3120 CONTINUE MAIN 545 DBAR=DBAR/SUMW MAIN 546 C DBAS IS EITHER 0 OR DBAR ACCORDING TO WHETHER MAIN 547 C STRESS FORMULA 1 OR 2 IS BEING USED. MAIN 548 DBAS = 0.0 MAIN 549 IF(LSCH.EQ.2) DBAS = DBAR MAIN 550 IF(IABS(SDSWIT).GE.10) MAIN 551 1 CALL FITM(MX,MY,LFITSW) MAIN 552 IF(IABS(SDSWIT).LT.10) MAIN 553 1 CALL FITP(MX,MY,SDSWIT,LFITRM) MAIN 554 C MAIN 555 C CALCULATE U, T, AND STRESS MAIN 556 C MAIN 557 3200 U=0.0 MAIN 558 T=0.0 MAIN 559 DO 3210 M=MX,MY MAIN 560 U=U+(DIST(M)-DHAT(M))**2*WW(M) MAIN 561 3210 T=T+(DIST(M)-DBAS)**2*WW(M) MAIN 562 3215 U=SQRT (U) MAIN 563 TEMP1=T MAIN 564 T=SQRT (T) MAIN 565 GRSTRS(NG) = U/T MAIN 566 STRESS = STRESS + GRSTRS(NG)**2 MAIN 567 IF(U.EQ.0.0) GO TO 3340 MAIN 568 3220 RUT=1.0/(U*T) MAIN 569 UOT3=U/(TEMP1*T) MAIN 570 C MAIN 571 C CALCULATE THE (NEGATIVE) GRADIENT MAIN 572 C MAIN 573 3300 DO 3330 M = MX,MY MAIN 574 IF(DIST(M).EQ.0.0) GO TO 3330 MAIN 575 LTEMP1=MOD(IJ(M),KPACK) MAIN 576 I=LTEMP1/IJPACK+1 MAIN 577 J=MOD(LTEMP1,IJPACK)+1 MAIN 578 FACTOR=UOT3*(DIST(M)-DBAS)-RUT*(DIST(M)-DHAT(M)) MAIN 579 FACTOR = (FACTOR/RM1POW(DIST(M)) ) * WW(M) MAIN 580 FACTOR = FACTOR * GRSTRS(NG) MAIN 581 DO 3320 L=1,LDIM MAIN 582 TEMP1 = FACTOR * RM1POW(X(I,L)-X(J,L)) MAIN 583 GR(I,L)=GR(I,L)+TEMP1 MAIN 584 3320 GR(J,L)=GR(J,L)-TEMP1 MAIN 585 3330 CONTINUE MAIN 586 IF(SDSWIT.GE.10) GO TO 3340 MAIN 587 C SAVE POLYNOMIAL COEFFICIENTS MAIN 588 M=MX MAIN 589 DO 3335 M1=1,SDSWIT MAIN 590 DHAT(M)=PCOEFF(M1) MAIN 591 3335 M=M+1 MAIN 592 IF(LFITRM.EQ.2) DHAT(MX)=0.0 MAIN 593 3340 CONTINUE MAIN 594 IF(STRESS .EQ. 0.0 ) GO TO 3700 MAIN 595 STRESS = SQRT( STRESS / FNGRPS ) MAIN 596 C MAIN 597 C AVOID MOVING FIXED POINTS MAIN 598 C MAIN 599 IF( LNFIXZ .EQ. 0) GO TO 3400 MAIN 600 DO 3360 I=1,LNFIXZ MAIN 601 DO 3360 L=1,LDIM MAIN 602 3360 GR(I,L) = 0.0 MAIN 603 C MAIN 604 C FIND SCALE FACTOR OF GRADIENT, ANGLE-COSINE BETWEEN GRADIENT MAIN 605 C AND PREVIOUS GRADIENT. MAIN 606 C MAIN 607 3400 SFGR=0.0 MAIN 608 CAGRGL=0.0 MAIN 609 DO 3410 I=1,N MAIN 610 DO 3410 L=1,LDIM MAIN 611 SFGR=SFGR+GR(I,L)**2 MAIN 612 3410 CAGRGL=CAGRGL+GR(I,L)*GL(I,L) MAIN 613 SFGR=SQRT (SFGR/FN) MAIN 614 C IF GRADIENT = 0.0, SKIP AHEAD. MAIN 615 IF(SFGR) 3420,3700,3420 MAIN 616 3420 TEMP1=SFGR*SFGL*FN MAIN 617 CAGRGL=CAGRGL/TEMP1 MAIN 618 C MAIN 619 IF(ITNO.EQ.0 .OR. NBAKUP.GE.4) GO TO 3500 MAIN 620 IF(CAGRGL.GT.(-0.95) .AND. STRESS/STRLST.LT. 1.2001 ) GOTO 3500MAIN 621 C MAIN 622 C BACK UP ALONG LAST GRADIENT MAIN 623 C MAIN 624 NBAKUP = NBAKUP + 1 MAIN 625 IF(NBAKUP.EQ.1) STBASC=1.0 MAIN 626 STBASC=STBASC*STBAMU MAIN 627 TEMP1 = STEP MAIN 628 STEP = STEP / 10.0 MAIN 629 WRITE (LPRINT,38) STRESS, CAGRGL, STEP MAIN 630 38 FORMAT(10X, F7.3, 14X, F7.3, 22X, F8.4 ) MAIN 631 TEMP1 = (TEMP1 - STEP) / SFGL MAIN 632 TEMP1=STBASC*TEMP1 MAIN 633 DO 3430 I = 1,N MAIN 634 DO 3430 L = 1,LDIM MAIN 635 X(I,L) = X(I,L) - TEMP1*GL(I,L) MAIN 636 3430 GR(I,L) = GL(I,L) MAIN 637 SFGR = SFGL MAIN 638 STRESS = STRLST MAIN 639 GO TO 3000 MAIN 640 C MAIN 641 C STEP SIZE CALCULATIONS MAIN 642 C MAIN 643 3500 IF(ITNO) 9999, 3510, 3520 MAIN 644 3510 SRAT=0.8 MAIN 645 GO TO 3530 MAIN 646 3520 SRAT=STRESS/STRLST MAIN 647 3530 CALL NEWSTP( STEP, ITNO, SFGR, STRESS, MAIN 648 1 CAGRGL, COSAV, ACSAV, COSAVW, ACSAVW, SRAT, SRATAV ) MAIN 649 NBAKUP = 0 MAIN 650 C MAIN 651 C PRINT CURRENT STATUS OF COMPUTATION MAIN 652 C MAIN 653 3700 IF(LHIPRT.EQ.2) WRITE (LPRINT,30) ITNO,STRESS,SRAT,SRATAV, MAIN 654 1 CAGRGL,COSAV,ACSAV,SFGR,STEP MAIN 655 30 FORMAT(I10,F7.3,F7.3,F7.3,F7.3,F7.3,F7.3,F8.4, F8.4) MAIN 656 C MAIN 657 C DECIDE WHETHER TO CONTINUE ITERATING MAIN 658 C MAIN 659 3800 IF(STRESS) 9999, 3840, 3810 MAIN 660 3810 IF(SFGR-SFGRMN) 3850, 3850, 3815 MAIN 661 3815 TEMP1 = 0.5 * (1.0+SRATST) MAIN 662 TEMP2 = 1.0 - TEMP1 MAIN 663 IF( ABS (SRAT-TEMP1) - TEMP2 ) 3816, 3816, 3820 MAIN 664 3816 IF( ABS (SRATAV-TEMP1) - TEMP2 ) 3850, 3850, 3820 MAIN 665 3820 IF(STRESS-STRMIN) 3860, 3860, 3830 MAIN 666 3830 IF(ITNO-NOIT) 3900, 3870, 9999 MAIN 667 3840 CONTINUE MAIN 668 WRITE (LPRINT, 31) MAIN 669 31 FORMAT(24H0ZERO STRESS WAS REACHED ) MAIN 670 GO TO 4000 MAIN 671 3850 CONTINUE MAIN 672 WRITE (LPRINT, 32) MAIN 673 32 FORMAT(21H0MINIMUM WAS ACHIEVED ) MAIN 674 GO TO 4000 MAIN 675 3860 CONTINUE MAIN 676 WRITE (LPRINT, 33) MAIN 677 33 FORMAT(32H0SATISFACTORY STRESS WAS REACHED ) MAIN 678 GO TO 4000 MAIN 679 3870 CONTINUE MAIN 680 WRITE (LPRINT, 34) MAIN 681 34 FORMAT(39H0MAXIMUM NUMBER OF ITERATIONS WERE USED ) MAIN 682 GO TO 4000 MAIN 683 C MAIN 684 C CONTINUE ITERATING MAIN 685 C MAIN 686 3900 ITNO=ITNO+1 MAIN 687 TEMP1=STEP/SFGR MAIN 688 DO 3910 I=1,N MAIN 689 DO 3910 L=1,LDIM MAIN 690 3910 X(I,L)=X(I,L)+TEMP1*GR(I,L) MAIN 691 GO TO 3000 MAIN 692 C MAIN 693 C STOP ITERATING *************************************MAIN 694 C MAIN 695 C ROTATE FINAL CONFIGURATION TO PRINCIPAL COMPONENTS MAIN 696 C MAIN 697 4000 IF(LCOORS.LT.2) GO TO 4002 MAIN 698 C COMPUTE X TRANSPOSE TIMES X AND STORE UPPER HALF IN RVEC MAIN 699 KK=0 MAIN 700 DO 4410 K=1,LDIM MAIN 701 DO 4410 J=1,K MAIN 702 SUM=0.0 MAIN 703 DO 4405 I=1,N MAIN 704 4405 SUM=SUM+X(I,J)*X(I,K) MAIN 705 KK=KK+1 MAIN 706 4410 RVEC(KK)=SUM MAIN 707 C COMPUTE MATRIX OF EIGENVECTORS, GL MAIN 708 CALL SGEV(LDIM,LDIM,GL) MAIN 709 C COMPUTE X TIMES GL MAIN 710 DO 4420 K=1,LDIM MAIN 711 DO 4420 J=1,N MAIN 712 SUM=0.0 MAIN 713 DO 4415 I=1,LDIM MAIN 714 4415 SUM=SUM+X(J,I)*GL(I,K) MAIN 715 4420 GR(J,K)=SUM MAIN 716 DO 4425 I=1,N MAIN 717 DO 4425 J=1,LDIM MAIN 718 4425 X(I,J)=GR(I,J) MAIN 719 WRITE(LPRINT,152) MAIN 720 152 FORMAT(66H0THE FINAL CONFIGURATION HAS BEEN ROTATED TO PRINCIPAL CMAIN 721 1OMPONENTS.) MAIN 722 GO TO 4004 MAIN 723 C MAIN 724 4002 IF(LCOORS.EQ.1) WRITE(LPRINT,153) MAIN 725 153 FORMAT(110H0THE CONFIGURATION HAS BEEN NORMALIZED DURING THE ITERAMAIN 726 1TIONS BUT HAS NOT BEEN ROTATED TO PRINCIPAL COMPONENTS.) MAIN 727 IF(LCOORS.EQ.0) WRITE(LPRINT,154) MAIN 728 154 FORMAT(70H0NO NORMALIZATION WAS DONE TO THE CONFIGURATION DURING TMAIN 729 1HE ITERATIONS.) MAIN 730 4004 WRITE (LPRINT, 40)N,LDIM,STRESS,LSCH,(L,L=1,LDIM) MAIN 731 40 FORMAT(27H0THE FINAL CONFIGURATION OF,I4, MAIN 732 1 10H POINTS IN,I3, 22H DIMENSIONS HAS STRESS,F7.3,8H FORMULA ,I2 MAIN 733 2 /1X/30H LABEL FOR CONFIGURATION PLOTS,10X,19HFINAL CONFIGURATION MAIN 734 3 / 39X,10I7) MAIN 735 IF(LPUNSW.EQ.2) GO TO 4005 MAIN 736 WRITE (LPUNCH,41) (TITLE(I),I=1,80),N,LDIM MAIN 737 41 FORMAT (14H CONFIGURATION/80A1/2I3) MAIN 738 WRITE (LPUNCH,52) MAIN 739 52 FORMAT (12H (2X,10F7.3)) MAIN 740 4005 DO 4010 I=1,N MAIN 741 WRITE(LPRINT,42) PTID(I),I,(X(I,L),L=1,LDIM) MAIN 742 IF(LPUNSW.EQ.2) GO TO 4010 MAIN 743 WRITE (LPUNCH, 43)I,(X(I,L),L=1,LDIM) MAIN 744 4010 CONTINUE MAIN 745 42 FORMAT(20X,A1,18X,I2,10F7.3) MAIN 746 43 FORMAT(I2,10F7.3) MAIN 747 WRITE (LPRINT,46) MAIN 748 DO 4020 NG=1,NOGRPS MAIN 749 MX=GRNO(NG) MAIN 750 MZ=GRNO(NG+1)-MX MAIN 751 SDSWIT=GRSDIS(NG) MAIN 752 LFITRM=SDSWIT/100 MAIN 753 SDSWIT=MOD(SDSWIT,100) MAIN 754 IF(SDSWIT-10) 4016,4013,4014 MAIN 755 4013 WRITE(LPRINT,60) NG,MZ,GRSTRS(NG) MAIN 756 60 FORMAT(1X,I5,2X,I5,F7.3,11H ASCENDING) MAIN 757 GO TO 4020 MAIN 758 4014 WRITE(LPRINT,62) NG,MZ,GRSTRS(NG) MAIN 759 62 FORMAT(1X,I5,2X,I5,F7.3,11H DESCENDING) MAIN 760 GO TO 4020 MAIN 761 4016 M1=SDSWIT+MX-1 MAIN 762 M=MX MAIN 763 DO 4017 I=1,SDSWIT MAIN 764 PCOEFF(I)=DHAT(M) MAIN 765 4017 M=M+1 MAIN 766 WRITE(LPRINT,61) NG,MZ,GRSTRS(NG),(PCOEFF(I),I=1,SDSWIT) MAIN 767 61 FORMAT(1X,I5,2X,I5,F7.3,12H POLYNOMIAL,5F15.7) MAIN 768 C RESTORE DHAT VALUES MAIN 769 DO 4019 M=MX,M1 MAIN 770 DA=DATA(M) MAIN 771 TEMP=0.0 MAIN 772 DO 4018 I=LFITRM,SDSWIT MAIN 773 4018 TEMP=TEMP+PCOEFF(I)*REGR(DA,I) MAIN 774 4019 DHAT(M)=TEMP MAIN 775 4020 CONTINUE MAIN 776 46 FORMAT(14H0DATA GROUP(S) /73H0SERIAL COUNT STRESS REGRESSION COEMAIN 777 .FFICIENTS (FROM DEGREE 0 TO MAX OF 4)) MAIN 778 4030 CONTINUE MAIN 779 IF(LDIPRT.EQ.2) MAIN 780 1 CALL DATAPR(GRNO,MM,NOGRPS,2) MAIN 781 C MAIN 782 C PLOTTING SECTION MAIN 783 C MAIN 784 IF(LPLSCT.EQ.0) GO TO 5100 MAIN 785 C MAIN 786 C PLOT SCATTER DIAGRAMS OF DIST AND DHAT VS. DATA MAIN 787 PTID(1)=DEE MAIN 788 PTID(2)=AIE MAIN 789 C MAIN 790 C PLOT UNSPLIT SCATTER DIAGRAMS MAIN 791 C IF NOGRPS=1 THERE IS ONLY ONE SPLIT MAIN 792 C IF LPLSCT=3 PRODUCE A JOINT PLOT OF ALL SPLITS MAIN 793 C MAIN 794 IF(NOGRPS.GT.1) GO TO 5091 MAIN 795 WRITE(LPRINT,5093) LDIM,LSCH,STRESS,TITLE MAIN 796 5093 FORMAT(52H1DIST(D) AND DHAT(-) (Y-AXIS) VS. DATA (X-AXIS), FOR,I3,MAIN 797 . 27H DIMENSIONS. STRESS,FORMULA,I2,2H,=,F8.4/ 30X,80A1) MAIN 798 LNYSCT=2 MAIN 799 CALL PLOT(DATA,PMAT,0.,0.,0.,0.,MM,LNYSCT,-1,1800,1) MAIN 800 GO TO 5094 MAIN 801 C MAIN 802 5091 IF(LPLSCT.LT.3) GO TO 5092 MAIN 803 WRITE(LPRINT,44) LDIM,LSCH,STRESS,TITLE MAIN 804 44 FORMAT(45H1PLOT OF DIST (Y-AXIS) VS. DATA (X-AXIS), FOR,I3, MAIN 805 . 27H DIMENSIONS. STRESS,FORMULA,I2,2H,=,F8.4/ 30X,80A1) MAIN 806 LNYSCT=-2 MAIN 807 CALL PLOT(DATA,PMAT,0.,0.,0.,0.,MM,LNYSCT,-1,1800,1) MAIN 808 GO TO 5094 MAIN 809 C MAIN 810 C PLOT SPLIT SCATTER DIAGRAMS MAIN 811 C IF LPLSCT=1 PLOT NOT MORE THAN THE FIRST FIVE SPLITS MAIN 812 C IF LPLSCT=2 PLOT ALL SPLITS MAIN 813 C MAIN 814 5092 KPLT=NOGRPS MAIN 815 IF(LPLSCT.EQ.1) KPLT=MIN0(KPLT,5) MAIN 816 DO 5060 NG=1,KPLT MAIN 817 MX=GRNO(NG) MAIN 818 MY=GRNO(NG+1)-1 MAIN 819 MZ=MY-MX+1 MAIN 820 WRITE(LPRINT,4062) NG,LSCH,GRSTRS(NG),LDIM,TITLE MAIN 821 4062 FORMAT(43H1DIST(D) AND DHAT(-) VS. DATA FOR GROUP NO.,I3,17H. STRMAIN 822 .ESS,FORMULA,I2,2H,=,F8.4,12H DIMENSION=,I3 / 30X,80A1) MAIN 823 LNYSCT=2 MAIN 824 5060 CALL PLOT(DATA,PMAT,0.,0.,0.,0.,MZ,LNYSCT,-1,1800,MX) MAIN 825 5094 PTID(1)=AAA MAIN 826 PTID(2)=BEE MAIN 827 C MAIN 828 C PLOT CONFIGURATION MAIN 829 C MAIN 830 5100 IF(LPLCON.EQ.0)GO TO 5220 MAIN 831 IF(LDIM.EQ.1)GO TO 5190 MAIN 832 DO 5138 J=1,N MAIN 833 5138 PTID2(J)=PTID(J) MAIN 834 5125 IF(LPLCON.EQ.1) GO TO 5160 MAIN 835 C MAIN 836 C PLOT ALL PAIRS OF DIMENSIONS MAIN 837 C MAIN 838 DO 5145 J=2,LDIM MAIN 839 JMIN1=J-1 MAIN 840 DO 5145 I=1,JMIN1 MAIN 841 DO 5150 K=1,N MAIN 842 X2(K)=X(K,I) MAIN 843 Y2(K)=X(K,J) MAIN 844 5150 PTID(K)=PTID2(K) MAIN 845 WRITE(LPRINT,5015) J,I,TITLE MAIN 846 5015 FORMAT(31H1CONFIGURATION PLOT. DIMENSION,I2,23H (Y-AXIS) VS. DIMEMAIN 847 .NSION,I3,9H (X-AXIS) / 30X,80A1) MAIN 848 CALL SORT(Y2,N,PTID,X2,DUMEST,2,1) MAIN 849 LNYSCT=1 MAIN 850 5145 CALL PLOT(X2,Y2,0.,0.,0.,0.,N,LNYSCT,2,100,1) MAIN 851 GO TO 5200 MAIN 852 C MAIN 853 C PLOT PAIRS OF DIMENSIONS SO THAT EACH DIMEN PLOTTED ONLY ONCE MAIN 854 C MAIN 855 5160 LNYSCT=1 MAIN 856 DO 5185 J=1,LDIM,2 MAIN 857 JJ1=J+1 MAIN 858 JJ2=J MAIN 859 IF((MOD(LDIM,2).EQ.0).OR.(J.EQ.1)) GO TO 5175 MAIN 860 JJ1=J MAIN 861 JJ2=J-1 MAIN 862 IF(J.EQ.3) JJ2=1 MAIN 863 5175 DO 5180 K=1,N MAIN 864 X2(K)=X(K,JJ2) MAIN 865 Y2(K)=X(K,JJ1) MAIN 866 5180 PTID(K)=PTID2(K) MAIN 867 WRITE(LPRINT,5015) JJ1,JJ2,TITLE MAIN 868 CALL SORT(Y2,N,PTID,X2,DUMEST,2,1) MAIN 869 5185 CALL PLOT(X2,Y2,0.,0.,0.,0.,N,LNYSCT,2,100,1) MAIN 870 GO TO 5200 MAIN 871 C MAIN 872 C ONE DIMENSIONAL CONFIGURATION PLOT MAIN 873 C MAIN 874 5190 WRITE(LPRINT,5195) TITLE MAIN 875 5195 FORMAT(1H1,20X,34HPLOT OF CONFIGURATION, DIMENSION 1/30X,80A1) MAIN 876 LNYSCT=1 MAIN 877 CALL PLOT(X(1,1),X(1,1),0.,0.,0.,0.,N,LNYSCT,3,100,1) MAIN 878 GO TO 5220 MAIN 879 C MAIN 880 C CHANGE DIMENSION MAIN 881 C MAIN 882 5200 DO 5210 I=1,N MAIN 883 5210 PTID(I)=PTID2(I) MAIN 884 5220 ISTC=ISTC+1 MAIN 885 STPL(ISTC)= STRESS MAIN 886 DIMSV(ISTC)=LDIM MAIN 887 4100 LDIM=LDIM-LDIMD MAIN 888 IF(LDIM-LDIMN)4110,4101,4101 MAIN 889 4101 WRITE(LPRINT,4102) MAIN 890 4102 FORMAT(1H1) MAIN 891 GO TO 2300 MAIN 892 C MAIN 893 C PLOT STRESS VERSUS DIMENSION IF MORE THAN TWO POINTS MAIN 894 C MAIN 895 4110 CONTINUE MAIN 896 IF(ISTC.LT.3) GO TO 100 MAIN 897 WRITE(LPRINT, 50) TITLE MAIN 898 50 FORMAT (1H1,10X,31HPLOT OF STRESS VERSUS DIMENSION / 30X,80A1) MAIN 899 C INVERT TO IMPROVE READING OF PLOT MAIN 900 NVT=ISTC/2 MAIN 901 DO 8050 I=1,NVT MAIN 902 IC=ISTC+1-I MAIN 903 TEMP = DIMSV(I) MAIN 904 SEMP = STPL(I) MAIN 905 DIMSV(I) = DIMSV(IC) MAIN 906 STPL(I) =STPL(IC) MAIN 907 DIMSV(IC) = TEMP MAIN 908 STPL(IC) = SEMP MAIN 909 8050 CONTINUE MAIN 910 C FIND MAX STRESS FOR Y-AXIS , AND FILL IN PLOTTING CHARACTERS MAIN 911 YMA=0.0 MAIN 912 DO 8231 I=1,ISTC MAIN 913 IDIM=DIMSV(I) MAIN 914 PTID(I)=XNUM(IDIM) MAIN 915 8231 IF(YMA .LE. STPL(I)) YMA=STPL(I) MAIN 916 CALL SORT(STPL,ISTC,PTID,DIMSV,DUMEST,2,1) MAIN 917 LNYSCT=1 MAIN 918 CALL PLOT(DIMSV,STPL,10.,YMA,0.,0.,ISTC,LNYSCT,-2,10,1) MAIN 919 PTID(1)=AAA MAIN 920 PTID(2)=BEE MAIN 921 PTID(3)=CEE MAIN 922 PTID(4)=DEE MAIN 923 PTID(5)=EEE MAIN 924 PTID(6)=FFF MAIN 925 C MAIN 926 C REINITIALIZE, AND RETURN FOR MORE INPUT. MAIN 927 GO TO 100 MAIN 928 C MAIN 929 C NORMAL TERMINATION, AFTER READING 'STOP' ONCONTROL CARD MAIN 930 9000 STOP MAIN 931 C MAIN 932 C TROUBLE EXIT MAIN 933 C MAIN 934 9999 WRITE (LPRINT, 99) MAIN 935 99 FORMAT(52H0KRUSKAL. IMPOSSIBLE BRANCH TAKEN FROM IF STATEMENT. ) MAIN 936 STOP MAIN 937 C MAIN 938 END MAIN 939 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 KYST VERSION 2 OCTOBER, 1976 BLOCK 6 C BLOCK 7 C THIS PROGRAM HOLDS THE TABLE OF WORDS WHICH CCACT CONSULTS BLOCK 8 C LOGICALLY IT CONSISTS OF SEVERAL TABLES BLOCK 9 C THE VALUES IN MTAB INDICATE WHICH ROWS START NEW TABLES BLOCK 10 C EACH ENTRY IN THE TABLE CONTAINS 4 ITEMS BLOCK 11 C THE FIRST IS THE CODED ALPHABETIC WORD BLOCK 12 C (ONE WORD MAY BE ENTERED SEVERAL TIMES ) BLOCK 13 C THE SECOND INDICATES THE NATURE OF THIS ENTRY BLOCK 14 C 1 INDICATES INTEGER PARAMETERS WHICH MUST BE EXPLIC9TLY BLOCK 15 C READ IN BLOCK 16 C 2 INDICATES REAL PARAMETERS WHICH MUST BE EXPLICITLY BLOCK 17 C READ IN BLOCK 18 C 3 INDICATES AN INTEGER PARAMETER WHICH IS IMPLICITLY BLOCK 19 C DEFINED BLOCK 20 C 4 INDICATES A REAL PARAMETER WHICH IS IMPLICITLY BLOCK 21 C DEFINED BLOCK 22 C 5 INDICATES A PARAMETER WHICH BELONGS TO CCACT ONLY BLOCK 23 C 6 INDICATES AN IMPLICITLY DEFINED INTEGER PARAMETER BLOCK 24 C EFFECTING PROGRAM EXECUTION UPON RETURN FROM CCACT BLOCK 25 C THE THIRD INDICATES WHICH PARAMTER IS INVOLVED BLOCK 26 C SEPARATE NUMBERING FOR MAIN PROGRAM PARAMETERS AND CCACT BLOCK 27 C THE FOURTH GIVES THE VALUE FOR ANY INTERNALLY DEFINED PARAMETERS BLOCK 28 C IF THE PARAMETER IS OF TYPE 3,4 OR 6.IF OF TYPE 1 OR 2, THIS BLOCK 29 C ENTRY GIVES THE NUMBER OF ITEMS TO BE READ IN. BLOCK 30 C BLOCK 31 INTEGER LPAR( 44 ) BLOCK 32 REAL PAR( 44 ) BLOCK 33 C BLOCK 34 INTEGER MTAB(13) BLOCK 35 INTEGER TAB1(4,18),TAB2(4,17),TAB3(4,16),TAB4(4,19),TAB5(4,19) BLOCK 36 INTEGER PTID(100),ITEM(101) BLOCK 37 C BLOCK 38 COMMON /ACCUR/ PRECSN,XMAG,XMAXN,FLPRIN,CHPRIN 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 EQUIVALENCE (LPAR(1),PAR(1)) BLOCK 44 C BLOCK 45 DATA PRECSN, XMAG, XMAXN, FLPRIN, CHPRIN BLOCK 46 . /1.5E-8, 1.0E-38, 1.0E38, 6.0, 10.0/, BLOCK 47 . LREAD,LPRINT,LPUNCH,LSCRAT /5,6,7,8/ BLOCK 48 DATA MTAB(1),MTAB(2),MTAB(3),MTAB(4),MTAB(5),MTAB(6),MTAB(7), BLOCK 49 . MTAB(8),MTAB(9),MTAB(10),MTAB(11),MTAB(12),MTAB(13) BLOCK 50 . /1,52,54,59,65,67,76,78,80,83,87,90,90/ BLOCK 51 C BLOCK 52 C TAB1 ENTRIES PERTAIN SUCCESSIVELY TO THESE CONTROL WORDS: BLOCK 53 C DIMMAX, DIMMIN, DIMDIF, CUTOFF, STRMIN, BLOCK 54 C SFGRMN, COSAVW, ACSAVW, DIAGON, MATRIX, BLOCK 55 C HALFMA, LOWERH, UPPERH, LOWERC, ARBITR, BLOCK 56 C TORSCA, DATAFI, SPREAD BLOCK 57 DATA TAB1(1,1),TAB1(2,1),TAB1(3,1),TAB1(4,1),TAB1(1,2),TAB1(2,2) BLOCK 58 .,TAB1(3,2),TAB1(4,2),TAB1(1,3),TAB1(2,3),TAB1(3,3),TAB1(4,3) BLOCK 59 .,TAB1(1,4),TAB1(2,4),TAB1(3,4),TAB1(4,4),TAB1(1,5),TAB1(2,5) BLOCK 60 .,TAB1(3,5),TAB1(4,5),TAB1(1,6),TAB1(2,6),TAB1(3,6),TAB1(4,6) BLOCK 61 .,TAB1(1,7),TAB1(2,7),TAB1(3,7),TAB1(4,7),TAB1(1,8),TAB1(2,8) BLOCK 62 .,TAB1(3,8),TAB1(4,8),TAB1(1,9),TAB1(2,9),TAB1(3,9),TAB1(4,9) BLOCK 63 .,TAB1(1,10),TAB1(2,10),TAB1(3,10),TAB1(4,10),TAB1(1,11) BLOCK 64 .,TAB1(2,11),TAB1(3,11),TAB1(4,11)/ BLOCK 65 1 6989, 1, 1, 1, BLOCK 66 2 5375, 1, 2, 1, BLOCK 67 3 6923, 1, 3, 1, BLOCK 68 4 13520, 2, 4, 1, BLOCK 69 5 390, 2, 5, 1, BLOCK 70 6 2553, 2, 6, 1, BLOCK 71 7 7303, 2, 7, 1, BLOCK 72 8 11226, 2, 8, 1, BLOCK 73 9 11136, 5, 1, 2, BLOCK 74 A 9277, 3, 10, 1, BLOCK 75 B 3527, 3, 10, 2/ BLOCK 76 DATA TAB1(1,12),TAB1(2,12),TAB1(3,12),TAB1(4,12),TAB1(1,13), BLOCK 77 . TAB1(2,13),TAB1(3,13),TAB1(4,13),TAB1(1,14),TAB1(2,14), BLOCK 78 . TAB1(3,14),TAB1(4,14),TAB1(1,15),TAB1(2,15),TAB1(3,15), BLOCK 79 . TAB1(4,15),TAB1(1,16),TAB1(2,16),TAB1(3,16),TAB1(4,16), BLOCK 80 . TAB1(1,17),TAB1(2,17),TAB1(3,17),TAB1(4,17),TAB1(1,18), BLOCK 81 . TAB1(2,18),TAB1(3,18),TAB1(4,18)/ BLOCK 82 C 6069, 3, 10, 2, BLOCK 83 D 8938, 3, 10, 3, BLOCK 84 E 1754, 3, 10, 4, BLOCK 85 F 10499, 3, 20, -99, BLOCK 86 G 2048, 3, 27, 4, BLOCK 87 H 3323, 1, 43, 1, BLOCK 88 I 6910, 2, 44, 1/ BLOCK 89 C BLOCK 90 C TAB2 ENTRIES PERTAIN SUCCESSIVELY TO THESE CONTROL WORDS: BLOCK 91 C UPPERC, BLOCKD, DISSIM, DISSIM, SIMILA, BLOCK 92 C SIMILA, CONFIG, COMPUT, PRIMAR, SECOND, BLOCK 93 C R , ITERAT, FIX , SRATST, STOP , BLOCK 94 C PRE-IT, COORDI BLOCK 95 DATA TAB2(1,1),TAB2(2,1),TAB2(3,1),TAB2(4,1),TAB2(1,2),TAB2(2,2) BLOCK 96 .,TAB2(3,2),TAB2(4,2),TAB2(1,3),TAB2(2,3),TAB2(3,3),TAB2(4,3) BLOCK 97 .,TAB2(1,4),TAB2(2,4),TAB2(3,4),TAB2(4,4),TAB2(1,5),TAB2(2,5) BLOCK 98 .,TAB2(3,5),TAB2(4,5),TAB2(1,6),TAB2(2,6),TAB2(3,6),TAB2(4,6) BLOCK 99 .,TAB2(1,7),TAB2(2,7),TAB2(3,7),TAB2(4,7),TAB2(1,8),TAB2(2,8) BLOCK100 .,TAB2(3,8),TAB2(4,8),TAB2(1,9),TAB2(2,9),TAB2(3,9),TAB2(4,9) BLOCK101 .,TAB2(1,10),TAB2(2,10),TAB2(3,10),TAB2(4,10),TAB2(1,11) BLOCK102 .,TAB2(2,11),TAB2(3,11),TAB2(4,11)/ BLOCK103 1 4623, 3, 10, 5, BLOCK104 2 8683, 5, 1, 5, BLOCK105 3 15096, 6, 12, 2, BLOCK106 4 15096, 3, 11, 10, BLOCK107 5 6868, 6, 12, 2, BLOCK108 6 6868, 3, 11, 11, BLOCK109 7 14846, 6, 12, 5, BLOCK110 8 11754, 6, 12, 6, BLOCK111 9 765, 3, 13, 1, BLOCK112 A 4119, 3, 13, 2, BLOCK113 B 16326, 2, 14, 1/ BLOCK114 DATA TAB2(1,12),TAB2(2,12),TAB2(3,12),TAB2(4,12),TAB2(1,13), BLOCK115 . TAB2(2,13),TAB2(3,13),TAB2(4,13),TAB2(1,14),TAB2(2,14), BLOCK116 . TAB2(3,14),TAB2(4,14),TAB2(1,15),TAB2(2,15),TAB2(3,15), BLOCK117 . TAB2(4,15),TAB2(1,16),TAB2(2,16),TAB2(3,16),TAB2(4,16), BLOCK118 . TAB2(1,17),TAB2(2,17),TAB2(3,17),TAB2(4,17)/ BLOCK119 C 15170, 1, 15, 1, BLOCK120 D 1855, 1, 28, 1, BLOCK121 E 3720, 2, 16, 1, BLOCK122 F 13171, 6, 12, 7, BLOCK123 G 10738, 1, 39, 1, BLOCK124 H 7261, 5, 1, 11/ BLOCK125 C BLOCK126 C TAB3 ENTRIES PERTAIN SUCCESSIVELY TO THESE CONTROL WORDS: BLOCK127 C WEIGHT, WFUNCT, SFORM1, SFORM2, CARDS, BLOCK128 C NOCARD, SPLIT , DATA , SAVE , RANDOM, BLOCK129 C PRINT , DFUNCT, REGRES, DCONST, WCONST, BLOCK130 C PLOT BLOCK131 DATA TAB3(1,1),TAB3(2,1),TAB3(3,1),TAB3(4,1),TAB3(1,2),TAB3(2,2) BLOCK132 .,TAB3(3,2),TAB3(4,2),TAB3(1,3),TAB3(2,3),TAB3(3,3),TAB3(4,3) BLOCK133 .,TAB3(1,4),TAB3(2,4),TAB3(3,4),TAB3(4,4),TAB3(1,5),TAB3(2,5) BLOCK134 .,TAB3(3,5),TAB3(4,5),TAB3(1,6),TAB3(2,6),TAB3(3,6),TAB3(4,6) BLOCK135 .,TAB3(1,7),TAB3(2,7),TAB3(3,7),TAB3(4,7),TAB3(1,8),TAB3(2,8) BLOCK136 .,TAB3(3,8),TAB3(4,8),TAB3(1,9),TAB3(2,9),TAB3(3,9),TAB3(4,9) BLOCK137 .,TAB3(1,10),TAB3(2,10),TAB3(3,10),TAB3(4,10),TAB3(1,11) BLOCK138 .,TAB3(2,11),TAB3(3,11),TAB3(4,11)/ BLOCK139 1 14543, 6, 12, 3, BLOCK140 2 11427, 6, 12, 4, BLOCK141 3 5318, 3, 17, 1, BLOCK142 4 6181, 3, 17, 2, BLOCK143 5 6927, 3, 18, 1, BLOCK144 6 16009, 3, 18, 2, BLOCK145 7 1966, 5, 1, 3, BLOCK146 8 6675, 6, 12, 2, BLOCK147 9 9189, 5, 1, 7, BLOCK148 A 8330, 1, 20, 1, BLOCK149 B 2775, 5, 1, 4/ BLOCK150 DATA TAB3(1,12),TAB3(2,12),TAB3(3,12),TAB3(4,12),TAB3(1,13), BLOCK151 . TAB3(2,13),TAB3(3,13),TAB3(4,13),TAB3(1,14),TAB3(2,14), BLOCK152 . TAB3(3,14),TAB3(4,14),TAB3(1,15),TAB3(2,15),TAB3(3,15), BLOCK153 . TAB3(4,15),TAB3(1,16),TAB3(2,16),TAB3(3,16),TAB3(4,16)/ BLOCK154 C 10575, 6, 12, 8, BLOCK155 D 14439, 5, 1, 6, BLOCK156 E 267, 2, 29, 5, BLOCK157 F 1119, 2, 34, 5, BLOCK158 G 6878, 5, 1, 8/ BLOCK159 C BLOCK160 C TAB4 ENTRIES PERTAIN SUCCESSIVELY TO THESE CONTROL WORDS: BLOCK161 C ABSENT, PRESEN, BYROWS, BYGROU, BYDECK, BLOCK162 C NOMORE, NOMORE, DATA , DISTAN, HISTOR, BLOCK163 C NODATA, NODIST, NOHIST, NO , YES , BLOCK164 C MONOTO, ASCEND, DESCEN, POLYNO BLOCK165 DATA TAB4(1,1),TAB4(2,1),TAB4(3,1),TAB4(4,1),TAB4(1,2),TAB4(2,2) BLOCK166 .,TAB4(3,2),TAB4(4,2),TAB4(1,3),TAB4(2,3),TAB4(3,3),TAB4(4,3) BLOCK167 .,TAB4(1,4),TAB4(2,4),TAB4(3,4),TAB4(4,4),TAB4(1,5),TAB4(2,5) BLOCK168 .,TAB4(3,5),TAB4(4,5),TAB4(1,6),TAB4(2,6),TAB4(3,6),TAB4(4,6) BLOCK169 .,TAB4(1,7),TAB4(2,7),TAB4(3,7),TAB4(4,7),TAB4(1,8),TAB4(2,8) BLOCK170 .,TAB4(3,8),TAB4(4,8),TAB4(1,9),TAB4(2,9),TAB4(3,9),TAB4(4,9) BLOCK171 .,TAB4(1,10),TAB4(2,10),TAB4(3,10),TAB4(4,10),TAB4(1,11) BLOCK172 .,TAB4(2,11),TAB4(3,11),TAB4(4,11)/ BLOCK173 1 4258, 3, 9, 2, BLOCK174 2 2575, 3, 9, 1, BLOCK175 3 7761, 3, 19, 100, BLOCK176 4 11782, 3, 19, 200, BLOCK177 5 11288, 3, 19, 300, BLOCK178 6 5274, 3, 19, 400, BLOCK179 7 5274, 3, 19, 2, BLOCK180 8 6675, 3, 21, 2, BLOCK181 9 9824, 3, 22, 2, BLOCK182 A 12801, 3, 24, 2, BLOCK183 B 16057, 3, 21, 1/ BLOCK184 DATA TAB4(1,12),TAB4(2,12),TAB4(3,12),TAB4(4,12),TAB4(1,13), BLOCK185 . TAB4(2,13),TAB4(3,13),TAB4(4,13),TAB4(1,14),TAB4(2,14), BLOCK186 . TAB4(3,14),TAB4(4,14),TAB4(1,15),TAB4(2,15),TAB4(3,15), BLOCK187 . TAB4(4,15),TAB4(1,16),TAB4(2,16),TAB4(3,16),TAB4(4,16), BLOCK188 . TAB4(1,17),TAB4(2,17),TAB4(3,17),TAB4(4,17),TAB4(1,18), BLOCK189 . TAB4(2,18),TAB4(3,18),TAB4(4,18),TAB4(1,19),TAB4(2,19), BLOCK190 . TAB4(3,19),TAB4(4,19)/ BLOCK191 C 5863, 3, 22, 1, BLOCK192 D 9395, 3, 24, 1, BLOCK193 E 9622, 3, 10, 100, BLOCK194 F 11125, 3, 10, 200, BLOCK195 G 15634, 3, 23, 0, BLOCK196 H 7782, 3, 23, 10, BLOCK197 I 11188, 3, 23, 11, BLOCK198 J 3756, 3, 26, 1/ BLOCK199 C BLOCK200 C TAB5 ENTRIES PERTAIN SUCCESSIVELY TO THESE CONTROL WORDS: BLOCK201 C POLYNO, CONSTA, NOCONS, MULTIV, MULTIV, BLOCK202 C DATA , CONFIG, CONFIG, SCATTE, ALL , BLOCK203 C SOME , NONE , ALL , SOME , NONE , BLOCK204 C JOINT , ROTATE, STANDA, AS-IS BLOCK205 DATA TAB5(1,1),TAB5(2,1),TAB5(3,1),TAB5(4,1),TAB5(1,2),TAB5(2,2) BLOCK206 .,TAB5(3,2),TAB5(4,2),TAB5(1,3),TAB5(2,3),TAB5(3,3),TAB5(4,3) BLOCK207 .,TAB5(1,4),TAB5(2,4),TAB5(3,4),TAB5(4,4),TAB5(1,5),TAB5(2,5) BLOCK208 .,TAB5(3,5),TAB5(4,5),TAB5(1,6),TAB5(2,6),TAB5(3,6),TAB5(4,6) BLOCK209 .,TAB5(1,7),TAB5(2,7),TAB5(3,7),TAB5(4,7),TAB5(1,8),TAB5(2,8) BLOCK210 .,TAB5(3,8),TAB5(4,8),TAB5(1,9),TAB5(2,9),TAB5(3,9),TAB5(4,9) BLOCK211 .,TAB5(1,10),TAB5(2,10),TAB5(3,10),TAB5(4,10),TAB5(1,11) BLOCK212 .,TAB5(2,11),TAB5(3,11),TAB5(4,11)/ BLOCK213 1 3756, 1, 23, 1, BLOCK214 2 14387, 3, 26, 100, BLOCK215 3 5018, 3, 26, 200, BLOCK216 4 3608, 3, 26, 0, BLOCK217 5 3608, 1, 23, 1, BLOCK218 6 6675, 3, 25, 2, BLOCK219 7 14846, 3, 27, 3, BLOCK220 8 14846, 5, 1, 9, BLOCK221 9 11109, 5, 1, 10, BLOCK222 A 5766, 3, 40, 2, BLOCK223 B 13660, 3, 40, 1/ BLOCK224 DATA TAB5(1,12),TAB5(2,12),TAB5(3,12),TAB5(4,12),TAB5(1,13), BLOCK225 . TAB5(2,13),TAB5(3,13),TAB5(4,13),TAB5(1,14),TAB5(2,14), BLOCK226 . TAB5(3,14),TAB5(4,14),TAB5(1,15),TAB5(2,15),TAB5(3,15), BLOCK227 . TAB5(4,15),TAB5(1,16),TAB5(2,16),TAB5(3,16),TAB5(4,16), BLOCK228 . TAB5(1,17),TAB5(2,17),TAB5(3,17),TAB5(4,17),TAB5(1,18), BLOCK229 . TAB5(2,18),TAB5(3,18),TAB5(4,18),TAB5(1,19),TAB5(2,19), BLOCK230 . TAB5(3,19),TAB5(4,19)/ BLOCK231 C 10008, 3, 40, 0, BLOCK232 D 5766, 3, 41, 2, BLOCK233 E 13660, 3, 41, 1, BLOCK234 F 10008, 3, 41, 0, BLOCK235 G 11053, 3, 41, 3, BLOCK236 H 4503, 3, 42, 2, BLOCK237 I 3418, 3, 42, 1, BLOCK238 J 9499, 3, 42, 0/ BLOCK239 C BLOCK240 DATA PTID(1),PTID(2),PTID(3),PTID(4),PTID(5),PTID(6),PTID(7), BLOCK241 . PTID(8),PTID(9),PTID(10),PTID(11),PTID(12),PTID(13),PTID(14), BLOCK242 . PTID(15),PTID(16),PTID(17),PTID(18),PTID(19),PTID(20),PTID(21), BLOCK243 . PTID(22),PTID(23),PTID(24),PTID(25),PTID(26),PTID(27),PTID(28), BLOCK244 . PTID(29),PTID(30),PTID(31),PTID(32),PTID(33),PTID(34),PTID(35), BLOCK245 . PTID(36),PTID(37),PTID(38),PTID(39),PTID(40),PTID(41),PTID(42), BLOCK246 . PTID(43),PTID(44),PTID(45),PTID(46),PTID(47),PTID(48),PTID(49), BLOCK247 . PTID(50)/ BLOCK248 . 1HA,1HB,1HC,1HD,1HE,1HF,1HG,1HH,1HI,1HJ,1HK,1HL,1HM, BLOCK249 1 1HN,1HO,1HP,1HQ,1HR,1HS,1HT,1HU,1HV,1HW,1HX,1HY, BLOCK250 2 1HZ,1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9,1HA,1HB,1HC, BLOCK251 3 1HD,1HE,1HF,1HG,1HH,1HI,1HJ,1HK,1HL,1HM,1HN,1HO/ BLOCK252 DATA PTID(51),PTID(52),PTID(53),PTID(54),PTID(55),PTID(56), BLOCK253 . PTID(57),PTID(58),PTID(59),PTID(60),PTID(61),PTID(62),PTID(63), BLOCK254 . PTID(64),PTID(65),PTID(66),PTID(67),PTID(68),PTID(69),PTID(70), BLOCK255 . PTID(71),PTID(72),PTID(73),PTID(74),PTID(75),PTID(76),PTID(77), BLOCK256 . PTID(78),PTID(79),PTID(80),PTID(81),PTID(82),PTID(83),PTID(84), BLOCK257 . PTID(85),PTID(86),PTID(87),PTID(88),PTID(89),PTID(90),PTID(91), BLOCK258 . PTID(92),PTID(93),PTID(94),PTID(95),PTID(96),PTID(97), BLOCK259 . PTID(98),PTID(99),PTID(100)/ BLOCK260 4 1HP,1HQ,1HR,1HS,1HT,1HU,1HV,1HW,1HX,1HY,1HZ,1H1,1H2, BLOCK261 5 1H3,1H4,1H5,1H6,1H7,1H8,1H9,1HA,1HB,1HC,1HD,1HE, BLOCK262 6 1HF,1HG,1HH,1HI,1HJ,1HK,1HL,1HM,1HN,1HO,1HP,1HQ, BLOCK263 7 1HR,1HS,1HT,1HU,1HV,1HW,1HX,1HY,1HZ,1H1,1H2,1H3,1H4/ BLOCK264 DATA ITEM(1),ITEM(2),ITEM(3),ITEM(4),ITEM(5),ITEM(6),ITEM(7), BLOCK265 . ITEM(8),ITEM(9),ITEM(10),ITEM(11),ITEM(12),ITEM(13),ITEM(14), BLOCK266 . ITEM(15),ITEM(16),ITEM(17),ITEM(18),ITEM(19),ITEM(20),ITEM(21), BLOCK267 . ITEM(22),ITEM(23),ITEM(24),ITEM(25),ITEM(26),ITEM(27),ITEM(28), BLOCK268 . ITEM(29),ITEM(30),ITEM(31),ITEM(32),ITEM(33),ITEM(34),ITEM(35), BLOCK269 . ITEM(36),ITEM(37),ITEM(38),ITEM(39),ITEM(40),ITEM(41),ITEM(42), BLOCK270 . ITEM(43),ITEM(44),ITEM(45),ITEM(46),ITEM(47),ITEM(48),ITEM(49)/ BLOCK271 1 1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H , BLOCK272 2 1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H , BLOCK273 3 1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H , BLOCK274 4 1H ,1H ,1H ,1H / BLOCK275 DATA BLOCK276 . ITEM(50),ITEM(51),ITEM(52),ITEM(53),ITEM(54),ITEM(55),ITEM(56), BLOCK277 . ITEM(57),ITEM(58),ITEM(59),ITEM(60),ITEM(61),ITEM(62),ITEM(63), BLOCK278 . ITEM(64),ITEM(65),ITEM(66),ITEM(67),ITEM(68),ITEM(69),ITEM(70), BLOCK279 . ITEM(71),ITEM(72),ITEM(73),ITEM(74),ITEM(75),ITEM(76),ITEM(77), BLOCK280 . ITEM(78),ITEM(79),ITEM(80),ITEM(81),ITEM(82),ITEM(83),ITEM(84), BLOCK281 . ITEM(85),ITEM(86),ITEM(87),ITEM(88),ITEM(89),ITEM(90),ITEM(91), BLOCK282 . ITEM(92),ITEM(93),ITEM(94),ITEM(95),ITEM(96),ITEM(97),ITEM(98), BLOCK283 . ITEM(99),ITEM(100),ITEM(101)/ BLOCK284 5 1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H , BLOCK285 6 1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H , BLOCK286 7 1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H ,1H , BLOCK287 8 1H ,1H ,1H ,1H ,1H ,1H ,1H / BLOCK288 END BLOCK289 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 C KYST VERSION 2 OCTOBER, 1976 BSEC1 7 DIMENSION A(1),B(1),W(1) BSEC1 8 REAL ETA,E,S BSEC1 9 C BSEC1 10 C W.KAHAN, ACCURATE EIGENVALUES OF A SYM- BSEC1 11 C METRIC TRIDIAGONAL MATRIX. TECH.REPORT BSEC1 12 C NO.CS41, JULY 22,1966, COMP.SC.DEPT., BSEC1 13 C STANFORD UNIVERSITY. BSEC1 14 C BSEC1 15 C BSEC1 IS CALLED BY SGEV, AND BSEC1 16 C USES IB BISECTION STEPS TO FIND THE ALGEBRAICALLY BSEC1 17 C GREATEST EIGENVALUE E OF THE SYMMETRIC BSEC1 18 C TRIDIAGONAL MATRIX WHOSE DIAGONAL ELEMENTS ARE BSEC1 19 C A(1) THROUGH A(N) BSEC1 20 C AND WHOSE SUB-AND-SUPERDIAGONAL ELEMENTS ARE BSEC1 21 C B(1) THROUGH B(N-1) BSEC1 22 C FURTHERMORE, BSEC1 23 C W(1) THROUGH W(N) BSEC1 24 C ARE SET EQUAL TO THE DIAGONAL ELEMENTS OF THE BSEC1 25 C MATRIX U IN THE LU-DECOMPOSITION REQUIRED BSEC1 26 C FOR SUBSEQUENT INVERSE ITERATION. A AND B BSEC1 27 C ARE RESCALED, THE CHANGE IF SCALE BEING BSEC1 28 C REFLECTED IN S. ETA IS LEAST POSITIVE MA- BSEC1 29 C CHINE NUMBER. BSEC1 30 C BSEC1 31 REAL BB,D,T,X BSEC1 32 IF(N.GT.1)GOTO 5 BSEC1 33 E=A(1) BSEC1 34 W(1)=ETA BSEC1 35 GOTO 70 BSEC1 36 5 B(N)=0.E0 BSEC1 37 T=0.E0 BSEC1 38 DO 10 I=1,N BSEC1 39 T=AMAX1(T,ABS(A(I))) BSEC1 40 10 T=AMAX1(T,ABS(B(I))) BSEC1 41 IF(T.NE.0.E0)GOTO 20 BSEC1 42 E=0.E0 BSEC1 43 DO 19 I=1,N BSEC1 44 19 W(I)=ETA BSEC1 45 GOTO 70 BSEC1 46 20 S=S*T BSEC1 47 DO 30 I=1,N BSEC1 48 A(I)=A(I)/T BSEC1 49 30 B(I)=B(I)/T BSEC1 50 E=3.E0 BSEC1 51 D=-E BSEC1 52 DO 50 K=1,IB BSEC1 53 X=(D+E)/2.E0 BSEC1 54 U=1.E0 BSEC1 55 NU=0 BSEC1 56 BB=0.E0 BSEC1 57 DO 40 I=1,N BSEC1 58 U=(A(I)-BB/U)-X BSEC1 59 IF(U.GT.ETA)GOTO 35 BSEC1 60 U=AMIN1(U,-ETA) BSEC1 61 NU=NU+1 BSEC1 62 35 BB=B(I)**2 BSEC1 63 40 W(I)=U BSEC1 64 IF(NU.LT.N)D=X BSEC1 65 IF(NU.EQ.N)E=X BSEC1 66 50 CONTINUE BSEC1 67 IF(NU.EQ.N)GOTO 70 BSEC1 68 U=1.E0 BSEC1 69 BB=0.E0 BSEC1 70 DO 60 I=1,N BSEC1 71 U=(A(I)-BB/U)-E BSEC1 72 IF(U.GT.ETA)GOTO 55 BSEC1 73 U=AMIN1(U,-ETA) BSEC1 74 55 BB=B(I)**2 BSEC1 75 60 W(I)=U BSEC1 76 70 RETURN BSEC1 77 END BSEC1 78 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 KYST VERSION 2 OCTOBER, 1976 CACT 6 C CCACT--CONTROL CARD ACTIVITY. READS AND INTERPRETS CONTROL CARDS.CACT 7 C CACT 8 REAL DUMMY(27),ATAB(4,89) CACT 9 INTEGER TAB(4,89),IFAC(6) CACT 10 INTEGER MTAB(13) CACT 11 C CACT 12 INTEGER LPAR( 44 ) CACT 13 REAL PAR( 44 ) CACT 14 C CACT 15 INTEGER IPARAM(5) CACT 16 C CACT 17 INTEGER DFLTSW CACT 18 INTEGER C(81),BLANK,DOT,EQUALS,COMMA,DOLLAR,MAP(40),X CACT 19 INTEGER CHTAB(13),WORD(18) CACT 20 INTEGER BLSW,NUMSW,DECSW,TYPE,XTYPE,T,TA, PARNO, TABNO CACT 21 C CACT 22 COMMON /MDSCL1/ LREAD, LPRINT, LPUNCH, LSCRAT CACT 23 COMMON /MDSCL2/ LPAR CACT 24 COMMON /MDSCL3/ MTAB, TAB CACT 25 EQUIVALENCE (LPAR(1),PAR(1)) CACT 26 EQUIVALENCE (MAP(1),DUMMY(1)), (MAP(28),CHTAB(1)) CACT 27 EQUIVALENCE (ATAB(1),TAB(1)) CACT 28 C CACT 29 DATA BLANK, DOT, EQUALS, COMMA, C(81) /1H ,1H.,1H=,1H,,1H / CACT 30 DATA DOLLAR/1H$/ CACT 31 DATA MODP,IFAC(1),IFAC(2),IFAC(3),IFAC(4),IFAC(5),IFAC(6) CACT 32 1 / 16381, 907, 887, 883, 881, 877, 863/ CACT 33 C CACT 34 DATA MAP(1),MAP(2),MAP(3),MAP(4),MAP(5),MAP(6),MAP(7),MAP(8), CACT 35 . MAP(9),MAP(10),MAP(11),MAP(12),MAP(13),MAP(14),MAP(15),MAP(16), CACT 36 . MAP(17),MAP(18),MAP(19),MAP(20),MAP(21),MAP(22),MAP(23),MAP(24), CACT 37 . MAP(25),MAP(26),MAP(27),MAP(28),MAP(29),MAP(30),MAP(31),MAP(32), CACT 38 . MAP(33),MAP(34),MAP(35),MAP(36),MAP(37),MAP(38),MAP(39),MAP(40)/ CACT 39 . 1H ,1HA,1HB,1HC,1HD,1HE,1HF,1HG,1HH,1HI,1HJ,1HK,1HL, CACT 40 . 1HM,1HN,1HO,1HP,1HQ,1HR,1HS,1HT,1HU,1HV,1HW,1HX,1HY, CACT 41 . 1HZ,1H0,1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9,1H-,1H+,1H./CACT 42 C CACT 43 C CACT 44 C CACT 45 1 FORMAT(80A1) CACT 46 11 FORMAT(1H0,80A1) CACT 47 2 FORMAT(1X,I1) CACT 48 3 FORMAT(1X,80A1) CACT 49 4 FORMAT(1X,I18) CACT 50 5 FORMAT(1X,F18.3) CACT 51 9 FORMAT(98H CCACT DIAGNOSTIC. A CONTROL CARD WAS EXPECTED ABOVE, BCACT 52 .UT THIS CARD IS NOT A PROPER CONTROL CARD.) CACT 53 99 FORMAT( 12H ITEM NUMBER, I5, 20H HAS EXPECTED TYPE , I5, CACT 54 1 18H AND ACTUAL TYPE , I5, 16H . THIS ITEM IS , 18A1) CACT 55 98 FORMAT(12H ITEM NUMBER,I5,39H HAS ILLEGAL CHARACTER. THIS ITEM ISCACT 56 . ,18A1) CACT 57 C CACT 58 C READ AND PRINT CONTROL CARD CACT 59 C CACT 60 100 READ (LREAD,1) (C(I),I=1,80) CACT 61 WRITE(LPRINT,11) (C(I),I=1,80) CACT 62 C ANALYZE CONTROL CARD INTO TOKENS CACT 63 C CACT 64 C EACH 'TOKEN' IS DELIMITED BY BLANKS CACT 65 C THERE ARE FOUR TYPES OF TOKENS. CACT 66 C ALPHABETIC, INTEGER, DECIMAL, AND DEFAULT CACT 67 C ALPHABETIC UNLESS FIRST CHARACTER IS DIGIT OR DEC POINT CACT 68 C OR PLUS OR MINUS OR DOLLAR SIGN CACT 69 C DECIMALS DISTINGUISHED BY DECIMAL POINT CACT 70 C DEFAULT = $ CACT 71 C CACT 72 REWIND LSCRAT CACT 73 BLSW=1 CACT 74 NUMSW=1 CACT 75 DECSW=1 CACT 76 DFLTSW=1 CACT 77 K=0 CACT 78 C CACT 79 DO 300 I=1,81 CACT 80 C CACT 81 X=C(I) CACT 82 GO TO (110,120), BLSW CACT 83 C CACT 84 110 IF(X.EQ.BLANK .OR. X.EQ.EQUALS .OR. X.EQ. COMMA ) GO TO 300CACT 85 BLSW=2 CACT 86 JA=I CACT 87 DO 115 KX=1,13 CACT 88 115 IF(X.EQ.CHTAB(KX))NUMSW=2 CACT 89 IF(X.EQ.DOLLAR) DFLTSW=2 CACT 90 IF(X.EQ.DOT) DECSW=2 CACT 91 GO TO 300 CACT 92 C CACT 93 120 IF(X.EQ.BLANK .OR. X.EQ.EQUALS .OR. X.EQ. COMMA ) GO TO 130CACT 94 IF(X.EQ.DOT) DECSW=2 CACT 95 GO TO 300 CACT 96 C CACT 97 130 K=K+1 CACT 98 JB = MIN0 (I-1,JA+16) CACT 99 JC=18-(JB-JA+1) CACT 100 TYPE=1 CACT 101 IF(NUMSW.EQ.2) TYPE=NUMSW+DECSW-1 CACT 102 IF(DFLTSW.EQ.2) TYPE=4 CACT 103 C CACT 104 WRITE (LSCRAT,2) TYPE CACT 105 IF(TYPE.EQ.4) GO TO 160 CACT 106 GO TO (140,150),NUMSW CACT 107 140 WRITE (LSCRAT,3) (C(J),J=JA,JB), (BLANK,J=1,JC) CACT 108 GO TO 160 CACT 109 150 WRITE (LSCRAT,3) (BLANK,J=1,JC), (C(J),J=JA,JB) CACT 110 160 BLSW=1 CACT 111 NUMSW=1 CACT 112 DECSW=1 CACT 113 DFLTSW=1 CACT 114 GO TO 300 CACT 115 C CACT 116 300 CONTINUE CACT 117 C CACT 118 C ANALYZE TOKENS AND SET PARAMETER VALUES ACCORDINGLY CACT 119 C CACT 120 KB=K CACT 121 IF(KB.EQ.0) RETURN CACT 122 REWIND LSCRAT CACT 123 XTYPE=1 CACT 124 IPARAM(1) = 1 CACT 125 NT=0 CACT 126 NTOKNS=0 CACT 127 NOPCC=0 CACT 128 C CACT 129 DO 1000 K=1,KB CACT 130 C CACT 131 READ (LSCRAT,2) TYPE CACT 132 IF( (XTYPE.EQ.1) .AND. (TYPE.NE.1) ) GO TO 995 CACT 133 IF( (XTYPE.NE.TYPE). AND. (TYPE.NE.4) ) GO TO 995 CACT 134 350 GO TO (400,410,420,430), TYPE CACT 135 C CACT 136 400 READ (LSCRAT,3) (WORD(L),L=1,18) CACT 137 C CONVERT FIRST SIX LETTERS OF ALPHABETIC WORD INTO CODE NUMBER. CACT 138 C LARGEST POSSIBLE CODE NUMBER IS LESS THAN 2**15. CACT 139 ICODE=0 CACT 140 DO 408 M1=1,6 CACT 141 X=WORD(M1) CACT 142 DO 402 M2=1,38 CACT 143 IF( X.EQ.MAP(M2) ) GO TO 404 CACT 144 402 CONTINUE CACT 145 C ILLEGAL CHARACTER CACT 146 WRITE (LPRINT,9) CACT 147 WRITE (LPRINT,98) K,(WORD(L),L=1,18) CACT 148 GO TO 999 CACT 149 404 M3=(M2-1)*IFAC(M1) CACT 150 IF( M3.GE.MODP ) M3=M3-MODP CACT 151 ICODE=ICODE+M3 CACT 152 IF( ICODE.GE.MODP ) ICODE=ICODE-MODP CACT 153 408 CONTINUE CACT 154 GO TO 510 CACT 155 C CACT 156 410 READ (LSCRAT,4) INTPAR CACT 157 ISUB=PARNO+NT CACT 158 LPAR(ISUB)=INTPAR CACT 159 GO TO 430 CACT 160 C CACT 161 420 READ (LSCRAT,5) DECPAR CACT 162 ISUB=PARNO+NT CACT 163 PAR(ISUB)=DECPAR CACT 164 C CACT 165 430 NT=NT+1 CACT 166 IF(NT.EQ.NTOKNS) XTYPE=1 CACT 167 GO TO 1000 CACT 168 C CACT 169 510 TABNO = IPARAM(1) CACT 170 MA = MTAB(TABNO) CACT 171 MB = MTAB(TABNO+1) - 1 CACT 172 C CACT 173 LMISSW = 0 CACT 174 DO 700 M=MA,MB CACT 175 IF( ICODE.NE.TAB(1,M) ) GO TO 700 CACT 176 LMISSW = 1 CACT 177 600 XTYPE=1 CACT 178 NT=0 CACT 179 IPARAM(1) = 1 CACT 180 PARNO=TAB(3,M) CACT 181 LTEMP=TAB(2,M) CACT 182 NTOKNS=TAB(4,M) CACT 183 IF(LTEMP.GT.2) NTOKNS=0 CACT 184 GO TO(610,620,630,640,650,660), LTEMP CACT 185 C CACT 186 C NAME OF INTEGER PARAMETER CACT 187 610 XTYPE=2 CACT 188 GO TO 700 CACT 189 C CACT 190 C NAME OF DECIMAL PARAMETER CACT 191 620 XTYPE=3 CACT 192 GO TO 700 CACT 193 C CACT 194 C IMPLICITLY SPECIFIED INTEGER PARAMETER CACT 195 C A SINGLE IMPLICITLY SPECIFIED PARAMTER CAN IN PRINCIPLE HOLD CACT 196 C AS MANY AS THREE PARAMTERS WHICH THE MAIN PROGRAM CONSIDERS CACT 197 C CONCEPTUALLY DISTINCT. ONE IS IN THE UNITS POSITION, ANOTHER IN CACT 198 C IN THE HUNDREDS POSITION, AND ANOTHER IN THE TEN-THOUSANDS CACT 199 C POSITION. CACT 200 630 LP = 1 CACT 201 LT=TAB(4,M) CACT 202 IF(LT.GE.100) LP = 100 CACT 203 IF(LT.GE.10000) LP = 10000 CACT 204 LQ = 100*LP CACT 205 LA = LPAR(PARNO) CACT 206 LPAR(PARNO)= LA-(MOD(LA,LQ)/LP)*LP +TAB(4,M) CACT 207 GO TO 700 CACT 208 C CACT 209 C IMPLICITLY SPECIFIED DECIMAL PARAMETER CACT 210 640 PAR(PARNO)=ATAB(4,M) CACT 211 GO TO 700 CACT 212 C CACT 213 C INTERNAL PARAMETER OF CCACT PROGRAM CACT 214 650 IPARAM(PARNO)=TAB(4,M) CACT 215 GO TO 700 CACT 216 C CACT 217 C IMPLICITLY SPECIFIED INTEGER PARAMETER CACT 218 C DETERMINES PROGRAM FLOW UPON RETURN FROM CCACT--SO ONLY ONE OF CACT 219 C THIS TYPE PER CONTROL CARD ALLOWED. CACT 220 660 IF(NOPCC)9999,665,996 CACT 221 665 NOPCC=1 CACT 222 LPAR(PARNO)=TAB(4,M) CACT 223 GO TO 700 CACT 224 C CACT 225 700 CONTINUE CACT 226 IF(LMISSW.EQ.0) GO TO 994 CACT 227 1000 CONTINUE CACT 228 IF(NT.LT.NTOKNS) GO TO 997 CACT 229 C CACT 230 1001 RETURN CACT 231 C CACT 232 994 WRITE(LPRINT,9) CACT 233 WRITE(LPRINT,998) (WORD(L),L=1,18) CACT 234 998 FORMAT(27H UNDEFINED CONTROL PHRASE ,18A1) CACT 235 GO TO 999 CACT 236 995 WRITE (LPRINT,9) CACT 237 IF(TYPE .LT.4) READ(LSCRAT,3) (WORD(L),L=1,18) CACT 238 IF (TYPE.NE.4) GO TO 9000 CACT 239 WORD(1)=DOLLAR CACT 240 DO 900 L=2,18 CACT 241 WORD(L)=BLANK CACT 242 900 CONTINUE CACT 243 9000 WRITE (LPRINT,99) K,XTYPE,TYPE, (WORD(L),L=1,18) CACT 244 WRITE(LPRINT,96) CACT 245 96 FORMAT(34H IS THIS ITEM SPELLED INCORRECTLY ) CACT 246 GO TO 999 CACT 247 996 WRITE(LPRINT,9) CACT 248 WRITE(LPRINT,6) CACT 249 6 FORMAT(69H0NO MORE THAN ONE CONTROL PHRASE OF TYPE 6 ALLOWED ON A CACT 250 .CONTROL CARD.) CACT 251 GO TO 999 CACT 252 997 WRITE(LPRINT,9) CACT 253 WRITE(LPRINT,97) CACT 254 97 FORMAT(47H CONTROL PHRASE NOT COMPLETED ON A SINGLE CARD.) CACT 255 999 STOP CACT 256 C CACT 257 9999 WRITE(LPRINT,7) CACT 258 7 FORMAT(23H0ERROR EXIT FROM CCACT.) CACT 259 GO TO 999 CACT 260 END CACT 261 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 KYST VERSION 2 OCTOBER, 1976 CONFIG 6 C ROUTINE TO OBTAIN A METRIC MULTIDIMENSIONAL SCALING SOLUTION CONFIG 7 C CONFIG 8 DIMENSION DATA(1800),WW(1800),IJ(1800),X(100,6) CONFIG 9 DIMENSION STORE(494),RWMEAN(100),ROOTS(6),XBEST(100,6), CONFIG10 . DATAIN(5430) CONFIG11 REAL MEAN CONFIG12 COMMON /KYST1/ DATA,WW,X,IJ CONFIG13 COMMON /KYST2/ STORE,RWMEAN,ROOTS,XBEST,DATAIN CONFIG14 C CONFIG15 C PARAMETERS CONFIG16 C X - CONFIGURATION AT OUTPUT CONFIG17 C N - NUMBER OF POINTS IN SCALING CONFIG18 C ND - NUMBER OF DIMENSIONS OF SOLUTION CONFIG19 C DATAIN - ON INPUT, DATA FROM INICON. DESTROYED BY CONFIG CONFIG20 C RWMEAN - ROW AVERAGES CONFIG21 C ROOTS - ON RETURN FROM SGEV, THE FIRST ND EIGENVALUES CONFIG22 C CONFIG23 ISUB(I1,I2)=((I2-1)*(I2-2))/2+I1 CONFIG24 MEAN=0.0 CONFIG25 FN=N CONFIG26 C CONFIG27 C COMPUTE ROW,COLUMN, AND OVERALL AVERAGES OF DATIN VALUES SQUARED CONFIG28 C CONFIG29 DO 130 J=1,N CONFIG30 RWMEAN(J)=0.0 CONFIG31 DO 129 I=1,N CONFIG32 IF(I-J) 126,129,127 CONFIG33 126 IJC=ISUB(I,J) CONFIG34 GO TO 128 CONFIG35 127 IJC=ISUB(J,I) CONFIG36 128 RWMEAN(J)=DATAIN(IJC)**2+RWMEAN(J) CONFIG37 129 CONTINUE CONFIG38 130 MEAN=RWMEAN(J)+MEAN CONFIG39 MEAN=MEAN/FN**2 CONFIG40 C CONFIG41 C CONVERT DATAIN VALUES TO SCALAR-PRODUCT-LIKE VALUES BY THE CONFIG42 C HOUSEHOLDER-TORGERSON PROCEDURE. WORK FROM LAST ROW AND CONFIG43 C COLUMN BACKWARDS. CONFIG44 C CONFIG45 DO 135 JJ=1,N CONFIG46 J=N-JJ+1 CONFIG47 JM1=J-1 CONFIG48 J1=(JM1*(J-2))/2 CONFIG49 J2=(J*JM1)/2 CONFIG50 DO 135 II=1,J CONFIG51 I=J-II+1 CONFIG52 IJ1=J1+I CONFIG53 IJ2=J2+I CONFIG54 IF(I.EQ.J) GO TO 133 CONFIG55 C HOUSEHOLDER-TORGERSON FORMULAS CONFIG56 DATAIN(IJ2)=((RWMEAN(I)+RWMEAN(J))/FN-MEAN-DATAIN(IJ1)**2)/2.0 CONFIG57 GO TO 135 CONFIG58 133 DATAIN(IJ2)=((2.0*RWMEAN(I))/FN-MEAN)/2.0 CONFIG59 135 CONTINUE CONFIG60 C CONFIG61 C OBTAIN EIGENROOTS AND EIGENVECTORS CONFIG62 CALL SGEV(N,ND,X) CONFIG63 C CONFIG64 C COMPUTE CONFIGURATION CONFIG65 C CONFIG66 DO 136 I=1,ND CONFIG67 136 ROOTS(I)=SQRT(ABS(ROOTS(I))) CONFIG68 DO 137 I=1,N CONFIG69 DO 137 J=1,ND CONFIG70 137 X(I,J)=X(I,J)*ROOTS(J) CONFIG71 RETURN CONFIG72 END CONFIG73 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 KYST VERSION 2 OCTOBER, 1976 DATAPR 5 C DATAPR 6 C DATA PRINT DATAPR 7 C DATAPR 8 DIMENSION DATA(1800),WW(1800),IJ(1800),X(100,6) DATAPR 9 DIMENSION GR(100,6),GL(100,6),STORE(1830),DIST(1800),DHAT(1800), DATAPR10 . P(50) DATAPR11 INTEGER Q(50),GRNO(1) DATAPR12 COMMON /KYST1/ DATA,WW,X,IJ DATAPR13 COMMON /KYST2/ GR,GL,STORE,DIST,DHAT DATAPR14 COMMON /MDSCL1/ LREAD, LPRINT, LPUNCH, LSCRAT DATAPR15 EQUIVALENCE(P(1),Q(1),STORE(1)) DATAPR16 DATA KPACK,IJPACK /10000,100/ DATAPR17 C DATAPR18 1 FORMAT( 1H1, 4(30H I J DATA WGHT GP NO ) ) DATAPR19 2 FORMAT(1H1, 3(40H I J DATA DIST DHAT WGHT GP NO ) ) DATAPR20 3 FORMAT(1X,4(2I4,F9.3,F5.2,I3,I4,1H,)) DATAPR21 4 FORMAT(1X,3(2I4,F9.3,3F5.2,I3,I4,1H,)) DATAPR22 C DATAPR23 100 NCOPA=3 DATAPR24 IF(LCODE.EQ.1) NCOPA=4 DATAPR25 C DATAPR26 MMLEFT=MM DATAPR27 200 MMNOW=MIN0( MMLEFT, 50*NCOPA ) DATAPR28 MMUSED=MM-MMLEFT DATAPR29 MMLEFT=MMLEFT-MMNOW DATAPR30 IF(LCODE.EQ.1) WRITE (LPRINT,1) DATAPR31 IF(LCODE.EQ.2) WRITE (LPRINT,2) DATAPR32 C DATAPR33 DO 300 LINENO=1,50 DATAPR34 NCONOW=(MMNOW-LINENO+50)/50 DATAPR35 IF(NCONOW.LE.0) GO TO 310 DATAPR36 C DATAPR37 L=0 DATAPR38 DO 270 LC=1,NCONOW DATAPR39 M = MMUSED + 50*(LC-1) + LINENO DATAPR40 C DATAPR41 IJM=MOD(IJ(M),KPACK) DATAPR42 I=IJM/IJPACK+1 DATAPR43 J=MOD(IJM,IJPACK)+1 DATAPR44 C DATAPR45 L=L+1 DATAPR46 Q(L)=I DATAPR47 L=L+1 DATAPR48 Q(L)=J DATAPR49 C DATAPR50 L=L+1 DATAPR51 P(L)=DATA(M) DATAPR52 C DATAPR53 IF(LCODE.EQ.1) GO TO 250 DATAPR54 C DATAPR55 L=L+1 DATAPR56 P(L)=DIST(M) DATAPR57 L=L+1 DATAPR58 P(L)=DHAT(M) DATAPR59 C DATAPR60 250 CONTINUE DATAPR61 C DATAPR62 L=L+1 DATAPR63 P(L)=WW(M) DATAPR64 C DATAPR65 DO 220 NG=1,NOGRPS DATAPR66 IF(M.GE.GRNO(NG)) NGM=NG DATAPR67 220 CONTINUE DATAPR68 L=L+1 DATAPR69 Q(L)=NGM DATAPR70 C DATAPR71 L=L+1 DATAPR72 Q(L)=M-GRNO(NGM)+1 DATAPR73 C DATAPR74 270 CONTINUE DATAPR75 C DATAPR76 LL=L DATAPR77 IF(LCODE.EQ.1) WRITE (LPRINT,3) (P(L),L=1,LL) DATAPR78 IF(LCODE.EQ.2) WRITE (LPRINT,4) (P(L),L=1,LL) DATAPR79 300 CONTINUE DATAPR80 C DATAPR81 310 IF(MMLEFT.GT.0) GO TO 200 DATAPR82 C DATAPR83 RETURN DATAPR84 C DATAPR85 END DATAPR86 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 KYST VERSION 2 OCTOBER, 1976 DFLT 7 C DFLT 8 C G.GOLUB AND W.KAHAN, CALCULATING THE DFLT 9 C SINGULAR VALUES AND PSEUDO-INVERSE OF A DFLT 10 C MATRIX. J.SIAM NUMER. ANAL. 2,205-224(1965). DFLT 11 C DFLT IS CALLED BY SGEV, AND DFLT 12 REAL A(1),B(1),E,X(1) DFLT 13 C DFLT 14 C DEFLATES THE SYMMETRIC TRIDIAGONAL MATRIX DFLT 15 C WITH DIAGONAL ELEMENTS A(1) THROUGH A(N) DFLT 16 C AND SUB- AND SUPER-DIAGONAL ELEMENTS B(1) DFLT 17 C THROUGH B(N-1) INTO A SYMMETRIC TRIDIAGONAL DFLT 18 C MATRIX OF ORDER N-1. E IS THE EIGENVALUE DFLT 19 C BEING REMOVED AND X AN ASSOCIATED EIGEN- DFLT 20 C VECTOR. THE TANGENTS OF HALF THE ROTATION DFLT 21 C ANGLES OVERWRITE THE FIRST (N-1) COMPONENTS DFLT 22 C OF X. DFLT 23 C DFLT 24 REAL H,W,F,V,T,Q,S,C,AINEW DFLT 25 H=A(1)-E DFLT 26 W=B(1) DFLT 27 N1=N-1 DFLT 28 IF(N1.LT.1)GOTO 30 DFLT 29 DO 20 I=1,N1 DFLT 30 F=X(I) DFLT 31 V=X(I+1) DFLT 32 IF(AMAX1(ABS(H),ABS(W)).LE.AMAX1(ABS(F),ABS(V)))GOTO 10 DFLT 33 F=W DFLT 34 V=-H DFLT 35 10 T=0.E0 DFLT 36 IF(F*F.EQ.0.E0)GOTO 20 DFLT 37 T=F DFLT 38 IF(V.GE.0.E0)T=-T DFLT 39 T=T/(SQRT(F*F+V*V)+ABS(V)) DFLT 40 S=(T+T)/(1.E0+T*T) DFLT 41 C=(1.E0-T)*(1.E0+T)/(1.E0+T*T) DFLT 42 IF(I.GT.1)B(I-1)=C*H+S*W DFLT 43 F=C*A(I)+S*B(I) DFLT 44 Q=C*B(I)+S*A(I+1) DFLT 45 W=S*B(I+1) DFLT 46 B(I+1)=C*B(I+1) DFLT 47 AINEW=F*C+Q*S DFLT 48 H=Q*C-F*S DFLT 49 A(I+1)=A(I)+A(I+1)-AINEW DFLT 50 A(I)=AINEW DFLT 51 X(I+1)=C*X(I+1)-S*X(I) DFLT 52 20 X(I)=T DFLT 53 B(N-1)=0.E0 DFLT 54 30 RETURN DFLT 55 END DFLT 56 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 KYST VERSION 2 OCTOBER, 1976 DTRAN 5 C DTRAN 6 C DTRAN--DATA TRANSFORMATION. DTRAN 7 DIMENSION DUMMY(28),DUM(6) DTRAN 8 COMMON /MDSCL2/ DUMMY,A,B,P,S,T,WCON1,WCON2,WCON3,WCON4,WCON5,DUM DTRAN 9 DTRAN=S+T*((A+B*XX)**P) DTRAN 10 RETURN DTRAN 11 END DTRAN 12 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 C KYST VERSION 2 OCTOBER, 1976 EQSOLV 5 CEQSOLV SOLVES SIMULTANEOUS LINEAR EQUATIONS EQSOLV 6 REAL A(5,5), B(5), C(5) EQSOLV 7 C EQSOLV 8 NTM1=NT-1 EQSOLV 9 IF(NT.EQ.NF) GO TO 60 EQSOLV10 DO 50 II=NF,NTM1 EQSOLV11 IIP1=II+1 EQSOLV12 DO 40 I=IIP1,NT EQSOLV13 E = A(I,II)/A(II,II) EQSOLV14 DO 30 J=NF,NT EQSOLV15 A(I,J)=A(I,J)-E*A(II,J) EQSOLV16 30 CONTINUE EQSOLV17 B(I)=B(I)-E*B(II) EQSOLV18 40 CONTINUE EQSOLV19 50 CONTINUE EQSOLV20 C EQSOLV21 60 DO 90 IC=NF,NT EQSOLV22 I=NT+NF-IC EQSOLV23 E=B(I) EQSOLV24 IF(I.EQ.NT) GO TO 80 EQSOLV25 IP1=I+1 EQSOLV26 DO 70 J=IP1,NT EQSOLV27 E=E-C(J)*A(I,J) EQSOLV28 70 CONTINUE EQSOLV29 80 C(I) = E / A(I,I) EQSOLV30 90 CONTINUE EQSOLV31 RETURN EQSOLV32 END EQSOLV33 CFITM FITM 1 SUBROUTINE FITM(ISTRT,MM,LFITSW) FITM 2 C FITM FOR KYST JANUARY,1973 FITM 3 C WRITTEN BY J.KRUSKAL FITM 4 C MODIFIED FOR KYST BY J.KRUSKAL AND J.SEERY JANUARY,1973 FITM 5 C KYST VERSION 2 OCTOBER, 1976 FITM 6 C FITM 7 C FITM PERFORMS WEIGHTED-LEAST-SQUARES MONOTONE REGRESSION FITM 8 C THIS ROUTINE FINDS THE VALUES FOR DHAT WHICH ARE MONOTONIC FITM 9 C AND WHICH BEST APPROXIMATE THE VALUES OF DIST, FITM 10 C IN THE SENSE THAT THE SUM OF THE SQUARED DEVIATIONS, FITM 11 C EACH ONE WEIGHTED BY THE CORRESPONDING VALUE IN WW , FITM 12 C IS A MINIMUM. FITM 13 C IT DIRECTLY USES SORT. FITM 14 C FITM 15 DIMENSION DATA(1800),IJ(1800),WW(1800),DIST(1800),DHAT(1800), FITM 16 . LBLOCK(1800),X(100,6),DUMMY(30),GR(100,6),GL(100,6) FITM 17 INTEGER SDSWIT FITM 18 COMMON /MDSCL1/ LREAD, LPRINT, LPUNCH, LSCRAT FITM 19 COMMON /KYST1/ DATA,WW,X,IJ FITM 20 COMMON /KYST2/ GR,GL,LBLOCK,DUMMY,DIST,DHAT FITM 21 C FITM 22 DATA KPACK/10000/ FITM 23 C FITM 24 C FORM FIRST APPROXIMATION TO CORRECT PARTITON FITM 25 C FITM 26 C IF LFITSW=1, USE PRIMARY APPROACH. THUS WE SORT EACH FITM 27 C BLOCK OF EQUAL VALUES OF DATA ACCORDING TO DIST VALUES. FITM 28 C THEN WE CREATE PARTITON INTO BLOCKS OF SIZE 1 TO START. FITM 29 C FITM 30 C IF LFITSW=2, USE SECONDARY APPROACH. THUS WE START WITHFITM 31 C PARTITION INTO BLOCKS OF EQUAL DATA VALUES. FITM 32 C FITM 33 C WITHIN EACH BLOCK OF WHATEVER SIZE, THE FIRST DHAT VALUEFITM 34 C GIVES THE WEIGHTED TOTAL OF THE DIST VALUES IN THAT BLOCK, WHILE FITM 35 C THE LAST DHAT VALUE GIVES THE SUM OF THE WEIGHTS FOR THATFITM 36 C BLOCK. SIMILARLY, WITHIN EACH BLOCK, THE FIRST LBLOCK FITM 37 C VALUE AND THE LAST LBLOCK VALUE BOTH GIVE THE SIZE OF THEFITM 38 C BLOCK. FITM 39 C BLOCKS OF SIZE 1 FORM A PARTIAL EXCEPTION. EVERYTHIN IS THE SAME FITM 40 C EXCEPT THAT THE SUM OF THE W VALUES IS NOT FOUND IN THE FITM 41 C LAST DHAT VALUE IN THE BLOCK. FITM 42 C FITM 43 100 MA=ISTRT FITM 44 110 KSUB=MA-1 FITM 45 115 KSUB=KSUB+1 FITM 46 KK=IJ(KSUB)-KPACK FITM 47 IF(KK) 115,120,120 FITM 48 120 MB=KSUB FITM 49 K=MB-MA+1 FITM 50 GO TO ( 200, 300 ), LFITSW FITM 51 C FITM 52 C PRIMARY APPROACH FITM 53 C FITM 54 200 IF(K-1) 9999, 220, 210 FITM 55 C IF K=1, SAVE SORTING TIME FITM 56 C FITM 57 210 IJ(MB)=IJ(MB)-KPACK FITM 58 CALL SORT( DIST(MA), K, IJ(MA), DATA(MA),WW(MA),3,+1) FITM 59 IJ(MB)=IJ(MB)+KPACK FITM 60 C 'SORT' WILL SORT THE K ELEMENTS OF 'DIST' IN ALGEBRAIC FITM 61 C ORDER, ASCENDING FITM 62 C THE ELEMENTS OF'IJ' AND 'DATA' WILL BE REARRANGED IN FITM 63 C EXACTLY THE SAME ORDER AS THOSE OF 'DIST'. FITM 64 C ALSO THE ELEMENTS OF 'WW'. FITM 65 C FITM 66 220 DO 230 M=MA,MB FITM 67 DHAT(M) = DIST(M) * WW(M) FITM 68 230 LBLOCK(M)=1 FITM 69 GO TO 400 FITM 70 C FITM 71 C SECONDARY APPROACH FITM 72 C FITM 73 C FITM 74 300 LBLOCK(MA)=K FITM 75 LBLOCK(MB)=K FITM 76 TEMP1=0.0 FITM 77 TEMP2 = 0.0 FITM 78 DO 310 M=MA,MB FITM 79 TEMP1 = TEMP1 + DIST(M) * WW(M) FITM 80 TEMP2 = TEMP2 + WW(M) FITM 81 310 CONTINUE FITM 82 DHAT(MA)=TEMP1 FITM 83 IF(K-1) 9999, 400, 320 FITM 84 320 DHAT(MB) = TEMP2 FITM 85 C PROCEED TO NEXT BOCK. QUERY. IS IT THE LAST FITM 86 C FITM 87 400 MA=MA+K FITM 88 IF(MA-MM-1) 110, 500, 9999 FITM 89 C FITM 90 C START MAIN COMPUTATION *************************************FITM 91 C FITM 92 500 MA=ISTRT FITM 93 510 LUD=2 FITM 94 NSATIS=0 FITM 95 520 K=LBLOCK(MA) FITM 96 MB=MA+K-1 FITM 97 IF(K-1) 9999, 530, 540 FITM 98 530 WT =WW(MB) FITM 99 GO TO 550 FITM 100 540 WT = DHAT(MB) FITM 101 550 DAV = DHAT(MA) / WT FITM 102 GO TO ( 600, 700 ), LUD FITM 103 C FITM 104 C IS BLOCK DOWN-SATISFIED. IF NOT, MERGE FITM 105 C FITM 106 600 IF(MA-ISTRT) 9999,630,610 FITM 107 610 MBD=MA-1 FITM 108 KD=LBLOCK(MBD) FITM 109 MAD=MBD-KD+1 FITM 110 IF(KD-1) 9999, 613, 615 FITM 111 613 WTD =WW(MBD) FITM 112 GO TO 617 FITM 113 615 WTD = DHAT(MBD) FITM 114 617 DAVD = DHAT(MAD) / WTD FITM 115 IF( DAVD-DAV ) 630, 620, 620 FITM 116 620 KNEW=K+KD FITM 117 LBLOCK(MAD)=KNEW FITM 118 LBLOCK(MB)=KNEW FITM 119 DTONEW = DHAT(MAD) + DHAT(MA) FITM 120 DHAT(MAD) = DTONEW FITM 121 DHAT(MB) = WT + WTD FITM 122 NSATIS=0 FITM 123 MA=MAD FITM 124 GO TO 800 FITM 125 630 NSATIS=NSATIS+1 FITM 126 GO TO 800 FITM 127 C FITM 128 C IS BLOCK UP-SATISFIED. IF NOT, MERGE FITM 129 C FITM 130 700 IF(MB-MM) 710, 730, 9999 FITM 131 710 MAU=MB+1 FITM 132 KU=LBLOCK(MAU) FITM 133 MBU=MAU+KU-1 FITM 134 IF(KU-1) 9999, 713, 715 FITM 135 713 WTU =WW(MBU) FITM 136 GO TO 717 FITM 137 715 WTU = DHAT(MBU) FITM 138 717 DAVU = DHAT(MAU) / WTU FITM 139 IF( DAV-DAVU ) 730, 720, 720 FITM 140 720 KNEW=K+KU FITM 141 LBLOCK(MA)=KNEW FITM 142 LBLOCK(MBU)=KNEW FITM 143 DTONEW = DHAT(MA) + DHAT(MAU) FITM 144 DHAT(MA) = DTONEW FITM 145 DHAT(MBU) = WT + WTU FITM 146 NSATIS=0 FITM 147 GO TO 800 FITM 148 730 NSATIS=NSATIS+1 FITM 149 GO TO 800 FITM 150 C FITM 151 C PROCEED TO NEXT BLOCK IF READY. FITM 152 C FITM 153 800 LUD = 3-LUD FITM 154 C QUERY. IS BLOCK BOTH UP AND DOWN SATISFIED. IF NOT, RETURFITM 155 IF(NSATIS-1) 520, 520, 810 FITM 156 C QUERY. IS THIS LAST BLOCK. IF NOT, GO ON TO NEXT BLOCK. FITM 157 810 IF(MB-MM) 820, 900, 9999 FITM 158 820 MA=MB+1 FITM 159 GO TO 510 FITM 160 C FITM 161 C MAIN COMPUTATION COMPLETE. PLACE ANSWERS IN DHAT. FITM 162 C FITM 163 900 MA=ISTRT FITM 164 910 K=LBLOCK(MA) FITM 165 MB=MA+K-1 FITM 166 IF(K-1) 9999, 940, 920 FITM 167 920 TEMP1 = DHAT(MA) / DHAT(MB) FITM 168 DO 930 M=MA,MB FITM 169 930 DHAT(M)=TEMP1 FITM 170 GO TO 945 FITM 171 940 DHAT(MA) = DIST(MA) FITM 172 945 MA = MB + 1 FITM 173 IF(MA-MM-1) 910, 950, 9999 FITM 174 950 RETURN FITM 175 C FITM 176 C TROUBLE EXIT FITM 177 C FITM 178 9999 WRITE (LPRINT, 99) FITM 179 99 FORMAT(50H0KRUSKAL. IMPOSSIBLE BRANCH TAKEN ON IF STATEMENT. ) FITM 180 STOP FITM 181 C FITM 182 END FITM 183 CFITMS FITMS 1 SUBROUTINE FITMS(ISTRT,MM) FITMS 2 C FITMS FOR KYST JANUARY,1973 FITMS 3 C WRITTEN FOR KYST BY J.SEERY JANUARY,1973 FITMS 4 C KYST VERSION 2 OCTOBER, 1976 FITMS 5 C FITMS 6 C FITMS PERFORMS LEAST-SQUARES MONOTONE REGRESSION FITMS 7 C THIS ROUTINE FINDS THE VALUES FOR DHATIN WHICH ARE MONOTONIC FITMS 8 C AND WHICH BEST APPROXIMATE THE VALUES OF DISTIN IN THE SENSE FITMS 9 C THAT THE SUM OF THE SQUARED DEVIATIONS IS A MINIMUM. FITMS 10 C CALLED ONLY FROM INICON FITMS 11 C ASSUMES NO TIES. DISTIN AND DHATIN SHARE THE SAME LOCATIONS. FITMS 12 C FITMS 13 DIMENSION GR(100,6),GL(100,6),LBLOCK(1800),DUMMY(30),DHATIN(1800),FITMS 14 . IJIN(1800),DISTIN(1800) FITMS 15 COMMON /MDSCL1/ LREAD, LPRINT, LPUNCH, LSCRAT FITMS 16 COMMON /KYST2/ GR,GL,LBLOCK,DUMMY,DISTIN,IJIN FITMS 17 EQUIVALENCE(DISTIN(1),DHATIN(1)) FITMS 18 C FITMS 19 C WITHIN EACH BLOCK, THE FIRST LBLOCK VALUE AND THE LAST LBLOCK FITMS 20 C VALUE BOTH GIVE THE SIZE OF THE BLOCK. FITMS 21 C FITMS 22 C FORM FIRST APPROXIMATION TO CORRECT PARTITON FITMS 23 C INITIALIZE LBLOCK TO 1 FITMS 24 DO 230 M=ISTRT,MM FITMS 25 230 LBLOCK(M)=1 FITMS 26 C FITMS 27 C START MAIN COMPUTATION *************************************FITMS 28 C FITMS 29 MA=ISTRT FITMS 30 510 LUD=2 FITMS 31 NSATIS=0 FITMS 32 520 K=LBLOCK(MA) FITMS 33 MB=MA+K-1 FITMS 34 WT=K FITMS 35 550 DAV=DHATIN(MA)/WT FITMS 36 GO TO ( 600, 700 ), LUD FITMS 37 C FITMS 38 C IS BLOCK DOWN-SATISFIED. IF NOT, MERGE FITMS 39 C FITMS 40 600 IF(MA-ISTRT)9999,630,610 FITMS 41 610 MBD=MA-1 FITMS 42 KD=LBLOCK(MBD) FITMS 43 MAD=MBD-KD+1 FITMS 44 WTD=KD FITMS 45 617 DAVD=DHATIN(MAD)/WTD FITMS 46 IF( DAVD-DAV ) 630, 620, 620 FITMS 47 620 KNEW=K+KD FITMS 48 LBLOCK(MAD)=KNEW FITMS 49 LBLOCK(MB)=KNEW FITMS 50 DTONEW=DHATIN(MAD)+DHATIN(MA) FITMS 51 DHATIN(MAD)=DTONEW FITMS 52 NSATIS=0 FITMS 53 MA=MAD FITMS 54 GO TO 800 FITMS 55 630 NSATIS=NSATIS+1 FITMS 56 GO TO 800 FITMS 57 C FITMS 58 C IS BLOCK UP-SATISFIED. IF NOT, MERGE FITMS 59 C FITMS 60 700 IF(MB-MM) 710, 730, 9999 FITMS 61 710 MAU=MB+1 FITMS 62 KU=LBLOCK(MAU) FITMS 63 MBU=MAU+KU-1 FITMS 64 WTU=KU FITMS 65 717 DAVU=DHATIN(MAU)/WTU FITMS 66 IF( DAV-DAVU ) 730, 720, 720 FITMS 67 720 KNEW=K+KU FITMS 68 LBLOCK(MA)=KNEW FITMS 69 LBLOCK(MBU)=KNEW FITMS 70 DTONEW=DHATIN(MA)+DHATIN(MAU) FITMS 71 DHATIN(MA)=DTONEW FITMS 72 NSATIS=0 FITMS 73 GO TO 800 FITMS 74 730 NSATIS=NSATIS+1 FITMS 75 GO TO 800 FITMS 76 C FITMS 77 C PROCEED TO NEXT BLOCK IF READY. FITMS 78 C FITMS 79 800 LUD = 3-LUD FITMS 80 C QUERY. IS BLOCK BOTH UP AND DOWN SATISFIED. IF NOT, RETURFITMS 81 IF(NSATIS-1) 520, 520, 810 FITMS 82 C QUERY. IS THIS LAST BLOCK. IF NOT, GO ON TO NEXT BLOCK. FITMS 83 810 IF(MB-MM) 820, 900, 9999 FITMS 84 820 MA=MB+1 FITMS 85 GO TO 510 FITMS 86 C FITMS 87 C MAIN COMPUTATION COMPLETE. PLACE ANSWERS IN DHATIN. FITMS 88 C FITMS 89 900 MA=ISTRT FITMS 90 910 K=LBLOCK(MA) FITMS 91 MB=MA+K-1 FITMS 92 IF(K-1) 9999, 940, 920 FITMS 93 920 XK=K FITMS 94 TEMP1=DHATIN(MA)/XK FITMS 95 DO 930 M=MA,MB FITMS 96 930 DHATIN(M)=TEMP1 FITMS 97 940 MA = MB + 1 FITMS 98 IF(MA-MM-1) 910, 950, 9999 FITMS 99 950 RETURN FITMS100 C FITMS101 C TROUBLE EXIT FITMS102 C FITMS103 9999 WRITE (LPRINT, 99) FITMS104 99 FORMAT(34H0ERROR EXIT FROM SUBROUTINE FITMS.) FITMS105 STOP FITMS106 END FITMS107 CFITP FITP 1 SUBROUTINE FITP(ISTRT,MM,NTI,K) FITP 2 C FITP FOR KYST JANUARY,1973 FITP 3 C WRITTEN BY J.KRUSKAL FITP 4 C KYST VERSION 2 OCTOBER, 1976 FITP 5 C FITP--LEAST SQUARES POLYNOMIAL REGRESSION FITP 6 C FITP 7 DIMENSION DATA(1800),WW(1800),IJ(1800),X(100,6) FITP 8 DIMENSION GR(100,6),GL(100,6),LBLOCK(1800),C(30),DIST(1800), FITP 9 . DHAT(1800),A(5,5),B(5),V(5) FITP 10 COMMON /KYST1/ DATA,WW,X,IJ FITP 11 COMMON /KYST2/ GR,GL,LBLOCK,C,DIST,DHAT FITP 12 C FITP 13 NT = MAX0( 1, MIN0( 5,NTI ) ) FITP 14 C FITP 15 SUMW=0.0 FITP 16 DO 1105 I=K,NT FITP 17 B(I) = 0.0 FITP 18 DO 1105 J=K,I FITP 19 A(I,J) = 0.0 FITP 20 1105 CONTINUE FITP 21 C FITP 22 DO 1150 M=ISTRT,MM FITP 23 DA=DATA(M) FITP 24 SUMW=SUMW+WW(M) FITP 25 C FITP 26 DO 1110 I=K,NT FITP 27 V(I)=REGR(DA,I) FITP 28 1110 CONTINUE FITP 29 C FITP 30 DO 1140 I=K,NT FITP 31 DO 1130 J=K,I FITP 32 A(I,J)=A(I,J)+V(I)*V(J)*WW(M) FITP 33 1130 CONTINUE FITP 34 B(I)=B(I)+V(I)*DIST(M)*WW(M) FITP 35 1140 CONTINUE FITP 36 C FITP 37 1150 CONTINUE FITP 38 C FITP 39 C FITP 40 DO 1155 I = K,NT FITP 41 DO 1155 J= K,I FITP 42 A(J,I) = A(I,J) FITP 43 1155 CONTINUE FITP 44 C FITP 45 CALL EQSOLV(A,B,C,NT,K) FITP 46 C FITP 47 DO 1170 M=ISTRT,MM FITP 48 DA = DATA(M) FITP 49 TEMP=0.0 FITP 50 DO 1160 I=K,NT FITP 51 TEMP=TEMP+C(I)*REGR(DA,I) FITP 52 1160 CONTINUE FITP 53 DHAT(M)=TEMP FITP 54 1170 CONTINUE FITP 55 C FITP 56 RETURN FITP 57 END FITP 58 CINICON INICON 1 SUBROUTINE INICON(N,ND,M1,NITER,LHIPRT,SDS,LSCH,SPREAD) INICON 2 C INICON FOR KYST JANUARY,1973 INICON 3 C WRITTEN BY F.YOUNG INICON 4 C MODIFIED FOR KYST BY J.KRUSKAL AND J.SEERY JANUARY,1973 INICON 5 C KYST VERSION 2 OCTOBER, 1976 INICON 6 C INICON 7 C WRITTEN BY FORREST W. YOUNG JUNE 1971 INICON 8 C ROUTINE TO COMPUTE INITIAL CONFIGURATION INICON 9 C USES THE TORGERSON-YOUNG PROCEDURE TO DERIVE THE BEST FITTING INICON10 C QUASI-NONMETRIC EUCLIDIAN CONFIGURATION INICON11 C INICON12 DIMENSION DATA(1800),WW(1800),IJ(1800),X(100,6) INICON13 DIMENSION GR(100,6),XBEST(100,6),DATAIN(5430),IJIN(1800), INICON14 . DHATIN(1800),DUMMY(1) INICON15 COMMON /ACCUR/ PRECSN,XMAG,XMAXN,FLPRIN,CHPRIN INICON16 COMMON /MDSCL1/ LREAD,LPRINT,LPUNCH,LSCRAT INICON17 COMMON /KYST1/ DATA,WW,X,IJ INICON18 COMMON /KYST2/ GR,XBEST,DATAIN INICON19 EQUIVALENCE (DHATIN(1),DATAIN(1831)) , (IJIN(1),DATAIN(3631)) INICON20 DATA NCOLS,KPACK,IJPACK /6,10000,100/ INICON21 ISUB(I1,I2)=((I2-1)*(I2-2))/2+I1 INICON22 C INICON23 C PARAMETERS INICON24 C INICON25 C CONSTANTS INICON26 C INICON27 NM1=N-1 INICON28 NM2=N-2 INICON29 N2=(N*NM1)/2 INICON30 FN2=N2 INICON31 MTIMES=NM1/NCOLS+1 INICON32 XMAGM=-XMAXN INICON33 STRBST = 1.0 INICON34 SPRDM1=SPREAD-1.0 INICON35 C INICON36 C CONVERT RAW DATA VALUES IN DATA TO AVERAGED VALUES, AND STORE INICON37 C THEM IN CORRESPONDING POSITIONS IN DATAIN. WORK WITH NCOLS (=6)INICON38 C ROWS AND COLUMNS AT A TIME. FLAG MISSING DATA VALUES. INICON39 C AT THE SAME TIME COMPUTE MAX,MIN, AND AVERAGE OF NON-MISSING INICON40 C DATAIN VALUES INICON41 C INICON42 AVG=0.0 INICON43 NUM=0 INICON44 DINMAX=XMAGM INICON45 DINMIN=1.0/XMAG INICON46 INSUB=0 INICON47 MEND=0 INICON48 DO 30 MT=1,MTIMES INICON49 MSTART=MEND+1 INICON50 MEND=MEND+NCOLS INICON51 MEND=MIN0(MEND,N) INICON52 MNDM1=MEND-1 INICON53 DO 5 L1=1,MNDM1 INICON54 DO 5 L2=1,NCOLS INICON55 X(L1,L2)=0. INICON56 5 XBEST(L1,L2)=0. INICON57 DO 15 L=1,M1 INICON58 IF(WW(L)) 15,15,8 INICON59 8 KK=MOD(IJ(L),KPACK) INICON60 I=KK/IJPACK+1 INICON61 J=MOD(KK,IJPACK)+1 INICON62 K=MIN0(I,J) INICON63 KK=MAX0(I,J) INICON64 IF((KK.LT.MSTART).OR.(KK.GT.MEND)) GO TO 15 INICON65 IF(K.EQ.KK) GO TO 15 INICON66 KK=MOD(KK-1,NCOLS)+1 INICON67 X(K,KK)=X(K,KK)+DATA(L)*WW(L) INICON68 XBEST(K,KK)=XBEST(K,KK)+WW(L) INICON69 15 CONTINUE INICON70 KK=0 INICON71 IF(MSTART.NE.1) GO TO 18 INICON72 MSTART=2 INICON73 KK=1 INICON74 18 DO 25 M=MSTART,MEND INICON75 MM1=M-1 INICON76 KK=KK+1 INICON77 DO 25 K=1,MM1 INICON78 INSUB=INSUB+1 INICON79 IF(XBEST(K,KK)) 20,20,22 INICON80 20 DATAIN(INSUB)=XMAGM INICON81 GO TO 25 INICON82 22 NUM=NUM+1 INICON83 DIN=X(K,KK)/XBEST(K,KK) INICON84 AVG=AVG+DIN INICON85 DATAIN(INSUB)=DIN INICON86 DINMAX=AMAX1(DINMAX,DIN) INICON87 DINMIN=AMIN1(DINMIN,DIN) INICON88 25 CONTINUE INICON89 30 CONTINUE INICON90 FNUM=NUM INICON91 AVG=AVG/FNUM INICON92 C INICON93 C REPLACE MISSING DATAIN VALUES WITH AVERAGE VALUE, AND INICON94 C CONVERT DATAIN TO DISTANCE-LIKE VALUES INICON95 C INICON96 IF(SDS) 35,99,38 INICON97 C HERE IF SIMILARITIES INICON98 35 XFAC=(SPREAD*DINMAX-DINMIN)/SPRDM1 INICON99 GO TO 40 INICON00 C HERE IF DISSIMILARITIES INICON01 38 XFAC=(DINMAX-SPREAD*DINMIN)/SPRDM1 INICON02 40 DO 60 L=1,N2 INICON03 IF(DATAIN(L).GT.XMAGM) GO TO 55 INICON04 DATAIN(L)=XFAC+SDS*AVG INICON05 GO TO 60 INICON06 55 DATAIN(L)=XFAC+SDS*DATAIN(L) INICON07 60 CONTINUE INICON08 C INICON09 C NORMALIZE THE CURRENT DATA INICON10 C INICON11 SSQ=0.0 INICON12 DO 70 I=1,N2 INICON13 70 SSQ=SSQ+DATAIN(I)*DATAIN(I) INICON14 SSQ=SQRT(SSQ/FN2) INICON15 DO 75 L=1,N2 INICON16 75 DATAIN(L)=DATAIN(L)/SSQ INICON17 C INICON18 WRITE(LPRINT,1000) N,ND INICON19 1000 FORMAT(46H0INITIAL CONFIGURATION COMPUTATION. NO. PTS.=,I4, INICON20 . 8H DIM=,I4) INICON21 C INICON22 C COMPUTE THE BEST FITTING METRIC EUCLIDIAN CONFIGURATION FOR THIS DATA INICON23 C INICON24 CALL CONFIG(N,ND) INICON25 C INICON26 C RECONVERT RAW DATA VALUES IN DATA TO AVERAGED VALUES: STORE THEM INICON27 C IN DHATIN, AND THEIR POSITIONS IN IJIN. IF PRE-ITERATIONS ARE INICON28 C TO BE PERFORMED, FLAG THE MISSING DATA AND POSITIONS. OTHERWISE, INICON29 C PACK THE AVERAGED VALUES TOGETHER. INICON30 C INICON31 ISTRT=1 INICON32 INSUB=0 INICON33 MEND=0 INICON34 DO 160 MT=1,MTIMES INICON35 MSTART=MEND+1 INICON36 MEND=MEND+NCOLS INICON37 MEND=MIN0(MEND,N) INICON38 MNDM1=MEND-1 INICON39 DO 130 L1=1,MNDM1 INICON40 DO 130 L2=1,NCOLS INICON41 GR(L1,L2)=0. INICON42 130 XBEST(L1,L2)=0. INICON43 DO 140 L=1,M1 INICON44 IF(WW(L)) 140,140,134 INICON45 134 KK=MOD(IJ(L),KPACK) INICON46 I=KK/IJPACK+1 INICON47 J=MOD(KK,IJPACK)+1 INICON48 K=MIN0(I,J) INICON49 KK=MAX0(I,J) INICON50 IF((KK.LT.MSTART).OR.(KK.GT.MEND)) GO TO 140 INICON51 IF(K.EQ.KK) GO TO 140 INICON52 KK=MOD(KK-1,NCOLS)+1 INICON53 GR(K,KK)=GR(K,KK)+DATA(L)*WW(L) INICON54 XBEST(K,KK)=XBEST(K,KK)+WW(L) INICON55 140 CONTINUE INICON56 KK=0 INICON57 IF(MSTART.NE.1) GO TO 144 INICON58 MSTART=2 INICON59 KK=1 INICON60 144 DO 155 M=MSTART,MEND INICON61 MM1=M-1 INICON62 KK=KK+1 INICON63 DO 155 K=1,MM1 INICON64 IF(XBEST(K,KK)) 146,146,152 INICON65 146 IF(NITER.EQ.0) GO TO 155 INICON66 ISTRT=ISTRT+1 INICON67 INSUB=INSUB+1 INICON68 DHATIN(INSUB)=XMAGM INICON69 IJIN(INSUB)=IJPACK*(K-1)+MM1 INICON70 GO TO 155 INICON71 152 INSUB=INSUB+1 INICON72 C NEGATE DATA VALUES IF SIMILARITIES INICON73 DHATIN(INSUB)=SDS*GR(K,KK)/XBEST(K,KK) INICON74 IJIN(INSUB)=IJPACK*(K-1)+MM1 INICON75 155 CONTINUE INICON76 160 CONTINUE INICON77 NELE=INSUB INICON78 NUM=INSUB-ISTRT+1 INICON79 FNUM=NUM INICON80 C INICON81 C SORT DHATIN INTO ASCENDING ORDER, AND REARRANGE IJIN. INICON82 200 CALL SORT(DHATIN,NELE,IJIN,DUMMY,DUMMY,1,1) INICON83 C INICON84 IF(NITER) 99,228,210 INICON85 C INICON86 C COMPUTE ITERATIVE BEST-FITTING EUCLIDIAN INITIAL CONFIGURATION INICON87 C INICON88 210 IF(LHIPRT.EQ.2) WRITE(LPRINT,1010) INICON89 1010 FORMAT(30H0 PRE-ITERATION STRESS/) INICON90 IT=-1 INICON91 C INICON92 C START OF ITERATIVE COMPUTATION INICON93 C INICON94 220 IT=IT+1 INICON95 C INICON96 C COMPUTE MINKOWSKI DISTANCES INICON97 C INICON98 228 DBAR=0.0 INICON99 DO 235 L=1,NELE INICON00 K=IJIN(L) INICON01 I=K/IJPACK+1 INICON02 J=MOD(K,IJPACK)+1 INICON03 SUM=0.0 INICON04 DO 230 K=1,ND INICON05 230 SUM=SUM+RPOWER(X(I,K)-X(J,K)) INICON06 DHATIN(L)=RROOT(SUM) INICON07 IF(L.LT.ISTRT) GO TO 235 INICON08 DBAR=DBAR+DHATIN(L) INICON09 235 CONTINUE INICON10 C INICON11 C COMPUTE THE BEST-FITTING PSEUDO-DISTANCES USING MONOTONE REGRESSION INICON12 C INICON13 CALL FITMS(ISTRT,NELE) INICON14 C INICON15 C COMPUTE STRESS OF THE CURRENT CONFIGURATION INICON16 C INICON17 U=0.0 INICON18 T=0.0 INICON19 C INICON20 C IF STRESS FORMULA 1 OPTED BY USER THEN DBAR=0 INICON21 C INICON22 DBAR=DBAR/FNUM INICON23 IF(LSCH.EQ.1) DBAR=0.0 INICON24 DO 240 L=ISTRT,NELE INICON25 K=IJIN(L) INICON26 I=K/IJPACK+1 INICON27 J=MOD(K,IJPACK)+1 INICON28 SUM=0.0 INICON29 DO 238 K=1,ND INICON30 238 SUM=SUM+RPOWER(X(I,K)-X(J,K)) INICON31 DISTL=RROOT(SUM) INICON32 U=U+(DISTL-DHATIN(L))**2 INICON33 240 T=T+(DISTL-DBAR)**2 INICON34 STRESS=SQRT(U/T) INICON35 C INICON36 IF(NITER) 99,875,245 INICON37 C INICON38 C WRITE RESULTS OF CURRENT ITERATION INICON39 C INICON40 245 IF(LHIPRT.EQ.2) WRITE(LPRINT,1015) IT,STRESS INICON41 1015 FORMAT(1H ,I15,F13.4) INICON42 C DETERMINE IF THIS IS THE BEST ITERATION SO FAR. INICON43 C INICON44 IF(STRESS) 99,247,248 INICON45 247 STRBST=STRESS INICON46 GO TO 850 INICON47 248 IF(STRESS-STRBST) 250,250,330 INICON48 C INICON49 C THIS IS THE BEST PRE-ITERATION SO FAR. IF THERE ARE MORE PRE- INICON50 C ITERATIONS TO GO, GO ON TO SEE IF A BETTER CONFIGURATION CAN BE INICON51 C FOUND. INICON52 C INICON53 250 STRBST=STRESS INICON54 IF(IT-NITER) 255,290,290 INICON55 C INICON56 C PUT DHATIN ENTRIES INTO CORRECT DATAIN POSITIONS INICON57 C FILL MISSING ENTRIES OF DATAIN WITH THE CORRESPONDING DISTANCES-- INICON58 C WHICH ARE STORED IN THE FIRST POSITIONS OF DHATIN INICON59 C INICON60 255 DO 260 L=1,NELE INICON61 K=IJIN(L) INICON62 I=K/IJPACK+1 INICON63 J=MOD(K,IJPACK)+1 INICON64 M=ISUB(I,J) INICON65 260 DATAIN(M)=DHATIN(L) INICON66 C INICON67 C SAVE THIS CONFIGURATION AS BEST SO FAR INICON68 C INICON69 DO 275 I=1,N INICON70 DO 275 J=1,ND INICON71 275 XBEST(I,J)=X(I,J) INICON72 C INICON73 CALL CONFIG(N,ND) INICON74 C INICON75 C GO TO TO NEXT PRE-ITERATION INICON76 GO TO 220 INICON77 C INICON78 290 WRITE(LPRINT,1020) NITER INICON79 1020 FORMAT(33H0MAXIMUM NUMBER OF PRE-ITERATIONS,I3,10H, REACHED.) INICON80 GO TO 850 INICON81 C INICON82 C IF THIS SOLUTION FITS NOT AS WELL AS THE PREVIOUS ONE, COME HERE INICON83 C TO RECOVER THE BEST SOLUTION. INICON84 C INICON85 330 IT=IT-1 INICON86 WRITE(LPRINT,1025) IT INICON87 1025 FORMAT(73H0STRESS STARTING TO INCREASE. BEST VALUE ACHIEVED ON PREINICON88 .-ITERATION NUMBER, I3) INICON89 DO 31 I=1,N INICON90 DO 31 J=1,ND INICON91 31 X(I,J)=XBEST(I,J) INICON92 C INICON93 850 WRITE(LPRINT,1030) N,ND,STRBST,LSCH INICON94 1030 FORMAT(34H0THE BEST INITIAL CONFIGURATION OF,I4,10H POINTS IN,I4, INICON95 . 22H DIMENSIONS HAS STRESS,F7.3,8H FORMULA,I2) INICON96 GO TO 900 INICON97 C INICON98 C NO PRE-ITERATIONS PERFORMED. INICON99 1005 FORMAT(38H0NO PRE-ITERATIONS PERFORMED. STRESS=,F8.4, INICON00 . 9H ,FORMULA,I2) INICON01 875 WRITE(LPRINT,1005) STRESS ,LSCH INICON02 900 RETURN INICON03 C INICON04 99 WRITE(LPRINT,98) INICON05 98 FORMAT(22H0ERROR EXIT FOR INICON) INICON06 STOP INICON07 END INICON08 CNEWSTP NEWSTP 1 SUBROUTINE NEWSTP( STEP, ITNO, SFGR, STRESS, NEWSTP 2 1 CAGRGL, COSAV, ACSAV, COSAVW, ACSAVW, SRAT, SRATAV ) NEWSTP 3 C NEWSTP FOR KYST JANUARY,1973 NEWSTP 4 C WRITTEN BY J.KRUSKAL NEWSTP 5 C KYST VERSION 2 OCTOBER, 1976 NEWSTP 6 C NEWSTP 7 C NEWSTP THIS SUBROUTINE COMPUTES THE STEP SIZE. NEWSTP 8 C NEWSTP 9 C THE MAIN PURPOSE OF THIS ROUTINE IS TO COMPUTE THE NEW NEWSTP10 C VALUE OF 'STEP'. NEWSTP11 C INCIDENTALLY, IT UPDATES 'COSAV', 'ACSAV', AND 'SRATAV'. NEWSTP12 C NEWSTP13 C UPDATE THREE AVERAGE QUANTITIES NEWSTP14 C NEWSTP15 COSAV = CAGRGL*COSAVW + COSAV*(1.0-COSAVW) NEWSTP16 ACSAV = ABS (CAGRGL)*ACSAVW + ACSAV*(1.0-ACSAVW) NEWSTP17 SRATAV = (SRAT**0.33334) * (SRATAV**0.66666) NEWSTP18 IF(ITNO) 100, 100, 200 NEWSTP19 C NEWSTP20 C GUESS INITIAL STEP SIZE NEWSTP21 C NEWSTP22 100 STEP= (25.0*STRESS) * SFGR NEWSTP23 RETURN NEWSTP24 C NEWSTP25 C FIND NEW STEP SIZE NEWSTP26 C NEWSTP27 200 ANG=4.0**COSAV NEWSTP28 TEMP1 = 1.0 + (AMIN1 (1.0,SRATAV) ) ** 5 NEWSTP29 TEMP2 = 1.0 + ( ACSAV - ABS (COSAV) ) NEWSTP30 BIAS = 1.6 / (TEMP1*TEMP2) NEWSTP31 GOODLK = SQRT (AMIN1 (1.0,SRAT) ) NEWSTP32 STEP = STEP * ANG * BIAS * GOODLK NEWSTP33 RETURN NEWSTP34 END NEWSTP35 CNRMLZ NRMLZ 1 SUBROUTINE NRMLZ(N,V) NRMLZ 2 C NRMLZ FOR KYST JANUARY,1973 NRMLZ 3 C WRITTEN FOR SGEV BY P.BUSINGER JANUARY,1972 NRMLZ 4 C KYST VERSION 2 OCTOBER, 1976 NRMLZ 5 C CALLED BY NVIT1 TO NORMALIZE A VECTOR OF LENGTH N NRMLZ 6 DIMENSION V(1) NRMLZ 7 C NRMLZ 8 REAL S NRMLZ 9 S=0.E0 NRMLZ 10 DO 10 K=1,N NRMLZ 11 10 S=AMAX1(S,ABS(V(K))) NRMLZ 12 DO 20 K=1,N NRMLZ 13 20 V(K)=V(K)/S NRMLZ 14 RETURN NRMLZ 15 END NRMLZ 16 CNVIT1 NVIT1 1 SUBROUTINE NVIT1(N,OMEGA,B,U,V) NVIT1 2 C NVIT1 FOR KYST JANUARY,1973 NVIT1 3 C WRITTEN AND TECHNIQUES ADVOCATED IN NVIT1 4 C REFERENCE IMPLEMENTED FOR SGEV NVIT1 5 C BY P.BUSINGER JANUARY,1972 NVIT1 6 C APRIL, 1973 VERSION NVIT1 7 C KYST VERSION 2 OCTOBER, 1976 NVIT1 8 DIMENSION B(1),U(1),V(1) NVIT1 9 C NVIT1 10 C G.GOLUB AND W.KAHAN, CALCULATING THE SINGULAR NVIT1 11 C VALUES AND PSEUDO-INVERSE OF A MATRIX. J.SIAM NVIT1 12 C NUMER.ANAL., SER.B, VOL.2 (1965), 205-224. NVIT1 13 C NVIT1 14 C NVIT1 IS CALLED BY SGEV, AND ASSUME THAT NVIT1 15 C GIVEN ARE THE DIAGONAL ELEMENTS NVIT1 16 C U(1) THROUGH U(N) NVIT1 17 C (ALL OF THE SAME SIGN) OF THE MATRIX NVIT1 18 C U IN THE LU-DECOMPOSITION OF A NEARLY NVIT1 19 C SINGULAR SYMMETRIC TRIDIAGONAL MATRIX NVIT1 20 C A. ALSO GIVEN ARE THE OFF-DIAGONAL NVIT1 21 C ELEMENTS NVIT1 22 C B(1) THROUGH B(N-1) NVIT1 23 C OF A AND OMEGA = LARGEST MACHINE NUMBER. NVIT1 24 C COMPUTED ARE THE COMPONENTS NVIT1 25 C V(1) THROUGH V(N) NVIT1 26 C OF THE EIGENVECTOR ASSOCIATED WITH THE NVIT1 27 C APPROXIMATE ZERO EIGENVALUE OF A. NVIT1 28 C NVIT1 29 C CALLS NRMLZ TO NORMALIZE A VECTOR NVIT1 30 REAL W,S,OMGA NVIT1 31 V(1)=1.E0 NVIT1 32 IF(N.LT.2)GOTO 40 NVIT1 33 S=1.E0 NVIT1 34 DO 1 I=1,N NVIT1 35 1 S=AMAX1(S,ABS(U(I))) NVIT1 36 OMGA=OMEGA/S NVIT1 37 DO 10 I=2,N NVIT1 38 I1=I-1 NVIT1 39 IF(ABS(B(I1)*V(I1)).GE.ABS(OMGA*U(I1)))CALL NRMLZ(I1,V) NVIT1 40 W=-(B(I1)*V(I1))/U(I1) NVIT1 41 10 V(I)=W+SIGN(V(1),W) NVIT1 42 CALL NRMLZ(N,V) NVIT1 43 V(N)=V(N)/U(N) NVIT1 44 DO 30 II=2,N NVIT1 45 I=N+1-II NVIT1 46 IF(ABS(V(I)-B(I)*V(I+1)).GE.ABS(OMGA*U(I)))CALL NRMLZ(N,V) NVIT1 47 30 V(I)=(V(I)-B(I)*V(I+1))/U(I) NVIT1 48 CALL NRMLZ(N,V) NVIT1 49 S=0.E0 NVIT1 50 DO 35 I=1,N NVIT1 51 35 S=S+V(I)**2 NVIT1 52 S=SQRT(S) NVIT1 53 DO 39 I=1,N NVIT1 54 39 V(I)=V(I)/S NVIT1 55 40 RETURN NVIT1 56 END NVIT1 57 CPLOT PLOT 1 SUBROUTINE PLOT(X,Y,XA,YA,XI,YI,NPOI,NY,ID,LONG,MX) PLOT 2 C KYST VERSION 2 OCTOBER, 1976 PLOT 3 C PLOT 4 C ROUTINE TO GENERATE A ONE PAGE PLOT OF A MATRIX Y VS SOME VECTOR X PLOT 5 C PLOT 6 C CALLING PARAMETERS PLOT 7 C Y--A MATRIX (LONG,1 OR 2) IT CONTAINS VALUES TO BE PLOTTED AGAINST PLOT 8 C X--A VECTOR (LONG); IT CONTAINS VALUES TO BE PLOTTED AGAINST Y PLOT 9 C XA-- X AXIS UPPER BOUND (MAY BE PASSED AS ZERO) PLOT 10 C YA-- Y AXIS UPPER BOUND (MAY BE PASSED AS ZERO) PLOT 11 C XI-- X AXIS LOWER BOUND (MAY BE PASSED AS ZERO) PLOT 12 C YI-- Y AXIS LOWER BOUND (MAY BE PASSED AS ZERO) PLOT 13 C NPOI-- NO. OF ELEMENTS IN X, NO. OF ROWS IN Y PLOT 14 C NY--PLOTTING CONTROL PARAMETER PLOT 15 C ID-- PRINT CONTROL PARAMETER, CONTROLS LABELING, AND AXIS GEN. PLOT 16 C LONG -- MAXIMUM NUMBER OF ELEMENTS IN X, MAX. NO. ROWS IN Y PLOT 17 C MX-- WORK WITH ELEMENTS MX THROUGH MX+NPOI-1 PLOT 18 C PLOT 19 C NOTES ABOUT CALLING PARAMETERS PLOT 20 C PLOT 21 C THE USER MAY PASS THE MIN AND MAX FOR X AND Y OR HE MAY LET THE PLOT 22 C PROGRAM FIND THEM. IF XA=XI THEN THE PROGRAM WILL FIND THE MAX AND MIPLOT 23 C FOR THE X VECTOR, IF YA=YI THE PROGRAM WILL FIND THE MAX AND MIN PLOT 24 C FOR THE Y MATRIX. IN EITHER CASE THE MAX AND MIN OF AN AXIS MAY BE PLOT 25 C ADJUSTED TO ENSURE THAT THE ASPECT RATIO IS CORRECT. PLOT 26 C PLOT 27 C IF -ID- IS POSITIVE, AXES WILL BE INCLUDED ON THE PLOT, PLOT 28 C AND THE ASPECT RATIO IS ADJUSTED. PLOT 29 C IF -ID- IS NEGATIVE, NO AXES WILL APPEAR ON THE PLOT, PLOT 30 C AND THE ASPECT RATIO IS NOT ADJUSTED. PLOT 31 C IF ID IS PLUS OR MINUS 1, THE POINTS WILL BE IDENTIFIED BY COLUMNS PLOT 32 C IF ID IS PLUS OR MINUS 2, THE POINTS WILL BE IDENTIFIED BY ROWS PLOT 33 C IF ID IS PLUS OR MINUS 3, A ONE DIMENSIONAL PLOT IS GENERATED. PLOT 34 C ONE DIMENSIONAL PLOTS ARE PLOTTED WITH NO AXES, ALL Y VALUES ARE PLOT 35 C FIXED AT (MAXY-MINY)/2, AND ALL POINTS ARE IDENTIFIED BY ROWS. PLOT 36 C PLOT 37 C WHEN NY=1 (ONE COLUMN IN Y) THEN BEFORE PLOT IS CALLED Y MUST BE PLOT 38 C SORTED INTO ASCENDING ORDER WITH X AND PTID REARRANGED ACCORDINGLY. PLOT 39 C WHEN NY=+2 (TWO COLUMNS IN Y, AND BOTH TO BE PLOTTED) THEN PLOT 40 C THE SECOND COLUMN MUST BE SORTED IN PLOT 41 C ASCENDING ORDER AND X AND THE FIRST COLUMN OF Y MUST BE REARRANGED. PLOT 42 C FURTHERMORE, WHEN THERE ARE TWO COLUMNS IN Y 'ID' MUST BE EITHER PLOT 43 C PLUS OR MINUS 1, AND THE FIRST TWO CHARACTERS OF PTID WILL BE USED PLOT 44 C TO LABEL THE COLUMNS OF Y. PLOT 45 C PLOT 46 C IF NY IS PLUS 1 -- CONFIGURATION IS PLOTTED. PLOT 47 C IF NY IS PLUS 2 -- SCATTER DIAGRAM IS PLOTTED. PLOT 48 C IF NY IS MINUS 2 -- JOINT SCATTER DIAGRAM (OMITTING FITTED PLOT 49 C FUNCTION) IS PLOTTED. PLOT 50 C PLOT 51 C PLOT 52 C NOTE THAT THE INFORMATION TO BE PLOTTED, PLOT 53 C WHICH IS CONTAINED IN X AND Y, IS DESTROYED. PLOT 54 C PLOT 55 C PLOT 56 C INTERNAL VARIABLE DEFINITIONS PLOT 57 C DELX-- X AXIS INCREMENT PLOT 58 C DELY-- Y AXIS INCREMENT PLOT 59 C FINC-- FLOATING INCREMENT FOR XLABEL PLOT 60 C IDAT-- COUNTER OVER ROWS OF X AND Y PLOT 61 C ILAB-- INTEGER LABEL FLAG PLOT 62 C ITEM(101)--ITEMS TO BE PRINTED INITIALLY BLANK PLOT 63 C LINE-- COUNTER OF CURRENT LINE PLOT 64 C NLPP-- NUMBER OF LINES PLOTTED PER PAGE PLOT 65 C PTID-- POINT IDENTIFICATION CHARACTERS PLOT 66 C SPANX-- RANGE OF X PLOT 67 C SPANY-- RANGE OF Y PLOT 68 C VALUE-- DATUM PRINTED ON Y AXIS PLOT 69 C XLABEL(11)--X-AXIS LABELS PLOT 70 C XNOW-- CURRENT VALUE OF X VECTOR NOW BEING WORKED ON PLOT 71 C YNOW-- CURRENT VALUE OF Y MATRIX NOW BEING WORKED ON PLOT 72 C PLOT 73 C FULL PAGE PLOTTER WRITTEN BY FORREST W. YOUNG, NOVEMBER, 1965 PLOT 74 C SINGLE LINE PLOTTER WRITTEN BY FORREST W. YOUNG & JAY R. LEVINSOHN, PLOT 75 C OCTOBER, 1971 PLOT 76 C ALTERED FOR KYST BY J.KRUSKAL AND J.SEERY JANUARY,1973 PLOT 77 C ALTERED FOR KYST-2 BY FORREST W. YOUNG, FEBRUARY, 1976 PLOT 78 C PLOT 79 INTEGER ITEM(101),PTID(100),MINUS,DD,STAR,BLANK,AIE,DUMX1(1) PLOT 80 REAL GR(100,6),GL(100,6),X2(5430),YMINV(2),YMAXV(2) PLOT 81 DIMENSION XLABL(11),X(LONG),XS(10),YS(10),IND(2),DUMX2(1) PLOT 82 C NOTE'Y IS DUMMY-DIMENSIONED; IT SHOULD HAVE 1 OR 2 COLUMNS PLOT 83 REAL Y(LONG,1) PLOT 84 COMMON /MDSCL1/ LREAD,LPRINT,LPUNCH,LSCRAT PLOT 85 COMMON /ACCUR/ PRECSN,XMAG,XMAXN,FLPRIN,CHPRIN PLOT 86 COMMON /PLTCHR/ PTID,ITEM PLOT 87 COMMON /KYST2/ GR,GL,X2 PLOT 88 DATA MINUS,DD,STAR,BLANK,AIE /1H-,1H.,1H*,1H ,1HI/ PLOT 89 DATA XS(1),XS(2),XS(3),XS(4),XS(5),XS(6),XS(7),XS(8),XS(9),XS(10) PLOT 90 . /0.,0.,0.,0.,0.,4.,2.,4.,1.,0./ PLOT 91 DATA YS(1),YS(2),YS(3),YS(4),YS(5),YS(6),YS(7),YS(8),YS(9),YS(10) PLOT 92 . /0.,0.,0.,0.,2.,4.,3.,4.,7.,2./ PLOT 93 C PLOT 94 C VARIABLE DEFINTIONS PLOT 95 C PLOT 96 JY=1 PLOT 97 NYSIGN=1 PLOT 98 IF (NY.LT.0) NYSIGN=-1 PLOT 99 NY=IABS(NY) PLOT 100 INDEX=0 PLOT 101 LINE=1 PLOT 102 IY0L=0 PLOT 103 IX0L=0 PLOT 104 IDAT=1 PLOT 105 NLPP=53 PLOT 106 XMAX=XA PLOT 107 YMAX=YA PLOT 108 YMIN=YI PLOT 109 XMIN=XI PLOT 110 LXFLG=0 PLOT 111 LYFLG=0 PLOT 112 IRPFLG=0 PLOT 113 NPOI2=NY*NPOI PLOT 114 IF(NYSIGN.LT.0) NPOI2=NPOI PLOT 115 NPOIRL=MX+NPOI-1 PLOT 116 IF(IABS(ID).GT.3.OR.ID.EQ.0)GO TO 9999 PLOT 117 IF(NY.LT.1.OR.NY.GT.2)GO TO 9999 PLOT 118 IF(NY.EQ.2.AND.IABS(ID).NE.1)GO TO 9999 PLOT 119 IF(XMAX.LT.XMIN)GO TO 9999 PLOT 120 IF(YMAX.LT.YMIN)GO TO 9999 PLOT 121 C PLOT 122 C DETERMINE MINIMUM AND MAXIMUM OF BOTH AXES, IF NECESSARY PLOT 123 C PLOT 124 DO 11 I=MX,NPOIRL PLOT 125 11 X2(I)=X(I) PLOT 126 IF(NY.LT.2)GO TO 118 PLOT 127 CALL SORT(Y(MX,1),NPOI,DUMX1,DUMX2,X2(MX),4,1) PLOT 128 118 IF (XMAX-XMIN)122,119,122 PLOT 129 119 XMAX=-XMAXN PLOT 130 XMIN=XMAXN PLOT 131 DO 121 I=MX,NPOIRL PLOT 132 IF(X(I)-XMAX) 125,123,123 PLOT 133 123 XMAX=X(I) PLOT 134 125 IF(X(I)-XMIN) 124,124,121 PLOT 135 124 XMIN=X(I) PLOT 136 121 CONTINUE PLOT 137 IF (IABS(ID)-3)122,305,122 PLOT 138 122 IF((NY.EQ.2).AND.(NYSIGN.LT.0)) GO TO 298 PLOT 139 DO 295 I=1,NY PLOT 140 YMAXV(I)=Y(NPOIRL,I) PLOT 141 YMINV(I)=Y(MX,I) PLOT 142 IF(NY.EQ.1)GO TO 300 PLOT 143 295 CONTINUE PLOT 144 YMAX=AMAX1(YMAXV(1),YMAXV(2)) PLOT 145 YMIN=AMIN1(YMINV(1),YMINV(2)) PLOT 146 IF(YMAX.NE.YMAXV(1)) JY=2 PLOT 147 YNOW=YMAX PLOT 148 YNOW1=AMIN1(YMAXV(1),YMAXV(2)) PLOT 149 GO TO 305 PLOT 150 298 YMAX=Y(NPOIRL,1) PLOT 151 YMIN=Y(MX,1) PLOT 152 YNOW=YMAX PLOT 153 YNOW1=-XMAXN PLOT 154 GO TO 305 PLOT 155 300 YNOW=Y(NPOIRL,1) PLOT 156 YNOW1=-XMAXN PLOT 157 IF(YMAX-YMIN) 305,302,305 PLOT 158 302 YMAX=Y(NPOIRL,1) PLOT 159 YMIN=Y(MX,1) PLOT 160 C PLOT 161 C DETERMINE RANGE AND INCREMENT OF X-AXIS PLOT 162 C PLOT 163 305 P=ALOG10(XMAX-XMIN+.0001) PLOT 164 IP=IFIX(P) PLOT 165 IF(P.LT.0.) IP=IP-1 PLOT 166 SPAN=(XMAX-XMIN)/(10.**IP) PLOT 167 IF(SPAN.LT.5.0) GO TO 310 PLOT 168 DEL=10.**IP PLOT 169 GO TO 330 PLOT 170 310 IF(SPAN.LT.2.0) GO TO 320 PLOT 171 DEL=(10.**IP)/2.0 PLOT 172 GO TO 330 PLOT 173 320 DEL=(10.**IP)/5.0 PLOT 174 330 DELAX=DEL PLOT 175 IF(XMAX.LT.0.) GO TO 336 PLOT 176 XMAX=FLOAT(IFIX(XMAX/DEL+.9999))*DEL PLOT 177 IF(XMIN) 334,332,332 PLOT 178 332 XMIN=FLOAT(IFIX(XMIN/DEL))*DEL PLOT 179 GO TO 338 PLOT 180 334 XMIN=FLOAT(IFIX(XMIN/DEL-.9999))*DEL PLOT 181 GO TO 338 PLOT 182 336 XMAX=FLOAT(IFIX(XMAX/DEL))*DEL PLOT 183 XMIN=FLOAT(IFIX(XMIN/DEL+.9999))*DEL PLOT 184 338 SPANX=XMAX-XMIN PLOT 185 NXUNIT=SPANX/DEL+.00001 PLOT 186 IF(NXUNIT.LE.10) GO TO 3338 PLOT 187 LXFLG=LXFLG+1 PLOT 188 IF(LXFLG-1) 9999,305,9999 PLOT 189 C PLOT 190 C DETERMINE RANGE AND INCREMENT OF Y-AXIS PLOT 191 C PLOT 192 3338 IF(IABS(ID).EQ.3)GO TO 3690 PLOT 193 IF(IRPFLG.EQ.1)GO TO 3668 PLOT 194 P=ALOG10(YMAX-YMIN+.0001) PLOT 195 IP=IFIX(P) PLOT 196 IF(P.LT.0.) IP=IP-1 PLOT 197 SPAN=(YMAX-YMIN)/(10.**IP) PLOT 198 IF(SPAN.LT.5.0) GO TO 340 PLOT 199 DEL=10.**IP PLOT 200 GO TO 360 PLOT 201 340 IF(SPAN.LT.2.0) GO TO 350 PLOT 202 DEL=(10.**IP)/2.0 PLOT 203 GO TO 360 PLOT 204 350 DEL=(10.**IP)/5.0 PLOT 205 360 DELAY=DEL PLOT 206 IF(YMAX.LT.0.) GO TO 366 PLOT 207 YMAX=FLOAT(IFIX(YMAX/DEL+.9999))*DEL PLOT 208 IF(YMIN) 364,362,362 PLOT 209 362 YMIN=FLOAT(IFIX(YMIN/DEL))*DEL PLOT 210 GO TO 368 PLOT 211 364 YMIN=FLOAT(IFIX(YMIN/DEL-.9999))*DEL PLOT 212 GO TO 368 PLOT 213 366 YMAX=FLOAT(IFIX(YMAX/DEL))*DEL PLOT 214 YMIN=FLOAT(IFIX(YMIN/DEL+.9999))*DEL PLOT 215 368 SPANY = YMAX-YMIN PLOT 216 NYUNIT=SPANY/DEL+.00001 PLOT 217 IF(NYUNIT.LE.10) GO TO 3668 PLOT 218 LYFLG=LYFLG+1 PLOT 219 IF(LYFLG-1)9999,3338,9999 PLOT 220 C PLOT 221 C ADJUST THE ASPECT RATIO PLOT 222 C PLOT 223 3668 IF(ID.LT.0)GO TO 3690 PLOT 224 DASPRT=(FLOAT(NLPP-1)/FLPRIN)/(100.0/CHPRIN) PLOT 225 ASPRT=SPANY/SPANX PLOT 226 IRPFLG=0 PLOT 227 IF(ASPRT-DASPRT)3670,3690,3680 PLOT 228 3670 ATEMP=NYUNIT PLOT 229 NYUNIT=ATEMP*DASPRT/ASPRT+0.5 PLOT 230 ATEMP=NYUNIT PLOT 231 ATEMP= ((ATEMP*DELAY)-(YMAX-YMIN))*0.5 PLOT 232 YMAX=YMAX+ATEMP PLOT 233 YMIN=YMIN-ATEMP PLOT 234 SPANY=YMAX-YMIN PLOT 235 IF(NYUNIT.LE.10)GO TO 3690 PLOT 236 LYFLG=0 PLOT 237 GO TO 3338 PLOT 238 3680 ATEMP=NXUNIT PLOT 239 NXUNIT=ATEMP*ASPRT/DASPRT+0.5 PLOT 240 ATEMP=NXUNIT PLOT 241 ATEMP= ((ATEMP*DELAX)-(XMAX-XMIN))*0.5 PLOT 242 XMAX=XMAX+ATEMP PLOT 243 XMIN=XMIN-ATEMP PLOT 244 SPANX=XMAX-XMIN PLOT 245 IF(NXUNIT.LE.10)GO TO 3690 PLOT 246 IRPFLG=1 PLOT 247 LXFLG=0 PLOT 248 GO TO 305 PLOT 249 C PLOT 250 C DETERMINE ORIGIN OF EACH AXIS PLOT 251 C PLOT 252 3690 NXP1=NXUNIT+1 PLOT 253 FINC=DELAX PLOT 254 DELX=SPANX/(100.-XS(NXUNIT)) PLOT 255 IXS=IFIX(XS(NXUNIT)/2.+.1) PLOT 256 XSF=IXS PLOT 257 IF(ID.LT.0) GO TO 339 PLOT 258 IF((XMIN.GE.0.).OR.(XMAX.LE.0.) ) GO TO 339 PLOT 259 IX0L=-XMIN/DELX+1.5 PLOT 260 IX0L=IX0L+IXS PLOT 261 339 IF(IABS(ID).EQ.3)GO TO 345 PLOT 262 DELY=SPANY/(FLOAT(NLPP-1)-YS(NYUNIT)) PLOT 263 DELTY=DELAY PLOT 264 VALUE=YMAX PLOT 265 YININC=(NLPP-1)/NYUNIT PLOT 266 LININC=IFIX(YININC) PLOT 267 IYLAB=LININC PLOT 268 IF(ID.LT.0) GO TO 369 PLOT 269 IF((YMAX.LE.0.).OR.(YMIN.GE.0.)) GO TO 369 PLOT 270 IY0L=YMAX/DELY+1.5 PLOT 271 GO TO 369 PLOT 272 345 VALUE=0.0 PLOT 273 DELTY=0.0 PLOT 274 LININC=27 PLOT 275 IYLAB=1 PLOT 276 NYUNIT=1 PLOT 277 C PLOT 278 C LABEL PLOT AT TOP PLOT 279 C PLOT 280 369 XLABL(1)=XMIN PLOT 281 DO 120 I=1,NXUNIT PLOT 282 120 XLABL(I+1)=XLABL(I)+FINC PLOT 283 GO TO (9999,9999,9999,365,370,375,380,385,390,395), NXUNIT PLOT 284 365 WRITE(LPRINT,3304) (XLABL(I),I=1,NXP1) PLOT 285 WRITE(LPRINT,3314) PLOT 286 3304 FORMAT(1H0,F20.4,4F25.4) PLOT 287 3314 FORMAT(1H ,14X,103H*.************************.********************PLOT 288 1****.************************.************************.*) PLOT 289 GO TO 400 PLOT 290 370 WRITE(LPRINT,3305) (XLABL(I),I=1,NXP1) PLOT 291 WRITE(LPRINT,3315) PLOT 292 3305 FORMAT(1H0,6F20.4) PLOT 293 3315 FORMAT(1H ,14X,103H*.*******************.*******************.*****PLOT 294 2**************.*******************.*******************.*) PLOT 295 GO TO 400 PLOT 296 375 WRITE(LPRINT,3306) (XLABL(I),I=1,NXP1) PLOT 297 WRITE(LPRINT,3316) PLOT 298 3306 FORMAT(1H0,6X,7F16.4) PLOT 299 3316 FORMAT(1H ,14X,103H***.***************.***************.***********PLOT 300 1****.***************.***************.***************.***) PLOT 301 GO TO 400 PLOT 302 380 WRITE(LPRINT,3307) (XLABL(I),I=1,NXP1) PLOT 303 WRITE(LPRINT,3317) PLOT 304 3307 FORMAT(1H0,7X,8F14.4) PLOT 305 3317 FORMAT(1H ,14X,103H**.*************.*************.*************.**PLOT 306 1***********.*************.*************.*************.**) PLOT 307 GO TO 400 PLOT 308 385 WRITE(LPRINT,3308) (XLABL(I),I=1,NXP1) PLOT 309 WRITE(LPRINT,3318) PLOT 310 3308 FORMAT(1H0,10X,9F12.4) PLOT 311 3318 FORMAT(1H ,14X,103H***.***********.***********.***********.*******PLOT 312 1****.***********.***********.***********.***********.***) PLOT 313 GO TO 400 PLOT 314 390 WRITE(LPRINT,3309) (XLABL(I),I=1,NXP1) PLOT 315 WRITE(LPRINT,3319) PLOT 316 3309 FORMAT(1H0,8X,10F11.3) PLOT 317 3319 FORMAT(1H ,14X,103H*.**********.**********.**********.**********.*PLOT 318 1*********.**********.**********.**********.**********.**) PLOT 319 GO TO 400 PLOT 320 395 WRITE(LPRINT,3310) (XLABL(I),I=1,NXP1) PLOT 321 WRITE(LPRINT,3320) PLOT 322 3310 FORMAT(1H0,9X,11F10.3) PLOT 323 3320 FORMAT(1H ,14X,103H*.*********.*********.*********.*********.*****PLOT 324 1****.*********.*********.*********.*********.*********.*) PLOT 325 400 CONTINUE PLOT 326 C PLOT 327 C PREPARE A SINGLE PLOTTING LINE PLOT 328 C PLOT 329 IND(1)=NPOIRL PLOT 330 IND(2)=IND(1) PLOT 331 4000 ITEMP=YS(NYUNIT)+.0001 PLOT 332 C IF BELOW BOTTOM OF PLOT, GO TO 2000 PLOT 333 IF(LINE-NLPP+ITEMP)130,130,2000 PLOT 334 130 IF (IDAT-1)700,701,700 PLOT 335 700 IF(IABS(ID)-3)706,701,706 PLOT 336 706 IF (LFLAG)702,702,701 PLOT 337 702 IF(IDAT-NPOI2) 703,703,701 PLOT 338 703 IND(JY)=IND(JY)-1 PLOT 339 INDEX=IND(JY) PLOT 340 IF(INDEX.LT.MX)GO TO 7704 PLOT 341 YNOW=Y(INDEX,JY) PLOT 342 IF(YNOW-YNOW1) 7702,701,701 PLOT 343 7704 IF(NYSIGN.LT.0) GO TO 701 PLOT 344 JY=NY-JY+1 PLOT 345 YNOW=YNOW1 PLOT 346 YNOW1=-XMAXN PLOT 347 GO TO 701 PLOT 348 7702 IF(NYSIGN.LT.0) GO TO 701 PLOT 349 JY=NY-JY+1 PLOT 350 TEMP=YNOW PLOT 351 YNOW=YNOW1 PLOT 352 YNOW1=TEMP PLOT 353 C IF ALL POINTS HAVE BEEN PLOTTED, GO TO 1000; ELSE, GO TO 806 PLOT 354 701 IF(IDAT-NPOI2) 806,806,1000 PLOT 355 806 IF (IABS(ID)-3)805,808,805 PLOT 356 805 IF (YMAX-YNOW)811,801,801 PLOT 357 801 IF (YNOW-YMIN)811,802,802 PLOT 358 802 I=(YMAX-YNOW)/DELY+1.50001 PLOT 359 GOTO 807 PLOT 360 808 IF (LINE-27)807,809,807 PLOT 361 809 DO 10 II=MX,NPOIRL PLOT 362 I=II PLOT 363 10 CALL XITEM(I,XMAX,XMIN,ID,DELX,X,JY,XSF) PLOT 364 GOTO 606 PLOT 365 807 IF (I-LINE)820,820,810 PLOT 366 810 LFLAG=1 PLOT 367 GOTO 1000 PLOT 368 811 LFLAG=0 PLOT 369 GOTO 606 PLOT 370 820 LFLAG=0 PLOT 371 INDEX=IND(JY) PLOT 372 IF((JY.EQ.2).AND.(NYSIGN.GT.0)) PLOT 373 *CALL XITEM(INDEX,XMAX,XMIN,ID,DELX,X,JY,XSF) PLOT 374 IF(JY.EQ.1) PLOT 375 *CALL XITEM(INDEX,XMAX,XMIN,ID,DELX,X2,JY,XSF) PLOT 376 606 IDAT=IDAT+1 PLOT 377 GOTO 4000 PLOT 378 C PLOT 379 C INSERT AXES INTO PLOTTING LINE IF DESIRED PLOT 380 C PLOT 381 1000 IF (ID)133,901,901 PLOT 382 901 IF(IABS(ID)-3) 902,133,902 PLOT 383 902 IF(IX0L) 907,907,903 PLOT 384 903 IF(ITEM(IX0L).NE.BLANK) GO TO 907 PLOT 385 906 ITEM(IX0L)=AIE PLOT 386 907 IF(IY0L) 133,133,909 PLOT 387 909 IF(LINE-IY0L) 133,904,133 PLOT 388 904 DO 905 I=1,101 PLOT 389 IF(ITEM(I).NE.BLANK) GO TO 905 PLOT 390 908 ITEM(I)=MINUS PLOT 391 905 CONTINUE PLOT 392 C PLOT 393 C PRINT ONE LINE OF PLOT PLOT 394 C PLOT 395 133 IF(IYLAB.NE.LININC) GO TO 915 PLOT 396 WRITE(LPRINT,3300) VALUE,DD,(ITEM(J),J=1,101),DD,VALUE PLOT 397 3300 FORMAT(1H ,F12.4,2X,103A1,F13.4) PLOT 398 VALUE=VALUE-DELTY PLOT 399 IYLAB=0 PLOT 400 GO TO 920 PLOT 401 915 WRITE(LPRINT,3410) STAR, (ITEM(J),J=1,101),STAR PLOT 402 3410 FORMAT(1H ,14X,103A1) PLOT 403 920 IYLAB=IYLAB+1 PLOT 404 DO 3350 I=1,101 PLOT 405 3350 ITEM(I)=BLANK PLOT 406 LINE=LINE+1 PLOT 407 GOTO 4000 PLOT 408 C PLOT 409 C LABEL PLOT ALONG THE BOTTOM PLOT 410 C PLOT 411 2000 GO TO (9999,9999,9999,965,970,975,980,985,990,995), NXUNIT PLOT 412 965 WRITE(LPRINT,3314) PLOT 413 WRITE(LPRINT,3304) (XLABL(I),I=1,NXP1) PLOT 414 GO TO 998 PLOT 415 970 WRITE(LPRINT,3315) PLOT 416 WRITE(LPRINT,3305) (XLABL(I),I=1,NXP1) PLOT 417 GO TO 998 PLOT 418 975 WRITE(LPRINT,3316) PLOT 419 WRITE(LPRINT,3306) (XLABL(I),I=1,NXP1) PLOT 420 GO TO 998 PLOT 421 980 WRITE(LPRINT,3317) PLOT 422 WRITE(LPRINT,3307) (XLABL(I),I=1,NXP1) PLOT 423 GO TO 998 PLOT 424 985 WRITE(LPRINT,3318) PLOT 425 WRITE(LPRINT,3308) (XLABL(I),I=1,NXP1) PLOT 426 GO TO 998 PLOT 427 990 WRITE(LPRINT,3319) PLOT 428 WRITE(LPRINT,3309) (XLABL(I),I=1,NXP1) PLOT 429 GO TO 998 PLOT 430 995 WRITE(LPRINT,3320) PLOT 431 WRITE(LPRINT,3310) (XLABL(I),I=1,NXP1) PLOT 432 998 CONTINUE PLOT 433 C PLOT 434 C RESTORE NY PLOT 435 NY=NYSIGN*NY PLOT 436 RETURN PLOT 437 C PLOT 438 9999 WRITE(LPRINT,9995) PLOT 439 9995 FORMAT(22H0ERROR EXIT FROM PLOT.) PLOT 440 STOP PLOT 441 C PLOT 442 END PLOT 443 CREGR REGR 1 FUNCTION REGR(DA,I) REGR 2 C REGR FOR KYST JANUARY,1973 REGR 3 C WRITTEN BY J.KRUSKAL REGR 4 C KYST VERSION 2 OCTOBER, 1976 REGR 5 CREGR CALCULATES VALUES OF FUNCTIONS REGRESSED OVER REGR 6 C REGR 7 J=MIN0(I,4) REGR 8 GO TO (10,20,30,40), J REGR 9 C REGR 10 10 REGR=1.0 REGR 11 RETURN REGR 12 C REGR 13 20 REGR=DA REGR 14 RETURN REGR 15 C REGR 16 30 REGR=DA*DA REGR 17 RETURN REGR 18 C REGR 19 40 REGR=DA**(I-1) REGR 20 RETURN REGR 21 C REGR 22 END REGR 23 CRM1POW RM1POW 1 FUNCTION RM1POW(YY) RM1POW 2 C RM1POW FOR KYST JANUARY,1973 RM1POW 3 C WRITTEN BY J.KRUSKAL RM1POW 4 C MODIFIED FOR KYST BY J.SEERY JANUARY,1973 RM1POW 5 C KYST VERSION 2 OCTOBER, 1976 RM1POW 6 C RM1POW--CALCULATES (R-1)-ST POWER RM1POW 7 C RM1POW 8 INTEGER RTYPE RM1POW 9 COMMON /METRIC/ RM1,RECR,R,RTYPE RM1POW10 Y=YY RM1POW11 GO TO (210,220,230), RTYPE RM1POW12 210 RM1POW=0.0 RM1POW13 IF(Y.NE.0.0) RM1POW=SIGN(1.0,Y) RM1POW14 RETURN RM1POW15 220 RM1POW=Y RM1POW16 RETURN RM1POW17 230 RM1POW=SIGN(ABS(Y)**RM1,Y) RM1POW18 RETURN RM1POW19 END RM1POW20 CRPOWER RPOWER 1 FUNCTION RPOWER(XX) RPOWER 2 C RPOWER FOR KYST JANUARY,1973 RPOWER 3 C WRITTEN BY J.KRUSKAL RPOWER 4 C MODIFIED FOR KYST BY J.SEERY JANUARY,1973 RPOWER 5 C KYST VERSION 2 OCTOBER, 1976 RPOWER 6 C RPOWER--CALCULATES R-TH POWER RPOWER 7 C RPOWER 8 INTEGER RTYPE RPOWER 9 COMMON /METRIC/ RM1,RECR,R,RTYPE RPOWER10 X=XX RPOWER11 GO TO (110,120,130), RTYPE RPOWER12 110 RPOWER=ABS(X) RPOWER13 RETURN RPOWER14 120 RPOWER=X*X RPOWER15 RETURN RPOWER16 130 RPOWER=ABS(X)**R RPOWER17 RETURN RPOWER18 END RPOWER19 CRROOT RROOT 1 FUNCTION RROOT(ZZ) RROOT 2 C RROOT FOR KYST JANUARY,1973 RROOT 3 C WRITTEN BY J.KRUSKAL RROOT 4 C MODIFIED FOR KYST BY J.SEERY JANUARY,1973 RROOT 5 C KYST VERSION 2 OCTOBER, 1976 RROOT 6 C RROOT--CALCULATES R-TH ROOT RROOT 7 C RROOT 8 INTEGER RTYPE RROOT 9 COMMON /METRIC/ RM1,RECR,R,RTYPE RROOT 10 Z=ZZ RROOT 11 GO TO (310,320,330), RTYPE RROOT 12 310 RROOT=Z RROOT 13 RETURN RROOT 14 320 RROOT=SQRT(Z) RROOT 15 RETURN RROOT 16 330 RROOT=Z**RECR RROOT 17 RETURN RROOT 18 END RROOT 19 CRUNIFV RUNIFV 1 FUNCTION RUNIFV(A) RUNIFV 2 C RUNIFV FOR KYST JANUARY,1973 RUNIFV 3 C WRITTEN BY J.KRUSKAL RUNIFV 4 C KYST VERSION 2 OCTOBER, 1976 RUNIFV 5 C RUNIFV 6 C THIS IS A SIMPLE ROUTINE FOR GENERATING RANDOM NUMBERS. RUNIFV 7 C THEY ARE UNIFORMLY DISTRIBUTED BETWEEN 0 AND 1. RUNIFV 8 C IT IS NOT FAST, NOR DOES IT PRODUCE VERY HIGH QUALITY NUMBERS. RUNIFV 9 C ITS MAIN VIRTUE IS THAT IT IS IN FORTRAN AND WILL WORK ON RUNIFV10 C ALMOST ANY MACHINE ON WHICH FORTRAN IS AVAILABLE. RUNIFV11 C RUNIFV12 DATA MODULO,FLMOD,K/8192,8192.0,1/ RUNIFV13 C MODULO = 2**13 RUNIFV14 C RUNIFV15 DO 10 I=1,15 RUNIFV16 10 K = MOD(5*K, MODULO) RUNIFV17 Z = FLOAT(K)/FLMOD RUNIFV18 RUNIFV = Z RUNIFV19 RETURN RUNIFV20 C RUNIFV21 END RUNIFV22 CSBK SBK 1 SUBROUTINE SBK(N,K,A,E,Z) SBK 2 C SBK FOR KYST JANUARY,1973 SBK 3 C TRANSLATED FROM ALGOL AND ADAPTED SBK 4 C FOR SGEV BY P.BUSINGER JANUARY,1972 SBK 5 C KYST VERSION 2 OCTOBER, 1976 SBK 6 REAL A(1),E(N),Z(100,K) SBK 7 C SBK 8 C R.S.MARTIN, C.REINSCH, AND J.H.WILKINSON, SBK 9 C HOUSEHOLDER-S TRIDIAGONALIZATION OF A SYM- SBK 10 C METRIC MATRIX. NUMER.MATH.11,181-195(1968). SBK 11 C SBK 12 C SBK IS CALLED BY SGEV, AND SBK 13 C USES THE INFORMATION PROVIDED BY THE SUB- SBK 14 C ROUTINE STO3 IN A AND E TO BACKTRANSFORM SBK 15 C THE EIGENVECTOR MATRIX Z. SBK 16 C SBK 17 REAL H,S SBK 18 IF(N.LT.2)GOTO 50 SBK 19 DO 40 I=2,N SBK 20 IF(E(I).EQ.0.E0)GOTO 40 SBK 21 L=I-1 SBK 22 MM=I*(I-1)/2+I-1 SBK 23 H=E(I)*A(MM) SBK 24 DO 30 J=1,K SBK 25 S=0.E0 SBK 26 DO 10 M=1,L SBK 27 MM=I*(I-1)/2+M SBK 28 10 S=S+A(MM)*Z(M,J) SBK 29 S=S/H SBK 30 DO 20 M=1,L SBK 31 MM=I*(I-1)/2+M SBK 32 20 Z(M,J)=Z(M,J)+S*A(MM) SBK 33 30 CONTINUE SBK 34 40 CONTINUE SBK 35 50 RETURN SBK 36 END SBK 37 CSERCH SERCH 1 SUBROUTINE SERCH(VEC,LEN,INDEX,KODE,JY,MX) SERCH 2 C SUBROUTINE SERCH IS NOT USED BY KYST-2A. SERCH 3 C THESE 6 CARDS CAN BE DELETED (ALTHOUGH NOT NECESSARY). SERCH 4 RETURN SERCH 5 END SERCH 6 CSGEV SGEV 1 SUBROUTINE SGEV(N,K,Z) SGEV 2 C SGEV FOR KYST JANUARY,1973 SGEV 3 C PROGRAMMED BY P.BUSINGER JANUARY,1972 SGEV 4 C KYST VERSION 2 OCTOBER, 1976 SGEV 5 C CALLS BSEC1,NVIT1,STO3,DFLT,SBK SGEV 6 DIMENSION B(100),C(100),D(100),U(100),STORE(194),E(6), SGEV 7 . XBEST(100,6),A(5430),Z(100,6) SGEV 8 COMMON /KYST2/ B,C,D,U,STORE,E,XBEST,A SGEV 9 COMMON /ACCUR/ EPS,ETA,OMEGA,FLPRIN,CHPRIN SGEV 10 C SGEV 11 C FINDS THE K ALGEBRAICALLY GREATEST EIGENVALUES SGEV 12 C E(1) .GE. E(2) .GE. ... .GE. E(K) SGEV 13 C OF THE N BY N REAL SYMMETRIC MATRIX WHOSE LOWER SGEV 14 C TRIANGULAR PART IS GIVEN ROW-COMPACTED IN A. SGEV 15 C ALSO PRODUCED IS SGEV 16 C AN ASSOCIATED ORTHOGONAL SET Z OF EIGENVECTORS. SGEV 17 C SGEV 18 C - CAUTION - SGEV 19 C NO A PRIORI ERROR BOUNDS ARE AVAILABLE AND NO SGEV 20 C A POSTERIORI CHECKS ARE PERFORMED ON THE RESULTS SGEV 21 C COMPUTED BY THIS SUBROUTINE. ONE CHECK WHICH A SGEV 22 C USER CAN PERFORM INEXPENSIVELY CONSISTS IN COMPU- SGEV 23 C TING THE RESIDUAL MATRICES AZ-ZE AND (ZT)Z-I SGEV 24 C WHOSE ELEMENTS SHOULD BE NEGLIGIBLE. SGEV 25 C SGEV 26 REAL S,SN,CS,T,Z1,Z2 SGEV 27 IF(K.LT.1)GOTO 99 SGEV 28 IF(N.LT.1)GOTO 99 SGEV 29 E(1)=A(1) SGEV 30 Z(1,1)=1.E0 SGEV 31 IF(N.EQ.1)GOTO 99 SGEV 32 CALL STO3(N,ETA/EPS,A,D,B) SGEV 33 DO 10 I=2,N SGEV 34 10 C(I-1)=B(I) SGEV 35 IB=5.E0-ALOG(EPS)/ALOG(2.E0) SGEV 36 S=1.E0 SGEV 37 DO 30 J=1,K SGEV 38 M=N+1-J SGEV 39 CALL BSEC1(IB,M,ETA,D,C,U,E(J),S) SGEV 40 DO 20 I=1,N SGEV 41 20 Z(I,J)=0.E0 SGEV 42 CALL NVIT1(M,OMEGA,C,U,Z(1,J)) SGEV 43 IF(J.NE.K)CALL DFLT(M,D,C,E(J),Z(1,J)) SGEV 44 30 E(J)=S*E(J) SGEV 45 IF(K.EQ.1)GOTO 80 SGEV 46 DO 70 JJ=2,K SGEV 47 J=K+1-JJ SGEV 48 NJ=N-K+JJ SGEV 49 Z(NJ,J)=1.E0 SGEV 50 NJ=N-J SGEV 51 DO 60 II=1,NJ SGEV 52 I=NJ+1-II SGEV 53 T=Z(I,J) SGEV 54 IF(T.EQ.0.E0)GOTO 60 SGEV 55 Z(I,J)=0.E0 SGEV 56 SN=(T+T)/(1.E0+T*T) SGEV 57 CS=(1.E0-T)*(1.E0+T)/(1.E0+T*T) SGEV 58 DO 50 M=J,K SGEV 59 Z1=Z(I,M) SGEV 60 Z2=Z(I+1,M) SGEV 61 Z(I,M)=CS*Z1-SN*Z2 SGEV 62 50 Z(I+1,M)=CS*Z2+SN*Z1 SGEV 63 60 CONTINUE SGEV 64 70 CONTINUE SGEV 65 80 CALL SBK(N,K,A,B,Z) SGEV 66 99 RETURN SGEV 67 END SGEV 68 CSORT SORT 1 SUBROUTINE SORT (A, N, B, C, D, K, SWITCH ) SORT 2 C SORT FOR KYST JANUARY,1973 SORT 3 C WRITTEN BY J.KRUSKAL SORT 4 C KYST VERSION 2 OCTOBER, 1976 SORT 5 C SORT 6 C THIS ROUTINE SORTS INPUT ARRAY 'A' AND REARRANGES, OPTIONALLY, SORT 7 C ARRAYS 'B', 'C', AND 'D', IN ORDER CORRESPONDING TO 'A'. SORT 8 C N = NUMBER OF ITEMS IN 'A' (AND 'B', 'C', 'D', IF USED) SORT 9 C K = 0--SORT 'A' ONLY, 1--REARRANGE 'B', 2--REARRANGE 'B' AND 'C', SORT 10 C 3--REARRANGE 'B', 'C', AND 'D' , 4--REARRANGE 'D' ONLY. SORT 11 C IF 'SWITCH' IS POSITIVE, SORT WILL BE IN ASCENDING ORDER, SORT 12 C IF ZERO OR NEGATIVE, IN DESCENDING ORDER. SORT 13 C ALGORITHM FROM CACM, JULY 1959, PAGE 30 BY D. L. SHELL SORT 14 C SORT 15 DIMENSION A(1),C(1),D(1) SORT 16 INTEGER SWITCH , B(1) SORT 17 105 KP1=K+1 SORT 18 IF(N.LE.1) GO TO 999 SORT 19 M = 1 SORT 20 106 M = M + M SORT 21 IF( M .LE. N ) GO TO 106 SORT 22 M = M - 1 SORT 23 994 M = M/2 SORT 24 IF(M.EQ.0) GO TO 999 SORT 25 KK = N-M SORT 26 J = 1 SORT 27 992 I = J SORT 28 996 IM = I + M SORT 29 IF(SWITCH) 810,810,800 SORT 30 800 IF(A(I).GT.A(IM)) GO TO 110 SORT 31 GO TO 995 SORT 32 810 IF(A(I).LT.A(IM)) GO TO 110 SORT 33 995 J = J+1 SORT 34 IF(J.GT.KK) GO TO 994 SORT 35 GO TO 992 SORT 36 110 TEMP=A(I) SORT 37 A(I) = A(IM) SORT 38 A(IM) = TEMP SORT 39 GO TO ( 140, 130, 120, 115, 112), KP1 SORT 40 112 TEMP=D(I) SORT 41 D(I)=D(IM) SORT 42 D(IM)=TEMP SORT 43 GO TO 140 SORT 44 115 TEMP = D(I) SORT 45 D(I) = D(IM) SORT 46 D(IM) = TEMP SORT 47 120 TEMP=C(I) SORT 48 C(I) = C(IM) SORT 49 C(IM) = TEMP SORT 50 130 ITEMP=B(I) SORT 51 B(I) = B(IM) SORT 52 B(IM) =ITEMP SORT 53 140 I = I-M SORT 54 IF(I.LT.1) GO TO 995 SORT 55 GO TO 996 SORT 56 999 RETURN SORT 57 END SORT 58 CSTO3 STO3 1 SUBROUTINE STO3(N,TOL,A,D,E) STO3 2 C STO3 FOR KYST JANUARY,1973 STO3 3 C TRANSLATED FROM ALGOL AND ADAPTED STO3 4 C FOR SGEV BY P.BUSINGER JANUARY,1972 STO3 5 C KYST VERSION 2 OCTOBER, 1976 STO3 6 REAL TOL,A(1),D(N),E(N) STO3 7 C STO3 8 C R.S.MARTIN, C.REINSCH, AND J.H.WILKINSON, STO3 9 C HOUSEHOLDER-S TRIDIAGONALIZATION OF A SYM- STO3 10 C METRIC MATRIX. NUMER.MATH.11,181-195(1968). STO3 11 C STO3 12 C STO3 IS CALLED BY SGEV, AND STO3 13 C USES HOUSEHOLDER TRANSFORMATIONS TO REDUCE STO3 14 C THE REAL SYMMETRIC N BY N MATRIX WHOSE STO3 15 C LOWER TRIANGULAR PART IS GIVEN ROW-COMPACTED STO3 16 C IN A TO TRIDIAGONAL FORM. THE DIAGONAL ELEMENTS STO3 17 C OF THE REDUCED MATRIX ARE STO3 18 C D(1) THROUGH D(N) STO3 19 C AND THE OFF-DIAGONAL ELEMENTS ARE STO3 20 C E(2) THROUGH E(N). STO3 21 C TOGETHER WITH THE E-S, THE INFORMATION STO3 22 C LEFT IN A ALLOWS THE SUBROUTINE SBAK TO STO3 23 C APPLY THE REDUCING TRANSFORMATION TO THE STO3 24 C EIGENVECTORS OF THE TRIDIAGONAL MATRIX. STO3 25 C TOL MUST BE GIVEN AS (LEAST POSITIVE MA- STO3 26 C CHINE NUMBER)/(RELATIVE MACHINE PRECISION). STO3 27 C STO3 28 REAL F,G,H STO3 29 DO 10 I=1,N STO3 30 M=I*(I-1)/2+I STO3 31 10 D(I)=A(M) STO3 32 DO 80 II=1,N STO3 33 I=N+1-II STO3 34 L=I-1 STO3 35 H=0.E0 STO3 36 IF(L.LT.1)GOTO 21 STO3 37 DO 20 K=1,L STO3 38 M=I*(I-1)/2+K STO3 39 20 H=H+A(M)**2 STO3 40 21 IF(H.GT.TOL)GOTO 30 STO3 41 E(I)=0.E0 STO3 42 GOTO 72 STO3 43 30 M=I*(I-1)/2+I-1 STO3 44 F=A(M) STO3 45 G=SQRT(H) STO3 46 IF(F.GE.0.E0)G=-G STO3 47 E(I)=G STO3 48 H=H-F*G STO3 49 A(M)=F-G STO3 50 F=0.E0 STO3 51 IF(L.LT.1)GOTO 61 STO3 52 DO 60 J=1,L STO3 53 G=0.E0 STO3 54 IF(J.LT.1)GOTO 41 STO3 55 DO 40 K=1,J STO3 56 M1=J*(J-1)/2+K STO3 57 M2=I*(I-1)/2+K STO3 58 40 G=G+A(M1)*A(M2) STO3 59 41 J1=J+1 STO3 60 IF(L.LT.J1)GOTO 51 STO3 61 DO 50 K=J1,L STO3 62 M1=K*(K-1)/2+J STO3 63 M2=I*(I-1)/2+K STO3 64 50 G=G+A(M1)*A(M2) STO3 65 51 G=G/H STO3 66 E(J)=G STO3 67 M=I*(I-1)/2+J STO3 68 60 F=F+G*A(M) STO3 69 61 H=F/(H+H) STO3 70 IF(L.LT.1)GOTO 72 STO3 71 DO 71 J=1,L STO3 72 M=I*(I-1)/2+J STO3 73 F=A(M) STO3 74 G=E(J)-H*F STO3 75 E(J)=G STO3 76 IF(J.LT.1)GOTO 71 STO3 77 DO 70 K=1,J STO3 78 M1=J*(J-1)/2+K STO3 79 M2=I*(I-1)/2+K STO3 80 70 A(M1)=A(M1)-F*E(K)-G*A(M2) STO3 81 71 CONTINUE STO3 82 72 H=D(I) STO3 83 M=I*(I-1)/2+I STO3 84 D(I)=A(M) STO3 85 80 A(M)=H STO3 86 RETURN STO3 87 END STO3 88 CWTRAN WTRAN 1 FUNCTION WTRAN(XX) WTRAN 2 C WTRAN FOR KYST JANUARY,1973 WTRAN 3 C WRITTEN FOR KYST BY J.SEERY JANUARY,1973 WTRAN 4 C KYST VERSION 2 OCTOBER, 1976 WTRAN 5 C WTRAN--WEIGHTS TRANSFORMATION. WTRAN 6 C WTRAN 7 DIMENSION DUMMY(28),DUM(6) WTRAN 8 COMMON /MDSCL2/ DUMMY,DCON1,DCON2,DCON3,DCON4,DCON5,A,B,P,S,T,DUM WTRAN 9 WTRAN=S+T*((A+B*XX)**P) WTRAN 10 RETURN WTRAN 11 END WTRAN 12 CXITEM XITEM 1 SUBROUTINE XITEM(INDEX,XMAX,XMIN,ID,DELX,X,JX,XSF) XITEM 2 C XITEM FOR KYST JANUARY,1973 XITEM 3 C JAY R LEVINSOHN OCTOBER 31,1971 XITEM 4 C ADAPTED FOR KYST BY J.SEERY JANUARY,1973 XITEM 5 C KYST VERSION 2 OCTOBER, 1976 XITEM 6 C XITEM 7 C ITS FUNCTION IS TO PLACE ELEMENTS OF THE X VECTOR IN THE XITEM 8 C PRINT LINE. XITEM 9 C XITEM 10 C XITEM 11 INTEGER ITEM(101),PTID(100) XITEM 12 REAL X(1) XITEM 13 COMMON /PLTCHR/ PTID,ITEM XITEM 14 C XITEM 15 XNOW=X(INDEX) XITEM 16 IF (XNOW-XMAX)610,610,605 XITEM 17 610 IF (XNOW-XMIN)605,615,615 XITEM 18 615 J=(XNOW-XMIN)/DELX + 1.5+XSF XITEM 19 C XITEM 20 C DETERMINE LABEL FOR OBSERVATION XITEM 21 C XITEM 22 IF(IABS(ID)-3)503,501,503 XITEM 23 503 IF(IABS(ID)-2)500,501,500 XITEM 24 500 IF(ITEM(J).NE.PTID(2)) ITEM(J)=PTID(JX) XITEM 25 GOTO 605 XITEM 26 501 ITEM(J)=PTID(INDEX) XITEM 27 605 CONTINUE XITEM 28 RETURN XITEM 29 END XITEM 30