1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAY 11, 1992 19:04:56 PAGE: 1 0REQUESTED OPTIONS (EXECUTE): OPT(2) XREF SXM 0OPTIONS IN EFFECT: NOLIST NOMAP XREF GOSTMT NODECK SOURCE TERM OBJECT FIXED TRMFLG SRCFLG NOSYM NORENT SDUMP(ISN) SXM NOVECTOR IL(DIM) NOTEST SC(*) NODC NOEC NOEMODE NOICA NODIRECTIVE NODBCS NOSAA NOPARALLEL NOSAVE NOTABS OPT(2) LANGLVL(77) NOFIPS FLAG(I) AUTODBL(NONE) LINECOUNT(60) CHARLEN(500) 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C C C CLUSPAC: Computer Programs for Mixture-Model Clustering C C C C COPYRIGHT (C) 1991, 1992 STANLEY L. SCLOVE C C C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C C C CMS FILE: MIXPDT CLUSPAC C C VERSION 2.4 11-MAY-92 C C VERSION 2.3 22-MAR-92 C C VERSION 2.2 15-NOV-89 C C C C PROGRAMMED BY: C C C C PROF. STANLEY L. SCLOVE, PH.D. 312/996-2681 C C INFORMATION & DECISION SCIENCES DEPT. M/C 294 C C COLLEGE OF BUSINESS ADMINISTRATION C C UNIVERSITY OF ILLINOIS AT CHICAGO C C BOX 4348, CHICAGO, IL 60680 C C C C C C C C C C CLUSPAC is a set of programs implementing clustering C C algorithms derived under the assumption of Gaussian C C class-conditional distributions. The ISDT* programs in C C CLUSPAC are based on the so-called "classification" C C likelihood. The MIX* programs are based on the mixture- C C model likelihood. C C C C Program MIXPDT in CLUSPAC is one of the C C mixture-model programs for clustering multivariate data. C C (For univariate data the "MIX1" programs may be used.) C C MIXPDT allows different covariance matrices; C C MIXPCM assumes a common covariance matrix across C C distributions. C C C C C C Input: C C ----- C C Number of clusters (K), initial values of means, prior C C probabilities, and covariance matrices. If desired, C C program ISDTPCM.CLUSPAC can be used to obtain these initial C C values. Use program MIXPDTA.CLUSPAC to try a range of C C numbers of clusters (values of K), with automatic setting C C of initial values. C 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAY 11, 1992 19:04:56 NAME:MAIN# PAGE: 2 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 C C C Program restrictions (can be modified): C C -------------------------------------- C C N, sample size, at most 1000; C C IP, number of variables, at most 20; C C K, number of clusters, at most 29; C C ITER, maximum number of iterations, 20. C C C C C C Subroutines called: C C MATEQ, which calls MATDT C C C C IV is a work array for subroutine MATEQ. C C C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C C C C C CONTROL CARDS: C C C C DATASET TITLE C C N, IN FORMAT (2X,I4) C C IP, IN FORMAT (3X,I2) C C DATFMT, IN FORMAT (18A4), E.G., (4F4.1) C C "DATFMT" WILL ALSO BE USED FOR OUTPUT, SO ALLOW AT LEAST C C ONE BLANK AT THE BEGINNING FOR CARRIAGE CONTROL. C C DATA, ONE CASE AT A TIME, IN FORMAT SPECIFIED BY DATFMT C C K, NUMBER OF CLUSTERS, IN FORMAT (2X,I2) C C MEANFT, in format (18A4) C C K INITIAL MEANS, IN FORMAT SPECIFIED BY MEANFT C C C C K INITIAL VALUES OF MIXING PROBABILITIES, C C IN FORMAT (5X,F3.2). C C C C COVFMT, in format (18A4). C C Initial value of SIGMA(1,JV,JW), JV,JW=1,IP, cov mx of Pop 1 C C Initial value of SIGMA(2,JV,JW), JV,JW=1,IP, cov mx of Pop 2 C C . C C . C C . C C Initial value of SIGMA(K,JV,JW), JV,JW=1,IP, cov mx of Pop K C C C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C 1 DIMENSION X(1000,20),SUM(29,20) 2 DIMENSION F(1000,29),PP(29,1000) 3 DIMENSION ICLUS(1000),ICLSOL(1000) 4 DIMENSION DENOM(1000),XMXPR(1000) 5 DIMENSION DSQ(29) 6 DIMENSION TITLE(18) 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAY 11, 1992 19:04:56 NAME:MAIN# PAGE: 3 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 7 DIMENSION NG(29),XBAR(29,20) 8 DIMENSION DATFMT(18) 9 DIMENSION MEANFT(18) 10 DIMENSION COVFMT(18) 11 DIMENSION SSD(29,20,20),SIGMA(29,20,20) 12 DIMENSION WGSS(20,20) 13 DIMENSION VARHAT(20,20) 14 DIMENSION XMIN(20),XMAX(20) 15 DIMENSION IV(20,20) 16 DIMENSION P(20,20),PREC(29,20,20) 17 DIMENSION XIDET(29),XLGDET(29),TRUDET(29) 18 DIMENSION PR(29) C 19 DOUBLE PRECISION SUM,SUMPXF 20 DOUBLE PRECISION WGSS,SSD 21 DOUBLE PRECISION VARHAT 22 DOUBLE PRECISION P 23 DOUBLE PRECISION DET,TRUDET 24 DOUBLE PRECISION DSQ 25 DOUBLE PRECISION XBAR 26 DOUBLE PRECISION DEVV,DEVW 27 DOUBLE PRECISION F C C C C FLOW OF PROGRAM: C C C READ DATA AND INITIAL PARAMETER ESTIMATES. C INVERT COVARIANCE MATRICES. C PRINT INITIAL PARAMETER ESTIMATES. C C ITERATION: C COMPUTE ESTIMATED VALUES OF PROBABILITY DENSITY FUNCTIONS. C COMPUTE LIKELIHOOD AND VALUES OF MODEL-SELECTION CRITERIA. C COMPUTE CLUSTER-MEMBERSHIP PROBABILITIES. C UPDATE PARAMETER ESTIMATES (INCLUDING C INVERSE COVARIANCE MATRICES). C CLUSTER BY MAXIMUM PROBABILITY OF CLUSTER MEMBERSHIP. C IF CLUSTERING HASN'T CHANGED, STOP AND PRINT RESULTS. C OTHERWISE DO ANOTHER ITERATION (UNLESS 20 HAVE C ALREADY BEEN DONE). C C C 28 DATA PI/3.141593/ C 29 READ (5,36000) TITLE C C WRITE PROGRAM INFORMATION. 30 WRITE (6,10000) 31 WRITE (6,40000) TITLE C C READ SAMPLE SIZE, N. 32 READ (5,12000) N 33 XN = N 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAY 11, 1992 19:04:56 NAME:MAIN# PAGE: 4 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 34 WRITE (6,15000) N C READ NUMBER OF VARIABLES, IP. 35 READ (5,54000) IP 36 WRITE (6,13000) IP C C READ DATA FORMAT. 37 READ (5,36000) DATFMT C C READ DATA. 38 DO 500 I = 1,N 1 39 READ (5,DATFMT) (X(I,JV), JV = 1,IP) 1 40 IF (I .EQ. 1) GO TO 100 1 41 GO TO 300 1 42 100 CONTINUE C COMPUTE MINIMA AND MAXIMA: 1 43 DO 200 JV = 1,IP 2 44 XMAX(JV) = X(1,JV) 2 45 XMIN(JV) = X(1,JV) 2 46 200 CONTINUE 1 47 300 CONTINUE 1 48 DO 400 JV = 1,IP 2 49 IF (X(I,JV) .LT. XMIN(JV)) XMIN(JV)=X(I,JV) 2 51 IF (X(I,JV) .GT. XMAX(JV)) XMAX(JV)=X(I,JV) 2 53 400 CONTINUE 1 54 500 CONTINUE C WRITE MINIMA AND MAXIMA: 55 WRITE (6,64000) 56 WRITE (6,DATFMT) (XMIN(JV),JV=1,IP) 57 WRITE (6,66000) 58 WRITE (6,DATFMT) (XMAX(JV),JV=1,IP) C C READ K, NUMBER OF CLUSTERS. 59 READ (5,11000) K 60 WRITE (6,26000) K C C READ INITIAL MEANS C READ INPUT FORMAT FOR MEANS: 61 READ (5,36000) MEANFT 62 DO 600 IC=1,K 1 63 READ (5,MEANFT) (XBAR(IC,JV), JV=1,IP) C 1 64 600 CONTINUE C READ INITIAL VALUES OF MIXING PROBABILITIES: 65 DO 650 IC=1,K 1 66 READ (5,25000) PR(IC) C 1 67 650 CONTINUE C C C READ INITIAL VALUES OF COVARIANCE MATRICES: C 68 READ(5,36000) COVFMT 69 DO 660 IC = 1,K 1 70 DO 660 JV = 1,IP 2 71 READ(5,COVFMT) (SIGMA(IC,JV,JW),JW=1,IP) 2 72 660 CONTINUE 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAY 11, 1992 19:04:56 NAME:MAIN# PAGE: 5 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 C C 73 WRITE(6,22000) 74 ITER = 1 75 700 CONTINUE 76 WRITE (6,32000) ITER C WRITE CURRENT ESTIMATES OF PARAMETERS: 77 WRITE (6,23000) C WRITE CURRENT ESTIMATES OF MEAN VECTORS:-- 78 WRITE (6,20000) 79 DO 710 IC = 1,K 1 80 WRITE (6,MEANFT) ( XBAR(IC,JV), JV=1,IP ) 1 81 710 CONTINUE 82 WRITE (6,35000) (PR(IC), IC = 1,K) C WRITE CURRENT ESTIMATES OF COVARIANCE MATRICES:- 83 WRITE(6,42000) 84 DO 670 IC = 1,K 1 85 WRITE(6,42500) IC 1 86 DO 670 JV = 1,IP 2 87 WRITE(6,COVFMT) (SIGMA(IC,JV,JW),JW=1,IP) 2 88 670 CONTINUE C 89 DO 680 IC = 1,K 1 90 DO 675 JV = 1,IP 2 91 DO 675 JW = 1,IP 3 92 VARHAT(JV,JW) = SIGMA(IC,JV,JW) 3 93 675 CONTINUE C C CALL SUBROUTINE TO COMPUTE INVERSE COVARIANCE MATRIX: C SET PARAMETERS OF SUBROUTINE CALL: 1 94 IDET = 1 1 95 NRS1 = 0 C C 1 96 CALL MATEQ(VARHAT,IP,20,JFLG,DET,IDET,IV,NRS1,P,20) C ON RETURN, P CONTAINS THE INVERSE MATRIX. C C C GENERAL FORM OF CALL IS: C CALL MATEQ(A,M,N,JFLG,DET,IDET,IV,NRS1,P,LL) C SEE SUBROUTINE LISTING FOR FULLER EXPLANATION. C DET(VARHAT) = DET*10.0**IDET C C ON RETURN, P CONTAINS THE INVERSE MATRIX. C C 1 97 DO 676 JV = 1,IP 2 98 DO 676 JW = 1,IP 3 99 PREC(IC,JV,JW) = P(JV,JW) 3 100 676 CONTINUE C 1 101 XIDET(IC) = IDET 1 102 XLGDET(IC) = DLOG(DET) + XIDET(IC)*ALOG(10.0) 1 103 TRUDET(IC) = DET*(10.0**IDET) 1 104 680 CONTINUE C 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAY 11, 1992 19:04:56 NAME:MAIN# PAGE: 6 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 C C C C COMMENCE COMPUTATION OF F(I,IC), I=1,...,N, IC=1,...,K: C C COMPUTE MAHALANOBIS D-SQUARE BETWEEN THE I-TH OBSERVATION AND C THE L-TH MEAN, L=1,2,...,K, I=1,2,...,N:-- C 105 DO 1200 I = 1,N C 1 106 DO 1000 L=1,K 2 107 DSQ(L) = 0.0 2 108 DO 1310 JV=1,IP 3 109 DEVV = XBAR(L,JV) - X(I,JV) 3 110 DO 1310 JW=1,IP 4 111 DEVW = XBAR(L,JW) - X(I,JW) 4 112 DSQ(L) = DSQ(L) + DEVV*PREC(L,JV,JW)*DEVW 4 113 1310 CONTINUE 2 114 ZSQ = DSQ(L) C IF D-SQ IS INORDINATELY LARGE, SET VALUE OF PDF TO ZERO C (IT IS EXTREMELY SMALL ANYWAY, AND THIS AVOIDS UNDERFLOW): 2 115 IF ( ZSQ/2.0 .LE. 174.673 ) GO TO 1090 2 116 F(I, L) = 0.0 2 117 GO TO 1100 C 2 118 1090 CONTINUE 2 119 F(I, L) = EXP(-ZSQ/2.0) 2 120 F(I, L) = F(I,L)/DSQRT(TRUDET(L)) 2 121 1100 CONTINUE 2 122 1000 CONTINUE 1 123 1200 CONTINUE C C COMPUTE LOG LIKELIHOOD AND VALUES OF MODEL-SELECTION CRITERIA: C 124 SUMLNF = 0.0 125 DO 2200 I=1,N 1 126 SUMPXF = 0.0 1 127 DO 2100 IC=1,K 2 128 SUMPXF = SUMPXF + PR(IC)*F(I,IC) 2 129 2100 CONTINUE C 1 130 IF ( SUMPXF .LE. 0.0 ) GO TO 2200 C 1 131 SUMLNF = SUMLNF + DLOG(SUMPXF) 1 132 2200 CONTINUE C 133 XMN2LL = -2.0*SUMLNF C 134 WRITE (6,30000) XMN2LL C C C COMPUTE MODEL SELECTION CRITERIA C C PARAMETERS: C K MEAN VECTORS OF DIMENSION P AND K P-BY-P COVARIANCE MATRICES, C WHERE P IS THE NUMBER OF VARIABLES, AND K-1 INDEPENDENT MIXING 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAY 11, 1992 19:04:56 NAME:MAIN# PAGE: 7 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 C PROBABILITIES. 135 NOPARM = K*IP + K*IP*(IP+1)/2 + (K-1) 136 WRITE (6,72000) NOPARM 137 AIC = XMN2LL + 2.0*NOPARM 138 SCH = XMN2LL + ALOG(XN)*NOPARM 139 WRITE (6,70000) AIC 140 WRITE (6,71000) SCH C C C COMPUTE POSTERIOR PROBABILITIES OF GROUP MEMBERSHIP, C PP(IC,I), THE CONDITIONAL PROBABILITY OF POP'N IC, GIVEN X(I), C AS PR(IC)*F(I,IC)/DENOM(I): 141 DO 1350 I = 1,N 1 142 DENOM(I) = 0.0 1 143 DO 1350 IC=1,K 2 144 DENOM(I) = DENOM(I) + PR(IC)*F(I,IC) 2 145 1350 CONTINUE 146 DO 1400 I = 1,N 1 147 DO 1400 IC=1,K 2 148 IF ( DENOM(I) .GT. 0.0 ) GO TO 1410 2 149 PP(IC,I) = 0.0 2 150 GO TO 1400 2 151 1410 CONTINUE 2 152 PP(IC,I)= PR(IC)*F(I,IC)/DENOM(I) 2 153 1400 CONTINUE C 154 WRITE(6,76000) 155 DO 1420 I = 1,4 1 156 WRITE(6,82000) ( PP(IC,I), IC = 1,K ) 1 157 1420 CONTINUE C C C UPDATE PARAMETER ESTIMATES: C C C UPDATE CLUSTER MIXING PROBABILITIES PR(IC):-- 158 DO 3800 IC = 1,K 1 159 PR(IC) = 0.0 1 160 DO 3800 I = 1,N 2 161 PR(IC) = PR(IC) + PP(IC,I) 2 162 3800 CONTINUE 163 DO 3900 IC = 1,K 1 164 PR(IC) = PR(IC)/N 1 165 3900 CONTINUE 166 DO 1750 IG = 1,K 1 167 DO 1750 JV = 1,IP 2 168 SUM(IG,JV) = 0.0 2 169 DO 1750 JW = 1,IP 3 170 SSD(IG,JV,JW) = 0.0 3 171 1750 CONTINUE C C UPDATE MEANS:-- 172 DO 1875 IC = 1,K 1 173 DO 1875 JV = 1,IP 2 174 DO 1875 I = 1,N 3 175 SUM(IC,JV) = SUM(IC,JV) + PP(IC,I)*X(I,JV) 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAY 11, 1992 19:04:56 NAME:MAIN# PAGE: 8 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 3 176 1875 CONTINUE 177 DO 1900 IG = 1,K 1 178 IF (N*PR(IG) .LT. 0.5) GO TO 2050 1 179 GO TO 2150 1 180 2050 WRITE (6,74000) IG 1 181 GO TO 4000 1 182 2150 CONTINUE 1 183 DO 1900 JV = 1,IP 2 184 XBAR(IG,JV) = SUM(IG,JV)/(N*PR(IG)) 2 185 1900 CONTINUE C C 186 DO 3600 IC=1,K 1 187 DO 3600 JV=1,IP 2 188 DO 3600 JW=JV,IP 3 189 DO 3600 I=1,N 4 190 DEVV=X(I,JV)-XBAR(IC,JV) 4 191 DEVW=X(I,JW)-XBAR(IC,JW) 4 192 TERM=PP(IC,I)*DEVV*DEVW 4 193 SSD(IC,JV,JW)=SSD(IC,JV,JW)+TERM 4 194 3600 CONTINUE C 195 DO 3700 IC=1,K 1 196 DO 3700 JV=2,IP 2 197 DO 3700 JW=1,JV-1 3 198 SSD(IC,JV,JW)=SSD(IC,JW,JV) 3 199 3700 CONTINUE C 200 DO 3710 IC=1,K 1 201 DO 3710 JV=1,IP 2 202 DO 3710 JW=1,IP 3 203 SIGMA(IC,JV,JW)=SSD(IC,JW,JV)/(N*PR(IC)) 3 204 3710 CONTINUE C C POOL: C 205 DO 1950 JV = 1,IP 1 206 DO 1950 JW = 1,IP 2 207 WGSS(JV,JW) = 0.0 2 208 1950 CONTINUE C 209 DO 2300 JV = 1,IP 1 210 DO 2300 JW = 1,IP 2 211 DO 2300 IC = 1,K 3 212 WGSS(JV,JW) = WGSS(JV,JW) + PR(IC)*SSD(IC,JV,JW) 3 213 2300 CONTINUE C COMPUTE VARHAT, MLE OF COMMON COVARIANCE MATRIX: 214 DO 2400 JV = 1,IP 1 215 DO 2400 JW = 1,IP 2 216 VARHAT(JV,JW) = WGSS(JV,JW)/N 2 217 2400 CONTINUE C (END PARAMETER-ESTIMATE UPDATE SEQUENCE) C 218 IF (ITER .EQ. 1) GO TO 900 C STORE OLD LABELS: 219 DO 800 I = 1,N 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAY 11, 1992 19:04:56 NAME:MAIN# PAGE: 9 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 1 220 ICLSOL(I) = ICLUS(I) 1 221 800 CONTINUE C 222 900 CONTINUE C C COMPUTE NEW LABELS BY MAX POSTERIOR PROBABILITY: 223 DO 1600 I = 1,N 1 224 XMXPR(I) = PP(1,I) 1 225 ICLUS(I) = 1 1 226 DO 1600 IC = 2,K 2 227 IF ( PP(IC,I) .GT. XMXPR(I) ) GO TO 1500 2 228 GO TO 1600 2 229 1500 XMXPR(I) = PP(IC,I) 2 230 ICLUS(I) = IC 2 231 1600 CONTINUE C C 232 IF (N .GE. 31) GO TO 250 233 WRITE (6,15500) C 234 WRITE (6,16000) 235 WRITE (6,18000) (I, ICLUS(I), I=1,N) 236 250 CONTINUE C C C C 237 IF (ITER .EQ. 1) GO TO 3000 238 DO 2900 I = 1,N 1 239 IF (ICLUS(I) .EQ. ICLSOL(I)) GO TO 2900 1 240 GO TO 3000 1 241 2900 CONTINUE 242 GO TO 3300 243 3000 CONTINUE 244 ITER = ITER + 1 245 IF (ITER.GE.21) GO TO 3100 246 GO TO 700 247 3100 WRITE (6,68000) 248 4000 STOP C C 249 3300 CONTINUE C C C COUNT NUMBERS IN CLUSTERS: 250 DO 3310 IC = 1,K 1 251 NG(IC) = 0 1 252 3310 CONTINUE 253 DO 3320 I = 1,N 1 254 IC = ICLUS(I) 1 255 NG(IC) = NG(IC) + 1 1 256 3320 CONTINUE C C C PRINT FINAL RESULTS C 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAY 11, 1992 19:04:56 NAME:MAIN# PAGE: 10 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 C 257 WRITE (6,34000) (NG(IC),IC=1,K) 258 WRITE (6,35000) (PR(IC),IC=1,K) 259 DO 2800 IC = 1,K 1 260 WRITE (6,28000) IC, (XBAR(IC,JV),JV=1,IP) 1 261 2800 CONTINUE C 262 WRITE (6,42000) 263 DO 2500 IC = 1,K 1 264 WRITE(6,42500) IC 1 265 DO 2500 JV=1,IP 2 266 WRITE (6,48000) (SIGMA(IC,JV,JW),JW=1,IP) 2 267 2500 CONTINUE C 268 WRITE (6,52000) 269 DO 3500 I = 1,N 1 270 WRITE (6,50000) I, ICLUS(I) 1 271 WRITE (6,DATFMT) (X(I,JV),JV=1,IP) 1 272 3500 CONTINUE C 273 WRITE (6,84000) C C 274 STOP C C C 275 10000 FORMAT('1','**********************************************'/// X' MIXPDT CLUSPAC - MIXTURE-MODEL CLUSTERING OF CASES', X' VARYING COVARIANCE MATRICES '// X' Developed and programmed by: '// X' Prof. Stanley L. Sclove, Ph.D. '/ X' Information & Decision Sciences Dept. M/C 294 '/ X' College of Business Administration '/ X' University of Illinois at Chicago '/ X' Box 4348, Chicago, IL 60680-4348 '// X' 312/996-2681 '// X' MIXPDT CLUSPAC VERSION 2.4 11-MAY-92 '// X' COPYRIGHT 1991 STANLEY LOUIS SCLOVE. ', X' ALL RIGHTS RESERVED.'//) 276 11000 FORMAT(2X,I2) 277 12000 FORMAT(2X,I4) 278 13000 FORMAT(1X,'NUMBER OF VARIABLES USED....................',I3/) 279 15000 FORMAT(1X,'NUMBER OF OBSERVATIONS (SAMPLE SIZE), N.....',I3/) 280 15500 FORMAT(//1X,'CLUSTERING:'/) 281 16000 FORMAT(/,1X,'CASES AND LABELS:--'/) 282 18000 FORMAT(15(I5,I3)) 283 20000 FORMAT(/1X,' MEANS: '/) 284 22000 FORMAT(' FIRST ITERATION USES THE INITIAL PARAMETER ESTIMATES', X' PROVIDED BY THE USER.'/) 285 25000 FORMAT(5X, F3.2) 286 23000 FORMAT(/1X, 'CURRENT ESTIMATES OF PARAMETERS: '//) 287 26000 FORMAT('1',//,1X,'K = ',I2,' CLUSTERS'///) 288 28000 FORMAT(1X,'MEAN VECTOR FOR CLUSTER ',I2,': ',(8F13.5/)) 289 30000 FORMAT(/1X,' MINUS 2 LOG LIKELIHOOD = ', F13.5/) 290 32000 FORMAT(//,1X,'ITERATION ', I2,//) 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAY 11, 1992 19:04:56 NAME:MAIN# PAGE: 11 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 291 34000 FORMAT(/,1X,'NUMBERS:',1X,9(I10,3X)/) 292 35000 FORMAT(/,1X,' MIXING PROBABILITIES: ',3X,9(F5.3,3X)/) 293 36000 FORMAT(18A4) 294 40000 FORMAT(//,1X,18A4//) 295 42000 FORMAT(/,' COVARIANCE MATRICES: ',/) 296 42500 FORMAT(/,' POPULATION ',I3/) 297 48000 FORMAT(1X,8F13.5/) 298 50000 FORMAT(1X,I4,1X,I2) 299 52000 FORMAT(//,1X,'CASE, LABEL / DATA'//) 300 54000 FORMAT(3X,I2) 301 58000 FORMAT(/1X,'CONVERGENCE: NO CASE CHANGED CLUSTERS AFTER ', X'ITERATION ',I2,'. RESULTS ARE PRINTED BELOW.'//) 302 60000 FORMAT(/,1X,'JFLG = ',I2,'. IF JFLG=0, COMPUTATION OF DET', X' WENT WELL; OTHERWISE, THERE WAS TROUBLE OR MATRIX WAS ', X'ILL-CONDITIONED.'//) 303 62000 FORMAT(/1X,'DET = ',F13.5,' IDET = ',I3,5X, X'ACTUAL DET. = DET*10**IDET',//) 304 64000 FORMAT(/1X,'MINIMUM FOR EACH VARIABLE: ',/) 305 66000 FORMAT(/1X,'MAXIMUM FOR EACH VARIABLE: ',/) 306 68000 FORMAT(1X,'PROGRAM HAS NOT CONVERGED IN 20 ITERATIONS. STOP') 307 70000 FORMAT(1X,'AIC = ', F15.5,/) 308 71000 FORMAT(1X,'SCHWARZ CRITERION = ', F15.5,/) 309 72000 FORMAT(/1X,'NUMBER OF PARAMETERS = ',I4/) 310 74000 FORMAT(1X,'EXPECTED NUMBER OF OBSERVATIONS IN CLUSTER ',I3, X' IS LESS THAN 0.5. STOP.'/) 311 76000 FORMAT(//' POSTERIOR PROBS OF GROUP MEMBERSHIP FOR 1ST 4 CASES:') 312 82000 FORMAT(1X,29F5.2) 313 84000 FORMAT(1X,'PROGRAM ENDED NORMALLY.') 314 END 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAY 11, 1992 19:04:56 NAME:MAIN# PAGE: 12 0SYMBOL CROSS REFERENCE DICTIONARY 0PROGRAM NAME: MAIN# 0TAGS: A-ARRAY I-INTRINSIC FUNCTION S-ASSIGNED C-COMMON K-NAMED CONSTANT T-EXPLICITLY TYPED D-DUMMY ARGUMENT N-ENTRY V-INITIAL VALUE E-EQUIVALENCED P-PROMOTED X-EXTERNAL SUBPROGRAM F-STATEMENT FUNCTION Q-PADDED Y-DYNAMIC COMMON G-GENERIC NAME R-SUBPROGRAM NAME Z-EXTENDED COMMON 0NAME TYPE TAG DECLARED REFS (F:REFD S:SET B:REFD/MAY BE SET) +______ ____ ______ ________ ________________________________________ 0AIC R*4 137S 139F ALOG R*4 I 102 138 COVFMT R*4 A 10 68S 71 87F DATFMT R*4 A 8 37S 39 56F 58F 271F DENOM R*4 A 4 142S 144F 144S 148F 152F DET R*8 T 23 96B 102F 103F DEVV R*8 T 26 109S 112F 190S 192F DEVW R*8 T 26 111S 112F 191S 192F DLOG R*8 I 102 131 DSQ R*8 AT 5 24 107S 112F 112S 114F DSQRT R*8 I 120 EXP GI 119 F R*8 AT 2 27 116S 119S 120F 120S 128F 144F 152F I I*4 38S 39 40F 49F 50F 51F 52B 105S 109F 111F 116F 119F 120F 120B 125S 128B 141S 142F 144F 144F 144B 146S 148F 149F 152F 152F 152B 155S 156B 160S 161B 174S 175F 175B 189S 190F 191F 192B 219S 220F 220B 223S 224F 224F 225F 227F 227F 229F 229F 230B 235S 235F 235B 238S 239F 239B 253S 254B 269S 270F 270F 271B IC I*4 62S 63B 65S 66B 69S 71B 79S 80B 82S 82S 84S 85F 87B 89S 92F 99F 101F 102F 102F 103B 127S 128F 128B 143S 144F 144B 147S 149F 152F 152F 152B 156S 156S 158S 159F 161F 161F 161B 163S 164F 164B 172S 175F 175F 175B 186S 190F 191F 192F 193F 193B 195S 198F 198B 200S 203F 203F 203B 211S 212F 212B 226S 227F 229F 230B 250S 251B 254S 255F 255F 257S 257S 258S 258S 259S 260F 260B 263S 264F 266B ICLSOL I*4 A 3 220S 239F 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAY 11, 1992 19:04:56 NAME:MAIN# PAGE: 13 0TAGS: A-ARRAY I-INTRINSIC FUNCTION S-ASSIGNED C-COMMON K-NAMED CONSTANT T-EXPLICITLY TYPED D-DUMMY ARGUMENT N-ENTRY V-INITIAL VALUE E-EQUIVALENCED P-PROMOTED X-EXTERNAL SUBPROGRAM F-STATEMENT FUNCTION Q-PADDED Y-DYNAMIC COMMON G-GENERIC NAME R-SUBPROGRAM NAME Z-EXTENDED COMMON 0NAME TYPE TAG DECLARED REFS (F:REFD S:SET B:REFD/MAY BE SET) +______ ____ ______ ________ ________________________________________ 0ICLUS I*4 A 3 220F 225S 230S 235F 239F 254F 270F IDET I*4 94S 96B 101F 103F IG I*4 166S 168F 170B 177S 178F 180F 184F 184F 184B IP I*4 35S 36F 39 43F 48F 56F 58F 63 70F 71 80F 86F 87F 90F 91F 96B 97F 98F 108F 110F 135F 135F 135F 167F 169F 173F 183F 187F 188F 196F 201F 202F 205F 206F 209F 210F 214F 215F 260F 265F 266F 271F ITER I*4 74S 76F 218F 237F 244F 244S 245F IV I*4 A 15 96B JFLG I*4 96B JV I*4 39S 39S 43S 44F 44F 45F 45B 48S 49F 49F 50F 50F 51F 51F 52F 52B 56S 56S 58S 58S 63S 63S 70S 71B 80S 80S 86S 87B 90S 92F 92B 97S 99F 99B 108S 109F 109F 112B 167S 168F 170B 173S 175F 175F 175B 183S 184F 184B 187S 188F 190F 190F 193F 193B 196S 197F 198F 198B 201S 203F 203B 205S 207B 209S 212F 212F 212B 214S 216F 216B 260S 260S 265S 266B 271S 271S JW I*4 71S 71S 87S 87S 91S 92F 92B 98S 99F 99B 110S 111F 111F 112B 169S 170B 188S 191F 191F 193F 193B 197S 198F 198B 202S 203F 203B 206S 207B 210S 212F 212F 212B 215S 216F 216B 266S 266S K I*4 59S 60F 62F 65F 69F 79F 82F 84F 89F 106F 127F 135F 135F 135F 143F 147F 156F 158F 163F 166F 172F 177F 186F 195F 200F 211F 226F 250F 257F 258F 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAY 11, 1992 19:04:56 NAME:MAIN# PAGE: 14 0TAGS: A-ARRAY I-INTRINSIC FUNCTION S-ASSIGNED C-COMMON K-NAMED CONSTANT T-EXPLICITLY TYPED D-DUMMY ARGUMENT N-ENTRY V-INITIAL VALUE E-EQUIVALENCED P-PROMOTED X-EXTERNAL SUBPROGRAM F-STATEMENT FUNCTION Q-PADDED Y-DYNAMIC COMMON G-GENERIC NAME R-SUBPROGRAM NAME Z-EXTENDED COMMON 0NAME TYPE TAG DECLARED REFS (F:REFD S:SET B:REFD/MAY BE SET) +______ ____ ______ ________ ________________________________________ 0 259F 263F L I*4 106S 107F 109F 111F 112F 112F 112F 114F 116F 119F 120F 120F 120B MATEQ X 96F MEANFT I*4 A 9 61S 63 80F N I*4 32S 33F 34F 38F 105F 125F 141F 146F 160F 164F 174F 178F 184F 189F 203F 216F 219F 223F 232F 235F 238F 253F 269F NG I*4 A 7 251S 255F 255S 257F NOPARM I*4 135S 136F 137F 138F NRS1 I*4 95S 96B P R*8 AT 16 22 96B 99F PI R*4 V 28 PP R*4 A 2 149S 152S 156F 161F 175F 192F 224F 227F 229F PR R*4 A 18 66S 82F 128F 144F 152F 159S 161F 161S 164F 164S 178F 184F 203F 212F 258F PREC R*4 A 16 99S 112F SCH R*4 138S 140F SIGMA R*4 A 11 71S 87F 92F 203S 266F SSD R*8 AT 11 20 170S 193F 193S 198F 198S 203F 212F SUM R*8 AT 1 19 168S 175F 175S 184F SUMLNF R*4 124S 131F 131S 133F SUMPXF R*8 T 19 126S 128F 128S 130F 131F TERM R*4 192S 193F TITLE R*4 A 6 29S 31F TRUDET R*8 AT 17 23 103S 120F VARHAT R*8 AT 13 21 92S 96B 216S WGSS R*8 AT 12 20 207S 212F 212S 216F X R*4 A 1 39S 44F 45F 49F 50F 51F 52F 109F 111F 175F 190F 191F 271F XBAR R*8 AT 7 25 63S 80F 109F 111F 184S 190F 191F 260F XIDET R*4 A 17 101S 102F XLGDET R*4 A 17 102S XMAX R*4 A 14 44S 51F 52S 58F XMIN R*4 A 14 45S 49F 50S 56F XMN2LL R*4 133S 134F 137F 138F XMXPR R*4 A 4 224S 227F 229S XN R*4 33S 138F ZSQ R*4 114S 115F 119F 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAY 11, 1992 19:04:56 NAME:MAIN# PAGE: 15 0 0VARIABLES REFERENCED BUT NOT SET. (* POSSIBLY SET AS ARGUMENT.) 0DET* IV* JFLG* P* 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAY 11, 1992 19:04:56 NAME:MAIN# PAGE: 16 0LABEL CROSS REFERENCE DICTIONARY 0TAGS: A-USED AS ARGUMENT F-FORMAT S-USED IN ASSIGN STATEMENT B-OBJECT OF BRANCH N-NON-EXECUTABLE 0 LABEL TAG DEFINED REFERENCED +____________ ___ _______ ___________________________________________ 0 100 B 42 40 200 46 43 250 B 236 232 300 B 47 41 400 53 48 500 54 38 600 64 62 650 67 65 660 72 69 70 670 88 84 86 675 93 90 91 676 100 97 98 680 104 89 700 B 75 246 710 81 79 800 221 219 900 B 222 218 1000 122 106 1090 B 118 115 1100 B 121 117 1200 123 105 1310 113 108 110 1350 145 141 143 1400 B 153 146 147 150 1410 B 151 148 1420 157 155 1500 B 229 227 1600 B 231 223 226 228 1750 171 166 167 169 1875 176 172 173 174 1900 185 177 183 1950 208 205 206 2050 B 180 178 2100 129 127 2150 B 182 179 2200 B 132 125 130 2300 213 209 210 211 2400 217 214 215 2500 267 263 265 2800 261 259 2900 B 241 238 239 3000 B 243 237 240 3100 B 247 245 3300 B 249 242 3310 252 250 3320 256 253 3500 272 269 3600 194 186 187 188 189 3700 199 195 196 197 3710 204 200 201 202 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAY 11, 1992 19:04:56 NAME:MAIN# PAGE: 17 0TAGS: A-USED AS ARGUMENT F-FORMAT S-USED IN ASSIGN STATEMENT B-OBJECT OF BRANCH N-NON-EXECUTABLE 0 LABEL TAG DEFINED REFERENCED +____________ ___ _______ ___________________________________________ 0 3800 162 158 160 3900 165 163 4000 B 248 181 10000 NF 275 30 11000 NF 276 59 12000 NF 277 32 13000 NF 278 36 15000 NF 279 34 15500 NF 280 233 16000 NF 281 234 18000 NF 282 235 20000 NF 283 78 22000 NF 284 73 23000 NF 286 77 25000 NF 285 66 26000 NF 287 60 28000 NF 288 260 30000 NF 289 134 32000 NF 290 76 34000 NF 291 257 35000 NF 292 82 258 36000 NF 293 29 37 61 68 40000 NF 294 31 42000 NF 295 83 262 42500 NF 296 85 264 48000 NF 297 266 50000 NF 298 270 52000 NF 299 268 54000 NF 300 35 58000 NF 301 UNREFERENCED 60000 NF 302 UNREFERENCED 62000 NF 303 UNREFERENCED 64000 NF 304 55 66000 NF 305 57 68000 NF 306 247 70000 NF 307 139 71000 NF 308 140 72000 NF 309 136 74000 NF 310 180 76000 NF 311 154 82000 NF 312 156 84000 NF 313 273 0*STATISTICS* SOURCE STATEMENTS: 312, PROGRAM SIZE: 663824 BYTES, PROGRAM NAME: MAIN#, PAGE: 1 *STATISTICS* NO DIAGNOSTICS GENERATED. **MAIN#** END OF COMPILATION 1 ****** TIME STAMP: 92.13219.04.56 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAY 11, 1992 19:04:57 PAGE: 18 0OPTIONS IN EFFECT: NOLIST NOMAP XREF GOSTMT NODECK SOURCE TERM OBJECT FIXED TRMFLG SRCFLG NOSYM NORENT SDUMP(ISN) SXM NOVECTOR IL(DIM) NOTEST SC(*) NODC NOEC NOEMODE NOICA NODIRECTIVE NODBCS NOSAA NOPARALLEL NOSAVE NOTABS OPT(2) LANGLVL(77) NOFIPS FLAG(I) AUTODBL(NONE) LINECOUNT(60) CHARLEN(500) 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 C C 1 SUBROUTINE MATEQ(A,M,N,JFLG,DET,IDET,IV,NRS1,P,LL) C SUBROUTINE MATEQ IS DMATEQ FROM THE UICC SUBROUTINE LIBRARY. C PAGE 1 C C SUBROUTINE DMATEQ C ***************** C THIS ROUTINE WILL SOLVE A REAL*8 SYSTEM OF LINEAR EQUATIONS,COMPUTE C THE DETERMINANT, WITHOUT UNDERFLOW OR OVERFLOW, OF A REAL*8 MATRIX, C AND/OR INVERT A REAL*8 MATRIX. C CALLING SEQUENCE: C CALL DMATEQ(A,N,IA,JFLG,DET,IDET,IV,NRS,P,IP) WHERE; C A (INPUT) - IS THE REAL*8 MATRIX ON WHICH THE ROUTINE IS C TO WORK. IN THE PROCESS OF COMPUTATION THE C CONTENTS OF THIS MATRIX ARE DESTROYED. C N (INPUT) - IS AN INTEGER*4 VARIABLE WHICH SPECIFIES THE C ORDER OF THE A MATRIX. C IA (INPUT) - IS AN INTEGER*4 VARIABLE WHICH SPECIFIES THE C ACTUAL ROW DIMENSION OF A AS DIMENSIONED IN C THE CALLING PROGRAM. IA MUST BE GREATER THAN C OR EQUAL TO N. C JFLG (OUTPUT) - IS AN INTEGER*4 RETURN CODE VARIABLE. UPON C RETURN FROM DMATEQ IF; C JFLG=0, ALL WENT WELL. C JFLG=1, THE A MATRIX WAS SINGULAR OR NEAR C SINGULAR AND THE COMPUTATIONS COULD NOT BE C COMPLETED. THE CONTENTS OF THE VARIABLES C A, DET, IDET AND P ARE MEANINGLESS. C DET (OUTPUT) - IS A REAL*8 VARIABLE WHICH CONTAINS THE C DETERMINANT OF A. (SEE IDET) C IDET (INPUT) - IS AN INTEGER*4 VARIABLE. ON INPUT IF; C IDET=0, NO DETERMINANT IS CALCULATED. C IDET NOT 0, THE DETERMINANT OF A IS COMPUTED. C ON OUTPUT IDET CONTAINS THE POWER OF 10 C THAT DET SHOULD BE MULTIPLIED BY TO GIVE THE C CORRECT VALUE OF THE DETERMINANT. I.E. C DET(A)=DET*10.0D0**IDET. C IF DET(A) CAN BE COMPUTED WITHOUT UNDER OR C OVERFLOW, THEN IDET=0 OTHERWISE IDET IS SET C TO THE PROPER VALUE SO THAT NO UNDER OR OVER- C FLOW WILL OCCUR IN COMPUTING DET. C IV (INPUT) - IS AN INTEGER*4 WORK ARRAY WHICH SHOULD BE C DIMENSIONED AT LEAST IV(N). C PAGE 2 C C NRS (INPUT) - IS AN INTEGER*4 VARIABLE WITH THE FOLLOWING C INTERPRETATION: C NRS>0, SOLVE A SYSTEM OF LINEAR EQUATIONS C WITH NRS RIGHT HAND SIDES. C NRS=0, INVERT THE A MATRIX. C NRS<0, ONLY COMPUTE THE DETERMINANT OF A. 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAY 11, 1992 19:04:57 NAME:MAIN# PAGE: 19 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 C IN THIS CASE IDET MUST BE DIFFERENT FROM 0. C P (INPUT) - IS A REAL*8 ARRAY WITH THE FOLLOWING INTER- C PRETATION: C IF NRS>0, THEN P CONTAINS THE NRS RIGHT HAND C SIDES STORED BY COLUMNS. IN THIS CASE P MUST C BE DIMENSIONED AT LEAST P(N,NRS). ON RETURN C THE COLUMNS OF P ARE REPLACED BY THE RESPEC- C TIVE SOLUTIONS. C IF NRS=0, THEN P MUST BE DIMENSIONED AT LEAST C P(N,N). ON RETURN P WILL CONTAIN THE INVERSE C OF A. C IF NRS<0,THEN P NEED ONLY BE A DUMMY VARIABLE C IN THIS CASE P IS NEVER ACCESSED BY DMATEQ. C IP (INPUT) - IS AN INTEGER*4 VARIABLE WHICH CONTAINS THE C ACTUAL ROW DIMENSION OF P AS DIMENSIONED IN C THE CALLING PROGRAM. IP MUST BE GREATER THAN C OR EQUAL TO N. C NOTE: IMMEDIATELY ON RETURN FROM DMATEQ THE CONDITION CODE FLAG, C JFLG, SHOULD BE INTERROGATED. IF JFLG=1, THEN THE ROUTINE C COULD NOT COMPUTE A SOLUTION. C METHOD - THE ALGORITHM USED IS GAUSSIAN ELIMINATION WITH PARTIAL C -1 C PIVOTING. IN ESSENCE THE ROUTINE GENERATES A MATRIX L SUCH C -1 C THAT L *A = U, WHERE U IS AN UPPER TRIANGULAR MATRIX. THEN IT C SOLVES THE SYSTEM A*X = P BY MEANS OF THE EQUIVALENT SYSTEM C -1 -1 C U*X = L *A*X = L *P BY BACK SUBSTITUTION. C -1 C THE L MATRIX CAN BE WRITTEN AS A PRODUCT OF THE FORM C -1 C L = L *P *....*L *P WHERE EACH P IS A PERMUTATION C N-1 N-1 1 1 K C MATRIX OBTAINED BY INTERCHANGING AT MOST TWO ROWS OF THE C IDENTITY MATRIX. ( THIS REPRESENTS THE INTERCHANGING OF TWO C ROWS). THE L MATRICES ARE ELIMINATION MATRICES WHICH ARE C K C CHOSEN TO INTRODUCE ZEROS IN THE LAST N-K ENTRIES OF THE K-TH C COLUMN OF THE MATRIX. C PAGE 3 C C -1 -1 C THE CALCULATIONS OF L *A AND L *P ARE DONE BY PERFORMING C THE PERMUTATIONS ON A AND P RESPECTIVELY. THE ACTUAL L AND P C K K C ARE NOT COMPUTED. C SUBROUTINES CALLED: DMATDT C REFERENCE: C G. W. STEWART, INTRODUCTION TO MATRIX COMPUTATIONS, C ACADEMIC PRESS, 1973. 2 REAL*8 A(N,1),DET,P(LL,1) 3 REAL*8 DNORM,DEN,DMULT,DSUM,DISIGN 4 DIMENSION IV(1) 5 NRS=NRS1 6 IF (NRS.EQ.0) IDET=1 8 DISIGN=1.0D+00 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAY 11, 1992 19:04:57 NAME:MATEQ PAGE: 20 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 9 DET=0.0D+00 10 JFLG=0 C C JFLG IS A TROUBLE FLAG.UPON EXIT IF JFLG=0 THEN THE MATRIX WAS PROCESS C WITHOUT TROUBLE.IF JFLG=1 EITHER THE MATRIX IS SINGULAR OR TROUBLE C OCCURED.ISIGN=-ISIGN EVERY TIME A ROW IS INTERCHANGED.THIS IS USED TO C INSURE THAT THE DETERMINANT HAS THE PROPER SIGN. C 11 M1=M-1 12 DO 100 I=1,M 1 13 100 IV(I)=I 14 IF (NRS) 500,200,500 15 200 DO 300 I=1,M 1 16 DO 300 J=1,M 2 17 300 P(I,J)=0.0D+00 18 DO 400 I=1,M 1 19 400 P(I,I)=1.0D+00 20 NRS=M C C INSTEAD OF ACTUALLY INTERCHANGING ROWS A POINTER ARRAY IS USED TO KEEP C TRACK OF THE ROW POSITIONS. C C BEGIN ELIMINATION LOOP. C 21 500 DO 1200 K=1,M1 1 22 ICOL=K 1 23 IPCOL=K C C SEARCHING FOR LARGEST ELEMENT IN ABSOLUTE VALUE IN COLUMN K. C 1 24 DNORM=A(IV(K),K) 1 25 IFLG=0 1 26 KK=K+1 1 27 DO 600 J=KK,M 2 28 IF (DABS(A(IV(J),K)).LE.DABS(DNORM)) GO TO 600 2 29 IFLG=1 2 30 IPCOL=IV(J) 2 31 DNORM=A(IPCOL,K) 2 32 600 CONTINUE C C IF IFLG=0 NO ROW INTERCHANGE TOOK PLACE.IF IFLG=1 A ROW INTERCHANGE C TOOK PLACE AND THE POINTER ARRAY IV MUST BE UPDATED. C 1 33 IF (IFLG.EQ.0) GO TO 800 1 34 ISAVE=IV(ICOL) 1 35 IV(ICOL)=IPCOL 1 36 ICOL1=ICOL+1 1 37 DO 700 L=ICOL1,M 2 38 IF (IV(L).EQ.IPCOL) IV(L)=ISAVE 2 40 700 CONTINUE 1 41 DISIGN=-DISIGN 1 42 800 IF (DNORM.EQ.0.0D+00) GO TO 1900 C C BEGIN ELIMINATION OF ROW BELOW IV(K).DEN IS THE PIVOT ELEMENT. C 1 43 K1=K+1 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAY 11, 1992 19:04:57 NAME:MATEQ PAGE: 21 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 1 44 DO 1100 IM=K1,M C C BEFORE ACTUALLY ELIMINATING WE CHECK TO SEE IF A(IV(IM),K) HAS C ALREADY BEEN ANNIHILATED. C 2 45 IF (A(IV(IM),K).EQ.0.0D+00) GO TO 1100 C C CALCULATE ELIMINATION FACTOR. C 2 46 DMULT=-A(IV(IM),K) C C WE NOW CALCULATE VALUE OF OTHER ELEMENTS IN ROW IV(IM). C 2 47 DO 900 NN=K1,M 3 48 900 A(IV(IM),NN)=(DMULT*A(IV(K),NN))/DNORM+A(IV(IM),NN) 2 49 IF (NRS.LE.0) GO TO 1100 2 50 DO 1000 IN=1,NRS 3 51 1000 P(IV(IM),IN)=(DMULT*P(IV(K),IN))/DNORM+P(IV(IM),IN) 2 52 1100 CONTINUE 1 53 1200 CONTINUE C C CALCULATE VALUE OF DETERMINANT. C 54 IF (A(IV(M),M).EQ.0.0D0) GO TO 1900 55 DET=DISIGN 56 IF (IDET.NE.0) CALL DMATDT(A,N,M,DET,IV,IDET) 58 IF (DET.EQ.0.0D+00) GO TO 1900 59 IF (NRS.LE.0) GO TO 2000 C C WE START SOLVING RIGHT HAND SIDES.THE SOLUTION REPLACES THE RIGHT HAND C VECTOR. C 60 1300 N1=M-1 61 DO 1600 JJ=1,NRS C C BEGIN BACK SUBSTITUTION. C 1 62 P(IV(M),JJ)=P(IV(M),JJ)/A(IV(M),M) 1 63 DO 1500 I=1,N1 2 64 DSUM=0.0D+00 2 65 DO 1400 J=1,I 3 66 1400 DSUM=DSUM-A(IV(M-I),M-J+1)*P(IV(M-J+1),JJ) 2 67 1500 P(IV(M-I),JJ)=(P(IV(M-I),JJ)+DSUM)/A(IV(M-I),M-I) 1 68 1600 CONTINUE 69 DO 1800 JJ=1,NRS 1 70 DO 1700 IND=1,M 2 71 1700 A(IND,1)=P(IV(IND),JJ) 1 72 DO 1800 IND=1,M 2 73 1800 P(IND,JJ)=A(IND,1) 74 RETURN 75 1900 JFLG=1 76 IDET=0 77 2000 RETURN 78 END 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAY 11, 1992 19:04:57 NAME:MATEQ PAGE: 22 0SYMBOL CROSS REFERENCE DICTIONARY 0PROGRAM NAME: MATEQ 0TAGS: A-ARRAY I-INTRINSIC FUNCTION S-ASSIGNED C-COMMON K-NAMED CONSTANT T-EXPLICITLY TYPED D-DUMMY ARGUMENT N-ENTRY V-INITIAL VALUE E-EQUIVALENCED P-PROMOTED X-EXTERNAL SUBPROGRAM F-STATEMENT FUNCTION Q-PADDED Y-DYNAMIC COMMON G-GENERIC NAME R-SUBPROGRAM NAME Z-EXTENDED COMMON 0NAME TYPE TAG DECLARED REFS (F:REFD S:SET B:REFD/MAY BE SET) +______ ____ ______ ________ ________________________________________ 0A R*8 ADT 1 2 24F 28F 31F 45F 46F 48F 48F 48S 54F 57B 62F 66F 67F 71S 73F DABS R*8 I 28 28 DEN R*8 T 3 UNREFERENCED DET R*8 DT 1 2 9S 55S 57B 58F DISIGN R*8 T 3 8S 41F 41S 55F DMATDT X 57F DMULT R*8 T 3 46S 48F 51F DNORM R*8 T 3 24S 28F 31S 42F 48F 51F DSUM R*8 T 3 64S 66F 66S 67F I I*4 12S 13F 13B 15S 17B 18S 19F 19B 63S 65F 66F 67F 67F 67F 67B ICOL I*4 22S 34F 35F 36F ICOL1 I*4 36S 37F IDET I*4 D 1 7S 56F 57B 76S IFLG I*4 25S 29S 33F IM I*4 44S 45F 46F 48F 48F 51F 51B IN I*4 50S 51F 51F 51B IND I*4 70S 71F 71B 72S 73F 73B IPCOL I*4 23S 30S 31F 35F 38F ISAVE I*4 34S 39F IV I*4 AD 1 4 13S 24F 28F 30F 34F 35S 38F 39S 45F 46F 48F 48F 48F 51F 51F 51F 54F 57B 62F 62F 62F 66F 66F 67F 67F 67F 71F J I*4 16S 17B 27S 28F 30B 65S 66F 66B JFLG I*4 D 1 10S 75S JJ I*4 61S 62F 62F 66F 67F 67B 69S 71F 73B K I*4 21S 22F 23F 24F 24F 26F 28F 31F 43F 45F 46F 48F 51B KK I*4 26S 27F K1 I*4 43S 44F 47F L I*4 37S 38F 39B 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAY 11, 1992 19:04:57 NAME:MATEQ PAGE: 23 0TAGS: A-ARRAY I-INTRINSIC FUNCTION S-ASSIGNED C-COMMON K-NAMED CONSTANT T-EXPLICITLY TYPED D-DUMMY ARGUMENT N-ENTRY V-INITIAL VALUE E-EQUIVALENCED P-PROMOTED X-EXTERNAL SUBPROGRAM F-STATEMENT FUNCTION Q-PADDED Y-DYNAMIC COMMON G-GENERIC NAME R-SUBPROGRAM NAME Z-EXTENDED COMMON 0NAME TYPE TAG DECLARED REFS (F:REFD S:SET B:REFD/MAY BE SET) +______ ____ ______ ________ ________________________________________ 0LL I*4 D 1 2 M I*4 D 1 11F 12F 15F 16F 18F 20F 27F 37F 44F 47F 54F 54F 57B 60F 62F 62F 62F 62F 66F 66F 66F 67F 67F 67F 67F 70F 72F MATEQ R 1 M1 I*4 11S 21F N I*4 D 1 2 57B NN I*4 47S 48F 48F 48B NRS I*4 5S 6F 14F 20S 49F 50F 59F 61F 69F NRS1 I*4 D 1 5F N1 I*4 60S 63F P R*8 ADT 1 2 17S 19S 51F 51F 51S 62F 62S 66F 67F 67S 71F 73S 0LABEL CROSS REFERENCE DICTIONARY 0TAGS: A-USED AS ARGUMENT F-FORMAT S-USED IN ASSIGN STATEMENT B-OBJECT OF BRANCH N-NON-EXECUTABLE 0 LABEL TAG DEFINED REFERENCED +____________ ___ _______ ___________________________________________ 0 100 13 12 200 B 15 14 300 17 15 16 400 19 18 500 B 21 14 14 600 B 32 27 28 700 40 37 800 B 42 33 900 48 47 1000 51 50 1100 B 52 44 45 49 1200 53 21 1300 60 UNREFERENCED 1400 66 65 1500 67 63 1600 68 61 1700 71 70 1800 73 69 72 1900 B 75 42 54 58 2000 B 77 59 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAY 11, 1992 19:04:57 NAME:MATEQ PAGE: 24 0*STATISTICS* SOURCE STATEMENTS: 75, PROGRAM SIZE: 3564 BYTES, PROGRAM NAME: MATEQ, PAGE: 18 *STATISTICS* NO DIAGNOSTICS GENERATED. **MATEQ** END OF COMPILATION 2 ****** TIME STAMP: 92.13219.04.57 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAY 11, 1992 19:04:58 PAGE: 25 0OPTIONS IN EFFECT: NOLIST NOMAP XREF GOSTMT NODECK SOURCE TERM OBJECT FIXED TRMFLG SRCFLG NOSYM NORENT SDUMP(ISN) SXM NOVECTOR IL(DIM) NOTEST SC(*) NODC NOEC NOEMODE NOICA NODIRECTIVE NODBCS NOSAA NOPARALLEL NOSAVE NOTABS OPT(2) LANGLVL(77) NOFIPS FLAG(I) AUTODBL(NONE) LINECOUNT(60) CHARLEN(500) 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 1 SUBROUTINE MATDT(A,IA,N,DET,IV,IDET) C SUBROUTINE MATDT IS DMATDT FROM THE UICC SUBROUTINE LIBRARY. 2 REAL*8 A(IA,1),DET,B,LOG16 3 INTEGER*4 IV(1),K 4 EQUIVALENCE (B,K) 5 NUM=16777216 6 LOG16=.120411998265592457D+01 7 IF (A(IV(N),N).EQ.0.0D+00) GO TO 300 8 L=0 9 DO 100 I=1,N 1 10 B=DABS(A(IV(I),I)) 1 11 K=K/NUM-64 1 12 L=L+K 1 13 100 DET=DET*(A(IV(I),I)/16.0D+00**K) 14 B=DABS(DET) 15 K=K/NUM-64 16 IW=L+K 17 IF ((IW.LT.-64).OR.(IW.GT.63)) GO TO 200 18 DET=DET*16.0D+00**L 19 IDET=0 20 GO TO 400 21 200 DET=DET*16.0D+00**(-K) 22 IDET=L+K 23 B=IDET*LOG16 24 IDET=B 25 B=B-DFLOAT(IDET) 26 DET=DET*1.0D+01**B 27 GO TO 400 28 300 DET=0.0D+00 29 IDET=0 30 400 RETURN 31 END 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAY 11, 1992 19:04:58 NAME:MATDT PAGE: 26 0SYMBOL CROSS REFERENCE DICTIONARY 0PROGRAM NAME: MATDT 0TAGS: A-ARRAY I-INTRINSIC FUNCTION S-ASSIGNED C-COMMON K-NAMED CONSTANT T-EXPLICITLY TYPED D-DUMMY ARGUMENT N-ENTRY V-INITIAL VALUE E-EQUIVALENCED P-PROMOTED X-EXTERNAL SUBPROGRAM F-STATEMENT FUNCTION Q-PADDED Y-DYNAMIC COMMON G-GENERIC NAME R-SUBPROGRAM NAME Z-EXTENDED COMMON 0NAME TYPE TAG DECLARED REFS (F:REFD S:SET B:REFD/MAY BE SET) +______ ____ ______ ________ ________________________________________ 0A R*8 ADT 1 2 7F 10F 13F B R*8 ET 2 4 10S 14S 23S 24F 25F 25S 26F DABS R*8 I 10 14 DET R*8 DT 1 2 13F 13S 14F 18F 18S 21F 21S 26F 26S 28S DFLOAT R*8 I 25 I I*4 9S 10F 10F 13F 13B IA I*4 D 1 2 IDET I*4 D 1 19S 22S 23F 24S 25F 29S IV I*4 ADT 1 3 7F 10F 13F IW I*4 16S 17F 17F K I*4 ET 3 4 11F 11S 12F 13F 15F 15S 16F 21F 22F L I*4 8S 12F 12S 16F 18F 22F LOG16 R*8 T 2 6S 23F MATDT R 1 N I*4 D 1 7F 7F 9F NUM I*4 5S 11F 15F 0LABEL CROSS REFERENCE DICTIONARY 0TAGS: A-USED AS ARGUMENT F-FORMAT S-USED IN ASSIGN STATEMENT B-OBJECT OF BRANCH N-NON-EXECUTABLE 0 LABEL TAG DEFINED REFERENCED +____________ ___ _______ ___________________________________________ 0 100 13 9 200 B 21 17 300 B 28 7 400 B 30 20 27 0*STATISTICS* SOURCE STATEMENTS: 31, PROGRAM SIZE: 1552 BYTES, PROGRAM NAME: MATDT, PAGE: 25 *STATISTICS* NO DIAGNOSTICS GENERATED. **MATDT** END OF COMPILATION 3 ****** TIME STAMP: 92.13219.04.58 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAY 11, 1992 19:04:58 PAGE: 27 0SUMMARY OF MESSAGES AND STATISTICS FOR ALL COMPILATIONS 0*STATISTICS* SOURCE STATEMENTS: 312, PROGRAM SIZE: 663824 BYTES, PROGRAM NAME: MAIN#, PAGE: 1 *STATISTICS* NO DIAGNOSTICS GENERATED. **MAIN#** END OF COMPILATION 1 ****** TIME STAMP: 92.13219.04.58 0*STATISTICS* SOURCE STATEMENTS: 75, PROGRAM SIZE: 3564 BYTES, PROGRAM NAME: MATEQ, PAGE: 18 *STATISTICS* NO DIAGNOSTICS GENERATED. **MATEQ** END OF COMPILATION 2 ****** TIME STAMP: 92.13219.04.58 0*STATISTICS* SOURCE STATEMENTS: 31, PROGRAM SIZE: 1552 BYTES, PROGRAM NAME: MATDT, PAGE: 25 *STATISTICS* NO DIAGNOSTICS GENERATED. **MATDT** END OF COMPILATION 3 ****** TIME STAMP: 92.13219.04.58 0******* SUMMARY STATISTICS ******* 0 DIAGNOSTICS GENERATED. HIGHEST SEVERITY CODE IS 0.