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
.