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 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 BOX 4348, CHICAGO, IL 60680 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 //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 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 C 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