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 PROGRAM ISDTPCM CLUSPAC C VERSION 5.4 19-SEP-91 C PREVIOUS UPDATE: VERSION 5.3 15-MAR-88 C C PROGRAMMED BY C DR. STANLEY L. SCLOVE 312/996-2681 C DEPARTMENT 312/996-2676 C OF INFORMATION AND DECISION SCIENCES M/C 294 C COLLEGE OF BUSINESS ADMINISTRATION C UNIVERSITY OF ILLINOIS AT CHICAGO C BOX 4348, CHICAGO, IL 60680 C C C C C PROGRAM ISDTPCM (ISoDaTa, P-variate data, CoMmon covariance C ** * * * * * C matrix) IN THE ISOPAC PACKAGE IS ONE OF THE "ISODATA" PROGRAMS C FOR CLUSTERING MULTIVARIATE DATA. (FOR UNIVARIATE DATA THE C "ISDT1" PROGRAMS MAY BE USED.) C USE ISDTPCM FOR DISTANCE IN THE METRIC OF THE COMMON C COVARIANCE MATRIX AND ISDTPDF FOR DIFFERENT COVARIANCE MATRICES. C ISDTPDT USES DIFFERENT COVARIANCE MATRICES, WITH ADJUSTMENT C BY THE DETERMINANTS, I.E., IT IS BASED ON THE GAUSSIAN LIKELIHOOD. C C MANUAL MODE: NUMBER OF CLUSTERS AND INITIAL MEANS ARE INPUT C USE PROGRAM ISODATA.COMMON.AUTO TO TRY A RANGE OF NUMBER OF C CLUSTERS, WITH AUTOMATIC SETTING OF INITIAL MEANS. C C C C C RESEARCH SUPPORTED IN PART BY C ONR CONTRACT N00014-80-C-0408, TASK NR042-443 C ARO CONTRACT DAAG29-82-K-0155 C C RESTRICTIONS (CAN BE MODIFIED): C N, SAMPLE SIZE, AT MOST 1000; C IP, NUMBER OF VARIABLES, AT MOST 20; C K, NUMBER OF CLUSTERS, AT MOST 29; C ITER, MAXIMUM NUMBER OF ITERATIONS, 20. C C SUBROUTINE(S) CALLED: C MATEQ, WHICH CALLS MATDT C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C DIMENSION X(1000,20),SUM(29,20) DIMENSION D(29),ICLUS(1000) DIMENSION TITLE(18) DIMENSION NG(29),XMEAN(29,20) DIMENSION FMT(18) DIMENSION SS(29,20,20),SSD(29,20,20) DIMENSION WGSS(20,20) DIMENSION VARHAT(20,20),WGMS(20,20) DIMENSION ICLSOL(1000) DIMENSION XMIN(20),XMAX(20) DIMENSION IV(20,20) DIMENSION P(20,20) DOUBLE PRECISION SS,SUM DOUBLE PRECISION WGSS,SSD DOUBLE PRECISION VARHAT DOUBLE PRECISION P DOUBLE PRECISION DET DOUBLE PRECISION D DOUBLE PRECISION XMEAN DOUBLE PRECISION TEMPIV,TEMPJV DOUBLE PRECISION F DOUBLE PRECISION CF C C IV IS A WORK ARRAY FOR SUBROUTINE MATEQ. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C CONTROL CARDS: C C DATASET TITLE C N, IN FORMAT (2X,I4) C IP, IN FORMAT (3X,I2) C FMT, IN FORMAT (18A4), E.G., (4F4.1) C "FMT" WILL ALSO BE USED FOR OUTPUT, SO ALLOW AT LEAST ONE BLANK C AT THE BEGINNING FOR CARRIAGE CONTROL. C DATA, ONE CASE AT A TIME, IN FORMAT SPECIFIED BY FMT C K, NUMBER OF CLUSTERS, IN FORMAT (2X,I2) C K INITIAL MEANS, IN FORMAT SPECIFIED BY FMT C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C READ (5,36000) TITLE C C WRITE PROGRAM INFORMATION. WRITE (6,24000) 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) FMT C C READ DATA. DO 500 I = 1,N READ (5,FMT) (X(I,IVAR), IVAR = 1,IP) IF (I .EQ. 1) GO TO 100 GO TO 300 100 CONTINUE DO 200 IVAR = 1,IP XMAX(IVAR) = X(1,IVAR) XMIN(IVAR) = X(1,IVAR) 200 CONTINUE 300 CONTINUE DO 400 IVAR = 1,IP IF (X(I,IVAR) .LT. XMIN(IVAR)) XMIN(IVAR)=X(I,IVAR) IF (X(I,IVAR) .GT. XMAX(IVAR)) XMAX(IVAR)=X(I,IVAR) 400 CONTINUE 500 CONTINUE WRITE (6,64000) WRITE (6,FMT) (XMIN(IVAR),IVAR=1,IP) WRITE (6,66000) WRITE (6,FMT) (XMAX(IVAR),IVAR=1,IP) C READ K, NUMBER OF CLUSTERS. READ (5,10000) K WRITE (6,26000) K C WRITE (6,20000) C READ INITIAL MEANS DO 600 IG=1,K READ (5,FMT) (XMEAN(IG,IVAR), IVAR=1,IP) C WRITE (6,22000) WRITE (6,FMT) ( XMEAN(IG,IVAR), IVAR=1,IP) 600 CONTINUE C SET CONSTANTS. PI = 3.1415927 ITER = 1 700 CONTINUE IF (ITER .EQ. 1) GO TO 900 DO 800 I = 1,N ICLSOL(I) = ICLUS(I) 800 CONTINUE C COMMENCE DISTANCE COMPUTATIONS. 900 CONTINUE DO 1600 I = 1,N DO 1000 L = 1,K D(L) = 0.0 1000 CONTINUE C FOR FIRST ITERATION, EUCLIDEAN DISTANCE IS USED BECAUSE C NO COVARIANCE MATRIX IS YET AVAILABLE. AFTER THE FIRST C ITERATION, DISTANCE WILL BE TAKEN IN THE METRIC OF THE C COVARIANCE MATRIX. C IF (ITER .GT. 1) GO TO 1200 DO 1100 L=1,K DO 1100 IVAR=1,IP D(L) = D(L) + ( XMEAN(L,IVAR) - X(I,IVAR) )**2 1100 CONTINUE GO TO 1400 1200 CONTINUE DO 1300 L=1,K DO 1300 IVAR=1,IP TEMPIV = XMEAN(L,IVAR) - X(I,IVAR) DO 1300 JV=1,IP TEMPJV = XMEAN(L,JV) - X(I,JV) D(L) = D(L) + TEMPIV*P(IVAR,JV)*TEMPJV 1300 CONTINUE 1400 CONTINUE F = D(1) ICLUS(I) = 1 DO 1600 L = 2,K IF ( D(L) - F ) 1500,1600,1600 1500 F = D(L) ICLUS(I) = L 1600 CONTINUE WRITE (6,14000) C WRITE (6,16000) WRITE (6,18000) (I, ICLUS(I), I=1,N) DO 1700 IG = 1,K NG(IG) = 0 DO 1700 IVAR = 1,IP SUM(IG,IVAR) = 0.0 DO 1700 JV = 1,IP SS(IG,IVAR,JV) = 0.0 SSD(IG,IVAR,JV) = 0.0 1700 CONTINUE DO 1800 I = 1,N IGROUP = ICLUS(I) NG(IGROUP) = NG(IGROUP) + 1 DO 1800 IVAR = 1,IP SUM(IGROUP,IVAR) = SUM(IGROUP,IVAR) + X(I,IVAR) DO 1800 JV = 1,IP SS(IGROUP,IVAR,JV) = SS(IGROUP,IVAR,JV) + X(I,IVAR) X *X(I,JV) 1800 CONTINUE DO 1900 IVAR = 1,IP DO 1900 JV = 1,IP WGSS(IVAR,JV) = 0.0 1900 CONTINUE C C DO 2200 IG = 1,K IF (NG(IG) .EQ. 0) GO TO 2000 GO TO 2100 2000 WRITE (6,74000) IG GO TO 3200 2100 CONTINUE DO 2200 IVAR = 1,IP XMEAN(IG,IVAR) = SUM(IG,IVAR)/NG(IG) DO 2200 JV = 1,IP CF = SUM(IG,IVAR)*SUM(IG,JV)/NG(IG) SSD(IG,IVAR,JV) = SS(IG,IVAR,JV) - CF 2200 CONTINUE C C C C POOL: C DO 2300 IG = 1,K DO 2300 IVAR = 1,IP DO 2300 JV = 1,IP WGSS(IVAR,JV) = WGSS(IVAR,JV) + SSD(IG,IVAR,JV) 2300 CONTINUE C C C COMPUTE VARHAT, MLE OF COMMON COVARIANCE MATRIX: DO 2400 IVAR = 1,IP DO 2400 JV = 1,IP VARHAT(IVAR,JV) = WGSS(IVAR,JV)/N 2400 CONTINUE WRITE (6,42000) DO 2500 IVAR=1,IP WRITE (6,48000) (VARHAT(IVAR,JV),JV=1,IP) 2500 CONTINUE IDET = 1 NRS1 = 0 CALL MATEQ(VARHAT,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 WRITE (6,62000) DET,IDET C WRITE (6,44000) DO 2600 IVAR=1,IP WRITE (6,48000) (P(IVAR,JV),JV=1,IP) 2600 CONTINUE XIDET = IDET XLGDET = DLOG(DET) + XIDET*ALOG(10.0) XMN2LL = N*(IP*ALOG(2.0*PI) + IP + XLGDET) C DO 2700 IVAR = 1,IP DO 2700 JV = 1,IP WGMS(IVAR,JV) = WGSS(IVAR,JV)/(N-K) 2700 CONTINUE WRITE (6,30000) XMN2LL C WRITE (6,32000) ITER DO 2800 IG = 1,K WRITE (6,28000) IG, (XMEAN(IG,IVAR),IVAR=1,IP) 2800 CONTINUE 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.21) GO TO 3100 GO TO 700 3100 WRITE (6,68000) 3200 STOP C C 3300 CONTINUE C WRITE (6,58000) ITER WRITE (6,34000) (NG(IG),IG=1,K) WRITE (6,46000) DO 3400 IVAR=1,IP WRITE (6,48000) (WGMS(IVAR,JV), JV=1,IP) 3400 CONTINUE C C C WRITE (6,52000) DO 3500 I = 1,N WRITE (6,50000) I, ICLUS(I) WRITE (6,FMT) (X(I,IVAR),IVAR=1,IP) 3500 CONTINUE 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 NOPARM = K*IP + IP*(IP+1)/2 WRITE (6,72000) NOPARM AIC = XMN2LL + 2.0*NOPARM SCH = XMN2LL + ALOG(XN)*NOPARM WRITE (6,70000) AIC WRITE (6,71000) SCH WRITE (6,80000) C C STOP 10000 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,'INITIAL MEANS'/) 22000 FORMAT(1X, (8F13.5/)//) 24000 FORMAT('1',' ', X//,1X,'PROGRAM ISDTPCM CLUSPAC'/ X,1X,'FOR CLUSTERING MULTIVARIATE DATA '/ X,1X,'USING DISTANCE IN THE METRIC OF THE COVARIANCE MATRIX'/ X//,1X,'DEVELOPED AND PROGRAMMED BY DR. STANLEY L. SCLOVE' X//,1X,'VERSION 5.4 19-SEP-91 '// Y//,1X,'PREVIOUS UPDATE: VERSION 5.3 15-MAR-88'/ Z' COPYRIGHT 1991 STANLEY LOUIS SCLOVE'//) 26000 FORMAT('1',//,1X,'K = ',I2,' CLUSTERS') 28000 FORMAT(1X,'MEAN VECTOR FOR CLUSTER ',I2,': ',(8F13.5/)) 30000 FORMAT(/1X,' MINUS 2 LOG LIKELIHOOD = ', F13.5//) 32000 FORMAT(/,1X,'ITERATION ', I2,//) 34000 FORMAT(/,1X,'NUMBERS:',3X,9(I10,3X)/) 36000 FORMAT(18A4) 38000 FORMAT(1X,'N = ',I3/) 40000 FORMAT(1X,18A4) 42000 FORMAT(///,1X,'COMMON COVARIANCE MATRIX (MLE):',//) 44000 FORMAT(///,1X,'INVERSE COVARIANCE MATRIX:',//) 46000 FORMAT(///,1X,'COMMON COVARIANCE MATRIX (UNBIASED ESTIMATE):',//) 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,'ISODATA HAS NOT CONVERGED IN 20 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,'NO OBSERVATIONS IN GROUP ',I3,'. STOP') 80000 FORMAT(1X,'PROGRAM ENDED NORMALLY.') END 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 ELIMINTING WE CHECK TO SEE IF A(IV(IM),K) HAS ALREADY C BEEN ANIHALATED. 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 END //GO.SYSIN DD * CMS DSN = RLESTATE DATAMLLR; Xs: STYLE;SFLA;GRADE;ASSESSED;MARKET N=0060 IP=05 (1X,F5.2,1X,F5.2,1X,F4.2,1X,F5.2,1X,F4.1) 1.25 10.11 0.90 39.00 57.6 CASE 01 1.00 7.21 0.90 32.40 49.2 CASE 02 1.00 7.60 0.90 39.90 53.7 CASE 03 1.00 6.62 0.90 29.40 51.9 CASE 04 1.25 9.31 0.90 36.90 57.0 CASE 05 1.00 7.22 0.75 11.10 38.4 CASE 06 1.75 16.02 1.00 46.50 71.1 CASE 07 1.00 6.77 0.90 34.20 54.6 CASE 08 1.25 10.47 0.90 45.00 65.4 CASE 09 1.00 8.38 0.90 38.40 53.4 CASE 10 1.25 9.18 0.90 48.30 61.8 CASE 11 1.00 7.87 0.90 33.60 58.8 CASE 12 1.75 8.58 0.75 33.60 53.7 CASE 13 1.25 11.71 0.90 48.30 72.6 CASE 14 1.00 8.99 0.90 34.80 56.4 CASE 15 1.00 5.77 0.75 33.30 39.3 CASE 16 1.00 5.38 0.75 42.30 29.1 CASE 17 1.25 10.56 0.90 38.70 62.7 CASE 18 1.25 8.83 0.90 38.70 59.4 CASE 19 1.25 10.33 0.90 37.80 60.6 CASE 20 1.25 10.75 0.90 45.60 64.2 CASE 21 1.25 10.70 0.90 34.20 65.4 CASE 22 1.00 10.60 0.90 36.00 57.6 CASE 23 1.25 11.38 0.90 38.40 58.5 CASE 24 1.00 5.44 0.90 34.80 37.8 CASE 25 1.00 8.25 0.90 42.30 57.0 CASE 26 2.00 14.35 0.90 32.10 58.2 CASE 27 1.00 8.34 0.75 27.00 51.9 CASE 28 1.00 7.80 0.75 33.90 46.5 CASE 29 1.00 7.67 0.75 33.30 50.4 CASE 30 1.25 10.11 0.90 42.00 65.4 CASE 31 1.25 10.24 0.90 40.50 63.3 CASE 32 1.00 6.91 0.75 33.90 49.5 CASE 33 1.00 9.66 0.90 45.30 66.0 CASE 34 1.25 11.06 0.90 47.40 64.2 CASE 35 1.00 11.00 0.90 45.00 62.4 CASE 36 1.00 8.90 0.90 30.30 52.5 CASE 37 1.00 9.20 0.75 31.20 46.8 CASE 38 2.00 18.04 0.90 45.90 68.1 CASE 39 1.12 9.65 0.90 21.90 55.8 CASE 40 1.25 6.61 0.75 35.70 46.5 CASE 41 1.75 11.64 1.00 44.10 62.7 CASE 42 1.25 10.79 0.90 36.30 60.9 CASE 43 1.00 8.15 0.90 21.60 51.6 CASE 44 1.75 10.60 1.00 43.80 67.2 CASE 45 2.00 12.98 1.00 35.40 67.8 CASE 46 2.00 12.49 0.90 33.00 70.8 CASE 47 1.00 8.14 0.90 22.20 52.2 CASE 48 1.00 5.21 0.75 11.70 39.0 CASE 49 1.50 10.40 0.90 33.60 60.6 CASE 50 1.25 10.51 0.90 39.60 62.1 CASE 51 1.00 7.12 0.90 31.80 63.6 CASE 52 1.00 6.94 0.90 42.00 56.1 CASE 53 1.50 10.52 0.90 30.30 59.4 CASE 54 1.00 9.23 0.75 36.90 45.0 CASE 55 1.00 9.67 0.90 39.00 66.3 CASE 56 1.00 7.43 0.90 39.30 52.2 CASE 57 1.50 12.37 0.90 25.50 59.7 CASE 58 1.00 9.26 0.90 27.33 56.1 CASE 59 1.00 8.02 0.90 38.10 54.0 CASE 60 K=03 1.00 7.00 0.50 15.00 30.0 1.10 9.00 0.75 30.00 50.0 1.20 11.00 1.00 45.00 70.0 /*