PDP-10 Archive: 43,50430/pca.f4 from decuslib10-06

Trailing-Edge - PDP-10 Archives - decuslib10-06 - 43,50430/pca.f4
There are no other files named pca.f4 in the archive.
	DIMENSION R(1),R1(1),R2(1),A(1),T(8),ST(13),IAA(1)
	COMMON KI,KO,INDSK,IODSK
	EQUIVALENCE (A,IAA)
	DIMENSION IPCA(4),IDISS(4),ICD(4)
	INTEGER RED1,RED2
C
C  PROGRAM TO DO PRINCIPAL COORDINATE AND COMPONENT ANALYSIS
C  USES SSP ROUTINES TO FIND EIGENVECTORS AND EIGENVALUES
C  OF THE SIMILARITY MATRICES CALCULATED BY BRACUR FOR COORD
C  AND CONDIT OUTPUT FOR COMPONENT ANALYSIS
C
	CALL INIT
	CALL GPCA(IDISS,IPCA,ICD)
	CALL IFILE(INDSK,'AVGE')
	READ(INDSK)I,J,K
	CALL RELEAS(INDSK)
C
C	CALCULATE MAXIMUM DIMENSIONS NEEDED FOR THE ARRAYS
C
	ICD(1)=MIN0(I,ICD(1))
	ICD(2)=MIN0(K,ICD(2))
	ICD(3)=MIN0(J,ICD(3))
	ICD(4)=MIN0(K,ICD(4))
	I=MAX0(I,J)
	I=MAX0(I,K)
	IRAMT=I*I
	IAAMT=(I*(I+1))/2
C
C	GET ENOUGH CORE FOR THE ARRAYS
C
	CALL MORCOR(R,IR,IRAMT)
	CALL MORCOR(A,IA,IAAMT)
	CALL MORCOR(R1,IR1,IRAMT)
	CALL MORCOR(R2,IR2,IRAMT)
C
C	NOW LOOP FOR AS MANY ARRAYS AS THERE ARE
C
	DO 5 J=1,4
	IF((IPCA(J).AND.1).EQ.0)GO TO 10
C
	ENCODE(5,1,FNAME)J
 1	FORMAT('DEN',I1,' ')
C
	CALL IFILE(INDSK,FNAME)
	READ(INDSK)NC
	READ(INDSK)T,ST
	CALL GETLT(R(IR),NC)
	CALL RELEAS(INDSK)
C
C	HAVING GOT THE SIMILARITY MATRIX WE NOW DO THE ANALYSIS
C
	CALL PCOORD(R(IR),A(IA),IAA(IA),NC,J,ST,ICD(J),IDISS)
C
 10	IF((IPCA(J).AND.2).EQ.0)GO TO 5
C
C  PRINCIPAL COMPONENT ANALYSIS
C
	ENCODE(5,12,FNAME)J
 12	FORMAT('MAT',I1,' ')
C
C  GET ARRAY AFTER CONDITIONING
C
	CALL IFILE(INDSK,FNAME)
	READ(INDSK)NR,NC
	CALL GETARR(R(IR),NC,NR)
	CALL RELEAS(INDSK)
C
C  NOW DO THE ANALYSIS
C
	CALL PCOMP(R(IR),A(IA),IAA(IA),R1(IR1),R2(IR2),NR,NC,J,ST)
C
 5	CONTINUE
C
	CALL EXIT
	END
	SUBROUTINE PCOORD(R,A,IA,NC,JJ,ST,NUMAX,IDISS)
C
C	SUBROUTINE TO DO PRINC COORD ANALYSIS
C
	DIMENSION R(NC,1),A(1),IA(1),T(8,4),ST(13),RT(4),IDISS(4)
	COMMON KI,KO,INDSK,IODSK
	DATA RT/'ENT1 ATTR ENT2 ATTR'/
	DATA (T(I,1),I=1,8)/'ENT1/ATTR ENT1 PRINC COORD ANALYSIS     '/
	DATA (T(I,2),I=1,8)/'ENT1/ATTR ATTR PRINC COORD ANALYSIS     '/
	DATA (T(I,3),I=1,8)/'ANT2/ATTR ENT2 PRINC COORD ANALYSIS     '/
	DATA (T(I,4),I=1,8)/'ENT2/ATTR ATTR PRINC COORD ANALYSIS     '/
C
C	SET UP THE DISSIMILARITY MATRIX
C
	DO 10 J=1,NC
	DO 10 K=1,J
	RTT=R(K,J)
	GO TO (5   ,   6   ,  5  ,  6  ,  5  ,  6  ),IDISS(JJ)
C             BC       MM     CM    DSQ   MTCH  RATIO
5	RTT=(100.0-RTT)/100.0
	GO TO 9
6	RTT=-0.5*RTT*RTT
	GO TO 9
9	R(K,J)=RTT
10	R(J,K)=RTT
C
C
C  NOW CALCULATE COLUMN AVERAGES AND OVERALL AVERAGE
C
	TSUM=0
	DO 15 I=1,NC
	SUM=0.0
	DO 17 J=1,NC
17	SUM=SUM+R(J,I)
	A(I)=SUM/NC
15	TSUM=TSUM+A(I)/NC
C
C	SCALE EACH ELEMENT OF R ACCORDINGLY
C
	DO 18 I=1,NC
	DO 18 J=1,I
18	R(J,I)=R(J,I)-A(I)-A(J)+TSUM
C
C	AND TRANSFER TO A FOR EIGEN TO OPERATE ON
C
	I=1
	DO 19 J=1,NC
	DO 19 K=1,J
	A(I)=R(K,J)
19	I=I+1
C
C	A NOW CONTAINS THE SIMILARITY MATRIX IN SYMMETRIC FORM
C	'COLUMN'-WISE
C	NOW CALCULATE THE EIGENVECTORS AND EIGENVALUES BASED ON IT
C
	CALL EIGEN(A,R,NC,0)
C
C	A CONTAINS EIGENVALUES IN DIAGONAL ELEMENTS
C	R CONTAINS EIGENVECTORS IN MY ROWS I.E. SSP'S COLS
C
C	COPY DIAGONAL OF A INTO FIRST NC LOCATIONS
C
	CALL DCPY(A,A,NC,1)
C
C	WRITE OUT THE EIGENVALUES
C
	WRITE(KO,20)(T(I,JJ),I=1,8),ST
 20	FORMAT(1H1,8A5,/,1X,13A5,/,' PRINCIPAL COORDINATE EIGENVALUES')
	WRITE(KO,25)(A(J),J=1,NC)
 25	FORMAT(12(1X,F9.3))
C
C	NOW NORMALIZE EIGENVECTORS SO THAT THEIR INDIVIDUAL SUM OF
C	SQUARES EQUALS THE CORRESPONDING EIGENVALUE
C
	DO 50 J=1,NC
	SSQ=0.0
	DO 30 I=1,NC
	AT=R(I,J)
 30	SSQ=SSQ+AT*AT
	IF(A(J).LT.0.AND.A(J).GT.-.001)A(J)=-A(J)
	RFACT=A(J)/SSQ
	RFACT=SQRT(RFACT)
	DO 35 I=1,NC
 35	R(I,J)=R(I,J)*RFACT
 50	CONTINUE
C
C	NOW R(I,J) CONTAINS THE ITH CO-ORDINATE W.R.T. NEW AXES
C	OF THE JTH ENTITY
C
	DO 60 I=1,NC
 60	IA(I)=I
C
	CALL MATWV(R,NC,NC,'COORD',RT(JJ),T(1,JJ),ST,IA,IA)
C
	RETURN
	END
	SUBROUTINE GETLT(R,NC)
C
C	A ROUTINE TO GET A LOWER TRIANGULAR ARRAY FROM INDSK INTO R
C	DIAGONAL ELEMENTS ARE SET TO 0.0
C
	COMMON KI,KO,INDSK
C
	DIMENSION R(NC,1)
	DO 5 I=2,NC
 5	READ(INDSK)(R(J,I),J=1,I-1)
	DO 10 I=1,NC
 10	R(I,I)=0.0
	RETURN
	END
	SUBROUTINE PCOMP(R,A,IA,R1,R2,NR,NC,JJ,ST)
C
C  PRINCIPAL COMPONENT ANALYSIS AFTER SEAL,BLACKITH+REYMENT
C	SOKAL+SNEATH AND UNCLE TOM COBBLY AND ALL
C
	DIMENSION R(NC,NR),R1(NR,NR),R2(NR,NR),A(1),ST(13),T(8,4),RT(4)
	DIMENSION IA(1)
	COMMON KI,KO,INDSK,IODSK
C
	DATA (T(I,1),I=1,8)/'ENT1/ATTR ENT1 PRINC COMP ANALYSIS      '/
	DATA (T(I,2),I=1,8)/'ENT1/ATTR ATTR PRINC COMP ANALYSIS      '/
	DATA (T(I,3),I=1,8)/'ENT2/ATTR ENT2 PRINC COMP ANALYSIS      '/
	DATA (T(I,4),I=1,8)/'ENT2/ATTR ATTR PRINC COMP ANALYSIS      '/
C
	DATA RT/'ENT1 ATTR ENT2 ATTR '/
C
C	STANDARDISE R SOTHAT IT HAS UNIT VARIANCE AND ZERO MEAN
C	FOR EACH ATTRIBUTE
C
	DO 55 I=1,NR
C
	NHERE=0
	SUM1=0
	SUM2=0
C
	DO 52 J=1,NC
	TR=R(J,I)
	IF(TR.EQ.99999)GO TO 52
	SUM1=SUM1+TR
	SUM2=SUM2+TR*TR
	NHERE=NHERE+1
52	CONTINUE
C
C   LEAVE A WHOLE ROW OF 99999'S AS IS
C
	IF(NHERE.GE.2)GO TO 51
	GO TO 55
C
51	RMEAN=SUM1/NHERE
	RDEV=(SUM2-(SUM1*SUM1)/NHERE)/(NHERE-1)
	RSD=SQRT(RDEV)
C
C   IF NO STD DEV MEAN WILL RESET ALL TO ZERO
C
	IF(RSD.EQ.0)RSD=1
C
	DO 53 J=1,NC
	IF(R(J,I).EQ.99999)GO TO 53
	R(J,I)=(R(J,I)-RMEAN)/RSD
53	CONTINUE
C
55	CONTINUE
C
C HAVING STANDARDISED R NOW FIND VARIANCE COVARIANCE MATRIX
C
	IACNT=1
	NEIG=MIN0(NR,NC-1)
C
	DO 20 I=1,NR
	DO 20 J=1,I
C
	NHERE=0
	SUMPRD=0
C
	DO 10 ICOL=1,NC
	EL1=R(ICOL,I)
	EL2=R(ICOL,J)
	IF(EL1.EQ.99999.OR.EL2.EQ.99999)GO TO 10
	NHERE=NHERE+1
	SUMPRD=SUMPRD+EL1*EL2
10	CONTINUE
	IF(NHERE.GE.2)GO TO 11
	A(IACNT)=0.0
	GO TO 20
C
11	A(IACNT)=SUMPRD/(NHERE-1)
20	IACNT=IACNT+1
C
C	NOW CALL EIGEN TO CALCULATE EVALS AND EVECTS OF A
C
	CALL EIGEN(A,R1,NR,0)
C
C	EVECTORS ARE IN R1 STORED IN COLS OF R1
C	EVALUES  ARE IN DIAGONAL OF A
C
	CALL DCPY(A,A,NR,1)
C
C	NOW MULTIPLY R1(NR*NR) INTO R(NC*NR) TO GET COORDINATES
C	OF ELEMENTS RELATIVE TO NEW  PRINCIPAL AXES
C	UNFORTUNATELY R1 IS STORED COL BY COL SO IT MUST BE RE-STORED
C	ROW BY ROW FOR MULTIPLICATION WITH R WHICH IS ROW BY ROW. SO
C	TRANSFER R1 TO R2 AND THEN PERFORM R1=R2*R
C
	CALL RCCON(R1,NR,NR,R2)
C
C	STANDARDISE EIGENVECTORS
C
	DO 27 J=1,NR
	SSQ=0
	DO 28 I=1,NR
28	SSQ=SSQ+R2(J,I)**2
	SSQ=SQRT(SSQ)
	DO 29 I=1,NR
29	R2(J,I)=R2(J,I)/SSQ
27	CONTINUE
C
C	NOW TRANSPOSE MATRIX R2
C	SO THAT EIGENVECTORS ARE IN ROWS
C
	DO 44 J=1,NR-1
	DO 44 I=J+1,NR
	TEMP=R2(I,J)
	R2(I,J)=R2(J,I)
44	R2(J,I)=TEMP
C
	CALL MULMAT(R2,R,R1,NR,NR,NC)
C
C	NOW HAVE
C		NEW COORDS IN R1
C		EIGENVECTORS IN R2
C		OLD COORDS IN R
C		EIGENVALS IN A
C
	WRITE(KO,30)(T(I,JJ),I=1,8),ST
30	FORMAT(1H1,8A5,/,1X,13A5,/,' PRINCIPAL COMPONENT',
	1	' ANALYSIS EIGENVALUES:',/)
	WRITE(KO,31)(A(J),J=1,NEIG)
31	FORMAT(12(1X,F9.4))
C
	DO 35 I=1,MAX0(NR,NC)
35	IA(I)=I
C
	CALL MATWV(R1,NEIG,NC,'COORD',RT(JJ),T(1,JJ),ST,IA,IA)
	CALL MATWV(R2,NEIG,NR,'NEWAX','WTOLD',T(1,JJ),ST,IA,IA)
C
	RETURN
	END