1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN FEB 15, 1994 23:20:25 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 CLUSTER ANALYSIS OF CASES BY MIXTURE MODEL C C C C PROGRAM MIXPDTA CLUSPAC: C C C C VARYING COVARIANCE MATRICES C C AUTOMATIC SETTING OF INITIAL PARAMETER ESTIMATES C C C C Programs are organized as C C programs within packages within libraries. C C PROGRAM: MIXPDTA PACKAGE: CLUSPAC LIBRARY: ISOPAC C C C C DOCUMENTATION (on my public disk for clustering programs): C C DIRETABL CLUSPAC, CLUSPAC DOC, MIXPDTA DOC C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C CMS FILE: MIXPDTA CLUSPAC C C VERSION 2.98 15-Feb-94 C C C C PROGRAMMED BY: C C C C DR. STANLEY L. SCLOVE 312/996-2681 C C INFORMATION AND 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-4348 C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C C C ISOPAC consists of the CLUSPAC clustering programs, C C the TSPAC time-series segmentation programs and the C C IMPAC image-segmentation programs. 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 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN FEB 15, 1994 23:20:25 NAME:MAIN# PAGE: 2 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 C model likelihood. C C C C Programs MIXP** in the CLUSPAC package are C C mixture-model programs for clustering multivariate data. C C (For univariate data the "MIX1**" programs may be used.) C C MIXPCM and MIXPCMA assume a common covariance matrix C C across distributions. MIXPDT and MIXPDTA allow different C C covariance matrices. C C MIXPCMA and MIXPDTA try a range of numbers of clusters, C C with automatic setting of initial values of the parameters. C C C C C C RESTRICTIONS (CAN BE MODIFIED): C C C C N, SAMPLE SIZE, AT MOST 9999; 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, 30. C C C C SUBROUTINE(S) CALLED: C C MATEQ, WHICH CALLS MATDT C C C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C 1 DIMENSION X(9999,20),SUM(29,20),ID(9999) C DIMENSION Y(N=IP+1,IP) 2 DIMENSION Y(21,20) 3 DIMENSION F(9999,29),PP(29,9999) 4 DIMENSION ICLUS(9999),ICLSOL(9999) 5 DIMENSION DENOM(9999),XMXPR(9999) 6 DIMENSION SSD(29,20,20),SIGMA(29,20,20) 7 DIMENSION PREC(29,20,20),CMPREC(20,20) 8 DIMENSION DSQ(29) 9 DIMENSION TITLE(18) 10 DIMENSION NG(29),XBAR(29,20) 11 DIMENSION DATFMT(18),INPFMT(18) C DIMENSION XIDET(29),XLGDET(29) 12 DIMENSION TRUDET(29) C 13 DIMENSION WGSS(20,20) 14 DIMENSION VARHAT(20,20),COVMX(20,20) 15 DIMENSION XMIN(20),XMAX(20) 16 DIMENSION XJ(9,9) 17 DIMENSION PR(29),XM2LL(9) 18 DIMENSION AICVEC(9),SCHVEC(9),AICPPH(9),SCHPPH(9),IPARM(9) 19 DIMENSION PPOLD(29,9999) 20 DIMENSION PRIRLN(9) 21 DIMENSION TOTAL(29) 22 DIMENSION XMEAN(29) 23 DIMENSION SUMSQS(20,20) 24 DIMENSION SSDEVS(20,20) 25 DIMENSION P(20,20),A(20,20) 26 DIMENSION IV(20,20) C 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN FEB 15, 1994 23:20:25 NAME:MAIN# PAGE: 3 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 C IV IS A WORK ARRAY FOR SUBROUTINE MATEQ. C 27 CHARACTER*6 ID C 28 DOUBLE PRECISION SUM 29 DOUBLE PRECISION WGSS,SSD,TERM,SSDEVS,SUMSQS,TOTAL 30 DOUBLE PRECISION DEVV,DEVW,PREC,CMPREC 31 DOUBLE PRECISION VARHAT,COVMX,SDJV,SDJW 32 DOUBLE PRECISION P,A,SIGMA 33 DOUBLE PRECISION DET C DOUBLE PRECISION XLGDET 34 DOUBLE PRECISION DSQ,ZSQ 35 DOUBLE PRECISION XBAR 36 DOUBLE PRECISION F,PR 37 DOUBLE PRECISION SUMPXF, SUMLNF 38 DOUBLE PRECISION SCHPPH,AICPPH,SPPHA,SPPHS,PP 39 DOUBLE PRECISION DENOM 40 DOUBLE PRECISION PPOLD 41 DOUBLE PRECISION PRIRLN 42 DOUBLE PRECISION TRUDET,COMDET 43 DOUBLE PRECISION TEST C DOUBLE PRECISION XLIKLY C 44 DOUBLE PRECISION PI C C C C ............................................................... C C CONTROL CARDS: C C DATASET TITLE C N, IN FORMAT (2X,I4) C IP, IN FORMAT (3X,I2) C INPFMT, IN FORMAT (18A4), E.G., (4F4.1,1X,A6) C DATFMT, IN FORMAT (18A4), E.G., (4F4.1) C "DATFMT" WILL ALSO BE USED FOR OUTPUT, SO ALLOW C AT LEAST ONE BLANK AT THE BEGINNING FOR CARRIAGE CONTROL. C DATA, ONE CASE AT A TIME, IN FORMAT SPECIFIED BY DATFMT C ............................................................... C C INPFMT allows for a last variable which is an alphanumeric C case label. C ............................................................... C C 45 READ (5,15000) TITLE C C WRITE PROGRAM INFORMATION. 46 WRITE (6,10050) 47 WRITE (6,10000) 48 WRITE (6,10050) 49 WRITE (6,10999) 50 WRITE (6,11001) TITLE C C THE FOLLOWING ARE IMSL DATE AND TIME FUNCTIONS: 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN FEB 15, 1994 23:20:25 NAME:MAIN# PAGE: 4 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 C CALL TDATE(IDAY,MONTH,IYEAR) C CALL DTIME(IHOUR,MINUTE,ISEC) C WRITE (6,11100) MONTH,IDAY,IYEAR,IHOUR,MINUTE,ISEC C C READ SAMPLE SIZE, N. 51 READ (5,12000) N C READ NUMBER OF VARIABLES, IP. 52 READ (5,54000) IP 53 WRITE (6,56000) IP 54 WRITE (6,13000) N C C READ INPUT FORMAT (INPUT INCLUDES DATA AND ID.) 55 READ (5,15000) INPFMT C READ DATA FORMAT. 56 READ (5,15000) DATFMT C C READ DATA. 57 DO 500 I = 1,N 1 58 READ (5,INPFMT) (X(I,JV), JV = 1,IP),ID(I) 1 59 IF (I .EQ. 1) GO TO 100 1 60 GO TO 300 1 61 100 CONTINUE 1 62 DO 200 JV = 1,IP 2 63 XMAX(JV) = X(1,JV) 2 64 XMIN(JV) = X(1,JV) 2 65 200 CONTINUE 1 66 300 CONTINUE 1 67 DO 400 JV = 1,IP 2 68 IF (X(I,JV) .LT. XMIN(JV)) XMIN(JV)=X(I,JV) 2 70 IF (X(I,JV) .GT. XMAX(JV)) XMAX(JV)=X(I,JV) 2 72 400 CONTINUE 1 73 500 CONTINUE 74 WRITE (6,19000) 75 WRITE (6,DATFMT) (XMIN(JV),JV=1,IP) 76 WRITE (6,21000) 77 WRITE (6,DATFMT) (XMAX(JV),JV=1,IP) C 78 WRITE (6,11111) 79 WRITE (6,11112) 80 DO 111 I=1,4 1 81 WRITE (6,INPFMT) (X(I,JV),JV=1,IP),ID(I) C 1 82 KUP = MIN(9,N) 1 83 KREACH = KUP 1 84 111 CONTINUE 85 WRITE (6,11096) KUP C C SET CONSTANTS: C 86 XN = N 87 PI = .314159265358979324D+01 88 NGFLAG = 0 C C SET INITIAL VALUES OF MIXING PROBABILITIES:-- C C XJ(IC,K), K=2 TO 9, IC = 1 TO K 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN FEB 15, 1994 23:20:25 NAME:MAIN# PAGE: 5 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 C OPTIMAL CLASS PROBABILITIES - NORMAL DISTRIBUTION C JOHARI & SCLOVE (1975). "PARTITIONING A DISTRIBUTION." C COMMUNICATIONS IN STATISTICS. C 89 XJ(1,2)=.5 90 XJ(2,2)=.5 91 XJ(1,3)=.2703 92 XJ(2,3)=.4594 93 XJ(3,3)=.2703 94 XJ(1,4)=.1631 95 XJ(2,4)=.3369 96 XJ(3,4)=.3369 97 XJ(4,4)=.1631 98 XJ(1,5)=.1068 99 XJ(2,5)=.2444 100 XJ(3,5)=.2976 101 XJ(4,5)=.2444 102 XJ(5,5)=.1068 103 XJ(1,6)=.0739 104 XJ(2,6)=.1810 105 XJ(3,6)=.2451 106 XJ(4,6)=.2451 107 XJ(5,6)=.1810 108 XJ(6,6)=.0739 109 XJ(1,7)=.0536 110 XJ(2,7)=.1375 111 XJ(3,7)=.1986 112 XJ(4,7)=.2106 113 XJ(5,7)=.1986 114 XJ(6,7)=.1375 115 XJ(7,7)=.0536 116 XJ(1,8)=.0402 117 XJ(2,8)=.1067 118 XJ(3,8)=.1613 119 XJ(4,8)=.1918 120 XJ(5,8)=.1918 121 XJ(6,8)=.1613 122 XJ(7,8)=.1067 123 XJ(8,8)=.0402 124 XJ(1,9)=.0310 125 XJ(2,9)=.0845 126 XJ(3,9)=.1324 127 XJ(4,9)=.1643 128 XJ(5,9)=.1756 129 XJ(6,9)=.1643 130 XJ(7,9)=.1324 131 XJ(8,9)=.0845 132 XJ(9,9)=.0310 C C EVALUATE MODEL-SELECTION CRITERIA FOR K=1 (NO CLUSTERING): C C COMPUTE SUM VECTOR AND SUM-OF-PRODUCTS MATRIX: 133 DO 105 JV = 1,IP 1 134 TOTAL(JV) = 0.0 1 135 DO 105 JW = 1,IP 2 136 SUMSQS(JV,JW) = 0.0 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN FEB 15, 1994 23:20:25 NAME:MAIN# PAGE: 6 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 2 137 105 CONTINUE 138 DO 110 JV = 1,IP 1 139 DO 110 I = 1,N 2 140 TOTAL(JV) = TOTAL(JV) + X(I,JV) 2 141 110 CONTINUE 142 DO 120 JV = 1,IP 1 143 DO 120 JW = 1,IP 2 144 DO 120 I = 1,N 3 145 SUMSQS(JV,JW) = SUMSQS(JV,JW) + X(I,JV)*X(I,JW) 3 146 120 CONTINUE 147 DO 130 JV = 1,IP 1 148 XMEAN(JV) = TOTAL(JV)/XN 1 149 130 CONTINUE 150 DO 140 JV = 1,IP 1 151 DO 140 JW = 1,IP 2 152 SSDEVS(JV,JW) = SUMSQS(JV,JW)- TOTAL(JV)*TOTAL(JW)/XN 2 153 VARHAT(JV,JW) = SSDEVS(JV,JW)/XN 2 154 140 CONTINUE C WRITE SUMMARY STATISTICS FOR WHOLE SAMPLE: 155 WRITE(6,10500) 156 WRITE(6,10100) ( XMEAN(JV), JV = 1,IP ) 157 WRITE(6,10200) 158 DO 150 JV = 1,IP 1 159 WRITE(6,48000) ( VARHAT(JV,JW), JW=1,IP ) 1 160 150 CONTINUE C C C CALL SUBROUTINE TO COMPUTE DETERMINANT OF SAMPLE COVARIANCE C MATRIX: C 161 IDET = 1 162 NRS1 = 0 C C COPY VARHAT INTO COVMX TO BE USED AS MATRIX INPUT TO MATEQ C BECAUSE MATRIX INPUT IS DESTROYED BY MATEQ AND VARHAT WILL C BE NEEDED LATER: 163 DO 2402 JV = 1,IP 1 164 DO 2402 JW = 1,IP 2 165 COVMX(JV,JW) = VARHAT(JV,JW) 2 166 2402 CONTINUE 167 CALL MATEQ(COVMX, IP,20,JFLG,DET,IDET,IV,NRS1,P,20) C C GENERAL FORM OF CALL IS: C CALL MATEQ(A,M,N,JFLG,DET,IDET,IV,NRS1,P,LL) C C THE INPUT MATRIX A IS DESTROYED IN THE COMPUTATION PROCESS. C ON RETURN, P CONTAINS THE INVERSE OF A. C SEE SUBROUTINE LISTING FOR FULLER EXPLANATION. C DET(VARHAT) = DET*10.0**IDET C 168 IF (JFLG .GT. 0) WRITE (6,60000) JFLG C C C XIDET(1) = IDET C XLGDET(1) = DLOG(DET) + XIDET(1)*DLOG(10.0D+00) 170 TRUDET(1) = DET*(10.0**IDET) 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN FEB 15, 1994 23:20:25 NAME:MAIN# PAGE: 7 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 C WRITE (6,62000) DET,IDET C C MAX LIKELIHOOD = C (2*PI)**(-N*P/2)*(DET OF COV MX)**(-N/2)*EXP(-N*P/2) C C COMPUTE LOG OF MAX LIKELIHOOD: C 171 XLL=-(XN*IP/2.0)*DLOG(2.0*PI)-(XN/2.0)*DLOG(TRUDET(1))-XN*IP/2.0 172 XMN2LL = (-2.0)*XLL 173 XM2LL(1) = XMN2LL 174 WRITE(6,10900) XMN2LL C NO. OF PARAMETERS FOR UNCLUSTERED SAMPLE IS: C 1 MEAN VECTOR + 1 COV MX. C 175 NOPARM = IP + IP*(IP+1)/2 176 IPARM(1) = NOPARM C 177 WRITE (6,72000) NOPARM 178 AIC = XMN2LL + 2.0*NOPARM 179 WRITE(6,70000) AIC 180 SCH = XMN2LL + ALOG(XN)*NOPARM 181 WRITE (6,78000) SCH 182 AICVEC(1) = AIC 183 SCHVEC(1) = SCH C 184 WRITE (6,10050) C C C PERFORM CLUSTERING FOR VARIOUS NUMBERS OF CLUSTERS, K: C 185 DO 850 K = 2,KUP 1 186 WRITE (6,26000) K C C C SET CONSTANTS FOR THIS VALUE OF K: C C COUNT NUMBER OF INDEPENDENT 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 C MIXING PROBABILITIES. C 1 187 NOPARM = K*IP + K*IP*(IP+1)/2 + (K-1) 1 188 IPARM(K) = NOPARM C 1 189 WRITE (6,72000) NOPARM C C SET INITIAL VALUES OF PARAMETER ESTIMATES FOR THIS VALUE OF K: C C SET INITIAL VALUES OF PRIOR PROBABILITIES EQUAL TO THE C OPTIMAL VALUES FOR A NORMAL DISTRIBUTION:-- C 1 190 DO 360 IC = 1,K 2 191 PR(IC) = XJ(IC,K) 2 192 360 CONTINUE C C TO SET INITIAL VALUES OF PRIOR PROBABILITIES EQUAL TO 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN FEB 15, 1994 23:20:25 NAME:MAIN# PAGE: 8 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 C BINOMIAL(X=IC-1;N=K-1,P=1/2). C FOR IC = 1,2,...,K, C REMOVE THE "C" FROM CC1 OF THE NEXT LINE: C PR(IC) = GAMMA(XK)/(GAMMA(XC)*GAMMA(XK-XC+1)*(2**K)) C C FOR EQUAL INITIAL VALUES OF PRIOR PROBABILITIES, C REMOVE THE "C" FROM CC1 OF THE NEXT LINE: C PR(IC) = 1.0/K . C C INITIALIZE POSTERIOR PROBABILITIES: 1 193 DO 549 I = 1,N 2 194 DO 549 IC = 1,K 3 195 PPOLD(IC,I) = 1.0/K 3 196 549 CONTINUE C 1 197 WRITE (6,20000) C C Set initial means equal to K of the observations:-- C 1 198 DO 550 IG = 1,K C 1st observation for IG=1, N-th observation for IG=K: 2 199 I = INT( 1 + (IG-1)*(N-1)/(K-1) ) C For more centered means, use instead: C I = INT( N/(2*K) ) + INT( (IG-1)*(N/K) ) 2 200 DO 550 JV = 1,IP 3 201 XBAR(IG,JV) = X(I,JV) 3 202 550 CONTINUE C C EXAMPLE: If N=150 and K=3, the initial means will be C cases 26, 76 AND 126. C C This initialization scheme may be good when the data are C ordered with respect to an important variable. C C C ALTERNATIVELY, HERE IS ANOTHER SCHEME FOR COMPUTING INITIAL MEANS: C (IN 2 DIMENSIONS THIS SCHEME RESULTS IN MEANS ORIENTED FROM C SOUTHWEST TO NORTHEAST): C DO 601 JV = 1,IP C DO 601 IG = 1,K C W = IG/(K+1.0) C XBAR(IG,JV) = (1-W)*XMIN(JV) + W*XMAX(JV) C 601 CONTINUE C C 1 203 DO 600 IG=1,K 2 204 WRITE (6,21999) 2 205 WRITE (6,22000) IG 2 206 WRITE (6,DATFMT) ( XBAR(IG,JV), JV=1,IP ) 2 207 600 CONTINUE C C SET INITIAL COVARIANCE MATRICES:-- C C First recall how the initial were C set equal to K of the observations:-- C 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN FEB 15, 1994 23:20:25 NAME:MAIN# PAGE: 9 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 C DO 550 IG = 1,K C I = INT( N/(2*K) ) + INT( (IG-1)*(N/K) ) C DO 550 JV = 1,IP C XBAR(IG,JV) = X(I,JV) C 550 CONTINUE C C The K initial covariance matrices will be set C using K blocks of P+1 observations. C Each block starts with one of the seed points. C C COMPUTE SUM VECTOR AND SUM-OF-PRODUCTS MATRIX: C C BY-PASS THIS METHOD OF INITIALIZING COV MCS: C 1 208 GO TO 113 C DO 142 IC = 1,K C C Get X's for this value of IC: C Store them as vectors Y(1),Y(2),...,Y(P),Y(P+1): C C Compute first subscript for this block of P+1: C ICN = INT( N/(2*K) ) + INT( (IC-1)*(N/K) ) C DO 107 I=1,IP+1 C DO 107 JV = 1,IP C Y(I,JV) = X(ICN+I-1,JV) C 107 CONTINUE C DO 106 JV = 1,IP C TOTAL(JV) = 0.0 C DO 106 JW = 1,IP C SUMSQS(JV,JW) = 0.0 C 106 CONTINUE C DO 112 JV = 1,IP C DO 112 I = 1,N C TOTAL(JV) = TOTAL(JV) + Y(I,JV) C 112 CONTINUE C DO 121 JV = 1,IP C DO 121 JW = 1,IP C DO 121 I = 1,N C SUMSQS(JV,JW) = SUMSQS(JV,JW) + Y(I,JV)*Y(I,JW) C 121 CONTINUE C DO 131 JV = 1,IP C XMEAN(JV) = TOTAL(JV)/XN C 131 CONTINUE C DO 141 JV = 1,IP C DO 141 JW = 1,IP C SSDEVS(JV,JW) = SUMSQS(JV,JW)- TOTAL(JV)*TOTAL(JW)/XN C VARHAT(JV,JW) = SSDEVS(JV,JW)/XN C SIGMA(IC,JV,JW) = VARHAT(JV,JW) C 141 CONTINUE C C 142 CONTINUE C C C C Here is an alternative way of setting initial cov. mcs.; C some correlation could be programmed in in stmt no. 114: 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN FEB 15, 1994 23:20:25 NAME:MAIN# PAGE: 10 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 C 1 209 113 CONTINUE C 1 210 DO 630 IC = 1,K 2 211 DO 610 JV=1,IP 3 212 DO 610 JW=1,IP 4 213 114 SIGMA(IC,JV,JW) = 0.00 4 214 610 CONTINUE 2 215 DO 620 JV = 1,IP 3 216 SIGMA(IC,JV,JV) = 1.0 3 217 620 CONTINUE 2 218 630 CONTINUE C SCALE UP THIS CORRELATION MATRIX TO A COVARIANCE MATRIX C USING THE WHOLE-SAMPLE VARIANCES: 1 219 DO 611 JV=1,IP 2 220 SDJV = DSQRT(VARHAT(JV,JV)) 2 221 DO 611 JW=1,IP 3 222 SDJW = DSQRT(VARHAT(JW,JW)) 3 223 DO 611 IC = 1,K 4 224 SIGMA(IC,JV,JW) = SDJV*SIGMA(IC,JV,JW)*SDJW 4 225 611 CONTINUE C C COMPUTE INITIAL INVERSE COVARIANCE MATRIX: C C 1 226 DO 640 IC=1,K C C 2 227 IDET = 1 2 228 NRS1 = 0 C C COPY SIGMA(IC,.,.) INTO COVMX TO BE USED AS MATRIX INPUT TO MATEQ C BECAUSE MATRIX INPUT IS DESTROYED BY MATEQ AND SIGMA(IC,.,.) C WILL BE NEEDED LATER: 2 229 DO 2404 JV = 1,IP 3 230 DO 2404 JW = 1,IP 4 231 COVMX(JV,JW) = SIGMA(IC,JV,JW) 4 232 2404 CONTINUE C 2 233 CALL MATEQ(COVMX, IP,20,JFLG,DET,IDET,IV,NRS1,P,20) 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 2 234 IF (JFLG .GT. 0) WRITE (6,60000) JFLG C C C SET INITIAL PRECISION MATRICES FOR THIS K: 2 236 DO 640 JV = 1,IP 3 237 DO 640 JW = 1,IP 4 238 PREC(IC,JV,JW) = P(JV,JW) 4 239 640 CONTINUE C WRITE(6,11500) DET,IDET 1 240 DO 680 IC = 1,K 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN FEB 15, 1994 23:20:25 NAME:MAIN# PAGE: 11 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 C XIDET(IC) = IDET C XLGDET(IC) = DLOG(DET) + XIDET(IC)*DLOG(10.0D+00) 2 241 TRUDET(IC) = DET*(10.0**IDET) 2 242 680 CONTINUE C END PARAMETER INITIALIZATION FOR THIS K; C BEGIN ITERATIONS: C C 1 243 ITER = 1 1 244 700 CONTINUE 1 245 WRITE(6,32000) ITER 1 246 IF (ITER .EQ. 1) GO TO 901 1 247 DO 800 I = 1,N 2 248 ICLSOL(I) = ICLUS(I) 2 249 800 CONTINUE C 1 250 901 CONTINUE 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 1 251 DO 1200 I = 1,N 2 252 DO 1000 L=1,K 3 253 DSQ(L) = 0.0D+00 3 254 DO 1310 JV=1,IP 4 255 DEVV = XBAR(L,JV) - X(I,JV) 4 256 DO 1310 JW=1,IP 5 257 DEVW = XBAR(L,JW) - X(I,JW) 5 258 DSQ(L) = DSQ(L) + DEVV*PREC(L,JV,JW)*DEVW 5 259 1310 CONTINUE 3 260 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): 3 261 IF ( ZSQ/2.0 .LE. 174.673D+00 ) GO TO 1090 3 262 F(I, L) = 0.0D+00 3 263 GO TO 1100 C 3 264 1090 CONTINUE C Compute C F(I,L) = DEXP(-ZSQ/2.0)/((2*PI)**(IP/2)*DSQRT(TRUDET(L))): C 3 265 XIP = IP 3 266 F(I,L) = (-ZSQ/2.0)-(XIP/2.0)*DLOG(2*PI)-0.5*DLOG(TRUDET(L)) 3 267 IF(F(I,L) .GT. 174.673) WRITE(6,87000) ZSQ,TRUDET(L) 3 269 F(I,L) = DEXP( F(I,L) ) 3 270 1100 CONTINUE 3 271 1000 CONTINUE 2 272 1200 CONTINUE C C COMPUTE LOG LIKELIHOOD AND VALUES OF MODEL-SELECTION CRITERIA: C C XLIKLY = 1.0 1 273 SUMLNF = 0.0 1 274 DO 2200 I=1,N 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN FEB 15, 1994 23:20:25 NAME:MAIN# PAGE: 12 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 2 275 SUMPXF = 0.0 2 276 DO 2100 IC=1,K 3 277 SUMPXF = SUMPXF + PR(IC)*F(I,IC) 3 278 2100 CONTINUE C C XLIKLY = XLIKLY*SUMPXF 2 279 IF ( SUMPXF .LE. 0.0 ) GO TO 2190 2 280 SUMLNF = SUMLNF + DLOG(SUMPXF) 2 281 2190 CONTINUE 2 282 2200 CONTINUE C C SUMLNF = DLOG(XLIKLY) 1 283 XMN2LL = -2.0*SUMLNF 1 284 XM2LL(K) = XMN2LL C 1 285 WRITE (6,10900) XMN2LL C C C COMPUTE MODEL SELECTION CRITERIA 1 286 AIC = XMN2LL + 2.0*NOPARM 1 287 SCH = XMN2LL + ALOG(XN)*NOPARM 1 288 WRITE (6,70000) AIC 1 289 WRITE (6,78000) SCH 1 290 AICVEC(K) = AIC 1 291 SCHVEC(K) = 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): 1 292 DO 1350 I = 1,N 2 293 DENOM(I) = 0.0 2 294 DO 1350 IC=1,K 3 295 DENOM(I) = DENOM(I) + PR(IC)*F(I,IC) 3 296 1350 CONTINUE C C IF DENOM(I)=0, PP(IC,I) IS SET EQUAL TO OLD PP(IC,I) SO THAT C THE I-TH CASE WILL STILL HAVE WEIGHT IN SUCCEEDING C COMPUTATIONS: C C 1 297 DO 1400 I = 1,N 2 298 IF ( DENOM(I) .GT. 0.0D+00 ) GO TO 1410 2 299 WRITE (6,14560) ID(I) 2 300 DO 1456 IC = 1,K 3 301 PP(IC,I) = PPOLD(IC,I) 3 302 1456 CONTINUE 2 303 GO TO 1401 C 2 304 1410 CONTINUE 2 305 DO 1402 IC=1,K 3 306 PP(IC,I) = PR(IC)*F(I,IC)/DENOM(I) 3 307 1402 CONTINUE C 2 308 1401 CONTINUE 2 309 1400 CONTINUE 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN FEB 15, 1994 23:20:25 NAME:MAIN# PAGE: 13 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 C 1 310 WRITE (6,55555) 1 311 IF ( N-15 ) 1421,1421,1422 1 312 1422 CONTINUE 1 313 WRITE(6,76000) 1 314 DO 1420 I = 1,4 2 315 WRITE(6,82000) ( PP(IC,I), IC = 1,K ) 2 316 1420 CONTINUE 1 317 GO TO 1423 1 318 1421 CONTINUE 1 319 DO 1424 I = 1,N 2 320 WRITE (6,25000) ID(I), (PP(IC,I),IC=1,K) 2 321 1424 CONTINUE 1 322 1423 CONTINUE C C C UPDATE PARAMETER ESTIMATES: C C C UPDATE CLUSTER PRIOR PROBABILITIES PR(IC):-- 1 323 DO 3800 IC = 1,K 2 324 PR(IC) = 0.0 2 325 DO 3800 I = 1,N 3 326 PR(IC) = PR(IC) + PP(IC,I) 3 327 3800 CONTINUE 1 328 DO 3900 IC = 1,K 2 329 PR(IC) = PR(IC)/N 2 330 3900 CONTINUE C IF MIXING PROBABILITY IS ZERO FOR ANY CLUSTER, STOP. 1 331 DO 3901 IC = 1,K 2 332 IF ( PR(IC) .GT. 0.0D+00 ) GO TO 3901 2 333 WRITE (6,40040) K,IC 2 334 40040 FORMAT(/' IN FITTING ',I1,' CLUSTERS, PR(',I1,')=0. STOP.'/) 2 335 GO TO 4004 2 336 3901 CONTINUE 1 337 DO 1750 IG = 1,K 2 338 DO 1750 JV = 1,IP 3 339 SUM(IG,JV) = 0.0 3 340 DO 1750 JW = 1,IP 4 341 SSD(IG,JV,JW) = 0.0 4 342 1750 CONTINUE C C UPDATE MEANS:-- 1 343 DO 1875 IC = 1,K 2 344 DO 1875 JV = 1,IP 3 345 DO 1875 I = 1,N 4 346 SUM(IC,JV) = SUM(IC,JV) + PP(IC,I)*X(I,JV) 4 347 1875 CONTINUE 1 348 DO 1900 IG = 1,K 2 349 TEST = XN*PR(IG) 2 350 IF ( TEST .EQ. 0.0D+00 ) GO TO 2050 C IF (N*PR(IG) .LT. 0.5) GO TO 2050 2 351 GO TO 2150 2 352 2050 WRITE (6,74000) IG 2 353 NGFLAG = 1 C IF EXP.NO. OF OBSERVATIONS FOR A CLUSTER IS LESS THAN 0.5, 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN FEB 15, 1994 23:20:25 NAME:MAIN# PAGE: 14 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 C END COMPUTATION FOR THIS NUMBER OF CLUSTERS C AND GO ON TO NEXT VALUE OF K: 2 354 GO TO 850 2 355 2150 CONTINUE 2 356 DO 1900 JV = 1,IP 3 357 XBAR(IG,JV) = SUM(IG,JV)/(XN*PR(IG)) 3 358 1900 CONTINUE C C 1 359 CALL XUFLOW(0) C 1 360 DO 3600 IC=1,K 2 361 DO 3600 JV=1,IP 3 362 DO 3600 JW=JV,IP 4 363 DO 3600 I=1,N 5 364 DEVV=X(I,JV)-XBAR(IC,JV) 5 365 DEVW=X(I,JW)-XBAR(IC,JW) 5 366 IF ( PP(IC,I) .EQ. 0.0 ) GO TO 3600 5 367 IF ( DEVV .EQ. 0.0 ) GO TO 3600 5 368 IF ( DEVW .EQ. 0.0 ) GO TO 3600 5 369 SSD(IC,JV,JW)=SSD(IC,JV,JW) + PP(IC,I)*DEVV*DEVW 5 370 3600 CONTINUE C 1 371 DO 3700 IC=1,K 2 372 DO 3700 JV=2,IP 3 373 DO 3700 JW=1,JV-1 4 374 SSD(IC,JV,JW)=SSD(IC,JW,JV) 4 375 3700 CONTINUE C 1 376 DO 3710 IC = 1,K 2 377 DO 3710 JV = 1,IP 3 378 DO 3710 JW = 1,IP 4 379 SIGMA(IC,JV,JW) = SSD(IC,JV,JW)/(N*PR(IC)) 4 380 3710 CONTINUE C C C POOL: C 1 381 DO 1950 JV = 1,IP 2 382 DO 1950 JW = 1,IP 3 383 WGSS(JV,JW) = 0.0 3 384 1950 CONTINUE C 1 385 DO 2300 JV = 1,IP 2 386 DO 2300 JW = 1,IP 3 387 DO 2300 IC = 1,K 4 388 WGSS(JV,JW) = WGSS(JV,JW) + PR(IC)*SSD(IC,JV,JW) 4 389 2300 CONTINUE C COMPUTE VARHAT, MLE OF COMMON COVARIANCE MATRIX: 1 390 DO 2400 JV = 1,IP 2 391 DO 2400 JW = 1,IP 3 392 VARHAT(JV,JW) = WGSS(JV,JW)/N 3 393 2400 CONTINUE C C UPDATE PRECISION MATRICES: C 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN FEB 15, 1994 23:20:25 NAME:MAIN# PAGE: 15 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 1 394 NRS1 = 0 1 395 IDET = 1 C C COPY VARHAT INTO COVMX TO BE USED AS MATRIX INPUT TO MATEQ C BECAUSE MATRIX INPUT IS DESTROYED BY MATEQ AND VARHAT WILL C BE NEEDED LATER: 1 396 DO 2406 JV = 1,IP 2 397 DO 2406 JW = 1,IP 3 398 COVMX(JV,JW) = VARHAT(JV,JW) 3 399 2406 CONTINUE C COMPUTE INVERSE OF POOLED COVARIANCE MATRIX:-- 1 400 CALL MATEQ(COVMX, IP,20,JFLG,DET,IDET,IV,NRS1,P,20) 1 401 DO 3813 JV = 1,IP 2 402 DO 3813 JW = 1,IP 3 403 CMPREC(JV, JW)=P(JV,JW) 3 404 3813 CONTINUE 1 405 IF (JFLG .GT. 0) WRITE (6,60000) JFLG 1 407 IF (JFLG .GT. 0) WRITE (6,62000) DET,IDET 1 409 COMDET = DET*(10.0**IDET) C C C UPDATE PRECISION MATRICES: C 1 410 DO 3711 IC = 1,K 2 411 DO 3712 JV = 1,IP 3 412 DO 3712 JW = 1,IP 4 413 A(JV,JW) = SIGMA(IC,JV,JW) 4 414 3712 CONTINUE C 2 415 IDET = 1 2 416 NRS1 = 0 C 2 417 CALL MATEQ(A,IP,20,JFLG,DET,IDET,IV,NRS1,P,20) C 2 418 DO 3713 JV = 1,IP 3 419 DO 3713 JW = 1,IP 4 420 PREC(IC,JV,JW)=P(JV,JW) 4 421 3713 CONTINUE C WRITE(6,11500) DET,IDET 2 422 IF (JFLG .GT. 0) WRITE (6,60000) JFLG 2 424 IF (JFLG .GT. 0) WRITE (6,62000) DET,IDET C XIDET(IC) = IDET C XLGDET(IC) = DLOG(DET) + XIDET(IC)*DLOG(10.0D+00) 2 426 TRUDET(IC) = DET*(10.0**IDET) C IF DET=0, REPLACE COV MX WITH COMMON COV MX: 2 427 IF ( TRUDET(IC) .GT. 0.0D+00 ) GO TO 3711 2 428 TRUDET(IC) = COMDET 2 429 WRITE (6, 17000) IC 2 430 DO 3814 JV = 1,IP 3 431 DO 3814 JW = 1,IP 4 432 PREC(IC,JV,JW) = CMPREC(JV,JW) 4 433 3814 CONTINUE 2 434 3711 CONTINUE C C C (END PARAMETER-ESTIMATE UPDATE SEQUENCE) 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN FEB 15, 1994 23:20:25 NAME:MAIN# PAGE: 16 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 C C 1 435 IF (ITER .EQ. 1) GO TO 900 C STORE OLD LABELS: 1 436 DO 825 I = 1,N 2 437 ICLSOL(I) = ICLUS(I) 2 438 825 CONTINUE C 1 439 900 CONTINUE C C COMPUTE NEW LABELS BY MAX POSTERIOR PROBABILITY: C 1 440 DO 1600 I = 1,N 2 441 XMXPR(I) = PP(1,I) 2 442 ICLUS(I) = 1 2 443 DO 1600 IC = 2,K 3 444 IF ( PP(IC,I) .GT. XMXPR(I) ) GO TO 1500 3 445 GO TO 1600 3 446 1500 XMXPR(I) = PP(IC,I) 3 447 ICLUS(I) = IC 3 448 1600 CONTINUE C 1 449 DO 1601 I = 1,N 2 450 DO 1601 IC = 1,K 3 451 PPOLD(IC,I) = 1.0/K 3 452 PPOLD(IC,I) = PP(IC,I) 3 453 1601 CONTINUE C 1 454 IF (N .GE. 31) GO TO 250 1 455 WRITE (6,14000) C 1 456 WRITE (6,16000) 1 457 WRITE (6,18000) (ID(I),ICLUS(I), I=1,N) 1 458 250 CONTINUE C C C 1 459 IF (ITER .EQ. 1) GO TO 3000 C IF NEW CLUSTERING IS SAME AS OLD, NO FURTHER ITERATION IS C NECESSARY AND RESULTS FOR THIS VALUE OF K WILL BE PRINTED. C OTHERWISE, FURTHER ITERATIONS WILL BE PERFORMED. C 1 460 DO 2900 I = 1,N 2 461 IF (ICLUS(I) .EQ. ICLSOL(I)) GO TO 2900 2 462 GO TO 3000 2 463 2900 CONTINUE 1 464 GO TO 3300 1 465 3000 CONTINUE C C COUNT NUMBERS IN CLUSTERS: 1 466 DO 3311 IC = 1,K 2 467 NG(IC) = 0 2 468 3311 CONTINUE 1 469 DO 3321 I = 1,N 2 470 IC = ICLUS(I) 2 471 NG(IC) = NG(IC) + 1 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN FEB 15, 1994 23:20:25 NAME:MAIN# PAGE: 17 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 2 472 3321 CONTINUE 1 473 WRITE (6,36000) (NG(IC),IC=1,K) C 1 474 ITER = ITER + 1 C 1 475 MAXITER = 30 1 476 IF ( ITER .GE. MAXITER+1 ) GO TO 3100 C C Unless the maximum number of iterations has already been C performed, go to stmt no. 700 and begin another: 1 477 GO TO 700 1 478 3100 WRITE (6,88000) C C 1 479 3300 CONTINUE C C 1 480 WRITE (6,40000) (PR(IC),IC=1,K) 1 481 DO 2800 IC = 1,K 2 482 WRITE (6,30000) IC, (XBAR(IC,JV),JV=1,IP) 2 483 2800 CONTINUE C 1 484 WRITE (6,42000) 1 485 DO 2500 JV=1,IP 2 486 WRITE (6,48000) (VARHAT(JV,JW),JW=1,IP) 2 487 2500 CONTINUE C 1 488 IF (N .GE. 31) GO TO 3510 C LIST CLUSTER MEMBERSHIP OF EACH CASE: 1 489 WRITE (6,52000) 1 490 DO 3500 I = 1,N 2 491 WRITE (6,50000) ID(I), ICLUS(I), ( X(I,JV),JV=1,IP ) 2 492 3500 CONTINUE 1 493 3510 CONTINUE C C C 1 494 WRITE (6,58000) ITER C C C C COMPUTE VALUES OF MODEL-SELECTION CRITERIA: C C 1 495 AIC = XMN2LL + 2.0*NOPARM 1 496 SCH = XMN2LL + ALOG(XN)*NOPARM C 1 497 IF ( N .GE. 301 ) GO TO 852 1 498 DO 851 IC=1,K 2 499 WRITE (6,44444) IC 2 500 DO 851 I = 1,N 3 501 IF ( ICLUS(I) .EQ. IC ) WRITE (6,44445) ID(I) 3 503 851 CONTINUE 1 504 852 CONTINUE C C 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN FEB 15, 1994 23:20:25 NAME:MAIN# PAGE: 18 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 C GO ON TO NEXT NUMBER OF CLUSTERS (NEXT VALUE OF K), C IF NOT FINISHED. C 1 505 WRITE (6,10050) C 1 506 850 CONTINUE C 507 4004 CONTINUE 508 KREACH = K-1 C C SCHPPH(K) = PR(H(K))* C CONST*(2*PI)**(IPARM(K)/2.0)*EXP(-SCHVEC(K)/2.0)*PRIOR PDF(K); C HERE PR(H(K)) IS TAKEN TO BE UNIFORM. C AICPPH(K) = PR(H(K))* C CONST*(2*PI)**(IPARM(K)/2.0)*EXP(-SCHVEC(K)/2.0)*PRIOR PDF(K); C HERE PR(H(K)) IS TAKEN TO BE Const*EXP(-NOPARM(K)). 509 SCHPPH(1) = (IPARM(1)/2.0)*DLOG(2.0*PI) - SCHVEC(1)/2.0 510 AICPPH(1) = (IPARM(1)/2.0)*DLOG(2.0*PI) - SCHVEC(1)/2.0 C Add in log of prior Pr[H(k)] for AIC: 511 AICPPH(1) = AICPPH(1) - IPARM(1) C Add in log of prior Pr[H(k)] for Schwarz criterion: 512 SCHPPH(1) = SCHPPH(1) + ALOG(1.0/KUP) C Approximate log of PRIOR PDF: C PRIRLN(1) = -(1*IP/2.0)*DLOG(2*PI) -IP/2.0 C Note that this assumes unit variance: It will have to be C corrected to use a variance more reflective of that in the C observations. C Meantime, set it equal to zero. 513 PRIRLN(1) = 0.0 C Add in log of approximate prior pdf: 514 SCHPPH(1) = SCHPPH(1) + PRIRLN(1) 515 AICPPH(1) = AICPPH(1) + PRIRLN(1) C C ADD CONSTANT SCHVEC(1)/2 OR AICVEC(1)/2 TO PREVENT UNDERFLOW: C (This is okay because the added term cancels out below.) C SCHPPH(1) = SCHPPH(1) + SCHVEC(1)/2 C AICPPH(1) = AICPPH(1) + AICVEC(1)/2 C 516 DO 6110 K = 2,KREACH 1 517 SCHPPH(K) = (IPARM(K)/2.0)*DLOG(2*PI) - SCHVEC(K)/2.0 1 518 AICPPH(K) = (IPARM(K)/2.0)*DLOG(2*PI) - SCHVEC(K)/2.0 C Add in log Pr[H(k)]: 1 519 AICPPH(K) = AICPPH(K) - IPARM(K) 1 520 SCHPPH(K) = SCHPPH(K) + ALOG(1.0/KUP) C Approximate log of PRIOR PDF: C PRIRLN(K) = -(K*IP/2.0)*DLOG(2*PI) -K*IP/2.0 C Note that this assumes unit variance: It will have to be C corrected to use a variance more reflective of that in the C observations. C Meantime, set it equal to zero. 1 521 PRIRLN(K) = 0.0 1 522 AICPPH(K) = AICPPH(K) + PRIRLN(K) 1 523 SCHPPH(K) = SCHPPH(K) + PRIRLN(K) C C ADD CONSTANT SCHVEC(1)/2 OR AICVEC(1)/2 TO PREVENT UNDERFLOW: C SCHPPH(K) = SCHPPH(K) + SCHVEC(1)/2 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN FEB 15, 1994 23:20:25 NAME:MAIN# PAGE: 19 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 C AICPPH(K) = AICPPH(K) + AICVEC(1)/2 C 1 524 6110 CONTINUE C AT THIS POINT, WE HAVE THE LOGS OF THE PPHS; WE NOW CONVERT C THEM TO THE PPHS: 525 SCHPPH(1) = DEXP(SCHPPH(1)) 526 AICPPH(1) = DEXP(AICPPH(1)) 527 DO 6119 K = 2,KREACH C If an argument is less or = -174.673, there will be underflow: 1 528 IF ( SCHPPH(K) .LE. -174. ) GO TO 6115 1 529 SCHPPH(K) = DEXP(SCHPPH(K)) 1 530 GO TO 6119 1 531 6115 SCHPPH(K) = 0.0 1 532 6119 CONTINUE 533 DO 6121 K = 2,KREACH 1 534 IF ( AICPPH(K) .LE. -174.673 ) GO TO 6126 1 535 AICPPH(K) = DEXP(AICPPH(K)) 1 536 GO TO 6121 1 537 6126 AICPPH(K) = 0.0 1 538 6121 CONTINUE 539 SPPHS = 0.0 540 SPPHA = 0.0 541 DO 6120 K = 1,KREACH 1 542 IF ( SCHPPH(K) .GT. 0.0D+00 ) SPPHS = SPPHS + SCHPPH(K) 1 544 IF ( AICPPH(K) .GT. 0.0D+00 ) SPPHA = SPPHA + AICPPH(K) 1 546 6120 CONTINUE 547 DO 6130 K = 1,KREACH 1 548 SCHPPH(K) = SCHPPH(K)/SPPHS 1 549 AICPPH(K) = AICPPH(K)/SPPHA 1 550 6130 CONTINUE C C C WRITE VALUES OF MSC'S FOR VARIOUS K: 551 WRITE (6,11000) TITLE 552 WRITE(6,34000) 553 DO 605 K = 1,KREACH 1 554 WRITE(6,46000) K,XM2LL(K),IPARM(K),AICVEC(K),AICPPH(K), 1 1SCHVEC(K),SCHPPH(K) 1 555 605 CONTINUE C C 556 WRITE (6,10050) 557 WRITE (6,11002) TITLE 558 WRITE (6,96666) 559 96666 FORMAT(' CONDITION CODE FOR THIS RUN: ') 560 IF ( NGFLAG .EQ. 0 ) GO TO 4002 561 WRITE (6,96000) 562 GO TO 4001 563 4002 CONTINUE 564 WRITE (6,98000) 565 4001 CONTINUE 566 WRITE (6,10050) 567 WRITE (6,11097) 568 WRITE (6,86000) 569 WRITE (6,11100) MONTH,IDAY,IYEAR,IHOUR,MINUTE,ISEC 570 CALL DTIME(IHOUR,MINUTE,ISEC) 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN FEB 15, 1994 23:20:25 NAME:MAIN# PAGE: 20 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 571 CALL TDATE(IDAY,MONTH,IYEAR) 572 WRITE (6,89000) 573 WRITE (6,11100) MONTH,IDAY,IYEAR,IHOUR,MINUTE,ISEC 574 STOP C C 575 10000 FORMAT('1','**********************************************'/// X' MIXPDTA CLUSPAC - MIXTURE-MODEL CLUSTERING OF CASES'// X' VARYING COVARIANCE MATRICES '/ X' AUTOMATIC SETTING OF INITIAL PARAMETER ESTIMATES '// X' Copyright 1991, 1992 Stanley Louis Sclove. '// X' Developed and programmed by: '// X19X,'Prof. Stanley L. Sclove, Ph.D. Phone (312) 996-2681'/ X19X,'Information & Decision Sciences Dept. M/C 294 '/ X19X,'College of Business Administration '/ X19X,'University of Illinois at Chicago '/ X19X,'Box 4348, Chicago, IL 60680-4348 '// X' Program Version: 2.98 15-Feb-94 '/ X' (VM/CMS) ') C 576 10050 FORMAT(/' ................................................', X'................................'//) 577 11112 FORMAT(' ----- ---- ---- ------'/) 578 19000 FORMAT( 1X,'MINIMUM OF EACH VARIABLE: ') 579 10100 FORMAT(' MEAN VECTOR: ',(9F8.2/)/) 580 10999 FORMAT( ' PROBLEM TITLE IS' ) 581 10200 FORMAT(/' COVARIANCE MATRIX: '/) 582 10500 FORMAT(/,' STATISTICS FOR WHOLE (UNCLUSTERED) SAMPLE: '/) 583 11001 FORMAT(1X,18A4/) 584 11002 FORMAT(/1X,18A4/) 585 11097 FORMAT(' MIXPDTA CLUSPAC - MIXTURE-MODEL CLUSTERING OF CASES'/) 586 11100 FORMAT(1X,I2,'/',I2,'/',I4,3X,'AT ',I2,':',I2,':',I2/) 587 25000 FORMAT(1X, A6, 9F8.3) 588 55555 FORMAT(/' CASE POSTERIOR PROBS OF CLUSTER MEMBERSHIP: '/) 589 20000 FORMAT(/ 1X,'INITIAL MEANS') 590 44444 FORMAT(' CLUSTER ',I2) 591 44445 FORMAT(1X,A6) 592 10900 FORMAT(/' MINUS 2 LOG LIKELIHOOD = ',F13.1) 593 11000 FORMAT(/,1X,18A4 ) 594 11500 FORMAT(' DET = ',E15.7, ' IDET = ',I2/) 595 12000 FORMAT(2X,I4) 596 13000 FORMAT( ' NUMBER OF OBSERVATIONS (SAMPLE SIZE), N . . ', I4/) 597 14000 FORMAT(1X, 'CLUSTERING:') 598 15000 FORMAT(18A4) 599 87000 FORMAT(1X,2E15.7) 600 86000 FORMAT(/' TIME BEGAN: ') 601 89000 FORMAT(/' TIME FINISHED: ') 602 16000 FORMAT(/,1X,'CASES AND LABELS:--'/) 603 17000 FORMAT(' PREC.MX. OF CLUSTER ',I2,' HAS BEEN SET TO POOLED ', X'PREC.MX.') 604 18000 FORMAT(1X, (7(1X,A6,': ',I2)/) ) 605 21999 FORMAT(' '//) 606 22000 FORMAT(1X,'MEAN VECTOR FOR CLUSTER ',I2,': ') 607 26000 FORMAT(/, 1X,'K = ',I2,' CLUSTERS'/) 608 30000 FORMAT(1X,'MEAN VECTOR FOR CLUSTER ',I2,': '/,(9F13.3/)) 609 32000 FORMAT(/,1X,'ITERATION ', I2) 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN FEB 15, 1994 23:20:25 NAME:MAIN# PAGE: 21 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 610 11096 FORMAT(/' NUMBER OF CLUSTERS TO REPORT. . . . . . . . 1 TO', XI2/) 611 11111 FORMAT(/' FIRST FOUR DATA POINTS') 612 14560 FORMAT(1X,'THE DENOMINATOR OF THE POSTERIOR PROBS FOR ', X1X,A6, ' IS 0, SO THE POSTERIOR PROBS PP(IC,I), IC=1,...,K,'/ X' HAVE BEEN SET EQUAL TO THOSE FROM THE ', X' PRECEDING ITERATION SO THAT THE I-TH CASE WILL STILL '/ X' HAVE WEIGHT IN SUCCEEDING COMPUTATIONS.') 613 21000 FORMAT(1X, 'MAXIMUM OF EACH VARIABLE: ') 614 36000 FORMAT(/,1X,'NUMBERS IN CLUSTERS:',3X,9(I5,3X)/) 615 40000 FORMAT(/,1X,'MIXING PROBABILITIES:',3X,9(F4.2,3X)/) 616 42000 FORMAT(/, 1X,'COMMON COVARIANCE MATRIX (MLE):',/ ) 617 34000 FORMAT(1X,'MODEL SELECTION CRITERIA'// X' MODEL-SELECTION CRITERIA '/ X' NUMBER MINUS 2 NUMBER (PPH=POST.PROB.OF MODEL):'/ X' OF TIMES OF AKAIKE SCHWARZ '/ X' CLUSTERS, LOG OF PARAM- ------------ ------------ '/ X' K LIKELIHOOD ETERS VALUE PPH VALUE PPH '/ X' --------- ---------- ----- ------- --- ------- ---'//) C 9 XXXX.X 143 XXXX.X .XX XXXX.X .XX C WRITE(6,46000) K,XM2LL(K),IPARM(K),AICVEC(K),AICPPH(K), C XSCHVEC(K),SCHPPH(K) 618 46000 FORMAT(5X,I1,6X,F8.1, 4X,I4,1X,F9.1,1X,F4.2,1X,F9.1,1X,F4.2) 619 48000 FORMAT(1X,5F13.3) 620 50000 FORMAT(1X,A6,': ',I2,(8F13.3/)) 621 52000 FORMAT(//,1X,'CLUSTER MEMBERSHIP OF EACH CASE:'/ X1X,'CASE LABEL DATA'/ ) 622 54000 FORMAT(3X,I2) 623 56000 FORMAT(/' NUMBER OF VARIABLES . . . . . . . . . . . . ',I3 ) 624 58000 FORMAT(1X, 'CONVERGENCE: NO CASE CHANGED CLUSTERS AFTER ', X'ITERATION ',I2,'.'/ ) 625 60000 FORMAT(/,1X,'JFLG = ',I2,'. IF JFLG=0, COMPUTATION OF DET', X' WENT WELL; OTHERWISE, THERE WAS TROUBLE OR MATRIX WAS ', X'ILL-CONDITIONED.'//) 626 62000 FORMAT(/1X,'DET = ',E13.5,' IDET = ',I3,5X, X'ACTUAL DET. = DET*10**IDET',//) 627 70000 FORMAT(1X,' AIC = ', F13.1) 628 78000 FORMAT(1X,' SCHWARZ CRITERION = ', F13.1) 629 72000 FORMAT( ' NUMBER OF PARAMETERS = ',I4) 630 74000 FORMAT(/1X,'EXPECTED NO. OF OBSERVATIONS IN GROUP ',I3,'IS ', X' LESS THAN 0.5. ', X' END COMPUTATION FOR THIS NUMBER OF CLUSTERS. '/) 631 76000 FORMAT(/' POSTERIOR PROB OF GROUP MEMBERSHIP FOR 1ST 4 CASES:') 632 82000 FORMAT(1X,29F5.2) 633 88000 FORMAT(' ALGORITHM HAS NOT CONVERGED IN 20 ITERATIONS. STOP') 634 96000 FORMAT(/' NGFLAG=1: FOR SOME K, AT LEAST ONE CLUSTER HAS '/ X' AN EXP.NO. OF OBSERVATIONS LESS TH AN 0.5 .') 635 98000 FORMAT(//1X,'PROGRAM ENDED NORMALLY.') 636 END 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN FEB 15, 1994 23:20:25 NAME:MAIN# PAGE: 22 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) +_______ ____ ______ ________ _______________________________________ 0A R*8 AT 25 32 413S 417B AIC R*4 178S 179F 182F 286S 288F 290F 495S AICPPH R*8 AT 18 38 510S 511F 511S 515F 515S 518S 519F 519S 522F 522S 526F 526S 534F 535F 535S 537S 544F 545F 549F 549S 554F AICVEC R*4 A 18 182S 290S 554F ALOG R*4 I 180 287 496 512 520 CMPREC R*8 AT 7 30 403S 432F COMDET R*8 T 42 409S 428F COVMX R*8 AT 14 31 165S 167B 231S 233B 398S 400B DATFMT R*4 A 11 56S 75F 77F 206F DENOM R*8 AT 5 39 293S 295F 295S 298F 306F DET R*8 T 33 167B 170F 233B 241F 400B 408F 409F 417B 425F 426F DEVV R*8 T 30 255S 258F 364S 367F 369F DEVW R*8 T 30 257S 258F 365S 368F 369F DEXP R*8 I 269 525 526 529 535 DLOG R*8 I 171 171 266 266 280 509 510 517 518 DSQ R*8 AT 8 34 253S 258F 258S 260F DSQRT R*8 I 220 222 DTIME X 570F F R*8 AT 3 36 262S 266S 267F 269F 269S 277F 295F 306F I I*4 57S 58 58 59F 68F 69F 70F 71B 80S 81F 81B 139S 140B 144S 145F 145B 193S 195B 199S 201F 247S 248F 248B 251S 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN FEB 15, 1994 23:20:25 NAME:MAIN# 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) +_______ ____ ______ ________ _______________________________________ 0 255F 257F 262F 266F 267F 269F 269B 274S 277B 292S 293F 295F 295F 295B 297S 298F 299F 301F 301F 306F 306F 306B 314S 315B 319S 320F 320B 325S 326B 345S 346F 346B 363S 364F 365F 366F 369B 436S 437F 437B 440S 441F 441F 442F 444F 444F 446F 446F 447B 449S 451F 452F 452B 457S 457F 457B 460S 461F 461B 469S 470B 490S 491F 491F 491B 500S 501F 502B IC I*4 190S 191F 191B 194S 195B 210S 213F 216B 223S 224F 224B 226S 231F 238B 240S 241B 276S 277F 277B 294S 295F 295B 300S 301F 301B 305S 306F 306F 306B 315S 315S 320S 320S 323S 324F 326F 326F 326B 328S 329F 329B 331S 332F 333B 343S 346F 346F 346B 360S 364F 365F 366F 369F 369F 369B 371S 374F 374B 376S 379F 379F 379B 387S 388F 388B 410S 413F 420F 426F 427F 428F 429F 432B 443S 444F 446F 447B 450S 451F 452F 452B 466S 467B 470S 471F 471F 473S 473S 480S 480S 481S 482F 482B 498S 499F 501B ICLSOL I*4 A 4 248S 437S 461F ICLUS I*4 A 4 248F 437F 442S 447S 457F 461F 470F 491F 501F ID CHAR AT 1 27 58S 81F 299F 320F 457F 491F 502F IDAY I*4 569F 571B 573F 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN FEB 15, 1994 23:20:25 NAME:MAIN# PAGE: 24 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) +_______ ____ ______ ________ _______________________________________ 0IDET I*4 161S 167B 170F 227S 233B 241F 395S 400B 408F 409F 415S 417B 425F 426F IG I*4 198S 199F 201B 203S 205F 206B 337S 339F 341B 348S 349F 352F 357F 357F 357B IHOUR I*4 569F 570B 573F INPFMT I*4 A 11 55S 58 81F INT GI 199 IP I*4 52S 53F 58 62F 67F 75F 77F 81F 133F 135F 138F 142F 143F 147F 150F 151F 156F 158F 159F 163F 164F 167B 171F 171F 175F 175F 175F 187F 187F 187F 200F 206F 211F 212F 215F 219F 221F 229F 230F 233B 236F 237F 254F 256F 265F 338F 340F 344F 356F 361F 362F 372F 377F 378F 381F 382F 385F 386F 390F 391F 396F 397F 400B 401F 402F 411F 412F 417B 418F 419F 430F 431F 482F 485F 486F 491F IPARM I*4 A 18 176S 188S 509F 510F 511F 517F 518F 519F 554F ISEC I*4 569F 570B 573F ITER I*4 243S 245F 246F 435F 459F 474F 474S 476F 494F IV I*4 A 26 167B 233B 400B 417B IYEAR I*4 569F 571B 573F JFLG I*4 167B 168F 169F 233B 234F 235F 400B 405F 406F 407F 417B 422F 423F 424F JV I*4 58S 58S 62S 63F 63F 64F 64B 67S 68F 68F 69F 69F 70F 70F 71F 71B 75S 75S 77S 77S 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN FEB 15, 1994 23:20:25 NAME:MAIN# PAGE: 25 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 81S 81S 133S 134F 136B 138S 140F 140F 140B 142S 145F 145F 145B 147S 148F 148B 150S 152F 152F 152F 153F 153B 156S 156S 158S 159B 163S 165F 165B 200S 201F 201B 206S 206S 211S 213B 215S 216F 216B 219S 220F 220F 224F 224B 229S 231F 231B 236S 238F 238B 254S 255F 255F 258B 338S 339F 341B 344S 346F 346F 346B 356S 357F 357B 361S 362F 364F 364F 369F 369B 372S 373F 374F 374B 377S 379F 379B 381S 383B 385S 388F 388F 388B 390S 392F 392B 396S 398F 398B 401S 403F 403B 411S 413F 413B 418S 420F 420B 430S 432F 432B 482S 482S 485S 486B 491S 491S JW I*4 135S 136B 143S 145F 145F 145B 151S 152F 152F 152F 153F 153B 159S 159S 164S 165F 165B 212S 213B 221S 222F 222F 224F 224B 230S 231F 231B 237S 238F 238B 256S 257F 257F 258B 340S 341B 362S 365F 365F 369F 369B 373S 374F 374B 378S 379F 379B 382S 383B 386S 388F 388F 388B 391S 392F 392B 397S 398F 398B 402S 403F 403B 412S 413F 413B 419S 420F 420B 431S 432F 432B 486S 486S K I*4 185S 186F 187F 187F 187F 188F 190F 191F 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN FEB 15, 1994 23:20:25 NAME:MAIN# PAGE: 26 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 194F 195F 198F 199F 203F 210F 223F 226F 240F 252F 276F 284F 290F 291F 294F 300F 305F 315F 320F 323F 328F 331F 333F 337F 343F 348F 360F 371F 376F 387F 410F 443F 450F 451F 466F 473F 480F 481F 498B 508F 516S 517F 517F 517F 518F 518F 518F 519F 519F 519F 520F 520F 521F 522F 522F 522F 523F 523F 523B 527S 528F 529F 529F 531B 533S 534F 535F 535F 537B 541S 542F 543F 544F 545B 547S 548F 548F 549F 549B 553S 554F 554F 554F 554F 554F 554F 554B KREACH I*4 83S 508S 516F 527F 533F 541F 547F 553F KUP I*4 82S 83F 85F 185F 512F 520F L I*4 252S 253F 255F 257F 258F 258F 258F 260F 262F 266F 266F 267F 268F 269F 269B MATEQ X 167F 233F 400F 417F MAXITER I*4 475S 476F MIN GI 82 MINUTE I*4 569F 570B 573F MONTH I*4 569F 571B 573F N I*4 51S 54F 57F 82F 86F 139F 144F 193F 199F 247F 251F 274F 292F 297F 311F 319F 325F 329F 345F 363F 379F 392F 436F 440F 449F 454F 457F 460F 469F 488F 490F 497F 500F NG I*4 A 10 467S 471F 471S 473F NGFLAG I*4 88S 353S 560F NOPARM I*4 175S 176F 177F 178F 180F 187S 188F 189F 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN FEB 15, 1994 23:20:25 NAME:MAIN# PAGE: 27 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 286F 287F 495F 496F NRS1 I*4 162S 167B 228S 233B 394S 400B 416S 417B P R*8 AT 25 32 167B 233B 238F 400B 403F 417B 420F PI R*8 T 44 87S 171F 266F 509F 510F 517F 518F PP R*8 AT 3 38 301S 306S 315F 320F 326F 346F 366F 369F 441F 444F 446F 452F PPOLD R*8 AT 19 40 195S 301F 451S 452S PR R*8 AT 17 36 191S 277F 295F 306F 324S 326F 326S 329F 329S 332F 349F 357F 379F 388F 480F PREC R*8 AT 7 30 238S 258F 420S 432S PRIRLN R*8 AT 20 41 513S 514F 515F 521S 522F 523F SCH R*4 180S 181F 183F 287S 289F 291F 496S SCHPPH R*8 AT 18 38 509S 512F 512S 514F 514S 517S 520F 520S 523F 523S 525F 525S 528F 529F 529S 531S 542F 543F 548F 548S 554F SCHVEC R*4 A 18 183S 291S 509F 510F 517F 518F 554F SDJV R*8 T 31 220S 224F SDJW R*8 T 31 222S 224F SIGMA R*8 AT 6 32 213S 216S 224F 224S 231F 379S 413F SPPHA R*8 T 38 540S 545F 545S 549F SPPHS R*8 T 38 539S 543F 543S 548F SSD R*8 AT 6 29 341S 369F 369S 374F 374S 379F 388F SSDEVS R*8 AT 24 29 152S 153F SUM R*8 AT 1 28 339S 346F 346S 357F SUMLNF R*8 T 37 273S 280F 280S 283F SUMPXF R*8 T 37 275S 277F 277S 279F 280F SUMSQS R*8 AT 23 29 136S 145F 145S 152F TDATE X 571F 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN FEB 15, 1994 23:20:25 NAME:MAIN# PAGE: 28 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) +_______ ____ ______ ________ _______________________________________ 0TERM R*8 T 29 UNREFERENCED TEST R*8 T 43 349S 350F TITLE R*4 A 9 45S 50F 551F 557F TOTAL R*8 AT 21 29 134S 140F 140S 148F 152F 152F TRUDET R*8 AT 12 42 170S 171F 241S 266F 268F 426S 427F 428S VARHAT R*8 AT 14 31 153S 159F 165F 220F 222F 392S 398F 486F WGSS R*8 AT 13 29 383S 388F 388S 392F X R*4 A 1 58S 63F 64F 68F 69F 70F 71F 81F 140F 145F 145F 201F 255F 257F 346F 364F 365F 491F XBAR R*8 AT 10 35 201S 206F 255F 257F 357S 364F 365F 482F XIP R*4 265S 266F XJ R*4 A 16 89S 90S 91S 92S 93S 94S 95S 96S 97S 98S 99S 100S 101S 102S 103S 104S 105S 106S 107S 108S 109S 110S 111S 112S 113S 114S 115S 116S 117S 118S 119S 120S 121S 122S 123S 124S 125S 126S 127S 128S 129S 130S 131S 132S 191F XLL R*4 171S 172F XMAX R*4 A 15 63S 70F 71S 77F XMEAN R*4 A 22 148S 156F XMIN R*4 A 15 64S 68F 69S 75F XMN2LL R*4 172S 173F 174F 178F 180F 283S 284F 285F 286F 287F 495F 496F XMXPR R*4 A 5 441S 444F 446S XM2LL R*4 A 17 173S 284S 554F XN R*4 86S 148F 152F 153F 171F 171F 171F 180F 287F 349F 357F 496F XUFLOW X 359F Y R*4 A 2 UNREFERENCED 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN FEB 15, 1994 23:20:25 NAME:MAIN# PAGE: 29 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) +_______ ____ ______ ________ _______________________________________ 0ZSQ R*8 T 34 260S 261F 266F 268F 0 0VARIABLES REFERENCED BUT NOT SET. (* POSSIBLY SET AS ARGUMENT.) 0DET* IDAY* IHOUR* ISEC* IV* IYEAR* JFLG* MINUTE* MONTH* P* 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN FEB 15, 1994 23:20:25 NAME:MAIN# PAGE: 30 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 61 59 105 137 133 135 110 141 138 139 111 84 80 113 B 209 208 114 213 UNREFERENCED 120 146 142 143 144 130 149 147 140 154 150 151 150 160 158 200 65 62 250 B 458 454 300 B 66 60 360 192 190 400 72 67 500 73 57 549 196 193 194 550 202 198 200 600 207 203 605 555 553 610 214 211 212 611 225 219 221 223 620 217 215 630 218 210 640 239 226 236 237 680 242 240 700 B 244 477 800 249 247 825 438 436 850 B 506 185 354 851 503 498 500 852 B 504 497 900 B 439 435 901 B 250 246 1000 271 252 1090 B 264 261 1100 B 270 263 1200 272 251 1310 259 254 256 1350 296 292 294 1400 309 297 1401 B 308 303 1402 307 305 1410 B 304 298 1420 316 314 1421 B 318 311 311 1422 B 312 311 1423 B 322 317 1424 321 319 1456 302 300 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN FEB 15, 1994 23:20:25 NAME:MAIN# PAGE: 31 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 1500 B 446 444 1600 B 448 440 443 445 1601 453 449 450 1750 342 337 338 340 1875 347 343 344 345 1900 358 348 356 1950 384 381 382 2050 B 352 350 2100 278 276 2150 B 355 351 2190 B 281 279 2200 282 274 2300 389 385 386 387 2400 393 390 391 2402 166 163 164 2404 232 229 230 2406 399 396 397 2500 487 485 2800 483 481 2900 B 463 460 461 3000 B 465 459 462 3100 B 478 476 3300 B 479 464 3311 468 466 3321 472 469 3500 492 490 3510 B 493 488 3600 B 370 360 361 362 363 366 367 368 3700 375 371 372 373 3710 380 376 377 378 3711 B 434 410 427 3712 414 411 412 3713 421 418 419 3800 327 323 325 3813 404 401 402 3814 433 430 431 3900 330 328 3901 B 336 331 332 4001 B 565 562 4002 B 563 560 4004 B 507 335 6110 524 516 6115 B 531 528 6119 B 532 527 530 6120 546 541 6121 B 538 533 536 6126 B 537 534 6130 550 547 10000 NF 575 47 10050 NF 576 46 48 184 505 556 566 10100 NF 579 156 10200 NF 581 157 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN FEB 15, 1994 23:20:25 NAME:MAIN# PAGE: 32 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 10500 NF 582 155 10900 NF 592 174 285 10999 NF 580 49 11000 NF 593 551 11001 NF 583 50 11002 NF 584 557 11096 NF 610 85 11097 NF 585 567 11100 NF 586 569 573 11111 NF 611 78 11112 NF 577 79 11500 NF 594 UNREFERENCED 12000 NF 595 51 13000 NF 596 54 14000 NF 597 455 14560 NF 612 299 15000 NF 598 45 55 56 16000 NF 602 456 17000 NF 603 429 18000 NF 604 457 19000 NF 578 74 20000 NF 589 197 21000 NF 613 76 21999 NF 605 204 22000 NF 606 205 25000 NF 587 320 26000 NF 607 186 30000 NF 608 482 32000 NF 609 245 34000 NF 617 552 36000 NF 614 473 40000 NF 615 480 40040 NF 334 333 42000 NF 616 484 44444 NF 590 499 44445 NF 591 502 46000 NF 618 554 48000 NF 619 159 486 50000 NF 620 491 52000 NF 621 489 54000 NF 622 52 55555 NF 588 310 56000 NF 623 53 58000 NF 624 494 60000 NF 625 169 235 406 423 62000 NF 626 408 425 70000 NF 627 179 288 72000 NF 629 177 189 74000 NF 630 352 76000 NF 631 313 78000 NF 628 181 289 82000 NF 632 315 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN FEB 15, 1994 23:20:25 NAME:MAIN# PAGE: 33 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 86000 NF 600 568 87000 NF 599 268 88000 NF 633 478 89000 NF 601 572 96000 NF 634 561 96666 NF 559 558 98000 NF 635 564 0 NUMBER MODULE LEVEL ISN VS FORTRAN ERROR MESSAGES + ______ ______ _____ ___ __ _______ _____ ________ 0ILX3470I VPCQ 0(I) "IDAY,IHOUR,ISEC,IYEAR,MINUTE,MONTH" MAY BE UNINITIALIZED WHEN USED 0*STATISTICS* SOURCE STATEMENTS: 624, PROGRAM SIZE: 8358440 BYTES, PROGRAM NAME: MAIN#, PAGE: 1 *STATISTICS* 1 DIAGNOSTIC GENERATED. SEVERITY CODE IS 0. **MAIN#** END OF COMPILATION 1 ****** TIME STAMP: 94.04623.20.25 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN FEB 15, 1994 23:20:34 PAGE: 34 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 C 1 SUBROUTINE MATEQ(A,M,N,JFLG,DET,IDET,IV,NRS1,P,LL) C 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. 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN FEB 15, 1994 23:20:34 NAME:MAIN# PAGE: 35 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 C NRS=0, INVERT THE A MATRIX. C NRS<0, ONLY COMPUTE THE DETERMINANT OF A. 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 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN FEB 15, 1994 23:20:34 NAME:MATEQ PAGE: 36 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 6 IF (NRS.EQ.0) IDET=1 8 DISIGN=1.0D+00 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. 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN FEB 15, 1994 23:20:34 NAME:MATEQ PAGE: 37 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 C 1 43 K1=K+1 1 44 DO 1100 IM=K1,M C C BEFORE ACTUALLY ELIMINATING WE CHECK TO SEE IF A(IV(IM),K) HAS ALREADY C BEEN ANIHILATED. 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 FEB 15, 1994 23:20:34 NAME:MATEQ PAGE: 38 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 FEB 15, 1994 23:20:34 NAME:MATEQ PAGE: 39 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 FEB 15, 1994 23:20:34 NAME:MATEQ PAGE: 40 0*STATISTICS* SOURCE STATEMENTS: 75, PROGRAM SIZE: 3564 BYTES, PROGRAM NAME: MATEQ, PAGE: 34 *STATISTICS* NO DIAGNOSTICS GENERATED. **MATEQ** END OF COMPILATION 2 ****** TIME STAMP: 94.04623.20.34 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN FEB 15, 1994 23:20:35 PAGE: 41 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 FEB 15, 1994 23:20:35 NAME:MATDT PAGE: 42 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: 41 *STATISTICS* NO DIAGNOSTICS GENERATED. **MATDT** END OF COMPILATION 3 ****** TIME STAMP: 94.04623.20.35 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN FEB 15, 1994 23:20:35 PAGE: 43 0SUMMARY OF MESSAGES AND STATISTICS FOR ALL COMPILATIONS 0 NUMBER MODULE LEVEL ISN VS FORTRAN ERROR MESSAGES + ______ ______ _____ ___ __ _______ _____ ________ 0ILX3470I VPCQ 0(I) "IDAY,IHOUR,ISEC,IYEAR,MINUTE,MONTH" MAY BE UNINITIALIZED WHEN USED 0*STATISTICS* SOURCE STATEMENTS: 624, PROGRAM SIZE: 8358440 BYTES, PROGRAM NAME: MAIN#, PAGE: 1 *STATISTICS* 1 DIAGNOSTIC GENERATED. SEVERITY CODE IS 0. **MAIN#** END OF COMPILATION 1 ****** TIME STAMP: 94.04623.20.35 0*STATISTICS* SOURCE STATEMENTS: 75, PROGRAM SIZE: 3564 BYTES, PROGRAM NAME: MATEQ, PAGE: 34 *STATISTICS* NO DIAGNOSTICS GENERATED. **MATEQ** END OF COMPILATION 2 ****** TIME STAMP: 94.04623.20.35 0*STATISTICS* SOURCE STATEMENTS: 31, PROGRAM SIZE: 1552 BYTES, PROGRAM NAME: MATDT, PAGE: 41 *STATISTICS* NO DIAGNOSTICS GENERATED. **MATDT** END OF COMPILATION 3 ****** TIME STAMP: 94.04623.20.35 0******* SUMMARY STATISTICS ******* 1 DIAGNOSTICS GENERATED. HIGHEST SEVERITY CODE IS 0.