/*JOBPARM T=(0,40),L=3,R=2048 // EXEC FORTVCLG //FORT.SYSIN DD * C 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 ISDTPDT C VERSION 1.7 23-MAR-92 C C PROGRAM ISDTPDT (ISoDaTa, P-variate data, using C DeTerminants of group covariance matrices) IS ONE OF THE C "ISODATA" (ISDTP)PROGRAMS FOR CLUSTERING MULTIVARIATE C DATA. (FOR UNIVARIATE DATA THE C "ISDT1" PROGRAMS MAY BE USED.) C ISODATA.EUCLID USES EUCLIDEAN DISTANCE (FOR RESEARCH C PURPOSES, NOT RECOMMENDED FOR DATA ANALYSIS). ISDTPCM C USES DISTANCE IN THE METRIC OF THE ESTIMATED COMMON C COVARIANCE MATRIX. ISDTPDF USES DIFFERENT COVARIANCE C MATRICES FOR THE C CLUSTERS. ISDTPDT USES DIFFERENT COVARIANCE MATRICES, WITH C ADJUSTMENT BY THE DETERMINANTS, I.E., IT USES THE ESTIMATED LOG C LIKELIHOOD FOR THE GAUSSIAN MODEL WITH DIFFERENT COVARIANCE C MATRICES. C MANUAL MODE: NUMBER OF CLUSTERS AND INITIAL MEANS ARE INPUT C USE PROGRAMS ISDT***A (A=AUTOMATIC) TO TRY A RANGE OF NUMBERS OF C C CLUSTERS, WITH AUTOMATIC SETTING OF INITIAL MEANS. C C C C C C PROGRAMMED BY C C DR. STANLEY L. SCLOVE 312/996-2681 C C DEPT. OF INFORMATION & DECISION SCIENCES 312/996-2676 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 PROGRAM ISOPAC:ISDTPDT C C COPYRIGHT 1991 STANLEY LOUIS SCLOVE. C C C C C C RESEARCH SUPPORTED IN PART BY: C C C C ONR CONTRACT N00014-80-C-0408, TASK NR042-443 C C ARO CONTRACT DAAG29-82-K-0155 C C C C RESTRICTIONS (CAN BE MODIFIED): 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 SUBROUTINE(S) CALLED: C C MATEQ, WHICH CALLS MATDT C C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 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) C DIMENSION A(20,20) C DIMENSION ET(29) DIMENSION PG(29,20,20) C 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 DOUBLE PRECISION A C DOUBLE PRECISION ET DOUBLE PRECISION PG C C IV IS A WORK ARRAY FOR SUBROUTINE MATEQ. 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 READ (5,38000) TITLE C C WRITE PROGRAM INFORMATION. WRITE (6,24000) WRITE (6,42000) TITLE C C READ SAMPLE SIZE, N. READ (5,12000) N XN = N WRITE (6,40000) N C READ NUMBER OF VARIABLES, IP. READ (5,56000) IP WRITE (6,58000) IP C C READ DATA FORMAT. READ (5,38000) 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,66000) WRITE (6,FMT) (XMIN(IVAR),IVAR=1,IP) WRITE (6,68000) 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 C PARAMETERS FOR MODEL WITH COMMON COVARIANCE MATRIX ARE C K MEAN VECTORS AND ONE COVARIANCE MATRIX. NOPARM = K*IP + IP*(IP+1)/2 C PARAMETERS FOR MODEL WITH DIFFERENT COVARIANCE MATRICES ARE C K MEAN VECTORS AND K COVARIANCE MATRICES. NPRMDF = K*IP + K*IP*(IP+1)/2 C 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 1700 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 1500 1200 CONTINUE DO 1400 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) C D(L) = D(L) + TEMPIV*PG(L,IVAR,JV)*TEMPJV C 1300 CONTINUE D(L) = D(L) + ET(L) 1400 CONTINUE 1500 CONTINUE F = D(1) ICLUS(I) = 1 DO 1700 L = 2,K IF ( D(L) - F ) 1600,1700,1700 1600 F = D(L) ICLUS(I) = L 1700 CONTINUE WRITE (6,14000) C IF ( N .GE. 61 ) GO TO 1701 WRITE (6,16000) WRITE (6,18000) (I, ICLUS(I), I=1,N) 1701 CONTINUE C DO 1800 IG = 1,K NG(IG) = 0 DO 1800 IVAR = 1,IP SUM(IG,IVAR) = 0.0 DO 1800 JV = 1,IP SS(IG,IVAR,JV) = 0.0 SSD(IG,IVAR,JV) = 0.0 1800 CONTINUE DO 1900 I = 1,N IGROUP = ICLUS(I) NG(IGROUP) = NG(IGROUP) + 1 DO 1900 IVAR = 1,IP SUM(IGROUP,IVAR) = SUM(IGROUP,IVAR) + X(I,IVAR) DO 1900 JV = 1,IP SS(IGROUP,IVAR,JV) = SS(IGROUP,IVAR,JV) + X(I,IVAR) X *X(I,JV) 1900 CONTINUE DO 2000 IVAR = 1,IP DO 2000 JV = 1,IP WGSS(IVAR,JV) = 0.0 2000 CONTINUE C C DO 2300 IG = 1,K IF (NG(IG) .EQ. 0) GO TO 2100 GO TO 2200 2100 WRITE (6,78000) IG GO TO 3500 2200 CONTINUE DO 2300 IVAR = 1,IP XMEAN(IG,IVAR) = SUM(IG,IVAR)/NG(IG) DO 2300 JV = 1,IP CF = SUM(IG,IVAR)*SUM(IG,JV)/NG(IG) SSD(IG,IVAR,JV) = SS(IG,IVAR,JV) - CF 2300 CONTINUE C C C C POOL: C DO 2400 IG = 1,K DO 2400 IVAR = 1,IP DO 2400 JV = 1,IP WGSS(IVAR,JV) = WGSS(IVAR,JV) + SSD(IG,IVAR,JV) 2400 CONTINUE C C C COMPUTE VARHAT, MLE OF COMMON COVARIANCE MATRIX: DO 2500 IVAR = 1,IP DO 2500 JV = 1,IP VARHAT(IVAR,JV) = WGSS(IVAR,JV)/N 2500 CONTINUE WRITE (6,44000) DO 2600 IVAR=1,IP WRITE (6,50000) (VARHAT(IVAR,JV),JV=1,IP) 2600 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,62000) JFLG WRITE (6,64000) DET,IDET C WRITE (6,46000) DO 2700 IVAR=1,IP WRITE (6,50000) (P(IVAR,JV),JV=1,IP) 2700 CONTINUE XIDET = IDET XLGDET = DLOG(DET) + XIDET*ALOG(10.0) XMN2LL = N*(IP*ALOG(2.0*PI) + IP + XLGDET) C WRITE (6,30000) XMN2LL C COMPUTE MODEL-SELECTION CRITERIA (AIC AND SCHWARZ' CRITERION): C PARAMETERS: C K MEAN VECTORS OF DIMENSION P AND A P-BY-P COVARIANCE MATRIX, C WHERE P IS THE NUMBER OF VARIABLES WRITE (6,76000) NOPARM AIC = XMN2LL + 2.0*NOPARM WRITE (6,72000) AIC SCH = XMN2LL + ALOG(XN)*NOPARM WRITE (6,73000) SCH TERM=0.0 DO 3000 L=1,K DO 2800 IVAR=1,IP DO 2800 JV = 1,IP A(IVAR,JV) = SSD(L,IVAR,JV)/NG(L) 2800 CONTINUE IDET = 1 NRS1 = 0 CALL MATEQ(A,IP,20,JFLG,DET,IDET,IV,NRS1,P,20) DO 2900 IVAR=1,IP DO 2900 JV=1,IP PG(L,IVAR,JV) = P(IVAR,JV) 2900 CONTINUE C ET(L) = DLOG(DET*10**IDET) TERM = TERM + NG(L)*ET(L) 3000 CONTINUE C COMPUTE MODEL-SELECTION CRITERIA BASED ON DIFFERENT COVARIANCE C MATRICES. C PARAMETERS: C K MEAN VECTORS OF DIMENSION P AND K P-BY-P COVARIANCE MATRICES, C WHERE P IS THE NUMBER OF VARIABLES XM2LLD = N*IP*ALOG(2*PI) + N*IP + TERM AICD = XM2LLD + 2.0*NPRMDF SCHD = XM2LLD + ALOG(XN)*NPRMDF WRITE (6,32000) XM2LLD WRITE (6,76000) NPRMDF WRITE (6,74000) AICD WRITE (6,75000) SCHD C C WRITE (6,34000) ITER DO 3100 IG = 1,K WRITE (6,28000) IG, (XMEAN(IG,IVAR),IVAR=1,IP) 3100 CONTINUE IF (ITER .EQ. 1) GO TO 3300 DO 3200 I = 1,N IF (ICLUS(I) .EQ. ICLSOL(I)) GO TO 3200 GO TO 3300 3200 CONTINUE GO TO 3600 3300 CONTINUE ITER = ITER + 1 IF (ITER.GE.21) GO TO 3400 GO TO 700 3400 WRITE (6,70000) 3500 STOP C C 3600 CONTINUE C WRITE (6,60000) ITER WRITE (6,36000) (NG(IG),IG=1,K) DO 3700 IVAR = 1,IP DO 3700 JV = 1,IP WGMS(IVAR,JV) = WGSS(IVAR,JV)/(N-K) 3700 CONTINUE WRITE (6,48000) DO 3800 IVAR=1,IP WRITE (6,50000) (WGMS(IVAR,JV), JV=1,IP) 3800 CONTINUE C WRITE (6,82000) DO 3900 L=1,K DO 3900 IVAR=1,IP DO 3900 JV=1,IP C3900 SSD(L,IVAR,JV) = SSD(L,IVAR,JV)/(NG(L)-1) 3900 SSD(L,IVAR,JV) = SSD(L,IVAR,JV)/NG(L) DO 4000 L=1,K WRITE (6,86000) L DO 4000 IVAR=1,IP 4000 WRITE (6,84000) (SSD(L,IVAR,JV),JV=1,IP) C C C IF ( N .GE. 61) GO TO 4001 C WRITE (6,54000) DO 4100 I = 1,N WRITE (6,52000) I, ICLUS(I) WRITE (6,FMT) (X(I,IVAR),IVAR=1,IP) 4100 CONTINUE 4001 CONTINUE C 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 ISDTPDT '/ X,1X,'FOR CLUSTERING MULTIVARIATE DATA '/ X,1X,'USING DISTANCE IN THE METRICS OF THE COVARIANCE MATRICES'/ X,1X,'ADJUSTED BY THE DETERMINANTS '/ X//,1X,'DEVELOPED AND PROGRAMMED BY DR. STANLEY L. SCLOVE' X//,1X,'VERSION 1.7 23-MAR-92 '//, Y//,1X,'PROGRAM ISOPAC:ISDTPDT'/ 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 FOR MODEL WITH COMMON', X' COVARIANCE MATRIX= ', F13.5//) 32000 FORMAT(/,1X,'MINUS 2 LOG LIKELIHOOD FOR MODEL WITH DIFFERENT ', X'COVARIANCE MATRICES = ', F13.5//) 34000 FORMAT(///,1X,'ITERATION ', I2,//) 36000 FORMAT(/,1X,'NUMBERS:',3X,9(I10,3X)/) 38000 FORMAT(18A4) 40000 FORMAT(1X,'N = ',I3/) 42000 FORMAT(1X,18A4) 44000 FORMAT(///,1X,'COMMON COVARIANCE MATRIX (MLE):',//) 46000 FORMAT(//,1X,'INVERSE COVARIANCE MATRIX:',//) 48000 FORMAT(///,1X,'COMMON COVARIANCE MATRIX (UNBIASED ESTIMATE):',//) 50000 FORMAT(1X,8F13.5/) 52000 FORMAT(1X,I4,1X,I2) 54000 FORMAT(//,1X,'CASE, LABEL / DATA'//) 56000 FORMAT(3X,I2) 58000 FORMAT(1X,'NUMBER OF VARIABLES = ',I2/) 60000 FORMAT(/1X,'CONVERGENCE: NO CASE CHANGED CLUSTERS AFTER ', X'ITERATION ',I2,'. RESULTS ARE PRINTED BELOW.'//) 62000 FORMAT(/,1X,'JFLG = ',I2,'. IF JFLG=0, COMPUTATION OF DET', X' WENT WELL; OTHERWISE, THERE WAS TROUBLE OR MATRIX WAS ', X'ILL-CONDITIONED.'//) 64000 FORMAT(/1X,'DET = ',F13.5,' IDET = ',I3,5X, X'ACTUAL DET. = DET*10**IDET',//) 66000 FORMAT(/1X,'MINIMUM FOR EACH VARIABLE: ',/) 68000 FORMAT(/1X,'MAXIMUM FOR EACH VARIABLE: ',/) 70000 FORMAT(1X,'ROUTINE HAS NOT CONVERGED IN 20 ITERATIONS. STOP') 72000 FORMAT(1X,'AIC FOR MODEL WITH COMMON COVARIANCE MATRIX = ',F15.5/) 73000 FORMAT(1X,'SCHWARZ CRITERION ', X'FOR MODEL WITH COMMON COVARIANCE MATRIX = ',F15.5/) 74000 FORMAT(1X,'AIC FOR MODEL WITH DIFFERENT COVARIANCE MATRICES = ', XF15.5/) 75000 FORMAT(1X,'SCHWARZ CRITERION ', X'FOR MODEL WITH DIFFERENT COVARIANCE MATRICES = ', XF15.5/) 76000 FORMAT(/,1X,'NUMBER OF PARAMETERS = ',I4//) 78000 FORMAT(1X,'NO OBSERVATIONS IN GROUP ',I3,'. STOP') 80000 FORMAT(//,1X,'PROGRAM ENDED NORMALLY.') 82000 FORMAT(//,1X,'COVARIANCE MATRICES (UNBIASED ESTIMATES)'//) 84000 FORMAT(/1X,15F5.2/(3X,15F5.2/)) 86000 FORMAT(//,1X,'COVARIANCE MATRIX FOR CLUSTER ', I2/) END SUBROUTINE MATEQ(A,M,N,JFLG,DET,IDET,IV,NRS1,P,LL) C SUBROUTINE MATEQ IS DMATEQ FROM THE UICC SUBROUTINE LIBRARY. 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 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 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 = REALEST DATAORIG; 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 0.99 7.22 0.74 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.65 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.01 5.21 0.76 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 6.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 /*