C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C C C CLUSPAC: Computer Programs for Mixture-Model Clustering C C C C COPYRIGHT (C) 1991, 1992 STANLEY LOUIS SCLOVE C C C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C C C CMS NAME OF PROGRAM: MIXPCM CLUSPAC C C VERSION 7.7 7-APR-92 C C C C 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 601 S. MORGAN ST. C CHICAGO, IL 60607-7124 C C C C C C COPYRIGHT (C) 1991 STANLEY LOUIS SCLOVE C C C C C C ISOPAC 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 ISOPAC 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 MIXPCM (MIXture model, P-variate data, CoMmon C C *** * * * C C covariance matrix) in the ISOPAC 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 assumes a common covariance matrix across C C distributions. MIXPDT allows different covariance C C matrices. C C C C Input: C C ----- C C Number of clusters (K), initial values of means, prior C C probabilities, and common covariance matrix. If desired, C C program ISDTPCM.ISOPAC can be used to obtain these initial C C values. Use program MIXPCMA.ISOPAC to try a range of C C numbers of clusters (values of K), with automatic setting C C of initial values. C C C C Program restrictions (can be modified): C C -------------------------------------- C C N, sample size, at most 1000; C C IP, number of variables, at most 20; C C K, number of clusters, at most 29; C C ITER, maximum number of iterations, 20. C C C C C C Subroutines called: C C MATEQ, which calls MATDT C C C C IV is a work array for subroutine MATEQ. C C C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C C C C C CONTROL CARDS: C C C C DATASET TITLE C C N, IN FORMAT (2X,I4) C C IP, IN FORMAT (3X,I2) C C DATFMT, IN FORMAT (18A4), E.G., (4F4.1) C C "DATFMT" WILL ALSO BE USED FOR OUTPUT, SO ALLOW AT LEAST C C ONE BLANK AT THE BEGINNING FOR CARRIAGE CONTROL. C C DATA, ONE CASE AT A TIME, IN FORMAT SPECIFIED BY DATFMT C C K, NUMBER OF CLUSTERS, IN FORMAT (2X,I2) C C MEANFT, in format (18A4) C C K INITIAL MEANS, IN FORMAT SPECIFIED BY MEANFT C C C C K INITIAL VALUES OF PRIOR PROBABILITIES, C C IN FORMAT (5X,F3.2). C C C C COVFMT, in format (18A4). C C Initial value of VARHAT, the common covariance matrix. C C C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C DIMENSION X(1000,20),SUM(29,20) DIMENSION F(1000,29),PP(29,1000) DIMENSION ICLUS(1000),ICLSOL(1000) DIMENSION DENOM(1000),XMXPR(1000) DIMENSION DSQ(29) DIMENSION TITLE(18) DIMENSION NG(29),XBAR(29,20) DIMENSION DATFMT(18) DIMENSION MEANFT(18) DIMENSION COVFMT(18) DIMENSION SSD(29,20,20) DIMENSION WGSS(20,20) DIMENSION VARHAT(20,20) DIMENSION XMIN(20),XMAX(20) DIMENSION IV(20,20) DIMENSION P(20,20) DIMENSION PR(29) C DOUBLE PRECISION SUM DOUBLE PRECISION WGSS,SSD DOUBLE PRECISION VARHAT DOUBLE PRECISION P DOUBLE PRECISION DET,TRUDET DOUBLE PRECISION DSQ DOUBLE PRECISION XBAR DOUBLE PRECISION DEVV,DEVW DOUBLE PRECISION F C C C C FLOW OF PROGRAM: C C C READ DATA AND INITIAL PARAMETER ESTIMATES. C INVERT COVARIANCE MATRIX. C PRINT INITIAL PARAMETER ESTIMATES. C C ITERATION: C COMPUTE ESTIMATED VALUES OF PROBABILITY DENSITY FUNCTIONS. C COMPUTE LIKELIHOOD AND VALUES OF MODEL-SELECTION CRITERIA. C COMPUTE CLUSTER-MEMBERSHIP PROBABILITIES. C UPDATE PARAMETER ESTIMATES (INCLUDING C INVERSE COVARIANCE MATRIX). C CLUSTER BY MAXIMUM PROBABILITY OF CLUSTER MEMBERSHIP. C IF CLUSTERING HASN'T CHANGED, STOP AND PRINT RESULTS. C OTHERWISE DO ANOTHER ITERATION (UNLESS 20 HAVE C ALREADY BEEN DONE). C C C DATA PI/3.141593/ C READ (5,36000) TITLE C C WRITE PROGRAM INFORMATION. WRITE (6,10000) WRITE (6,10050) WRITE (6,40000) TITLE C C READ SAMPLE SIZE, N. READ (5,12000) N XN = N WRITE (6,38000) N C READ NUMBER OF VARIABLES, IP. READ (5,54000) IP WRITE (6,56000) IP C C READ DATA FORMAT. READ (5,36000) DATFMT C C READ DATA. DO 500 I = 1,N READ (5,DATFMT) (X(I,JV), JV = 1,IP) IF (I .EQ. 1) GO TO 100 GO TO 300 100 CONTINUE C COMPUTE MINIMA AND MAXIMA: 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 C WRITE MINIMA AND MAXIMA: WRITE (6,64000) WRITE (6,DATFMT) (XMIN(JV),JV=1,IP) WRITE (6,66000) WRITE (6,DATFMT) (XMAX(JV),JV=1,IP) C C READ K, NUMBER OF CLUSTERS. READ (5,11000) K WRITE (6,26000) K C C READ INITIAL MEANS C READ INPUT FORMAT FOR MEANS: READ (5,36000) MEANFT DO 600 IC=1,K READ (5,MEANFT) (XBAR(IC,JV), JV=1,IP) C 600 CONTINUE C READ INITIAL VALUES OF PRIOR PROBABILITIES: DO 650 IC=1,K READ (5,25000) PR(IC) C 650 CONTINUE C C C READ INITIAL VALUE OF COVARIANCE MATRIX, VARHAT: C READ(5,36000) COVFMT DO 660 JV = 1,IP READ(5,COVFMT) (VARHAT(JV,JW),JW=1,IP) 660 CONTINUE C C WRITE(6,22000) ITER = 1 700 CONTINUE WRITE (6,32000) ITER C WRITE CURRENT ESTIMATES OF PARAMETERS: C WRITE CURRENT ESTIMATES OF MEAN VECTORS:-- WRITE (6,20000) DO 710 IC = 1,K WRITE (6,MEANFT) ( XBAR(IC,JV), JV=1,IP ) 710 CONTINUE WRITE (6,23000) WRITE (6,35000) (PR(IC), IC = 1,K) C WRITE CURRENT ESTIMATE OF COVARIANCE MATRIX, VARHAT:-- WRITE(6,42000) DO 670 JV = 1,IP WRITE(6,COVFMT) (VARHAT(JV,JW),JW=1,IP) 670 CONTINUE C C CALL SUBROUTINE TO COMPUTE INVERSE COVARIANCE MATRIX: C SET PARAMETERS OF SUBROUTINE CALL: IDET = 1 NRS1 = 0 C C CALL MATEQ(VARHAT,IP,20,JFLG,DET,IDET,IV,NRS1,P,20) C ON RETURN, P CONTAINS THE INVERSE MATRIX. C C C GENERAL FORM OF CALL IS: C CALL MATEQ(A,M,N,JFLG,DET,IDET,IV,NRS1,P,LL) C SEE SUBROUTINE LISTING FOR FULLER EXPLANATION. C DET(VARHAT) = DET*10.0**IDET C C ON RETURN, P CONTAINS THE INVERSE MATRIX. C C C C WRITE VALUE OF INVERSE COVARIANCE MATRIX: C WRITE(6,44000) C DO 680 JV = 1,IP C WRITE(6,COVFMT) (P(JV,JW),JW=1,IP) C 680 CONTINUE C C If determinant is not positive, replace it by an approximation C based on diagonal version of the covariance matrix: IF ( DET .LE. 0.0 ) GO TO 601 GO TO 602 601 CONTINUE DET = 1.0 DO 603 IVAR=1,IP DET = DET*VARHAT(IVAR,IVAR) 603 CONTINUE TRUDET = DET GO TO 604 602 CONTINUE C XIDET = IDET XLGDET = DLOG(DET) + XIDET*ALOG(10.0) TRUDET = DET*(10.0**IDET) 604 CONTINUE C C C COMMENCE COMPUTATION OF F(I,IC), I=1,...,N, IC=1,...,K: C C COMPUTE MAHALANOBIS D-SQUARE BETWEEN THE I-TH OBSERVATION AND C THE L-TH MEAN, L=1,2,...,K, I=1,2,...,N:-- C DO 1200 I = 1,N C DO 1000 L=1,K DSQ(L) = 0.0 DO 1310 JV=1,IP DEVV = X(I,JV) - XBAR(L,JV) DO 1310 JW=1,IP DEVW = X(I,JW) - XBAR(L,JW) DSQ(L) = DSQ(L) + DEVV*P(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): C IF ( ZSQ/2.0 .LE. 174.673 ) GO TO 1090 IF ( ZSQ/2.0 .LE. 100.0 ) GO TO 1090 F(I, L) = 0.0 GO TO 1100 C 1090 CONTINUE F(I, L) = EXP(-ZSQ/2.0)/((2*PI)**(IP/2)) F(I, L) = F(I,L)/DSQRT(TRUDET) 1100 CONTINUE 1000 CONTINUE 1200 CONTINUE C C COMPUTE LOG LIKELIHOOD AND VALUES OF MODEL-SELECTION CRITERIA: C 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 IF ( SUMPXF .EQ. 0.0 ) GO TO 2200 SUMLNF = SUMLNF + ALOG(SUMPXF) 2200 CONTINUE C XMN2LL = -2.0*SUMLNF C WRITE (6,30000) 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. NOPARM = K*IP + IP*(IP+1)/2 + (K-1) WRITE (6,72000) NOPARM AIC = XMN2LL + 2.0*NOPARM SCH = XMN2LL + ALOG(XN)*NOPARM WRITE (6,70000) AIC WRITE (6,71000) SCH C C C COMPUTE POSTERIOR PROBABILITIES OF GROUP MEMBERSHIP, C PP(IC,I), THE CONDITIONAL PROBABILITY OF POP'N IC, GIVEN X(I), C AS PR(IC)*F(I,IC)/DENOM(I): 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 DO 1400 I = 1,N DO 1400 IC=1,K IF ( DENOM(I) .GT. 0.0 ) GO TO 1410 PP(IC,I) = 0.0 GO TO 1400 1410 CONTINUE PP(IC,I)= PR(IC)*F(I,IC)/DENOM(I) 1400 CONTINUE C WRITE(6,76000) DO 1420 I = 1,4 WRITE(6,82000) ( PP(IC,I), IC = 1,K ) 1420 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 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 IF (N*PR(IG) .LT. .05) GO TO 2050 GO TO 2150 2050 WRITE (6,74000) IG GO TO 4000 2150 CONTINUE DO 1900 JV = 1,IP XBAR(IG,JV) = SUM(IG,JV)/(N*PR(IG)) 1900 CONTINUE C 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) TERM=PP(IC,I)*DEVV*DEVW SSD(IC,JV,JW)=SSD(IC,JV,JW)+TERM 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 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 (END PARAMETER-ESTIMATE UPDATE SEQUENCE) C IF (ITER .EQ. 1) GO TO 900 C STORE OLD LABELS: DO 800 I = 1,N ICLSOL(I) = ICLUS(I) 800 CONTINUE C 900 CONTINUE C C COMPUTE NEW LABELS BY MAX POSTERIOR PROBABILITY: 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 C IF (N .GE. 31) GO TO 250 WRITE (6,14000) C WRITE (6,16000) WRITE (6,18000) (I, ICLUS(I), I=1,N) 250 CONTINUE C C C C IF (ITER .EQ. 1) GO TO 3000 DO 2900 I = 1,N IF (ICLUS(I) .EQ. ICLSOL(I)) GO TO 2900 GO TO 3000 2900 CONTINUE GO TO 3300 3000 CONTINUE ITER = ITER + 1 IF (ITER.GE.51) GO TO 3100 GO TO 700 3100 WRITE (6,68000) 4000 STOP C C 3300 CONTINUE C C C COUNT NUMBERS IN CLUSTERS: DO 3310 IC = 1,K NG(IC) = 0 3310 CONTINUE DO 3320 I = 1,N IC = ICLUS(I) NG(IC) = NG(IC) + 1 3320 CONTINUE C C C PRINT FINAL RESULTS C C WRITE (6,34000) (NG(IC),IC=1,K) WRITE (6,35000) (PR(IC),IC=1,K) DO 2800 IC = 1,K WRITE (6,28000) IC WRITE (6,MEANFT) (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 .GT. 60 ) GO TO 3510 WRITE (6,52000) DO 3500 I = 1,N WRITE (6,50000) I, ICLUS(I) WRITE (6,DATFMT) (X(I,JV),JV=1,IP) 3500 CONTINUE 3510 CONTINUE C WRITE (6,84000) C C STOP C C C 10000 FORMAT('1','****************************************'//// X' CLUSTER ANALYSIS OF CASES '/ X' USING COMMON COVARIANCE MATRIX'// X' CMS FILE: MIXPCM CLUSPAC '// X' DEVELOPED AND PROGRAMMED BY DR. STANLEY L. SCLOVE'// X' MIXPCM CLUSPAC VERSION 7.7 7-APR-92 '//' COPYRIGHT', Z' 1991, 1992 STANLEY L. SCLOVE ') 10050 FORMAT(' ************************************************', X'********************************'//) 11000 FORMAT(2X,I2) 12000 FORMAT(2X,I4) 14000 FORMAT(//1X,'CLUSTERING:'/) 16000 FORMAT(/,1X,'CASES AND LABELS:--'/) 18000 FORMAT(15(I5,I3)) 20000 FORMAT(//1X,'MEANS'/) 22000 FORMAT(' FIRST ITERATION USES THE INITIAL PARAMETER ESTIMATES'/ X' PROVIDED BY THE USER.'/) 25000 FORMAT(5X, F3.2) 23000 FORMAT(//1X, 'CURRENT ESTIMATES OF PRIOR PROBABILITIES:') 26000 FORMAT('1',//,1X,'K = ',I2,' CLUSTERS') 28000 FORMAT(1X,'MEAN VECTOR FOR CLUSTER ',I2,': ') 30000 FORMAT(/1X,' MINUS 2 LOG LIKELIHOOD = ', F13.5//) 32000 FORMAT(1X,'ITERATION ', I2) 34000 FORMAT(/,1X,'NUMBERS IN CLUSTERS:',3X,29(I3,3X)/) 35000 FORMAT(/,1X,'MIXING PROBABILITIES:',3X,10(F5.3,3X)/) 36000 FORMAT(18A4) 38000 FORMAT(1X,'NUMBER OF OBSERVATIONS (SAMPLE SIZE), N = ',I4/) 40000 FORMAT(1X,18A4//) 42000 FORMAT(//,1X,'COMMON COVARIANCE MATRIX (MLE):',/) C 44000 FORMAT(//,1X,'INVERSE COVARIANCE MATRIX:',/) C 48000 FORMAT(1X,8F13.5/) 50000 FORMAT(1X,I4,1X,I2) 52000 FORMAT(//,1X,'CASE, LABEL / DATA'//) 54000 FORMAT(3X,I2) 56000 FORMAT(1X,'NUMBER OF VARIABLES = ',I2/) 58000 FORMAT(/1X,'CONVERGENCE: NO CASE CHANGED CLUSTERS AFTER ', X'ITERATION ',I2,'. RESULTS ARE PRINTED BELOW.'//) 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 = ',F13.5,' IDET = ',I3,5X, X'ACTUAL DET. = DET*10**IDET',//) 64000 FORMAT(/1X,'MINIMUM FOR EACH VARIABLE: ',/) 66000 FORMAT(/1X,'MAXIMUM FOR EACH VARIABLE: ',/) 68000 FORMAT(1X,'PROGRAM HAS NOT CONVERGED IN 50 ITERATIONS: STOP'//) 70000 FORMAT(1X,'AIC = ', F15.5,/) 71000 FORMAT(1X,'SCHWARZ CRITERION = ', F15.5,/) 72000 FORMAT(//1X,'NUMBER OF PARAMETERS = ',I4//) 74000 FORMAT(1X,'EXPECTED NUMBER OF OBSERVATIONS IN CLUSTER ',I3, X' IS LESS THAN .05. STOP.'/) 76000 FORMAT(' POSTERIOR PROBS OF GROUP MEMBERSHIP FOR 1ST 4 CASES:') 82000 FORMAT(/,1X,29F5.2/) 84000 FORMAT(/,1X,'PROGRAM ENDED NORMALLY.') END C C SUBROUTINE MATEQ(A,M,N,JFLG,DET,IDET,IV,NRS1,P,LL) C SUBROUTINE MATEQ IS DMATEQ FROM THE UICC SUBROUTINE LIBRARY. C PAGE 1 C C SUBROUTINE DMATEQ C ***************** C THIS ROUTINE WILL SOLVE A REAL*8 SYSTEM OF LINEAR EQUATIONS,COMPUTE C THE DETERMINANT, WITHOUT UNDERFLOW OR OVERFLOW, OF A REAL*8 MATRIX, C AND/OR INVERT A REAL*8 MATRIX. C CALLING SEQUENCE: C CALL DMATEQ(A,N,IA,JFLG,DET,IDET,IV,NRS,P,IP) WHERE; C A (INPUT) - IS THE REAL*8 MATRIX ON WHICH THE ROUTINE IS C TO WORK. IN THE PROCESS OF COMPUTATION THE C CONTENTS OF THIS MATRIX ARE DESTROYED. C N (INPUT) - IS AN INTEGER*4 VARIABLE WHICH SPECIFIES THE C ORDER OF THE A MATRIX. C IA (INPUT) - IS AN INTEGER*4 VARIABLE WHICH SPECIFIES THE C ACTUAL ROW DIMENSION OF A AS DIMENSIONED IN C THE CALLING PROGRAM. IA MUST BE GREATER THAN C OR EQUAL TO N. C JFLG (OUTPUT) - IS AN INTEGER*4 RETURN CODE VARIABLE. UPON C RETURN FROM DMATEQ IF; C JFLG=0, ALL WENT WELL. C JFLG=1, THE A MATRIX WAS SINGULAR OR NEAR C SINGULAR AND THE COMPUTATIONS COULD NOT BE C COMPLETED. THE CONTENTS OF THE VARIABLES C A, DET, IDET AND P ARE MEANINGLESS. C DET (OUTPUT) - IS A REAL*8 VARIABLE WHICH CONTAINS THE C DETERMINANT OF A. (SEE IDET) C IDET (INPUT) - IS AN INTEGER*4 VARIABLE. ON INPUT IF; C IDET=0, NO DETERMINANT IS CALCULATED. C IDET NOT 0, THE DETERMINANT OF A IS COMPUTED. C ON OUTPUT IDET CONTAINS THE POWER OF 10 C THAT DET SHOULD BE MULTIPLIED BY TO GIVE THE C CORRECT VALUE OF THE DETERMINANT. I.E. C DET(A)=DET*10.0D0**IDET. C IF DET(A) CAN BE COMPUTED WITHOUT UNDER OR C OVERFLOW, THEN IDET=0 OTHERWISE IDET IS SET C TO THE PROPER VALUE SO THAT NO UNDER OR OVER- C FLOW WILL OCCUR IN COMPUTING DET. C IV (INPUT) - IS AN INTEGER*4 WORK ARRAY WHICH SHOULD BE C DIMENSIONED AT LEAST IV(N). C PAGE 2 C C NRS (INPUT) - IS AN INTEGER*4 VARIABLE WITH THE FOLLOWING C INTERPRETATION: C NRS>0, SOLVE A SYSTEM OF LINEAR EQUATIONS C WITH NRS RIGHT HAND SIDES. C NRS=0, INVERT THE A MATRIX. C NRS<0, ONLY COMPUTE THE DETERMINANT OF A. 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 C ALREADY BEEN ANNIHILATED. 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