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 ISDTPDTA C VERSION 3.1 19-SEP-82 C VERSION 3.1 CALCULATES BOTH AIC AND BIC (SCHWARZ'S C CRITERION) FOR MODEL SELECTION C C THE "ISDTP"(ISODATA PROCEDURE, P VARIABLES) PROGRAMS ARE C FOR CLUSTERING MULTIVARIATE DATA. (FOR UNIVARIATE DATA THE C "ISDT1" PROGRAMS MAY BE USED.) C ISODATA.EUCLID USES EUCLIDEAN DISTANCE. ISDTPCM C USES DISTANCE IN THE METRIC OF THE ESTIMATED COMMON COVARIANCE C MATRIX. ISDTPDF USES DIFFERENT COVARIANCE MATRICES. C ISDTPDT USES DIFFERENT COVARIANCE MATRICES, WITH ADJUSTMENT C BY THE DETERMINANTS, I.E., IT USES THE ESTIMATED LOG C LIKELIHOOD FOR THE GAUSSIAN MODEL WITH DIFFERENT COVARIANCE C MATRICES. C AUTOMATIC MODE: C A RANGE OF NUMBERS OF CLUSTERS IS TRIED, WITH AUTOMATIC SETTING OF C INITIAL MEANS. C C C PROGRAMMED BY C DR. STANLEY L. SCLOVE 312/996-2681 C QUANTITATIVE METHODS DEPARTMENT 312/996-2676 C UNIVERSITY OF ILLINOIS AT CHICAGO C BOX 4348, CHICAGO, IL 60680 C C RESEARCH SUPPORTED IN PART BY: C 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 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) 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 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 WRITE (6,40000) N C READ NUMBER OF VARIABLES, IP. READ (5,56000) IP WRITE (6,58000) IP C Read min and max no. of clusters: READ (5,13000) KLO READ (5,13000) KHI 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 C SET CONSTANTS. PI = 3.1415927 XN = N C SET VALUE OF PER-PARAMETER PENALTY IN MODEL SELECTION CRITERION C ( LOG(N) FOR SCHWARZ', 2 FOR AKAIKE'S ): XM = ALOG(XN) C C PERFORM CLUSTERING FOR VARIOUS NUMBERS OF CLUSTERS, K: DO 850 K = KLO,KHI WRITE (6,26000) K C WRITE (6,20000) C COMPUTE INITIAL MEANS DO 601 IVAR = 1,IP DO 601 IG = 1,K W = IG/(K+1.0) XMEAN(IG,IVAR) = (1-W)*XMIN(IVAR) + W*XMAX(IVAR) 601 CONTINUE C DO 600 IG=1,K WRITE (6,22000) WRITE (6,FMT) ( XMEAN(IG,IVAR), IVAR=1,IP) 600 CONTINUE C SET CONSTANTS FOR THIS VALUE OF K. 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 C ADJUST DISTANCE BY ADDING LOG OF DETERMINANT OF COVARIANCE MATRIX: 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 WRITE (6,16000) WRITE (6,18000) (I, ICLUS(I), I=1,N) 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 AIC. 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 C XM IS THE PER-PARAMETER PENALTY IN THE MODEL SELECTION CRITERION C (VALUE WAS SET ABOVE). AIC = XMN2LL + 2.0*NOPARM AICS = XMN2LL + XM*NOPARM WRITE (6,72000) AICS WRITE (6,72005) AIC 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 AIC BASED ON DIFFERENT COVARIANCE 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 AICDS = XM2LLD + XM*NPRMDF WRITE (6,32000) XM2LLD WRITE (6,76000) NPRMDF WRITE (6,74000) AICDS WRITE (6,74005) AICD 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 3900 SSD(L,IVAR,JV) = SSD(L,IVAR,JV)/(NG(L)-1) 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 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 C WRITE (6,80000) C C C GO ON TO NEXT NUMBER OF CLUSTERS, IF NOT FINISHED. 850 CONTINUE STOP 10000 FORMAT(2X,I2) 12000 FORMAT(2X,I4) 13000 FORMAT(4X,I2) 14000 FORMAT(//1X,'CLUSTERING:'/) 16000 FORMAT(/,1X,'CASES AND LABELS:--'/) 18000 FORMAT(15(I5,I3)) 20000 FORMAT(//1X,'INITIAL MEANS'/) 22000 FORMAT(1X, (8E13.5/)//) 24000 FORMAT('1','****************************************', X//,1X,'PROGRAM ISDT.DET '/ 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 3.2 4-AUG-83 '//, Y//,1X,'PROGRAM ISOPAC:ISDTPDTA'/ Z' COPYRIGHT 1991 STANLEY LOUIS SCLOVE. '//) 26000 FORMAT('1',//,1X,'K = ',I2,' CLUSTERS') 28000 FORMAT(1X,'MEAN VECTOR FOR CLUSTER ',I2,': ',(8E13.5/)) 30000 FORMAT(/,1X,'MINUS 2 LOG LIKELIHOOD FOR MODEL WITH COMMON', X' COVARIANCE MATRIX= ', E13.5//) 32000 FORMAT(/,1X,'MINUS 2 LOG LIKELIHOOD FOR MODEL WITH DIFFERENT ', X'COVARIANCE MATRICES = ', E13.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,8E13.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 = ',E13.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,'PROCEDURE HAS NOT CONVERGED IN 20 ITERATIONS. STOP') 72000 FORMAT(1X,'VALUE OF MODEL SELECTION CRITERION ', X'FOR MODEL WITH COMMON COVARIANCE MATRIX = ',F16.5/) 72005 FORMAT(1X,'VALUE OF AKAIKE INFORMATION CRITERION ', X'FOR MODEL WITH COMMON COVARIANCE MATRIX = ',F16.5/) 74000 FORMAT(1X,'VALUE OF MODEL SELECTION CRITERION ', X'FOR MODEL WITH DIFFERENT COVARIANCE MATRICES = ', XF16.5/) 74005 FORMAT(1X,'VALUE OF AKAIKE INFORMATION CRITERION ', X'FOR MODEL WITH DIFFERENT COVARIANCE MATRICES = ', XF16.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