C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C C C CLUSPAC: Computer Programs for Mixture-Model Clustering C C C C COPYRIGHT (C) 1991, 1992, 1993 STANLEY L. SCLOVE C C C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C C C CMS DSN = EVALUATE CLUSPAC C C VERSION 1.4 30-JAN-90 C C C C PROGRAMMED BY: C C C C DR. STANLEY L. SCLOVE 312/996-2681 C C INFORMATION AND DECISION SCIENCE DEPT. M/C 294 C C COLLEGE OF BUSINESS ADMINISTRATION C C UNIVERSITY OF ILLINOIS AT CHICAGO C C 601 S. MORGAN ST. C C CHICAGO, IL 60607-7124 C C C C C C COPYRIGHT 1991 STANLEY LOUIS SCLOVE. C C C C ..............................................................C C C C C C Program EVALUATE CLUSPAC evaluates the likelihood C C of classified (labeled) data under the assumption of C C Gaussian class-conditional distributions. The evaluation C C is made both for a common covariance matrix and varying C C covariance matrices. Values of model-selection criteria C C are computed. These may be compared with those obtained C C by clustering. C C C C ..............................................................C C C C Input: C C ----- C C Number of groups (K), data and labels. C C C C Program restrictions (can be modified): C C -------------------------------------- C C N, sample sizes, at most 1000; C C IP, number of variables, at most 20; C C K, number of groups, at most 29; 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 C ..............................................................C C C C C C C C CONTROL CARDS: C C C C FISHER IRIS DATA C C K=03 C C IP=04 C C (4(1X,F3.1),1X,I1) C C N=0050 C C C C THESE ARE: C C C C DATASET TITLE C C K, NUMBER OF GROUPS, IN FORMAT (2X,I2) C C NO. OF VARIABLES, 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 N FOR FIRST SAMPLE, IN FORMAT (2X,I4) C DATA FOR FIRST SAMPLE C C (DATA, ONE CASE AT A TIME, IN FORMAT SPECIFIED BY DATFMT) C C N FOR SECOND SAMPLE, IN FORMAT (2X,I4) C DATA FOR SECOND SAMPLE C . C . C . C C N FOR K-TH SAMPLE, IN FORMAT (2X,I4) C DATA FOR K-TH SAMPLE C C C C C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C DIMENSION X(1000,20),SUM(29,20),XMEAN(20) C ELEMENTS OF X ARE X(I,JV), I=1,2,...,N, JV=1,2,...,IP. DIMENSION TITLE(18) DIMENSION DATFMT(18) DIMENSION NG(29),TOTAL(29),SUMSQS(20,20) DIMENSION SSD(29,20,20),SIGMA(29,20,20),SSDEVS(20,20) DIMENSION WGSS(20,20) DIMENSION VARHAT(20,20),COMSIG(20,20),COVMX(20,20) DIMENSION XMIN(20),XMAX(20) DIMENSION IV(20,20) DIMENSION XLLIKE(29) C DOUBLE PRECISION SUM DOUBLE PRECISION WGSS,SSD DOUBLE PRECISION VARHAT,COVMX,TOTAL,SUMSQS,SSDEVS DOUBLE PRECISION P DOUBLE PRECISION DET,TRUDET DOUBLE PRECISION DSQ DOUBLE PRECISION DEVV,DEVW DOUBLE PRECISION F DOUBLE PRECISION XLLIKE C C C C FLOW OF PROGRAM: C C C FOR EACH SAMPLE:-- C READ DATA. C COMPUTE MEAN VECTOR AND COVARIANCE MATRIX. C EVALUATE LIKELIHOOD. C L = EXP(-NP/2)/( (2*PI)**NP/2*DET(COV MX)**N/2 ) C EVALUATE LIKELIHOOD WITH VARYING COVARIANCE MATRICES BY C MULTIPLYING LIKELIHOODS FOR THE SEPARATE SAMPLES. C COMPUTE MODEL-SELECTION CRITERIA. C C POOL TO OBTAIN COMMON COVARIANCE MATRIX. C EVALUATE LIKELIHOOD WITH COMMON COVARIANCE MATRIX. C COMPUTE MODEL-SELECTION CRITERIA. C PRINT RESULTS. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C PI = .314159265358979324D+01 C READ (5,36000) TITLE C C WRITE PROGRAM INFORMATION. WRITE (6,10000) WRITE (6,10050) WRITE (6,40000) TITLE C C READ K, NUMBER OF GROUPS. READ (5,11000) K WRITE (6,26000) K C READ NUMBER OF VARIABLES, IP. READ (5,54000) IP WRITE (6,13000) IP WRITE (6,10050) C C READ DATA FORMAT. READ (5,36000) DATFMT C C CALL XUFLOW(0) C C DO FOR EACH SAMPLE: C C DO 4100 ISAMPL = 1,K WRITE (6,70000) ISAMPL C READ SAMPLE SIZE, N. READ (5,12000) N XN = N NG(ISAMPL) = N WRITE (6,15000) N 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 C C COMPUTE SUM VECTOR AND SUM-OF-PRODUCTS MATRIX: DO 105 JV = 1,IP TOTAL(JV) = 0.0D+00 DO 105 JW = 1,IP SUMSQS(JV,JW) = 0.0D+00 SSDEVS(JV,JW) = 0.0D+00 105 CONTINUE DO 110 JV = 1,IP DO 110 I = 1,N TOTAL(JV) = TOTAL(JV) + X(I,JV) 110 CONTINUE C DO 130 JV = 1,IP XMEAN(JV) = TOTAL(JV)/XN 130 CONTINUE DO 120 JV = 1,IP DO 120 JW = 1,IP DO 120 I = 1,N DEVV=X(I,JV)-XMEAN(JV) DEVW=X(I,JW)-XMEAN(JW) SSDEVS(JV,JW) = SSDEVS(JV,JW) + DEVV*DEVW 120 CONTINUE DO 140 JV = 1,IP DO 140 JW = 1,IP VARHAT(JV,JW) = SSDEVS(JV,JW)/XN SIGMA(ISAMPL,JV,JW) = SSDEVS(JV,JW) 140 CONTINUE C WRITE SUMMARY STATISTICS FOR THIS SAMPLE: WRITE(6,21000) WRITE(6,22000) 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 DO 2404 JV = 1,IP DO 2404 JW = 1,IP COVMX(JV,JW) = VARHAT(JV,JW) 2404 CONTINUE WRITE (6,11500) WRITE (6,11700) COVMX(1,1) CALL MATEQ(COVMX,IP,20,JFLG,DET,IDET,IV,NRS1,P,20) WRITE (6,11600) 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(A) = DET*10.0**IDET C WRITE (6,60000) JFLG C C XIDET = IDET C XLGDET = DLOG(DET) + XIDET*ALOG(10.0) TRUDET = DET*(10.0**IDET) WRITE (6,62000) DET,IDET C C COMPUTE LIKELIHOOD ACCORDING TO 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: XLL=-(XN*IP/2.0)*DLOG(2.0D+00*PI)-(XN/2.0)*DLOG(TRUDET)-XN*IP/2.0 XLLIKE(ISAMPL) = XLL C C C C GO ON TO NEXT SAMPLE (I.E., NEXT VALUE OF ISAMPL): C 4100 CONTINUE C C EVALUATE LIKELIHOOD WITH VARYING COVARIANCE MATRICES BY C MULTIPLYING LIKELIHOODS FOR THE SEPARATE SAMPLES, C I.E., BY ADDING THEIR LOGARITHMS: XLOGL = 0.0 DO 4150 ISAMPL=1,K XLOGL = XLOGL + XLLIKE(ISAMPL) 4150 CONTINUE C COMPUTE MODEL-SELECTION CRITERIA: WRITE (6,69000) C COMPUTE MODEL-SELECTION CRITERIA FOR MODEL WITH C VARYING COVARIANCE MATRICES:-- WRITE (6,69100) C C PARAMETERS: C K MEAN VECTORS OF DIMENSION P AND K P-BY-P COVARIANCE MATRICES, C WHERE P IS THE NUMBER OF VARIABLES:-- C NOPARM = K*IP + K*IP*(IP+1)/2 C WRITE (6,72000) NOPARM XMN2LL = -2.0*XLOGL WRITE (6,76000) XMN2LL AIC = XMN2LL + 2.0*NOPARM WRITE (6,78000) AIC SCH = XMN2LL + ALOG(XN)*NOPARM WRITE (6,80000) SCH C COMPUTE MODEL-SELECTION CRITERIA FOR MODEL WITH C COMMON COVARIANCE MATRIX:-- C C BEGIN BY COMPUTING POOLED ESTIMATE OF COVARIANCE MATRIX: C C POOL: C DO 1950 JV = 1,IP DO 1950 JW = 1,IP WGSS(JV,JW) = 0.0D+00 1950 CONTINUE C DO 2300 JV = 1,IP DO 2300 JW = 1,IP DO 2300 IC = 1,K WGSS(JV,JW) = WGSS(JV,JW) + SIGMA(IC,JV,JW) 2300 CONTINUE C N = 0 DO 4160 ISAMPL=1,K N = N + NG(ISAMPL) 4160 CONTINUE XN = N C C COMPUTE COMSIG, MLE OF COMMON COVARIANCE MATRIX: DO 2400 JV = 1,IP DO 2400 JW = 1,IP COMSIG(JV,JW) = WGSS(JV,JW)/XN 2400 CONTINUE C WRITE (6,48500) DO 2410 JV = 1,IP WRITE (6,48000) ( COMSIG(JV,JW),JW=1,IP ) 2410 CONTINUE C C C C CALL SUBROUTINE TO COMPUTE DETERMINANT OF SAMPLE COVARIANCE C MATRIX: C IDET = 1 NRS1 = 0 DO 2405 JV = 1,IP DO 2405 JW = 1,IP COVMX(JV,JW) = COMSIG(JV,JW) 2405 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 SEE SUBROUTINE LISTING FOR FULLER EXPLANATION. C DET(VARHAT) = DET*10.0**IDET C WRITE (6,60000) JFLG C C XIDET = IDET C XLGDET = DLOG(DET) + XIDET*ALOG(10.0) TRUDET = DET*(10.0**IDET) IF (JFLG .GT. 0) WRITE (6,62000) DET,IDET C C COMPUTE LIKELIHOOD OF MODEL WITH COMMON COVARIANCE MATRIX C ACCORDING TO 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: XLL=-(XN*IP/2.0)*DLOG(2.0D+00*PI)-(XN/2.0)*DLOG(TRUDET)-XN*IP/2.0 C XMN2LL = -2.0*XLL C WRITE (6,68000) NOPARM = K*IP + IP*(IP+1)/2 WRITE (6,72000) NOPARM WRITE (6,76000) XMN2LL AIC = XMN2LL + 2.0*NOPARM WRITE (6,78000) AIC SCH = XMN2LL + ALOG(XN)*NOPARM WRITE (6,80000) SCH C C WRITE (6,84000) C C STOP C C C 10000 FORMAT('1','**********************************************'/// X' EVALUATION OF LIKELIHOOD IN ONE-WAY CLASSIFICATION '/ X' WITH COMMON AND VARYING COVARIANCE MATRICES '// X' CMS FILE: EVALUATE CLUSPAC '/// X' DEVELOPED AND PROGRAMMED BY: '// X' DR. STANLEY L. SCLOVE '/ X' INFORMATION & DECISION SCIENCES DEPT. M/C 294 '/ X' COLLEGE OF BUSINESS ADMINISTRATION '/ X' UNIVERSITY OF ILLINOIS AT CHICAGO '/ X' BOX 4348, CHICAGO, IL 60680 '// X' 312/996-2681 '// X' EVALUATE CLUSPAC VERSION 1.4 30-JAN-90 '// X' COPYRIGHT 1991 STANLEY LOUIS SCLOVE. ', X' '//) 10050 FORMAT(/' ................................................', X'................................'//) 11000 FORMAT(2X,I2) 11500 FORMAT(2X,'ENTER SUBROUTINES.') 11700 FORMAT(2X,'(1,1) ELEMENT OF INPUT MATRIX = ',E15.7) 11600 FORMAT(2X,'LEAVE SUBROUTINES.') 12000 FORMAT(2X,I4) 13000 FORMAT(1X,'NUMBER OF VARIABLES USED .................',I3/) 26000 FORMAT(1X,'NUMBER OF GROUPS, K ......................',I3/) 15000 FORMAT(1X,'NUMBER OF OBSERVATIONS (SAMPLE SIZE), N...',I3/) 18000 FORMAT(15(I5,I3)) 20000 FORMAT(/1X,' MEANS: '/) 25000 FORMAT(5X, F3.2) 21000 FORMAT(/ ' STATISTICS FOR THIS SAMPLE: ' ) 22000 FORMAT( ' ---------- --- ---- ------ ' /) 10100 FORMAT(' MEAN VECTOR: ',(9F8.2/)//) 10200 FORMAT(/' COVARIANCE MATRIX: '/) 28000 FORMAT(1X,'MEAN VECTOR FOR GROUP ',I2,': ',(8F13.5/)) 30000 FORMAT(/1X,' MINUS 2 LOG LIKELIHOOD = ', F13.5/) 36000 FORMAT(18A4) 40000 FORMAT(//,1X,18A4//) 42000 FORMAT(/,' COVARIANCE MATRICES: ',/) 42500 FORMAT(1X,' GROUP ',I2/) 48000 FORMAT(1X,8F13.5/) 48500 FORMAT(1X,'MLE OF COMMON COVARIANCE MATRIX: '/) 50000 FORMAT(1X,I4,1X,I2) 54000 FORMAT(3X,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 = ',E18.7,' IDET = ',I3,5X, X'ACTUAL DET. = DET*10**IDET',//) 64000 FORMAT(/1X,'MINIMUM FOR EACH VARIABLE: ',/) 66000 FORMAT(/1X,'MAXIMUM FOR EACH VARIABLE: ',/) 69000 FORMAT(1X,'VALUES OF MODEL SELECTION CRITERIA: '/) 69100 FORMAT(1X,'FOR MODEL WITH VARYING COVARIANCE MATRICES:--') 68000 FORMAT(1X,'FOR MODEL WITH COMMON COVARIANCE MATRIX:--') 70000 FORMAT(1X,'RESULTS FOR SAMPLE ',I2/) 72000 FORMAT(1X,'NUMBER OF PARAMETERS = ',5X,I4/) 76000 FORMAT(1X,' MINUS 2 LOG L = ', F15.5 /) 78000 FORMAT(1X,' AIC = ', F15.5 /) 80000 FORMAT(1X,' SCHWARZ CRITERION = ', F15.5 /) 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 ANIHILATED. C IF (A(IV(IM),K).EQ.0.0D+00) GO TO 1100 C C CACULATE 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 C CONTROL CARDS FOR PROGRAM EVALUATE CLUSPAC: C C C C //GO.SYSIN DD * C FISHER IRIS DATA C K=03 C IP=04 C (4(1X,F3.1),1X,I1) C N=0050 C C C THESE ARE: C JCL "GO"-STEP CARD C DATASET TITLE C C K, NUMBER OF GROUPS, IN FORMAT (2X,I2) C C NO. OF VARIABLES, IP, IN FORMAT (3X,I2) 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 N FOR FIRST SAMPLE, IN FORMAT (2X,I4) C DATA FOR FIRST SAMPLE C C (DATA, ONE CASE AT A TIME, IN FORMAT SPECIFIED BY DATFMT) C C N FOR SECOND SAMPLE, IN FORMAT (2X,I4) C DATA FOR SECOND SAMPLE C . C . C . C C N FOR K-TH SAMPLE, IN FORMAT (2X,I4) C DATA FOR K-TH SAMPLE C END //GO.SYSIN DD * FISHER IRIS DATA K=03 IP=04 (4(1X,F3.1)) N=0050 5.1 3.5 1.4 0.2 1 001 4.9 3.0 1.4 0.2 1 002 4.7 3.2 1.3 0.2 1 003 4.6 3.1 1.5 0.2 1 004 5.0 3.6 1.4 0.2 1 005 5.4 3.9 1.7 0.4 1 006 4.6 3.4 1.4 0.3 1 007 5.0 3.4 1.5 0.2 1 008 4.4 2.9 1.4 0.2 1 009 4.9 3.1 1.5 0.1 1 010 5.4 3.7 1.5 0.2 1 011 4.8 3.4 1.6 0.2 1 012 4.8 3.0 1.4 0.1 1 013 4.3 3.0 1.1 0.1 1 014 5.8 4.0 1.2 0.2 1 015 5.7 4.4 1.5 0.4 1 016 5.4 3.9 1.3 0.4 1 017 5.1 3.5 1.4 0.3 1 018 5.7 3.8 1.7 0.3 1 019 5.1 3.8 1.5 0.3 1 020 5.4 3.4 1.7 0.2 1 021 5.1 3.7 1.5 0.4 1 022 4.6 3.6 1.0 0.2 1 023 5.1 3.3 1.7 0.5 1 024 4.8 3.4 1.9 0.2 1 025 5.0 3.0 1.6 0.2 1 026 5.0 3.4 1.6 0.4 1 027 5.2 3.5 1.5 0.2 1 028 5.2 3.4 1.4 0.2 1 029 4.7 3.2 1.6 0.2 1 030 4.8 3.1 1.6 0.2 1 031 5.4 3.4 1.5 0.4 1 032 5.2 4.1 1.5 0.1 1 033 5.5 4.2 1.4 0.2 1 034 4.9 3.1 1.5 0.2 1 035 5.0 3.2 1.2 0.2 1 036 5.5 3.5 1.3 0.2 1 037 4.9 3.6 1.4 0.1 1 038 4.4 3.0 1.3 0.2 1 039 5.1 3.4 1.5 0.2 1 040 5.0 3.5 1.3 0.3 1 041 4.5 2.3 1.3 0.3 1 042 4.4 3.2 1.3 0.2 1 043 5.0 3.5 1.6 0.6 1 044 5.1 3.8 1.9 0.4 1 045 4.8 3.0 1.4 0.3 1 046 5.1 3.8 1.6 0.2 1 047 4.6 3.2 1.4 0.2 1 048 5.3 3.7 1.5 0.2 1 049 5.0 3.3 1.4 0.2 1 050 N=0050 7.0 3.2 4.7 1.4 2 051 6.4 3.2 4.5 1.5 2 052 6.9 3.1 4.9 1.5 2 053 5.5 2.3 4.0 1.3 2 054 6.5 2.8 4.6 1.5 2 055 5.7 2.8 4.5 1.3 2 056 6.3 3.3 4.7 1.6 2 057 4.9 2.4 3.3 1.0 2 058 6.6 2.9 4.6 1.3 2 059 5.2 2.7 3.9 1.4 2 060 5.0 2.0 3.5 1.0 2 061 5.9 3.0 4.2 1.5 2 062 6.0 2.2 4.0 1.0 2 063 6.1 2.9 4.7 1.4 2 064 5.6 2.9 3.6 1.3 2 065 6.7 3.1 4.4 1.4 2 066 5.6 3.0 4.5 1.5 2 067 5.8 2.7 4.1 1.0 2 068 6.2 2.2 4.5 1.5 2 069 5.6 2.5 3.9 1.1 2 070 5.9 3.2 4.8 1.8 2 071 6.1 2.8 4.0 1.3 2 072 6.3 2.5 4.9 1.5 2 073 6.1 2.8 4.7 1.2 2 074 6.4 2.9 4.3 1.3 2 075 6.6 3.0 4.4 1.4 2 076 6.8 2.8 4.8 1.4 2 077 6.7 3.0 5.0 1.7 2 078 6.0 2.9 4.5 1.5 2 079 5.7 2.6 3.5 1.0 2 080 5.5 2.4 3.8 1.1 2 081 5.5 2.4 3.7 1.0 2 082 5.8 2.7 3.9 1.2 2 083 6.0 2.7 5.1 1.6 2 084 5.4 3.0 4.5 1.5 2 085 6.0 3.4 4.5 1.6 2 086 6.7 3.1 4.7 1.5 2 087 6.3 2.3 4.4 1.3 2 088 5.6 3.0 4.1 1.3 2 089 5.5 2.5 4.0 1.3 2 090 5.5 2.6 4.4 1.2 2 091 6.1 3.0 4.6 1.4 2 092 5.8 2.6 4.0 1.2 2 093 5.0 2.3 3.3 1.0 2 094 5.6 2.7 4.2 1.3 2 095 5.7 3.0 4.2 1.2 2 096 5.7 2.9 4.2 1.3 2 097 6.2 2.9 4.3 1.3 2 098 5.1 2.5 3.0 1.1 2 099 5.7 2.8 4.1 1.3 2 100 N=0050 6.3 3.3 6.0 2.5 3 101 5.8 2.7 5.1 1.9 3 102 7.1 3.0 5.9 2.1 3 103 6.3 2.9 5.6 1.8 3 104 6.5 3.0 5.8 2.2 3 105 7.6 3.0 6.6 2.1 3 106 4.9 2.5 4.5 1.7 3 107 7.3 2.9 6.3 1.8 3 108 6.7 2.5 5.8 1.8 3 109 7.2 3.6 6.1 2.5 3 110 6.5 3.2 5.1 2.0 3 111 6.4 2.7 5.3 1.9 3 112 6.8 3.0 5.5 2.1 3 113 5.7 2.5 5.0 2.0 3 114 5.8 2.8 5.1 2.4 3 115 6.4 3.2 5.3 2.3 3 116 6.5 3.0 5.5 1.8 3 117 7.7 3.8 6.7 2.2 3 118 7.7 2.6 6.9 2.3 3 119 6.0 2.2 5.0 1.5 3 120 6.9 3.2 5.7 2.3 3 121 5.6 2.8 4.9 2.0 3 122 7.7 2.8 6.7 2.0 3 123 6.3 2.7 4.9 1.8 3 124 6.7 3.3 5.7 2.1 3 125 7.2 3.2 6.0 1.8 3 126 6.2 2.8 4.8 1.8 3 127 6.1 3.0 4.9 1.8 3 128 6.4 2.8 5.6 2.1 3 129 7.2 3.0 5.8 1.6 3 130 7.4 2.8 6.1 1.9 3 131 7.9 3.8 6.4 2.0 3 132 6.4 2.8 5.6 2.2 3 133 6.3 2.8 5.1 1.5 3 134 6.1 2.6 5.6 1.4 3 135 7.7 3.0 6.1 2.3 3 136 6.3 3.4 5.6 2.4 3 137 6.4 3.1 5.5 1.8 3 138 6.0 3.0 4.8 1.8 3 139 6.9 3.1 5.4 2.1 3 140 6.7 3.1 5.6 2.4 3 141 6.9 3.1 5.1 2.3 3 142 5.8 2.7 5.1 1.9 3 143 6.8 3.2 5.9 2.3 3 144 6.7 3.3 5.7 2.5 3 145 6.7 3.0 5.2 2.3 3 146 6.3 2.5 5.0 1.9 3 147 6.5 3.0 5.2 2.0 3 148 6.2 3.4 5.4 2.3 3 149 5.9 3.0 5.1 1.8 3 150 /* (1X,F3.1,1X,F3.1,1X,F3.1,1X,F3.1) 5.1 3.5 1.4 0.2 1 001 7.0 3.2 4.7 1.4 2 051 6.3 3.3 6.0 2.5 3 101 .33 = P(1) .33 = P(2) .34 = P(3) (1X, 5(4X,F9.5)) 1.00000 0.50000 0.50000 0.50000 0.50000 1.00000 0.50000 0.50000 0.50000 0.50000 1.00000 0.50000 0.50000 0.50000 0.50000 1.00000 1.00000 0.50000 0.50000 0.50000 0.50000 1.00000 0.50000 0.50000 0.50000 0.50000 1.00000 0.50000 0.50000 0.50000 0.50000 1.00000 1.00000 0.50000 0.50000 0.50000 0.50000 1.00000 0.50000 0.50000 0.50000 0.50000 1.00000 0.50000 0.50000 0.50000 0.50000 1.00000 /* (4(1X,F3.1),1X,I1) /*