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 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 DIMENSION X(9999,20),SUM(29,20),ID(9999) C DIMENSION Y(N=IP+1,IP) DIMENSION Y(21,20) DIMENSION F(9999,29),PP(29,9999) DIMENSION ICLUS(9999),ICLSOL(9999) DIMENSION DENOM(9999),XMXPR(9999) DIMENSION SSD(29,20,20),SIGMA(29,20,20) DIMENSION PREC(29,20,20),CMPREC(20,20) DIMENSION DSQ(29) DIMENSION TITLE(18) DIMENSION NG(29),XBAR(29,20) DIMENSION DATFMT(18),INPFMT(18) C DIMENSION XIDET(29),XLGDET(29) DIMENSION TRUDET(29) C DIMENSION WGSS(20,20) DIMENSION VARHAT(20,20),COVMX(20,20) DIMENSION XMIN(20),XMAX(20) DIMENSION XJ(9,9) DIMENSION PR(29),XM2LL(9) DIMENSION AICVEC(9),SCHVEC(9),AICPPH(9),SCHPPH(9),IPARM(9) DIMENSION PPOLD(29,9999) DIMENSION PRIRLN(9) DIMENSION TOTAL(29) DIMENSION XMEAN(29) DIMENSION SUMSQS(20,20) DIMENSION SSDEVS(20,20) DIMENSION P(20,20),A(20,20) DIMENSION IV(20,20) C C IV IS A WORK ARRAY FOR SUBROUTINE MATEQ. C CHARACTER*6 ID C DOUBLE PRECISION SUM DOUBLE PRECISION WGSS,SSD,TERM,SSDEVS,SUMSQS,TOTAL DOUBLE PRECISION DEVV,DEVW,PREC,CMPREC DOUBLE PRECISION VARHAT,COVMX,SDJV,SDJW DOUBLE PRECISION P,A,SIGMA DOUBLE PRECISION DET C DOUBLE PRECISION XLGDET DOUBLE PRECISION DSQ,ZSQ DOUBLE PRECISION XBAR DOUBLE PRECISION F,PR DOUBLE PRECISION SUMPXF, SUMLNF DOUBLE PRECISION SCHPPH,AICPPH,SPPHA,SPPHS,PP DOUBLE PRECISION DENOM DOUBLE PRECISION PPOLD DOUBLE PRECISION PRIRLN DOUBLE PRECISION TRUDET,COMDET DOUBLE PRECISION TEST C DOUBLE PRECISION XLIKLY C 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 READ (5,15000) TITLE C C WRITE PROGRAM INFORMATION. WRITE (6,10050) WRITE (6,10000) WRITE (6,10050) WRITE (6,10999) WRITE (6,11001) TITLE C C THE FOLLOWING ARE IMSL DATE AND TIME FUNCTIONS: 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. READ (5,12000) N C READ NUMBER OF VARIABLES, IP. READ (5,54000) IP WRITE (6,56000) IP WRITE (6,13000) N C C READ INPUT FORMAT (INPUT INCLUDES DATA AND ID.) READ (5,15000) INPFMT C READ DATA FORMAT. READ (5,15000) DATFMT C C READ DATA. DO 500 I = 1,N READ (5,INPFMT) (X(I,JV), JV = 1,IP),ID(I) IF (I .EQ. 1) GO TO 100 GO TO 300 100 CONTINUE DO 200 JV = 1,IP XMAX(JV) = X(1,JV) XMIN(JV) = X(1,JV) 200 CONTINUE 300 CONTINUE DO 400 JV = 1,IP IF (X(I,JV) .LT. XMIN(JV)) XMIN(JV)=X(I,JV) IF (X(I,JV) .GT. XMAX(JV)) XMAX(JV)=X(I,JV) 400 CONTINUE 500 CONTINUE WRITE (6,19000) WRITE (6,DATFMT) (XMIN(JV),JV=1,IP) WRITE (6,21000) WRITE (6,DATFMT) (XMAX(JV),JV=1,IP) C WRITE (6,11111) WRITE (6,11112) DO 111 I=1,4 WRITE (6,INPFMT) (X(I,JV),JV=1,IP),ID(I) C KUP = MIN(9,N) KREACH = KUP 111 CONTINUE WRITE (6,11096) KUP C C SET CONSTANTS: C XN = N PI = .314159265358979324D+01 NGFLAG = 0 C C SET INITIAL VALUES OF MIXING PROBABILITIES:-- C C XJ(IC,K), K=2 TO 9, IC = 1 TO K C OPTIMAL CLASS PROBABILITIES - NORMAL DISTRIBUTION C JOHARI & SCLOVE (1975). "PARTITIONING A DISTRIBUTION." C COMMUNICATIONS IN STATISTICS. C XJ(1,2)=.5 XJ(2,2)=.5 XJ(1,3)=.2703 XJ(2,3)=.4594 XJ(3,3)=.2703 XJ(1,4)=.1631 XJ(2,4)=.3369 XJ(3,4)=.3369 XJ(4,4)=.1631 XJ(1,5)=.1068 XJ(2,5)=.2444 XJ(3,5)=.2976 XJ(4,5)=.2444 XJ(5,5)=.1068 XJ(1,6)=.0739 XJ(2,6)=.1810 XJ(3,6)=.2451 XJ(4,6)=.2451 XJ(5,6)=.1810 XJ(6,6)=.0739 XJ(1,7)=.0536 XJ(2,7)=.1375 XJ(3,7)=.1986 XJ(4,7)=.2106 XJ(5,7)=.1986 XJ(6,7)=.1375 XJ(7,7)=.0536 XJ(1,8)=.0402 XJ(2,8)=.1067 XJ(3,8)=.1613 XJ(4,8)=.1918 XJ(5,8)=.1918 XJ(6,8)=.1613 XJ(7,8)=.1067 XJ(8,8)=.0402 XJ(1,9)=.0310 XJ(2,9)=.0845 XJ(3,9)=.1324 XJ(4,9)=.1643 XJ(5,9)=.1756 XJ(6,9)=.1643 XJ(7,9)=.1324 XJ(8,9)=.0845 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: DO 105 JV = 1,IP TOTAL(JV) = 0.0 DO 105 JW = 1,IP SUMSQS(JV,JW) = 0.0 105 CONTINUE DO 110 JV = 1,IP DO 110 I = 1,N TOTAL(JV) = TOTAL(JV) + X(I,JV) 110 CONTINUE DO 120 JV = 1,IP DO 120 JW = 1,IP DO 120 I = 1,N SUMSQS(JV,JW) = SUMSQS(JV,JW) + X(I,JV)*X(I,JW) 120 CONTINUE DO 130 JV = 1,IP XMEAN(JV) = TOTAL(JV)/XN 130 CONTINUE DO 140 JV = 1,IP DO 140 JW = 1,IP SSDEVS(JV,JW) = SUMSQS(JV,JW)- TOTAL(JV)*TOTAL(JW)/XN VARHAT(JV,JW) = SSDEVS(JV,JW)/XN 140 CONTINUE C WRITE SUMMARY STATISTICS FOR WHOLE SAMPLE: WRITE(6,10500) WRITE(6,10100) ( XMEAN(JV), JV = 1,IP ) WRITE(6,10200) DO 150 JV = 1,IP WRITE(6,48000) ( VARHAT(JV,JW), JW=1,IP ) 150 CONTINUE C C C CALL SUBROUTINE TO COMPUTE DETERMINANT OF SAMPLE COVARIANCE C MATRIX: C IDET = 1 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: DO 2402 JV = 1,IP DO 2402 JW = 1,IP COVMX(JV,JW) = VARHAT(JV,JW) 2402 CONTINUE 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 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) TRUDET(1) = DET*(10.0**IDET) 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 XLL=-(XN*IP/2.0)*DLOG(2.0*PI)-(XN/2.0)*DLOG(TRUDET(1))-XN*IP/2.0 XMN2LL = (-2.0)*XLL XM2LL(1) = XMN2LL WRITE(6,10900) XMN2LL C NO. OF PARAMETERS FOR UNCLUSTERED SAMPLE IS: C 1 MEAN VECTOR + 1 COV MX. C NOPARM = IP + IP*(IP+1)/2 IPARM(1) = NOPARM C WRITE (6,72000) NOPARM AIC = XMN2LL + 2.0*NOPARM WRITE(6,70000) AIC SCH = XMN2LL + ALOG(XN)*NOPARM WRITE (6,78000) SCH AICVEC(1) = AIC SCHVEC(1) = SCH C WRITE (6,10050) C C C PERFORM CLUSTERING FOR VARIOUS NUMBERS OF CLUSTERS, K: C DO 850 K = 2,KUP 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 NOPARM = K*IP + K*IP*(IP+1)/2 + (K-1) IPARM(K) = NOPARM C 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 DO 360 IC = 1,K PR(IC) = XJ(IC,K) 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: DO 549 I = 1,N DO 549 IC = 1,K PPOLD(IC,I) = 1.0/K 549 CONTINUE C WRITE (6,20000) C C Set initial means equal to K of the observations:-- C DO 550 IG = 1,K C 1st observation for IG=1, N-th observation for IG=K: 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) ) DO 550 JV = 1,IP XBAR(IG,JV) = X(I,JV) 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 DO 600 IG=1,K WRITE (6,21999) WRITE (6,22000) IG WRITE (6,DATFMT) ( XBAR(IG,JV), JV=1,IP ) 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 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 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: C 113 CONTINUE C DO 630 IC = 1,K DO 610 JV=1,IP DO 610 JW=1,IP 114 SIGMA(IC,JV,JW) = 0.00 610 CONTINUE DO 620 JV = 1,IP SIGMA(IC,JV,JV) = 1.0 620 CONTINUE 630 CONTINUE C SCALE UP THIS CORRELATION MATRIX TO A COVARIANCE MATRIX C USING THE WHOLE-SAMPLE VARIANCES: DO 611 JV=1,IP SDJV = DSQRT(VARHAT(JV,JV)) DO 611 JW=1,IP SDJW = DSQRT(VARHAT(JW,JW)) DO 611 IC = 1,K SIGMA(IC,JV,JW) = SDJV*SIGMA(IC,JV,JW)*SDJW 611 CONTINUE C C COMPUTE INITIAL INVERSE COVARIANCE MATRIX: C C DO 640 IC=1,K C C IDET = 1 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: DO 2404 JV = 1,IP DO 2404 JW = 1,IP COVMX(JV,JW) = SIGMA(IC,JV,JW) 2404 CONTINUE C 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 IF (JFLG .GT. 0) WRITE (6,60000) JFLG C C C SET INITIAL PRECISION MATRICES FOR THIS K: DO 640 JV = 1,IP DO 640 JW = 1,IP PREC(IC,JV,JW) = P(JV,JW) 640 CONTINUE C WRITE(6,11500) DET,IDET DO 680 IC = 1,K C XIDET(IC) = IDET C XLGDET(IC) = DLOG(DET) + XIDET(IC)*DLOG(10.0D+00) TRUDET(IC) = DET*(10.0**IDET) 680 CONTINUE C END PARAMETER INITIALIZATION FOR THIS K; C BEGIN ITERATIONS: C C ITER = 1 700 CONTINUE WRITE(6,32000) ITER IF (ITER .EQ. 1) GO TO 901 DO 800 I = 1,N ICLSOL(I) = ICLUS(I) 800 CONTINUE C 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 DO 1200 I = 1,N DO 1000 L=1,K DSQ(L) = 0.0D+00 DO 1310 JV=1,IP DEVV = XBAR(L,JV) - X(I,JV) DO 1310 JW=1,IP DEVW = XBAR(L,JW) - X(I,JW) DSQ(L) = DSQ(L) + DEVV*PREC(L,JV,JW)*DEVW 1310 CONTINUE 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): IF ( ZSQ/2.0 .LE. 174.673D+00 ) GO TO 1090 F(I, L) = 0.0D+00 GO TO 1100 C 1090 CONTINUE C Compute C F(I,L) = DEXP(-ZSQ/2.0)/((2*PI)**(IP/2)*DSQRT(TRUDET(L))): C XIP = IP F(I,L) = (-ZSQ/2.0)-(XIP/2.0)*DLOG(2*PI)-0.5*DLOG(TRUDET(L)) IF(F(I,L) .GT. 174.673) WRITE(6,87000) ZSQ,TRUDET(L) F(I,L) = DEXP( F(I,L) ) 1100 CONTINUE 1000 CONTINUE 1200 CONTINUE C C COMPUTE LOG LIKELIHOOD AND VALUES OF MODEL-SELECTION CRITERIA: C C XLIKLY = 1.0 SUMLNF = 0.0 DO 2200 I=1,N SUMPXF = 0.0 DO 2100 IC=1,K SUMPXF = SUMPXF + PR(IC)*F(I,IC) 2100 CONTINUE C C XLIKLY = XLIKLY*SUMPXF IF ( SUMPXF .LE. 0.0 ) GO TO 2190 SUMLNF = SUMLNF + DLOG(SUMPXF) 2190 CONTINUE 2200 CONTINUE C C SUMLNF = DLOG(XLIKLY) XMN2LL = -2.0*SUMLNF XM2LL(K) = XMN2LL C WRITE (6,10900) XMN2LL C C C COMPUTE MODEL SELECTION CRITERIA AIC = XMN2LL + 2.0*NOPARM SCH = XMN2LL + ALOG(XN)*NOPARM WRITE (6,70000) AIC WRITE (6,78000) SCH AICVEC(K) = AIC 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): DO 1350 I = 1,N DENOM(I) = 0.0 DO 1350 IC=1,K DENOM(I) = DENOM(I) + PR(IC)*F(I,IC) 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 DO 1400 I = 1,N IF ( DENOM(I) .GT. 0.0D+00 ) GO TO 1410 WRITE (6,14560) ID(I) DO 1456 IC = 1,K PP(IC,I) = PPOLD(IC,I) 1456 CONTINUE GO TO 1401 C 1410 CONTINUE DO 1402 IC=1,K PP(IC,I) = PR(IC)*F(I,IC)/DENOM(I) 1402 CONTINUE C 1401 CONTINUE 1400 CONTINUE C WRITE (6,55555) IF ( N-15 ) 1421,1421,1422 1422 CONTINUE WRITE(6,76000) DO 1420 I = 1,4 WRITE(6,82000) ( PP(IC,I), IC = 1,K ) 1420 CONTINUE GO TO 1423 1421 CONTINUE DO 1424 I = 1,N WRITE (6,25000) ID(I), (PP(IC,I),IC=1,K) 1424 CONTINUE 1423 CONTINUE C C C UPDATE PARAMETER ESTIMATES: C C C UPDATE CLUSTER PRIOR PROBABILITIES PR(IC):-- DO 3800 IC = 1,K PR(IC) = 0.0 DO 3800 I = 1,N PR(IC) = PR(IC) + PP(IC,I) 3800 CONTINUE DO 3900 IC = 1,K PR(IC) = PR(IC)/N 3900 CONTINUE C IF MIXING PROBABILITY IS ZERO FOR ANY CLUSTER, STOP. DO 3901 IC = 1,K IF ( PR(IC) .GT. 0.0D+00 ) GO TO 3901 WRITE (6,40040) K,IC 40040 FORMAT(/' IN FITTING ',I1,' CLUSTERS, PR(',I1,')=0. STOP.'/) GO TO 4004 3901 CONTINUE DO 1750 IG = 1,K DO 1750 JV = 1,IP SUM(IG,JV) = 0.0 DO 1750 JW = 1,IP SSD(IG,JV,JW) = 0.0 1750 CONTINUE C C UPDATE MEANS:-- DO 1875 IC = 1,K DO 1875 JV = 1,IP DO 1875 I = 1,N SUM(IC,JV) = SUM(IC,JV) + PP(IC,I)*X(I,JV) 1875 CONTINUE DO 1900 IG = 1,K TEST = XN*PR(IG) IF ( TEST .EQ. 0.0D+00 ) GO TO 2050 C IF (N*PR(IG) .LT. 0.5) GO TO 2050 GO TO 2150 2050 WRITE (6,74000) IG NGFLAG = 1 C IF EXP.NO. OF OBSERVATIONS FOR A CLUSTER IS LESS THAN 0.5, C END COMPUTATION FOR THIS NUMBER OF CLUSTERS C AND GO ON TO NEXT VALUE OF K: GO TO 850 2150 CONTINUE DO 1900 JV = 1,IP XBAR(IG,JV) = SUM(IG,JV)/(XN*PR(IG)) 1900 CONTINUE C C CALL XUFLOW(0) C DO 3600 IC=1,K DO 3600 JV=1,IP DO 3600 JW=JV,IP DO 3600 I=1,N DEVV=X(I,JV)-XBAR(IC,JV) DEVW=X(I,JW)-XBAR(IC,JW) IF ( PP(IC,I) .EQ. 0.0 ) GO TO 3600 IF ( DEVV .EQ. 0.0 ) GO TO 3600 IF ( DEVW .EQ. 0.0 ) GO TO 3600 SSD(IC,JV,JW)=SSD(IC,JV,JW) + PP(IC,I)*DEVV*DEVW 3600 CONTINUE C DO 3700 IC=1,K DO 3700 JV=2,IP DO 3700 JW=1,JV-1 SSD(IC,JV,JW)=SSD(IC,JW,JV) 3700 CONTINUE C DO 3710 IC = 1,K DO 3710 JV = 1,IP DO 3710 JW = 1,IP SIGMA(IC,JV,JW) = SSD(IC,JV,JW)/(N*PR(IC)) 3710 CONTINUE C C C POOL: C DO 1950 JV = 1,IP DO 1950 JW = 1,IP WGSS(JV,JW) = 0.0 1950 CONTINUE C DO 2300 JV = 1,IP DO 2300 JW = 1,IP DO 2300 IC = 1,K WGSS(JV,JW) = WGSS(JV,JW) + PR(IC)*SSD(IC,JV,JW) 2300 CONTINUE C COMPUTE VARHAT, MLE OF COMMON COVARIANCE MATRIX: DO 2400 JV = 1,IP DO 2400 JW = 1,IP VARHAT(JV,JW) = WGSS(JV,JW)/N 2400 CONTINUE C C UPDATE PRECISION MATRICES: C NRS1 = 0 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: DO 2406 JV = 1,IP DO 2406 JW = 1,IP COVMX(JV,JW) = VARHAT(JV,JW) 2406 CONTINUE C COMPUTE INVERSE OF POOLED COVARIANCE MATRIX:-- CALL MATEQ(COVMX, IP,20,JFLG,DET,IDET,IV,NRS1,P,20) DO 3813 JV = 1,IP DO 3813 JW = 1,IP CMPREC(JV, JW)=P(JV,JW) 3813 CONTINUE IF (JFLG .GT. 0) WRITE (6,60000) JFLG IF (JFLG .GT. 0) WRITE (6,62000) DET,IDET COMDET = DET*(10.0**IDET) C C C UPDATE PRECISION MATRICES: C DO 3711 IC = 1,K DO 3712 JV = 1,IP DO 3712 JW = 1,IP A(JV,JW) = SIGMA(IC,JV,JW) 3712 CONTINUE C IDET = 1 NRS1 = 0 C CALL MATEQ(A,IP,20,JFLG,DET,IDET,IV,NRS1,P,20) C DO 3713 JV = 1,IP DO 3713 JW = 1,IP PREC(IC,JV,JW)=P(JV,JW) 3713 CONTINUE C WRITE(6,11500) DET,IDET IF (JFLG .GT. 0) WRITE (6,60000) JFLG IF (JFLG .GT. 0) WRITE (6,62000) DET,IDET C XIDET(IC) = IDET C XLGDET(IC) = DLOG(DET) + XIDET(IC)*DLOG(10.0D+00) TRUDET(IC) = DET*(10.0**IDET) C IF DET=0, REPLACE COV MX WITH COMMON COV MX: IF ( TRUDET(IC) .GT. 0.0D+00 ) GO TO 3711 TRUDET(IC) = COMDET WRITE (6, 17000) IC DO 3814 JV = 1,IP DO 3814 JW = 1,IP PREC(IC,JV,JW) = CMPREC(JV,JW) 3814 CONTINUE 3711 CONTINUE C C C (END PARAMETER-ESTIMATE UPDATE SEQUENCE) C C IF (ITER .EQ. 1) GO TO 900 C STORE OLD LABELS: DO 825 I = 1,N ICLSOL(I) = ICLUS(I) 825 CONTINUE C 900 CONTINUE C C COMPUTE NEW LABELS BY MAX POSTERIOR PROBABILITY: C DO 1600 I = 1,N XMXPR(I) = PP(1,I) ICLUS(I) = 1 DO 1600 IC = 2,K IF ( PP(IC,I) .GT. XMXPR(I) ) GO TO 1500 GO TO 1600 1500 XMXPR(I) = PP(IC,I) ICLUS(I) = IC 1600 CONTINUE C DO 1601 I = 1,N DO 1601 IC = 1,K PPOLD(IC,I) = 1.0/K PPOLD(IC,I) = PP(IC,I) 1601 CONTINUE C IF (N .GE. 31) GO TO 250 WRITE (6,14000) C WRITE (6,16000) WRITE (6,18000) (ID(I),ICLUS(I), I=1,N) 250 CONTINUE C C C 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 DO 2900 I = 1,N IF (ICLUS(I) .EQ. ICLSOL(I)) GO TO 2900 GO TO 3000 2900 CONTINUE GO TO 3300 3000 CONTINUE C C COUNT NUMBERS IN CLUSTERS: DO 3311 IC = 1,K NG(IC) = 0 3311 CONTINUE DO 3321 I = 1,N IC = ICLUS(I) NG(IC) = NG(IC) + 1 3321 CONTINUE WRITE (6,36000) (NG(IC),IC=1,K) C ITER = ITER + 1 C MAXITER = 30 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: GO TO 700 3100 WRITE (6,88000) C C 3300 CONTINUE C C WRITE (6,40000) (PR(IC),IC=1,K) DO 2800 IC = 1,K WRITE (6,30000) IC, (XBAR(IC,JV),JV=1,IP) 2800 CONTINUE C WRITE (6,42000) DO 2500 JV=1,IP WRITE (6,48000) (VARHAT(JV,JW),JW=1,IP) 2500 CONTINUE C IF (N .GE. 31) GO TO 3510 C LIST CLUSTER MEMBERSHIP OF EACH CASE: WRITE (6,52000) DO 3500 I = 1,N WRITE (6,50000) ID(I), ICLUS(I), ( X(I,JV),JV=1,IP ) 3500 CONTINUE 3510 CONTINUE C C C WRITE (6,58000) ITER C C C C COMPUTE VALUES OF MODEL-SELECTION CRITERIA: C C AIC = XMN2LL + 2.0*NOPARM SCH = XMN2LL + ALOG(XN)*NOPARM C IF ( N .GE. 301 ) GO TO 852 DO 851 IC=1,K WRITE (6,44444) IC DO 851 I = 1,N IF ( ICLUS(I) .EQ. IC ) WRITE (6,44445) ID(I) 851 CONTINUE 852 CONTINUE C C C GO ON TO NEXT NUMBER OF CLUSTERS (NEXT VALUE OF K), C IF NOT FINISHED. C WRITE (6,10050) C 850 CONTINUE C 4004 CONTINUE 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)). SCHPPH(1) = (IPARM(1)/2.0)*DLOG(2.0*PI) - SCHVEC(1)/2.0 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: AICPPH(1) = AICPPH(1) - IPARM(1) C Add in log of prior Pr[H(k)] for Schwarz criterion: 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. PRIRLN(1) = 0.0 C Add in log of approximate prior pdf: SCHPPH(1) = SCHPPH(1) + PRIRLN(1) 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 DO 6110 K = 2,KREACH SCHPPH(K) = (IPARM(K)/2.0)*DLOG(2*PI) - SCHVEC(K)/2.0 AICPPH(K) = (IPARM(K)/2.0)*DLOG(2*PI) - SCHVEC(K)/2.0 C Add in log Pr[H(k)]: AICPPH(K) = AICPPH(K) - IPARM(K) 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. PRIRLN(K) = 0.0 AICPPH(K) = AICPPH(K) + PRIRLN(K) 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 C AICPPH(K) = AICPPH(K) + AICVEC(1)/2 C 6110 CONTINUE C AT THIS POINT, WE HAVE THE LOGS OF THE PPHS; WE NOW CONVERT C THEM TO THE PPHS: SCHPPH(1) = DEXP(SCHPPH(1)) AICPPH(1) = DEXP(AICPPH(1)) DO 6119 K = 2,KREACH C If an argument is less or = -174.673, there will be underflow: IF ( SCHPPH(K) .LE. -174. ) GO TO 6115 SCHPPH(K) = DEXP(SCHPPH(K)) GO TO 6119 6115 SCHPPH(K) = 0.0 6119 CONTINUE DO 6121 K = 2,KREACH IF ( AICPPH(K) .LE. -174.673 ) GO TO 6126 AICPPH(K) = DEXP(AICPPH(K)) GO TO 6121 6126 AICPPH(K) = 0.0 6121 CONTINUE SPPHS = 0.0 SPPHA = 0.0 DO 6120 K = 1,KREACH IF ( SCHPPH(K) .GT. 0.0D+00 ) SPPHS = SPPHS + SCHPPH(K) IF ( AICPPH(K) .GT. 0.0D+00 ) SPPHA = SPPHA + AICPPH(K) 6120 CONTINUE DO 6130 K = 1,KREACH SCHPPH(K) = SCHPPH(K)/SPPHS AICPPH(K) = AICPPH(K)/SPPHA 6130 CONTINUE C C C WRITE VALUES OF MSC'S FOR VARIOUS K: WRITE (6,11000) TITLE WRITE(6,34000) DO 605 K = 1,KREACH WRITE(6,46000) K,XM2LL(K),IPARM(K),AICVEC(K),AICPPH(K), 1SCHVEC(K),SCHPPH(K) 605 CONTINUE C C WRITE (6,10050) WRITE (6,11002) TITLE WRITE (6,96666) 96666 FORMAT(' CONDITION CODE FOR THIS RUN: ') IF ( NGFLAG .EQ. 0 ) GO TO 4002 WRITE (6,96000) GO TO 4001 4002 CONTINUE WRITE (6,98000) 4001 CONTINUE WRITE (6,10050) WRITE (6,11097) WRITE (6,86000) WRITE (6,11100) MONTH,IDAY,IYEAR,IHOUR,MINUTE,ISEC CALL DTIME(IHOUR,MINUTE,ISEC) CALL TDATE(IDAY,MONTH,IYEAR) WRITE (6,89000) WRITE (6,11100) MONTH,IDAY,IYEAR,IHOUR,MINUTE,ISEC STOP C C 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 10050 FORMAT(/' ................................................', X'................................'//) 11112 FORMAT(' ----- ---- ---- ------'/) 19000 FORMAT( 1X,'MINIMUM OF EACH VARIABLE: ') 10100 FORMAT(' MEAN VECTOR: ',(9F8.2/)/) 10999 FORMAT( ' PROBLEM TITLE IS' ) 10200 FORMAT(/' COVARIANCE MATRIX: '/) 10500 FORMAT(/,' STATISTICS FOR WHOLE (UNCLUSTERED) SAMPLE: '/) 11001 FORMAT(1X,18A4/) 11002 FORMAT(/1X,18A4/) 11097 FORMAT(' MIXPDTA CLUSPAC - MIXTURE-MODEL CLUSTERING OF CASES'/) 11100 FORMAT(1X,I2,'/',I2,'/',I4,3X,'AT ',I2,':',I2,':',I2/) 25000 FORMAT(1X, A6, 9F8.3) 55555 FORMAT(/' CASE POSTERIOR PROBS OF CLUSTER MEMBERSHIP: '/) 20000 FORMAT(/ 1X,'INITIAL MEANS') 44444 FORMAT(' CLUSTER ',I2) 44445 FORMAT(1X,A6) 10900 FORMAT(/' MINUS 2 LOG LIKELIHOOD = ',F13.1) 11000 FORMAT(/,1X,18A4 ) 11500 FORMAT(' DET = ',E15.7, ' IDET = ',I2/) 12000 FORMAT(2X,I4) 13000 FORMAT( ' NUMBER OF OBSERVATIONS (SAMPLE SIZE), N . . ', I4/) 14000 FORMAT(1X, 'CLUSTERING:') 15000 FORMAT(18A4) 87000 FORMAT(1X,2E15.7) 86000 FORMAT(/' TIME BEGAN: ') 89000 FORMAT(/' TIME FINISHED: ') 16000 FORMAT(/,1X,'CASES AND LABELS:--'/) 17000 FORMAT(' PREC.MX. OF CLUSTER ',I2,' HAS BEEN SET TO POOLED ', X'PREC.MX.') 18000 FORMAT(1X, (7(1X,A6,': ',I2)/) ) 21999 FORMAT(' '//) 22000 FORMAT(1X,'MEAN VECTOR FOR CLUSTER ',I2,': ') 26000 FORMAT(/, 1X,'K = ',I2,' CLUSTERS'/) 30000 FORMAT(1X,'MEAN VECTOR FOR CLUSTER ',I2,': '/,(9F13.3/)) 32000 FORMAT(/,1X,'ITERATION ', I2) 11096 FORMAT(/' NUMBER OF CLUSTERS TO REPORT. . . . . . . . 1 TO', XI2/) 11111 FORMAT(/' FIRST FOUR DATA POINTS') 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.') 21000 FORMAT(1X, 'MAXIMUM OF EACH VARIABLE: ') 36000 FORMAT(/,1X,'NUMBERS IN CLUSTERS:',3X,9(I5,3X)/) 40000 FORMAT(/,1X,'MIXING PROBABILITIES:',3X,9(F4.2,3X)/) 42000 FORMAT(/, 1X,'COMMON COVARIANCE MATRIX (MLE):',/ ) 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) 46000 FORMAT(5X,I1,6X,F8.1, 4X,I4,1X,F9.1,1X,F4.2,1X,F9.1,1X,F4.2) 48000 FORMAT(1X,5F13.3) 50000 FORMAT(1X,A6,': ',I2,(8F13.3/)) 52000 FORMAT(//,1X,'CLUSTER MEMBERSHIP OF EACH CASE:'/ X1X,'CASE LABEL DATA'/ ) 54000 FORMAT(3X,I2) 56000 FORMAT(/' NUMBER OF VARIABLES . . . . . . . . . . . . ',I3 ) 58000 FORMAT(1X, 'CONVERGENCE: NO CASE CHANGED CLUSTERS AFTER ', X'ITERATION ',I2,'.'/ ) 60000 FORMAT(/,1X,'JFLG = ',I2,'. IF JFLG=0, COMPUTATION OF DET', X' WENT WELL; OTHERWISE, THERE WAS TROUBLE OR MATRIX WAS ', X'ILL-CONDITIONED.'//) 62000 FORMAT(/1X,'DET = ',E13.5,' IDET = ',I3,5X, X'ACTUAL DET. = DET*10**IDET',//) 70000 FORMAT(1X,' AIC = ', F13.1) 78000 FORMAT(1X,' SCHWARZ CRITERION = ', F13.1) 72000 FORMAT( ' NUMBER OF PARAMETERS = ',I4) 74000 FORMAT(/1X,'EXPECTED NO. OF OBSERVATIONS IN GROUP ',I3,'IS ', X' LESS THAN 0.5. ', X' END COMPUTATION FOR THIS NUMBER OF CLUSTERS. '/) 76000 FORMAT(/' POSTERIOR PROB OF GROUP MEMBERSHIP FOR 1ST 4 CASES:') 82000 FORMAT(1X,29F5.2) 88000 FORMAT(' ALGORITHM HAS NOT CONVERGED IN 20 ITERATIONS. STOP') 96000 FORMAT(/' NGFLAG=1: FOR SOME K, AT LEAST ONE CLUSTER HAS '/ X' AN EXP.NO. OF OBSERVATIONS LESS TH AN 0.5 .') 98000 FORMAT(//1X,'PROGRAM ENDED NORMALLY.') END C C C 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. 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. REAL*8 A(N,1),DET,P(LL,1) REAL*8 DNORM,DEN,DMULT,DSUM,DISIGN DIMENSION IV(1) NRS=NRS1 IF (NRS.EQ.0) IDET=1 DISIGN=1.0D+00 DET=0.0D+00 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 M1=M-1 DO 100 I=1,M 100 IV(I)=I IF (NRS) 500,200,500 200 DO 300 I=1,M DO 300 J=1,M 300 P(I,J)=0.0D+00 DO 400 I=1,M 400 P(I,I)=1.0D+00 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 500 DO 1200 K=1,M1 ICOL=K IPCOL=K C C SEARCHING FOR LARGEST ELEMENT IN ABSOLUTE VALUE IN COLUMN K. C DNORM=A(IV(K),K) IFLG=0 KK=K+1 DO 600 J=KK,M IF (DABS(A(IV(J),K)).LE.DABS(DNORM)) GO TO 600 IFLG=1 IPCOL=IV(J) DNORM=A(IPCOL,K) 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 IF (IFLG.EQ.0) GO TO 800 ISAVE=IV(ICOL) IV(ICOL)=IPCOL ICOL1=ICOL+1 DO 700 L=ICOL1,M IF (IV(L).EQ.IPCOL) IV(L)=ISAVE 700 CONTINUE DISIGN=-DISIGN 800 IF (DNORM.EQ.0.0D+00) GO TO 1900 C C BEGIN ELIMINATION OF ROW BELOW IV(K).DEN IS THE PIVOT ELEMENT. C K1=K+1 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 IF (A(IV(IM),K).EQ.0.0D+00) GO TO 1100 C C CALCULATE ELIMINATION FACTOR. C DMULT=-A(IV(IM),K) C C WE NOW CALCULATE VALUE OF OTHER ELEMENTS IN ROW IV(IM). C DO 900 NN=K1,M 900 A(IV(IM),NN)=(DMULT*A(IV(K),NN))/DNORM+A(IV(IM),NN) IF (NRS.LE.0) GO TO 1100 DO 1000 IN=1,NRS 1000 P(IV(IM),IN)=(DMULT*P(IV(K),IN))/DNORM+P(IV(IM),IN) 1100 CONTINUE 1200 CONTINUE C C CALCULATE VALUE OF DETERMINANT. C IF (A(IV(M),M).EQ.0.0D0) GO TO 1900 DET=DISIGN IF (IDET.NE.0) CALL DMATDT(A,N,M,DET,IV,IDET) IF (DET.EQ.0.0D+00) GO TO 1900 IF (NRS.LE.0) GO TO 2000 C C WE START SOLVING RIGHT HAND SIDES.THE SOLUTION REPLACES THE RIGHT HAND C VECTOR. C 1300 N1=M-1 DO 1600 JJ=1,NRS C C BEGIN BACK SUBSTITUTION. C P(IV(M),JJ)=P(IV(M),JJ)/A(IV(M),M) DO 1500 I=1,N1 DSUM=0.0D+00 DO 1400 J=1,I 1400 DSUM=DSUM-A(IV(M-I),M-J+1)*P(IV(M-J+1),JJ) 1500 P(IV(M-I),JJ)=(P(IV(M-I),JJ)+DSUM)/A(IV(M-I),M-I) 1600 CONTINUE DO 1800 JJ=1,NRS DO 1700 IND=1,M 1700 A(IND,1)=P(IV(IND),JJ) DO 1800 IND=1,M 1800 P(IND,JJ)=A(IND,1) RETURN 1900 JFLG=1 IDET=0 2000 RETURN END SUBROUTINE MATDT(A,IA,N,DET,IV,IDET) C SUBROUTINE MATDT IS DMATDT FROM THE UICC SUBROUTINE LIBRARY. REAL*8 A(IA,1),DET,B,LOG16 INTEGER*4 IV(1),K EQUIVALENCE (B,K) NUM=16777216 LOG16=.120411998265592457D+01 IF (A(IV(N),N).EQ.0.0D+00) GO TO 300 L=0 DO 100 I=1,N B=DABS(A(IV(I),I)) K=K/NUM-64 L=L+K 100 DET=DET*(A(IV(I),I)/16.0D+00**K) B=DABS(DET) K=K/NUM-64 IW=L+K IF ((IW.LT.-64).OR.(IW.GT.63)) GO TO 200 DET=DET*16.0D+00**L IDET=0 GO TO 400 200 DET=DET*16.0D+00**(-K) IDET=L+K B=IDET*LOG16 IDET=B B=B-DFLOAT(IDET) DET=DET*1.0D+01**B GO TO 400 300 DET=0.0D+00 IDET=0 400 RETURN END