1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAR 03, 1994 21:51:41 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 1991 STANLEY LOUIS SCLOVE. C C C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C C C C C CMS FILE NAME OF PROGRAM: MIXPCMA CLUSPAC C C VERSION 4.6 3-MAR-94 C C C C C C DEVELOPED AND PROGRAMMED BY: C C C C DR. STANLEY L. SCLOVE 312/996-2681 C C DEPARTMENT C C OF INFORMATION AND DECISION SCIENCES 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 C C C COPYRIGHT 1991 STANLEY LOUIS SCLOVE. C 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 C model likelihood. C C C C Program MIXPCMA in the CLUSPAC package 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 MIXPCM and MIXPCMA assume a common covariance matrix C C across distributions. MIXPDT and MIXPDTA allow different C C covariance matrices. C C MIXPCMA tries a range of numbers of clusters, with C C automatic setting of initial values of the parameters. C C C C C 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAR 03, 1994 21:51:41 NAME:MAIN# PAGE: 2 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 C ..............................................................C C C C C C C RESTRICTIONS (CAN BE MODIFIED): C N, SAMPLE SIZE, AT MOST 1000; C IP, NUMBER OF VARIABLES, AT MOST 20; C K, NUMBER OF CLUSTERS, AT MOST 29; C ITER, MAXIMUM NUMBER OF ITERATIONS, 20. C C SUBROUTINE(S) CALLED: C MATEQ, WHICH CALLS MATDT C C C C C ..............................................................C C C 1 DIMENSION X(1000,20),SUM(29,20),ID(1000) 2 DIMENSION F(1000,29),PP(29,1000),PPOLD(29,1000) 3 DIMENSION ICLUS(1000),ICLSOL(1000) 4 DIMENSION DENOM(1000),XMXPR(1000) 5 DIMENSION DSQ(29) 6 DIMENSION TITLE(18) 7 DIMENSION NG(29),XBAR(29,20) 8 DIMENSION DATFMT(18),INPFMT(18) 9 DIMENSION SSD(29,20,20) 10 DIMENSION WGSS(20,20) 11 DIMENSION VARHAT(20,20),COVMX(20,20) 12 DIMENSION XMIN(20),XMAX(20) 13 DIMENSION IV(20,20) 14 DIMENSION P(20,20) 15 DIMENSION XJ(9,9) 16 DIMENSION PR(29),IPARM(29),XM2LL(29) 17 DIMENSION AICVEC(9),SCHVEC(9),SCHPPH(9),AICPPH(9),PRIRLN(9) 18 DIMENSION TOTAL(29) 19 DIMENSION XMEAN(29) 20 DIMENSION SUMSQS(20,20) 21 DIMENSION SSDEVS(20,20) C 22 CHARACTER*6 ID C 23 DOUBLE PRECISION SUM 24 DOUBLE PRECISION WGSS,SSD 25 DOUBLE PRECISION PI 26 DOUBLE PRECISION VARHAT,COVMX 27 DOUBLE PRECISION SSDEVS,SUMSQS,TOTAL 28 DOUBLE PRECISION SUMLNF,SUMPXF 29 DOUBLE PRECISION PR,DENOM,PP,AICPPH,SCHPPH,PRIRLN,PPOLD 30 DOUBLE PRECISION P 31 DOUBLE PRECISION DET,ZSQ,TRUDET 32 DOUBLE PRECISION DSQ 33 DOUBLE PRECISION XBAR 34 DOUBLE PRECISION DEVV,DEVW 35 DOUBLE PRECISION F 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAR 03, 1994 21:51:41 NAME:MAIN# PAGE: 3 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 C C C C C IV IS A WORK ARRAY FOR SUBROUTINE MATEQ. 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 AT LEAST ONE C BLANK AT THE BEGINNING FOR CARRIAGE CONTROL. C DATA, ONE CASE AT A TIME, IN FORMAT SPECIFIED BY DATFMT C C 36 READ (5,15000) TITLE C C WRITE PROGRAM INFORMATION. 37 WRITE (6,10050) 38 WRITE (6,10000) 39 WRITE (6,10050) 40 WRITE (6,10999) 41 WRITE (6,11000) TITLE C THE FOLLOWING ARE IMSL DATE AND TIME FUNCTIONS: 42 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. 43 READ (5,12000) N C READ NUMBER OF VARIABLES, IP. 44 READ (5,17000) IP 45 WRITE (6,19000) IP 46 WRITE (6,13000) N C C READ DATA FORMAT. 47 READ (5,15000) INPFMT 48 READ (5,15000) DATFMT C C READ DATA. 49 DO 500 I = 1,N 1 50 READ (5,INPFMT) (X(I,JV), JV = 1,IP),ID(I) 1 51 IF (I .EQ. 1) GO TO 100 1 52 GO TO 300 1 53 100 CONTINUE 1 54 DO 200 JV = 1,IP 2 55 XMAX(JV) = X(1,JV) 2 56 XMIN(JV) = X(1,JV) 2 57 200 CONTINUE 1 58 300 CONTINUE 1 59 DO 400 JV = 1,IP 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAR 03, 1994 21:51:41 NAME:MAIN# PAGE: 4 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 2 60 IF (X(I,JV) .LT. XMIN(JV)) XMIN(JV)=X(I,JV) 2 62 IF (X(I,JV) .GT. XMAX(JV)) XMAX(JV)=X(I,JV) 2 64 400 CONTINUE 1 65 500 CONTINUE 66 WRITE (6,16000) 67 WRITE (6,DATFMT) (XMIN(JV),JV=1,IP) 68 WRITE (6,15500) 69 WRITE (6,DATFMT) (XMAX(JV),JV=1,IP) C 70 WRITE (6,11111) 71 WRITE (6,11112) 72 DO 111 I=1,4 1 73 WRITE (6,INPFMT) (X(I,JV),JV=1,IP),ID(I) 1 74 111 CONTINUE 75 KUP = MIN(9,N) 76 WRITE (6,11096) KUP 77 11096 FORMAT(/' NUMBER OF CLUSTERS TO REPORT. . . . . . . . 1 TO', XI2/) C C SET CONSTANTS. 78 XN = N 79 PI = .314159265358979324D+01 C C SET INITIAL VALUES OF MIXING PROBABILITIES:-- C C SET CONSTANTS. C XJ(IC,K), K=1 TO 9, IC = 1 TO K C OPTIMAL CLASS PROBABILITIES - NORMAL DISTRIBUTION C JOHARI & SCLOVE (1975). "PARTITIONING A DISTRIBUTION." C COMMUNICATIONS IN STATISTICS. 80 XJ(1,2)=.5 81 XJ(2,2)=.5 82 XJ(1,3)=.2703 83 XJ(2,3)=.4594 84 XJ(3,3)=.2703 85 XJ(1,4)=.1631 86 XJ(2,4)=.3369 87 XJ(3,4)=.3369 88 XJ(4,4)=.1631 89 XJ(1,5)=.1068 90 XJ(2,5)=.2444 91 XJ(3,5)=.2976 92 XJ(4,5)=.2444 93 XJ(5,5)=.1068 94 XJ(1,6)=.0739 95 XJ(2,6)=.1810 96 XJ(3,6)=.2451 97 XJ(4,6)=.2451 98 XJ(5,6)=.1810 99 XJ(6,6)=.0739 100 XJ(1,7)=.0536 101 XJ(2,7)=.1375 102 XJ(3,7)=.1986 103 XJ(4,7)=.2106 104 XJ(5,7)=.1986 105 XJ(6,7)=.1375 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAR 03, 1994 21:51:41 NAME:MAIN# PAGE: 5 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 106 XJ(7,7)=.0536 107 XJ(1,8)=.0402 108 XJ(2,8)=.1067 109 XJ(3,8)=.1613 110 XJ(4,8)=.1918 111 XJ(5,8)=.1918 112 XJ(6,8)=.1613 113 XJ(7,8)=.1067 114 XJ(8,8)=.0402 115 XJ(1,9)=.0310 116 XJ(2,9)=.0845 117 XJ(3,9)=.1324 118 XJ(4,9)=.1643 119 XJ(5,9)=.1756 120 XJ(6,9)=.1643 121 XJ(7,9)=.1324 122 XJ(8,9)=.0845 123 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: 124 DO 105 JV = 1,IP 1 125 TOTAL(JV) = 0.0D+00 1 126 DO 105 JW = 1,IP 2 127 SUMSQS(JV,JW) = 0.0D+00 2 128 105 CONTINUE 129 DO 110 JV = 1,IP 1 130 DO 110 I = 1,N 2 131 TOTAL(JV) = TOTAL(JV) + X(I,JV) 2 132 110 CONTINUE 133 DO 120 JV = 1,IP 1 134 DO 120 JW = 1,IP 2 135 DO 120 I = 1,N 3 136 SUMSQS(JV,JW) = SUMSQS(JV,JW) + X(I,JV)*X(I,JW) 3 137 120 CONTINUE 138 DO 130 JV = 1,IP 1 139 XMEAN(JV) = TOTAL(JV)/XN 1 140 130 CONTINUE 141 DO 140 JV = 1,IP 1 142 DO 140 JW = 1,IP 2 143 SSDEVS(JV,JW) = SUMSQS(JV,JW)- TOTAL(JV)*TOTAL(JW)/XN 2 144 VARHAT(JV,JW) = SSDEVS(JV,JW)/XN 2 145 140 CONTINUE C WRITE SUMMARY STATISTICS FOR WHOLE SAMPLE: 146 WRITE (6,10050) 147 WRITE(6,21000) 148 WRITE(6,22000) 149 WRITE(6,10100) ( XMEAN(JV), JV = 1,IP ) 150 WRITE(6,10200) 151 DO 150 JV = 1,IP 1 152 WRITE(6,48000) ( VARHAT(JV,JW), JW=1,IP ) 1 153 150 CONTINUE C C C CALL SUBROUTINE TO COMPUTE DETERMINANT OF SAMPLE COVARIANCE 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAR 03, 1994 21:51:41 NAME:MAIN# PAGE: 6 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 C MATRIX: C 154 IDET = 1 155 NRS1 = 0 156 DO 2404 JV = 1,IP 1 157 DO 2404 JW = 1,IP 2 158 COVMX(JV,JW) = VARHAT(JV,JW) 2 159 2404 CONTINUE 160 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 161 IF (JFLG .GT. 0) WRITE (6,60000) JFLG C C 163 XIDET = IDET C XLGDET = DLOG(DET) + XIDET*ALOG(10.0) 164 TRUDET = DET*(10.0**IDET) 165 IF (JFLG .GT. 0) 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 COMPUTE LOG OF MAX LIKELIHOOD: 167 XLL=-(XN*IP/2.0)*DLOG(2.0*PI)-(XN/2.0)*DLOG(TRUDET)-XN*IP/2.0 168 XMN2LL = (-2.0)*XLL 169 XM2LL(1) = XMN2LL 170 WRITE(6,37000) XMN2LL C NO. OF PARAMETERS FOR UNCLUSTERED SAMPLE IS: C 1 MEAN VECTOR + 1 COV MX. 171 NOPARM = IP + IP*(IP+1)/2 172 IPARM(1) = NOPARM C 173 AIC = XMN2LL + 2.0*NOPARM 174 WRITE(6,70000) AIC 175 SCH = XMN2LL + ALOG(XN)*NOPARM 176 WRITE (6,71000) SCH 177 AICVEC(1) = AIC 178 SCHVEC(1) = SCH 179 WRITE (6,10050) C C C PERFORM CLUSTERING FOR VARIOUS NUMBERS OF CLUSTERS, K: C 180 DO 850 K = 2,KUP 1 181 WRITE (6,26000) K C C C SET CONSTANTS FOR THIS VALUE OF K. C PARAMETERS FOR MODEL WITH COMMON COVARIANCE MATRIX ARE C K MEAN VECTORS AND ONE COVARIANCE MATRIX AND K-1 INDEPENDENT C MIXING PROBABILITIES. C 1 182 NOPARM = K*IP + IP*(IP+1)/2 + (K-1) 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAR 03, 1994 21:51:41 NAME:MAIN# PAGE: 7 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 C C 1 183 WRITE (6,72000) NOPARM 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 184 DO 360 IC = 1,K 2 185 PR(IC) = XJ(IC,K) 2 186 360 CONTINUE C C TO SET INITIAL VALUES OF PRIOR PROBABILITIES EQUAL TO 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 187 DO 549 I = 1,N 2 188 DO 549 IC = 1,K 3 189 PPOLD(IC,I) = 1.0/K 3 190 549 CONTINUE C 1 191 WRITE (6,20000) C C SET INITIAL MEANS EQUAL TO K OF THE OBSERVATIONS 1 192 DO 550 JV = 1,IP 2 193 XBAR(1,JV) = X(1,JV) 2 194 550 CONTINUE 1 195 DO 555 IG = 2,K 2 196 I = INT( (IG-1)*N/(K-1) ) 2 197 DO 555 JV = 1,IP 3 198 XBAR(IG,JV) = X(I,JV) 3 199 555 CONTINUE 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 200 DO 600 IG=1,K 2 201 WRITE (6,DATFMT) ( XBAR(IG,JV), JV=1,IP) 2 202 600 CONTINUE C C SET INITIAL COVARIANCE MATRIX: C (HERE SOME SLIGHT POSITIVE CORRELATION IN PROGRAMMED IN. 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAR 03, 1994 21:51:41 NAME:MAIN# PAGE: 8 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 C THIS COULD BE CHANGED.) 1 203 DO 610 JV=1,IP 2 204 DO 610 JW=1,IP 3 205 VARHAT(JV,JW)=0.40 3 206 610 CONTINUE 1 207 DO 620 JV = 1,IP 2 208 VARHAT(JV,JV) = 1.0 2 209 620 CONTINUE C C COMPUTE INITIAL INVERSE COVARIANCE MATRIX: C 1 210 IDET = 1 1 211 NRS1 = 0 1 212 DO 2406 JV = 1,IP 2 213 DO 2406 JW = 1,IP 3 214 COVMX(JV,JW) = VARHAT(JV,JW) 3 215 2406 CONTINUE 1 216 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 1 217 IF (JFLG .GT. 0) WRITE (6,60000) JFLG C C 1 219 XIDET = IDET C XLGDET = DLOG(DET) + XIDET*ALOG(10.0) 1 220 TRUDET = DET*(10.0**IDET) C C 1 221 ITER = 1 1 222 700 CONTINUE 1 223 WRITE(6,32000) ITER 1 224 IF (ITER .EQ. 1) GO TO 910 1 225 DO 800 I = 1,N 2 226 ICLSOL(I) = ICLUS(I) 2 227 800 CONTINUE C 1 228 910 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 229 DO 1200 I = 1,N C 2 230 DO 1000 L=1,K 3 231 DSQ(L) = 0.0D+00 3 232 DO 1310 JV=1,IP 4 233 DEVV = XBAR(L,JV) - X(I,JV) 4 234 DO 1310 JW=1,IP 5 235 DEVW = XBAR(L,JW) - X(I,JW) 5 236 DSQ(L) = DSQ(L) + DEVV*P(JV,JW)*DEVW 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAR 03, 1994 21:51:41 NAME:MAIN# PAGE: 9 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 5 237 1310 CONTINUE 3 238 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 239 IF ( ZSQ/2.0 .LE. 174.673 ) GO TO 1090 3 240 F(I, L) = 0.0D+00 3 241 GO TO 1100 C 3 242 1090 CONTINUE 3 243 F(I, L) =DEXP(-ZSQ/2.0)/((2.0*PI)**(IP/2.0)*DSQRT(TRUDET)) 3 244 1100 CONTINUE 3 245 1000 CONTINUE 2 246 1200 CONTINUE C C COMPUTE LOG LIKELIHOOD AND VALUES OF MODEL-SELECTION CRITERIA: C 1 247 SUMLNF = 0.0D+00 1 248 DO 2200 I=1,N 2 249 SUMPXF = 0.0D+00 2 250 DO 2100 IC=1,K 3 251 SUMPXF = SUMPXF + PR(IC)*F(I,IC) 3 252 2100 CONTINUE C 2 253 IF ( SUMPXF .EQ. 0.0 ) GO TO 2200 2 254 SUMLNF = SUMLNF + DLOG(SUMPXF) 2 255 2200 CONTINUE C 1 256 XMN2LL = -2.0*SUMLNF 1 257 XM2LL(K) = XMN2LL C 1 258 WRITE (6,69000) XMN2LL C C C COMPUTE MODEL SELECTION CRITERIA C C PARAMETERS: C K MEAN VECTORS OF DIMENSION P AND A P-BY-P COVARIANCE MATRIX, C WHERE P IS THE NUMBER OF VARIABLES, AND K-1 INDEPENDENT C MIXING PROBABILITIES. C 1 259 NOPARM = K*IP + IP*(IP+1)/2 + (K-1) 1 260 IPARM(K) = NOPARM C 1 261 AIC = XMN2LL + 2.0*NOPARM 1 262 SCH = XMN2LL + ALOG(XN)*NOPARM 1 263 WRITE (6,70000) AIC 1 264 WRITE (6,71000) SCH 1 265 AICVEC(K) = AIC 1 266 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 267 DO 1350 I = 1,N 2 268 DENOM(I) = 0.0D+00 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAR 03, 1994 21:51:41 NAME:MAIN# PAGE: 10 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 2 269 DO 1350 IC=1,K 3 270 DENOM(I) = DENOM(I) + PR(IC)*F(I,IC) 3 271 1350 CONTINUE C 1 272 DO 1400 I = 1,N 2 273 IF ( DENOM(I) .GT. 0.0D+00 ) GO TO 1410 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 2 274 WRITE (6,14560) ID(I) 2 275 14560 FORMAT(1X,'THE DENOMINATOR OF THE POSTERIOR PROBS FOR '/ 2 X1X,A6, ' IS 0, SO THE POSTERIOR PROBS PP(IC,I), IC=1,...,K, '/ 2 X' HAVE BEEN SET EQUAL TO THOSE FROM THE PRECEDING ITERATION '/ 2 X' SO THAT THE I-TH CASE WILL STILL ', 2 X'HAVE WEIGHT IN SUCCEEDING COMPUTATIONS.') C 2 276 DO 1456 IC = 1,K 3 277 PP(IC,I) = PPOLD(IC,I) 3 278 1456 CONTINUE 2 279 GO TO 1401 C 2 280 1410 CONTINUE 2 281 DO 1402 IC=1,K 3 282 PP(IC,I)= PR(IC)*F(I,IC)/DENOM(I) 3 283 1402 CONTINUE C 2 284 1401 CONTINUE 2 285 1400 CONTINUE C 1 286 IF ( N-15 ) 1421,1421,1422 1 287 1422 CONTINUE 1 288 WRITE(6,76000) 1 289 DO 1420 I = 1,4 2 290 WRITE(6,96000) ( PP(IC,I), IC = 1,K ) 2 291 1420 CONTINUE 1 292 GO TO 1423 1 293 1421 CONTINUE 1 294 DO 1424 I = 1,N 2 295 WRITE (6,25000) I, (PP(IC,I),IC=1,K) 2 296 1424 CONTINUE 1 297 1423 CONTINUE C C C UPDATE PARAMETER ESTIMATES: C C C UPDATE CLUSTER PRIOR PROBABILITIES PR(IC):-- 1 298 DO 3800 IC = 1,K 2 299 PR(IC) = 0.0D+00 2 300 DO 3800 I = 1,N 3 301 PR(IC) = PR(IC) + PP(IC,I) 3 302 3800 CONTINUE 1 303 DO 3900 IC = 1,K 2 304 PR(IC) = PR(IC)/N 2 305 3900 CONTINUE 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAR 03, 1994 21:51:41 NAME:MAIN# PAGE: 11 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 1 306 DO 1750 IG = 1,K 2 307 DO 1750 JV = 1,IP 3 308 SUM(IG,JV) = 0.0D+00 3 309 DO 1750 JW = 1,IP 4 310 SSD(IG,JV,JW) = 0.0D+00 4 311 1750 CONTINUE C C UPDATE MEANS:-- 1 312 DO 1875 IC = 1,K 2 313 DO 1875 JV = 1,IP 3 314 DO 1875 I = 1,N 4 315 SUM(IC,JV) = SUM(IC,JV) + PP(IC,I)*X(I,JV) 4 316 1875 CONTINUE 1 317 DO 1900 IG = 1,K 2 318 IF (N*PR(IG) .LT. 0.5) GO TO 2050 2 319 GO TO 2150 2 320 2050 WRITE (6,74000) IG 2 321 WRITE (6,10050) C IF INSUFFICIENT WEIGHT IN ANY CLUSTER, GO ON TO NEXT K: 2 322 GO TO 850 2 323 2150 CONTINUE 2 324 DO 1900 JV = 1,IP 3 325 XBAR(IG,JV) = SUM(IG,JV)/(N*PR(IG)) 3 326 1900 CONTINUE C C C 1 327 CALL XUFLOW(0) 1 328 DO 3601 IC=1,K 2 329 DO 3601 JV=1,IP 3 330 DO 3601 JW=1,IP 4 331 DO 3601 I=1,N 5 332 DEVV=X(I,JV)-XBAR(IC,JV) 5 333 DEVW=X(I,JW)-XBAR(IC,JW) 5 334 IF ( PP(IC,I) .EQ. 0.0D+00 ) GO TO 3601 5 335 IF ( DEVV .EQ. 0.0D+00 ) GO TO 3601 5 336 IF ( DEVW .EQ. 0.0D+00 ) GO TO 3601 5 337 SSD(IC,JV,JW)=SSD(IC,JV,JW) + PP(IC,I)*DEVV*DEVW 5 338 3601 CONTINUE 1 339 DO 3600 IC=1,K 2 340 DO 3600 JV=1,IP 3 341 DO 3600 I=1,N 4 342 DEVV=X(I,JV)-XBAR(IC,JV) 4 343 IF ( PP(IC,I) .EQ. 0.0D+00 ) GO TO 3600 4 344 IF ( DEVV .EQ. 0.0D+00 ) GO TO 3600 4 345 SSD(IC,JV,JV)=SSD(IC,JV,JV) + PP(IC,I)*(DEVV**2) 4 346 3600 CONTINUE C C C C C POOL: C 1 347 DO 1950 JV = 1,IP 2 348 DO 1950 JW = 1,IP 3 349 WGSS(JV,JW) = 0.0D+00 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAR 03, 1994 21:51:41 NAME:MAIN# PAGE: 12 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 3 350 1950 CONTINUE C 1 351 DO 2300 JV = 1,IP 2 352 DO 2300 JW = 1,IP 3 353 DO 2300 IC = 1,K 4 354 WGSS(JV,JW) = WGSS(JV,JW) + PR(IC)*SSD(IC,JV,JW) 4 355 2300 CONTINUE C COMPUTE VARHAT, MLE OF COMMON COVARIANCE MATRIX: 1 356 DO 2400 JV = 1,IP 2 357 DO 2400 JW = 1,IP 3 358 VARHAT(JV,JW) = WGSS(JV,JW)/XN 3 359 2400 CONTINUE C C UPDATE PRECISION MATRIX AND DETERMINANT OF COVARIANCE MATRIX:-- C 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 360 DO 2402 JV = 1,IP 2 361 DO 2402 JW = 1,IP 3 362 COVMX(JV,JW) = VARHAT(JV,JW) 3 363 2402 CONTINUE 1 364 NRS1 = 0 1 365 IDET = 1 1 366 CALL MATEQ(COVMX, IP,20,JFLG,DET,IDET,IV,NRS1,P,20) C 1 367 IF (JFLG .GT. 0) WRITE (6,60000) JFLG 1 369 IF (JFLG .GT. 0) WRITE (6,62000) DET,IDET 1 371 TRUDET = DET*(10.0**IDET) C C (END PARAMETER-ESTIMATE UPDATE SEQUENCE) C C 1 372 IF (ITER .EQ. 1) GO TO 900 C STORE OLD LABELS: 1 373 DO 825 I = 1,N 2 374 ICLSOL(I) = ICLUS(I) 2 375 825 CONTINUE C 1 376 900 CONTINUE C C COMPUTE NEW LABELS BY MAX POSTERIOR PROBABILITY: 1 377 DO 1600 I = 1,N 2 378 XMXPR(I) = PP(1,I) 2 379 ICLUS(I) = 1 2 380 DO 1600 IC = 2,K 3 381 IF ( PP(IC,I) .GT. XMXPR(I) ) GO TO 1500 3 382 GO TO 1600 3 383 1500 XMXPR(I) = PP(IC,I) 3 384 ICLUS(I) = IC 3 385 1600 CONTINUE C 1 386 DO 1601 I = 1,N 2 387 DO 1601 IC = 1,K 3 388 PPOLD(IC,I) = PP(IC,I) 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAR 03, 1994 21:51:41 NAME:MAIN# PAGE: 13 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 3 389 1601 CONTINUE C 1 390 IF (N .GE. 31) GO TO 250 1 391 WRITE (6,14000) C 1 392 IF ( N .GE. 31 ) GO TO 250 1 393 WRITE (6,24000) 1 394 WRITE (6,18000) (ID(I),ICLUS(I), I=1,N) 1 395 250 CONTINUE C C C C 1 396 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. 1 397 DO 2900 I = 1,N 2 398 IF (ICLUS(I) .EQ. ICLSOL(I)) GO TO 2900 2 399 GO TO 3000 2 400 2900 CONTINUE 1 401 GO TO 3300 1 402 3000 CONTINUE 1 403 ITER = ITER + 1 1 404 IF (ITER.GE.21) GO TO 3100 C COUNT NUMBERS IN CLUSTERS: 1 405 DO 3311 IC = 1,K 2 406 NG(IC) = 0 2 407 3311 CONTINUE 1 408 DO 3321 I = 1,N 2 409 IC = ICLUS(I) 2 410 NG(IC) = NG(IC) + 1 2 411 3321 CONTINUE 1 412 WRITE (6,34000) (NG(IC),IC=1,K) C 1 413 GO TO 700 1 414 3100 WRITE (6,68000) C IF MORE THAN 20 ITERATIONS, GO ON TO NEXT K: 1 415 GO TO 850 C C C 1 416 3300 CONTINUE C C C COUNT NUMBERS IN CLUSTERS: 1 417 DO 3310 IC = 1,K 2 418 NG(IC) = 0 2 419 3310 CONTINUE 1 420 DO 3320 I = 1,N 2 421 IC = ICLUS(I) 2 422 NG(IC) = NG(IC) + 1 2 423 3320 CONTINUE C C C PRINT FINAL RESULTS FOR THIS VALUE OF K: 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAR 03, 1994 21:51:41 NAME:MAIN# PAGE: 14 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 C C 1 424 WRITE (6,34000) (NG(IC),IC=1,K) 1 425 WRITE (6,36000) (PR(IC),IC=1,K) 1 426 DO 2800 IC = 1,K 2 427 WRITE (6,28000) IC, (XBAR(IC,JV),JV=1,IP) 2 428 2800 CONTINUE C 1 429 WRITE (6,42000) 1 430 DO 2500 JV=1,IP 2 431 WRITE (6,48000) (VARHAT(JV,JW),JW=1,IP) 2 432 2500 CONTINUE C 1 433 IF ( N .GE. 31 ) GO TO 3650 1 434 WRITE (6,52000) 1 435 DO 3500 I = 1,N 2 436 WRITE (6,50000) ID(I), ICLUS(I), ( X(I,JV),JV=1,IP ) 2 437 3500 CONTINUE 1 438 3650 CONTINUE C C C C C C C C 1 439 WRITE (6,58000) ITER 1 440 WRITE (6,34000) (NG(IG),IG=1,K) 1 441 DO 3400 JV=1,IP 2 442 3400 CONTINUE C C C C COMPUTE VALUES OF MODEL-SELECTION CRITERIA: C C 1 443 WRITE (6,72000) NOPARM 1 444 AIC = XMN2LL + 2.0*NOPARM 1 445 SCH = XMN2LL + ALOG(XN)*NOPARM 1 446 WRITE (6,71000) SCH 1 447 WRITE (6,70000) AIC C C 1 448 WRITE (6,44000) 1 449 DO 851 IC=1,K 2 450 WRITE (6,45000) IC 2 451 DO 851 I = 1,N 3 452 IF ( ICLUS(I) .EQ. IC ) WRITE (6,46000) ID(I) 3 454 851 CONTINUE C C C GO ON TO NEXT NUMBER OF CLUSTERS, IF NOT FINISHED. C 1 455 WRITE (6,10050) 1 456 850 CONTINUE 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAR 03, 1994 21:51:41 NAME:MAIN# PAGE: 15 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 C C C COMPUTE VALUES OF MODEL-SELECTION CRITERIA AND CORRESPONDING C POSTERIOR PROBABILITIES OF THE MODELS: C C C SCHPPH(K) = PR(H(K))* C CONST*(2*PI)**(IPARM(K)/2)*EXP(-SCHVEC(K)/2)*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)*EXP(-SCHVEC(K)/2)*PRIOR PDF(K); C HERE PR(H(K)) IS TAKEN TO BE Const*EXP(-NOPARM(K)). C [SEE KASHYAP (1981).] C APPROXIMATE PRIOR PDF: 457 PRIRLN(1) = -(1*IP/2.0)*DLOG(2.0*PI) -IP/2.0 458 SCHPPH(1) = (IPARM(1)/2.0)*DLOG(2.0*PI) - SCHVEC(1)/2.0 459 SCHPPH(1) = SCHPPH(1) + PRIRLN(1) 460 AICPPH(1) = (IPARM(1)/2.0)*DLOG(2.0*PI) - SCHVEC(1)/2.0 461 AICPPH(1) = AICPPH(1) + PRIRLN(1) 462 AICPPH(1) = AICPPH(1) - IPARM(1) C ADD CONSTANT SCHVEC(1)/2 OR AICVEC(1)/2 TO PREVENT UNDERFLOW: 463 SCHPPH(1) = SCHPPH(1) + SCHVEC(1)/2 464 AICPPH(1) = AICPPH(1) + AICVEC(1)/2 465 DO 6110 K = 2,KUP C C APPROXIMATE PRIOR PDF: 1 466 PRIRLN(K) = -(K*IP/2.0)*DLOG(2.0*PI) -K*IP/2.0 C This is the prior on the k distribution parameters. C It is suggested by a Gaussian prior with cov mx = I, C maximized. C 1 467 SCHPPH(K) = (IPARM(K)/2.0)*DLOG(2.0*PI) - SCHVEC(K)/2.0 1 468 SCHPPH(K) = SCHPPH(K) + PRIRLN(K) 1 469 AICPPH(K) = (IPARM(K)/2.0)*DLOG(2.0*PI) - SCHVEC(K)/2.0 1 470 AICPPH(K) = AICPPH(K) + PRIRLN(K) 1 471 AICPPH(K) = AICPPH(K) - IPARM(K) C ADD CONSTANT SCHVEC(1)/2 OR AICVEC(1)/2 TO PREVENT UNDERFLOW: 1 472 SCHPPH(K) = SCHPPH(K) + SCHVEC(1)/2 1 473 AICPPH(K) = AICPPH(K) + AICVEC(1)/2 1 474 6110 CONTINUE C AT THIS POINT, WE HAVE THE LOGS OF THE PPHS; WE NOW CONVERT C THEM TO THE PPHS: 475 SCHPPH(1) = DEXP(SCHPPH(1)) 476 AICPPH(1) = DEXP(AICPPH(1)) 477 DO 6119 K = 2,KUP 1 478 IF ( SCHPPH(K) .LE. -174.673 ) GO TO 6115 1 479 SCHPPH(K) = DEXP(SCHPPH(K)) 1 480 GO TO 6119 1 481 6115 SCHPPH(K) = 0.0D+00 1 482 6119 CONTINUE 483 DO 6121 K = 2,KUP 1 484 IF ( AICPPH(K) .LE. -174.673 ) GO TO 6126 1 485 AICPPH(K) = DEXP(AICPPH(K)) 1 486 GO TO 6121 1 487 6126 AICPPH(K) = 0.0D+00 1 488 6121 CONTINUE 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAR 03, 1994 21:51:41 NAME:MAIN# PAGE: 16 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 489 SPPHS = 0.0D+00 490 SPPHA = 0.0D+00 491 DO 6120 K = 1,KUP 1 492 IF ( AICPPH(K) .GT. 0.0D+00 ) SPPHA = SPPHA + AICPPH(K) 1 494 6120 CONTINUE 495 DO 6127 K = 1,KUP 1 496 IF ( SCHPPH(K) .GT. 0.0D+00 ) SPPHS = SPPHS + SCHPPH(K) 1 498 6127 CONTINUE 499 DO 6130 K = 1,KUP 1 500 SCHPPH(K) = SCHPPH(K)/SPPHS 1 501 AICPPH(K) = AICPPH(K)/SPPHA 1 502 6130 CONTINUE C C C WRITE VALUES OF MSC'S FOR VARIOUS K: C 503 WRITE (6,11001) TITLE 504 WRITE(6,38000) 505 DO 605 K = 1,KUP 1 506 WRITE(6,33000) K,XM2LL(K),IPARM(K),AICVEC(K),AICPPH(K), 1 XSCHVEC(K),SCHPPH(K) 1 507 605 CONTINUE C C C C 508 WRITE (6,10050) 509 WRITE (6,11001) TITLE 510 WRITE (6,11097) C C WRITE (6,86000) C WRITE (6,11100) MONTH,IDAY,IYEAR,IHOUR,MINUTE,ISEC C CALL DTIME(IHOUR,MINUTE,ISEC) 511 CALL TDATE(IDAY,MONTH,IYEAR) C WRITE (6,89000) C WRITE (6,11100) MONTH,IDAY,IYEAR,IHOUR,MINUTE,ISEC 512 WRITE (6,98000) 513 STOP C C C ..............................................................C C C 514 10050 FORMAT(/' ................................................', X'................................'//) 515 11111 FORMAT(/' FIRST FOUR DATA POINTS') 516 11112 FORMAT(' ----------------------') 517 13000 FORMAT(/' NUMBER OF OBSERVATIONS (SAMPLE SIZE), N . . ',I3/) 518 19000 FORMAT(/' NUMBER OF VARIABLES . . . . . . . . . . . . ',I3 ) 519 10000 FORMAT('1','**********************************************'/// X' MIXPCMA CLUSPAC - MIXTURE-MODEL CLUSTERING OF CASES'// X' COMMON COVARIANCE MATRIX '/ X' AUTOMATIC SETTING OF INITIAL PARAMETER VALUES'// X' Copyright 1991 Stanley Louis Sclove. '/// X' Developed and programmed by: '// X19X,'Prof. Stanley L. Sclove, Ph.D. Phone (312) 996-2681'/ 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAR 03, 1994 21:51:41 NAME:MAIN# PAGE: 17 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 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: 4.6 3-Mar-94 '/ X' (VM/CMS) ') 520 10100 FORMAT(' MEAN VECTOR: ',(9F8.2/)//) 521 11097 FORMAT(/' MIXPCMA CLUSPAC - MIXTURE-MODEL CLUSTERING OF CASES'/) 522 10200 FORMAT(/' COVARIANCE MATRIX: '/) 523 21000 FORMAT(//' STATISTICS FOR WHOLE (UNCLUSTERED) SAMPLE: ') 524 22000 FORMAT(' ---------- --- ----- ------------- ------ '//) 525 37000 FORMAT(/' MINUS 2 LOG LIKELIHOOD =',F13.1 ) 526 10999 FORMAT(/' PROBLEM TITLE IS' ) 527 11000 FORMAT(1X,18A4/) 528 11001 FORMAT(/1X,18A4/) 529 15500 FORMAT(/1X,'MAXIMUM OF EACH VARIABLE: ',/) 530 16000 FORMAT(/1X,'MINIMUM OF EACH VARIABLE: ',/) 531 12000 FORMAT(2X,I4) 532 11100 FORMAT(1X,I2,'/',I2,'/',I4,3X,'AT ',I2,':',I2,':',I2/) 533 17000 FORMAT(3X,I2) 534 14000 FORMAT(//1X,'CLUSTERING:'/) 535 15000 FORMAT(18A4) 536 18000 FORMAT(1X, (7(A6,':',I2,2X)/) ) 537 24000 FORMAT(/,1X,'CASES AND LABELS:--'/) 538 20000 FORMAT(//1X,'INITIAL MEANS'/) 539 25000 FORMAT(1X, 'I=',I4,' PPS: ', 9F8.3) 540 26000 FORMAT(/,1X,'K = ',I2,' CLUSTERS'/) 541 28000 FORMAT(1X,'MEAN VECTOR FOR CLUSTER ',I2,': ',(8F13.5/)) 542 32000 FORMAT(/,1X,'ITERATION ', I2,/ ) C WRITE(6,33000) K,XM2LL(K),IPARM(K),AICVEC(K),AICPPH(K), C XSCHVEC(K),SCHPPH(K) 543 33000 FORMAT(5X,I1,6X,F8.1, 4X,I4,3X,F7.1,1X,F4.2,3X,F7.1,1X,F4.2) 544 38000 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 545 34000 FORMAT(/,1X,'NUMBERS IN CLUSTERS:',2X,9(I5,1X)/) 546 36000 FORMAT(/,1X,'MIXING PROBABILITIES:',3X,9(F4.2,1X)//) 547 42000 FORMAT(//, 1X,'COMMON COVARIANCE MATRIX (MLE):',//) 548 44000 FORMAT(/' LIST OF CASES, BY CLUSTER'//) 549 45000 FORMAT(' CLUSTER ',I2) 550 46000 FORMAT('+',1X,A6) 551 48000 FORMAT(1X,8F13.5/) 552 50000 FORMAT(1X,A6,':',I2,(8F13.3/)) 553 52000 FORMAT(//,1X,'CASE, LABEL / DATA'//) 554 58000 FORMAT(/1X,'CONVERGENCE: NO CASE CHANGED CLUSTERS AFTER ', X'ITERATION ',I2,'. RESULTS ARE PRINTED BELOW.'//) 555 60000 FORMAT(/,1X,'JFLG = ',I2,'. IF JFLG=0, COMPUTATION OF DET', X' WENT WELL; OTHERWISE, THERE WAS TROUBLE OR MATRIX WAS ', X'ILL-CONDITIONED.'//) 556 62000 FORMAT(/1X,'DET = ',E13.5,' IDET = ',I3,5X, 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAR 03, 1994 21:51:41 NAME:MAIN# PAGE: 18 0 IF DO ISN *....*...1.........2.........3.........4.........5.........6.........7.*.......8 0 X'ACTUAL DET. = DET*10**IDET',//) 557 68000 FORMAT(' ALGORITHM HAS NOT CONVERGED IN 20 ITERATIONS. STOP') 558 69000 FORMAT(/' MINUS 2 LOG LIKELIHOOD =',F13.1 ) 559 70000 FORMAT( ' AIC = ', F13.1 ) 560 71000 FORMAT( ' SCHWARZ CRITERION = ', F13.1/) 561 72000 FORMAT(1X, ' NUMBER OF PARAMETERS = ',I4 ) 562 74000 FORMAT(1X,'NO OBSERVATIONS IN GROUP ',I3,'. END PROCESSING ', X'FOR THIS NUMBER OF CLUSTERS.'/) 563 76000 FORMAT(' POSTERIOR PROBS OF GROUP MEMBERSHIP FOR 1ST 4 CASES:') C6000 FORMAT(/' TIME BEGAN: ') C9000 FORMAT(/' TIME FINISHED: ') 564 96000 FORMAT(1X,29F5.2) 565 98000 FORMAT(/1X,'PROGRAM ENDED NORMALLY.') 566 END 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAR 03, 1994 21:51:41 NAME:MAIN# PAGE: 19 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 173S 174F 177F 261S 263F 265F 444S 447F AICPPH R*8 AT 17 29 460S 461F 461S 462F 462S 464F 464S 469S 470F 470S 471F 471S 473F 473S 476F 476S 484F 485F 485S 487S 492F 493F 501F 501S 506F AICVEC R*4 A 17 177S 265S 464F 473F 506F ALOG R*4 I 175 262 445 COVMX R*8 AT 11 26 158S 160B 214S 216B 362S 366B DATFMT R*4 A 8 48S 67F 69F 201F DENOM R*8 AT 4 29 268S 270F 270S 273F 282F DET R*8 T 31 160B 164F 166F 216B 220F 366B 370F 371F DEVV R*8 T 34 233S 236F 332S 335F 337F 342S 344F 345F DEVW R*8 T 34 235S 236F 333S 336F 337F DEXP R*8 I 243 475 476 479 485 DLOG R*8 I 167 167 254 457 458 460 466 467 469 DSQ R*8 AT 5 32 231S 236F 236S 238F DSQRT R*8 I 243 F R*8 AT 2 35 240S 243S 251F 270F 282F I I*4 49S 50 50 51F 60F 61F 62F 63B 72S 73F 73B 130S 131B 135S 136F 136B 187S 189B 196S 198F 225S 226F 226B 229S 233F 235F 240F 243B 248S 251B 267S 268F 270F 270F 270B 272S 273F 274F 277F 277F 282F 282F 282B 289S 290B 294S 295F 295B 300S 301B 314S 315F 315B 331S 332F 333F 334F 337B 341S 342F 343F 345B 373S 374F 374B 377S 378F 378F 379F 381F 381F 383F 383F 384B 386S 388F 388B 394S 394F 394B 397S 398F 398B 408S 409B 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAR 03, 1994 21:51:41 NAME:MAIN# PAGE: 20 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 420S 421B 435S 436F 436F 436B 451S 452F 453B IC I*4 184S 185F 185B 188S 189B 250S 251F 251B 269S 270F 270B 276S 277F 277B 281S 282F 282F 282B 290S 290S 295S 295S 298S 299F 301F 301F 301B 303S 304F 304B 312S 315F 315F 315B 328S 332F 333F 334F 337F 337F 337B 339S 342F 343F 345F 345F 345B 353S 354F 354B 380S 381F 383F 384B 387S 388F 388B 405S 406B 409S 410F 410F 412S 412S 417S 418B 421S 422F 422F 424S 424S 425S 425S 426S 427F 427B 449S 450F 452B ICLSOL I*4 A 3 226S 374S 398F ICLUS I*4 A 3 226F 374F 379S 384S 394F 398F 409F 421F 436F 452F ID CHAR AT 1 22 50S 73F 274F 394F 436F 453F IDAY I*4 42B 511B IDET I*4 154S 160B 163F 164F 166F 210S 216B 219F 220F 365S 366B 370F 371F IG I*4 195S 196F 198B 200S 201B 306S 308F 310B 317S 318F 320F 325F 325F 325B 440S 440S INPFMT I*4 A 8 47S 50 73F INT GI 196 IP I*4 44S 45F 50 54F 59F 67F 69F 73F 124F 126F 129F 133F 134F 138F 141F 142F 149F 151F 152F 156F 157F 160B 167F 167F 171F 171F 171F 182F 182F 182F 192F 197F 201F 203F 204F 207F 212F 213F 216B 232F 234F 243F 259F 259F 259F 307F 309F 313F 324F 329F 330F 340F 347F 348F 351F 352F 356F 357F 360F 361F 366B 427F 430F 431F 436F 441F 457F 457F 466F 466F IPARM I*4 A 16 172S 260S 458F 460F 462F 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAR 03, 1994 21:51:41 NAME:MAIN# PAGE: 21 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 467F 469F 471F 506F ITER I*4 221S 223F 224F 372F 396F 403F 403S 404F 439F IV I*4 A 13 160B 216B 366B IYEAR I*4 42B 511B JFLG I*4 160B 161F 162F 165F 216B 217F 218F 366B 367F 368F 369F JV I*4 50S 50S 54S 55F 55F 56F 56B 59S 60F 60F 61F 61F 62F 62F 63F 63B 67S 67S 69S 69S 73S 73S 124S 125F 127B 129S 131F 131F 131B 133S 136F 136F 136B 138S 139F 139B 141S 143F 143F 143F 144F 144B 149S 149S 151S 152B 156S 158F 158B 192S 193F 193B 197S 198F 198B 201S 201S 203S 205B 207S 208F 208B 212S 214F 214B 232S 233F 233F 236B 307S 308F 310B 313S 315F 315F 315B 324S 325F 325B 329S 332F 332F 337F 337B 340S 342F 342F 345F 345F 345F 345B 347S 349B 351S 354F 354F 354B 356S 358F 358B 360S 362F 362B 427S 427S 430S 431B 436S 436S 441B JW I*4 126S 127B 134S 136F 136F 136B 142S 143F 143F 143F 144F 144B 152S 152S 157S 158F 158B 204S 205B 213S 214F 214B 234S 235F 235F 236B 309S 310B 330S 333F 333F 337F 337B 348S 349B 352S 354F 354F 354B 357S 358F 358B 361S 362F 362B 431S 431S K I*4 180S 181F 182F 182F 184F 185F 188F 189F 195F 196F 200F 230F 250F 257F 259F 259F 260F 265F 266F 269F 276F 281F 290F 295F 298F 303F 306F 312F 317F 328F 339F 353F 380F 387F 405F 412F 417F 424F 425F 426F 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAR 03, 1994 21:51:41 NAME:MAIN# PAGE: 22 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 440F 449B 465S 466F 466F 466F 467F 467F 467F 468F 468F 468F 469F 469F 469F 470F 470F 470F 471F 471F 471F 472F 472F 473F 473B 477S 478F 479F 479F 481B 483S 484F 485F 485F 487B 491S 492F 493B 495S 496F 497B 499S 500F 500F 501F 501B 505S 506F 506F 506F 506F 506F 506F 506B KUP I*4 75S 76F 180F 465F 477F 483F 491F 495F 499F 505F L I*4 230S 231F 233F 235F 236F 236F 238F 240F 243B MATEQ X 160F 216F 366F MIN GI 75 MONTH I*4 42B 511B N I*4 43S 46F 49F 75F 78F 130F 135F 187F 196F 225F 229F 248F 267F 272F 286F 294F 300F 304F 314F 318F 325F 331F 341F 373F 377F 386F 390F 392F 394F 397F 408F 420F 433F 435F 451F NG I*4 A 7 406S 410F 410S 412F 418S 422F 422S 424F 440F NOPARM I*4 171S 172F 173F 175F 182S 183F 259S 260F 261F 262F 443F 444F 445F NRS1 I*4 155S 160B 211S 216B 364S 366B P R*8 AT 14 30 160B 216B 236F 366B PI R*8 T 25 79S 167F 243F 457F 458F 460F 466F 467F 469F PP R*8 AT 2 29 277S 282S 290F 295F 301F 315F 334F 337F 343F 345F 378F 381F 383F 388F PPOLD R*8 AT 2 29 189S 277F 388S PR R*8 AT 16 29 185S 251F 270F 282F 299S 301F 301S 304F 304S 318F 325F 354F 425F PRIRLN R*8 AT 17 29 457S 459F 461F 466S 468F 470F SCH R*4 175S 176F 178F 262S 264F 266F 445S 446F SCHPPH R*8 AT 17 29 458S 459F 459S 463F 463S 467S 468F 468S 472F 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAR 03, 1994 21:51:41 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 472S 475F 475S 478F 479F 479S 481S 496F 497F 500F 500S 506F SCHVEC R*4 A 17 178S 266S 458F 460F 463F 467F 469F 472F 506F SPPHA R*4 490S 493F 493S 501F SPPHS R*4 489S 497F 497S 500F SSD R*8 AT 9 24 310S 337F 337S 345F 345S 354F SSDEVS R*8 AT 21 27 143S 144F SUM R*8 AT 1 23 308S 315F 315S 325F SUMLNF R*8 T 28 247S 254F 254S 256F SUMPXF R*8 T 28 249S 251F 251S 253F 254F SUMSQS R*8 AT 20 27 127S 136F 136S 143F TDATE X 42F 511F TITLE R*4 A 6 36S 41F 503F 509F TOTAL R*8 AT 18 27 125S 131F 131S 139F 143F 143F TRUDET R*8 T 31 164S 167F 220S 243F 371S VARHAT R*8 AT 11 26 144S 152F 158F 205S 208S 214F 358S 362F 431F WGSS R*8 AT 10 24 349S 354F 354S 358F X R*4 A 1 50S 55F 56F 60F 61F 62F 63F 73F 131F 136F 136F 193F 198F 233F 235F 315F 332F 333F 342F 436F XBAR R*8 AT 7 33 193S 198S 201F 233F 235F 325S 332F 333F 342F 427F XIDET R*4 163S 219S XJ R*4 A 15 80S 81S 82S 83S 84S 85S 86S 87S 88S 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 185F XLL R*4 167S 168F XMAX R*4 A 12 55S 62F 63S 69F XMEAN R*4 A 19 139S 149F XMIN R*4 A 12 56S 60F 61S 67F XMN2LL R*4 168S 169F 170F 173F 175F 256S 257F 258F 261F 262F 444F 445F XMXPR R*4 A 4 378S 381F 383S XM2LL R*4 A 16 169S 257S 506F 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAR 03, 1994 21:51:41 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) +______ ____ ______ ________ ________________________________________ 0XN R*4 78S 139F 143F 144F 167F 167F 167F 175F 262F 358F 445F XUFLOW X 327F ZSQ R*8 T 31 238S 239F 243F 0 0VARIABLES REFERENCED BUT NOT SET. (* POSSIBLY SET AS ARGUMENT.) 0DET* IDAY* IV* IYEAR* JFLG* MONTH* P* 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAR 03, 1994 21:51:41 NAME:MAIN# PAGE: 25 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 53 51 105 128 124 126 110 132 129 130 111 74 72 120 137 133 134 135 130 140 138 140 145 141 142 150 153 151 200 57 54 250 B 395 390 392 300 B 58 52 360 186 184 400 64 59 500 65 49 549 190 187 188 550 194 192 555 199 195 197 600 202 200 605 507 505 610 206 203 204 620 209 207 700 B 222 413 800 227 225 825 375 373 850 B 456 180 322 415 851 454 449 451 900 B 376 372 910 B 228 224 1000 245 230 1090 B 242 239 1100 B 244 241 1200 246 229 1310 237 232 234 1350 271 267 269 1400 285 272 1401 B 284 279 1402 283 281 1410 B 280 273 1420 291 289 1421 B 293 286 286 1422 B 287 286 1423 B 297 292 1424 296 294 1456 278 276 1500 B 383 381 1600 B 385 377 380 382 1601 389 386 387 1750 311 306 307 309 1875 316 312 313 314 1900 326 317 324 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAR 03, 1994 21:51:41 NAME:MAIN# PAGE: 26 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 1950 350 347 348 2050 B 320 318 2100 252 250 2150 B 323 319 2200 B 255 248 253 2300 355 351 352 353 2400 359 356 357 2402 363 360 361 2404 159 156 157 2406 215 212 213 2500 432 430 2800 428 426 2900 B 400 397 398 3000 B 402 396 399 3100 B 414 404 3300 B 416 401 3310 419 417 3311 407 405 3320 423 420 3321 411 408 3400 442 441 3500 437 435 3600 B 346 339 340 341 343 344 3601 B 338 328 329 330 331 334 335 336 3650 B 438 433 3800 302 298 300 3900 305 303 6110 474 465 6115 B 481 478 6119 B 482 477 480 6120 494 491 6121 B 488 483 486 6126 B 487 484 6127 498 495 6130 502 499 10000 NF 519 38 10050 NF 514 37 39 146 179 321 455 508 10100 NF 520 149 10200 NF 522 150 10999 NF 526 40 11000 NF 527 41 11001 NF 528 503 509 11096 NF 77 76 11097 NF 521 510 11100 NF 532 UNREFERENCED 11111 NF 515 70 11112 NF 516 71 12000 NF 531 43 13000 NF 517 46 14000 NF 534 391 14560 NF 275 274 15000 NF 535 36 47 48 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAR 03, 1994 21:51:41 NAME:MAIN# PAGE: 27 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 15500 NF 529 68 16000 NF 530 66 17000 NF 533 44 18000 NF 536 394 19000 NF 518 45 20000 NF 538 191 21000 NF 523 147 22000 NF 524 148 24000 NF 537 393 25000 NF 539 295 26000 NF 540 181 28000 NF 541 427 32000 NF 542 223 33000 NF 543 506 34000 NF 545 412 424 440 36000 NF 546 425 37000 NF 525 170 38000 NF 544 504 42000 NF 547 429 44000 NF 548 448 45000 NF 549 450 46000 NF 550 453 48000 NF 551 152 431 50000 NF 552 436 52000 NF 553 434 58000 NF 554 439 60000 NF 555 162 218 368 62000 NF 556 166 370 68000 NF 557 414 69000 NF 558 258 70000 NF 559 174 263 447 71000 NF 560 176 264 446 72000 NF 561 183 443 74000 NF 562 320 76000 NF 563 288 96000 NF 564 290 98000 NF 565 512 0*STATISTICS* SOURCE STATEMENTS: 556, PROGRAM SIZE: 946332 BYTES, PROGRAM NAME: MAIN#, PAGE: 1 *STATISTICS* NO DIAGNOSTICS GENERATED. **MAIN#** END OF COMPILATION 1 ****** TIME STAMP: 94.06221.51.41 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAR 03, 1994 21:51:45 PAGE: 28 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 MAR 03, 1994 21:51:45 NAME:MAIN# PAGE: 29 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 MAR 03, 1994 21:51:45 NAME:MATEQ PAGE: 30 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 MAR 03, 1994 21:51:45 NAME:MATEQ PAGE: 31 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 ANIHALATED. C 2 45 IF (A(IV(IM),K).EQ.0.0D+00) GO TO 1100 C C CACULATE 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 MAR 03, 1994 21:51:45 NAME:MATEQ PAGE: 32 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 MAR 03, 1994 21:51:45 NAME:MATEQ PAGE: 33 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 MAR 03, 1994 21:51:45 NAME:MATEQ PAGE: 34 0*STATISTICS* SOURCE STATEMENTS: 75, PROGRAM SIZE: 3564 BYTES, PROGRAM NAME: MATEQ, PAGE: 28 *STATISTICS* NO DIAGNOSTICS GENERATED. **MATEQ** END OF COMPILATION 2 ****** TIME STAMP: 94.06221.51.45 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAR 03, 1994 21:51:46 PAGE: 35 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 MAR 03, 1994 21:51:46 NAME:MATDT PAGE: 36 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: 35 *STATISTICS* NO DIAGNOSTICS GENERATED. **MATDT** END OF COMPILATION 3 ****** TIME STAMP: 94.06221.51.46 1LEVEL 2.5.0 (JUNE 1991) VS FORTRAN MAR 03, 1994 21:51:46 PAGE: 37 0SUMMARY OF MESSAGES AND STATISTICS FOR ALL COMPILATIONS 0*STATISTICS* SOURCE STATEMENTS: 556, PROGRAM SIZE: 946332 BYTES, PROGRAM NAME: MAIN#, PAGE: 1 *STATISTICS* NO DIAGNOSTICS GENERATED. **MAIN#** END OF COMPILATION 1 ****** TIME STAMP: 94.06221.51.46 0*STATISTICS* SOURCE STATEMENTS: 75, PROGRAM SIZE: 3564 BYTES, PROGRAM NAME: MATEQ, PAGE: 28 *STATISTICS* NO DIAGNOSTICS GENERATED. **MATEQ** END OF COMPILATION 2 ****** TIME STAMP: 94.06221.51.46 0*STATISTICS* SOURCE STATEMENTS: 31, PROGRAM SIZE: 1552 BYTES, PROGRAM NAME: MATDT, PAGE: 35 *STATISTICS* NO DIAGNOSTICS GENERATED. **MATDT** END OF COMPILATION 3 ****** TIME STAMP: 94.06221.51.46 0******* SUMMARY STATISTICS ******* 0 DIAGNOSTICS GENERATED. HIGHEST SEVERITY CODE IS 0.