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