C cancor.f C The authors of this software are Jih-Jie Chang and J D Carroll. C Copyright (c) 1993 by AT&T. C Permission to use, copy, modify, and distribute this software for any C purpose without fee is hereby granted, provided that this entire C notice is included in all copies of any software which is or includes C a copy or modification of this software and in all copies of the C supporting documentation for such software. C THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMP- C LIED WARRANTY. IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY C REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY C OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE. C This software comes from the FIRST MDS Package of AT&T Bell Laboratories C For explanation of the method this software implements, see C Carroll, J.D. (1968), "Generalization of canonical correlation C analysis to three or more sets of variables" in Proceedings of the C 76th Annual Convention of the American Psychological Association, 3, C 227-228. C----+----@----+----@----+----@----+----@----+----@----+----@----+----@ $ FORTRAN $ INCODE IBMF C CANCOR - CANONICAL CORRELATION ANALYSIS C OF THREE OR MORE SETS OF VARIABLES C CHANG FOR CARROLL, FEBRUARY 1969 C NSET - NUMBER OF SETS (MAX=100) C NS - NUMBER OF STIMULI (MAX=100). C NCS - NUMBER OF COMPONENTS WITHIN EACH SET. C NC - NUMBER OF COMPONENTS TO BE COMPUTED C WF - WEIGHTING FACTOR, ONE FOR EACH SET. C IRX - =1 THE SET IS PUNCHED NCS BY NS C =0 THE SET IS PUNCHED NS BY NCS C ICOR - =1 COMPUTE CORRELATIONS BETWEEN SETS C =0 DO NOT COMPUTE CORRELATIONS BETWEEN SETS C X - DATA, STORES ONE SET OF VARIABLE. SCRATCH DISK IS USED TO C STORE ALL SETS. C TIT - TITLE C SNAME - SUB-TITLE, ONE FOR EACH SET. C FMT - VARIABLE FORMAT FOR READING DATA C CORL - CORRELATIONS BETWEEN EACH SET AND THE CANONICAL COMPONENTS DIMENSION X(100,10),Q(100,100),Z(100,100),Y(100,10),C(100,10) DIMENSION NCS(100),SNAME(12,100), FMT(12),TIT(12) DIMENSION CORL(10,100),SUMZ(10),SSQZ(10) EQUIVALENCE(C,Z) ND1=100 ND2=10 NDD=1000 NDH=11 500 READ100,NSET,NS,NC,ICOR IF(NSET.EQ.0) GOTO999 XNS=NS READ101,(TIT(I),I=1,12) PRINT102,(TIT(I),I=1,12) DO5I=1,100 DO5J=1,100 5 Q(I,J)=0. PRINT113 113 FORMAT(13H0I. DATA SET) DO10K=1,NSET READ101,(SNAME(I,K),I=1,12) READ112,NCS(K),IRX,WF M=NCS(K) CALL READX(X,NS,M,IRX,ND1) C X SUBTRACT XMEAN DO15I=1,M SUM=0. DO16J=1,NS 16 SUM=SUM&X(J,I) XMEAN=SUM/XNS DO17J=1,NS 17 X(J,I)=X(J,I)-XMEAN 15 CONTINUE C PUT X ON DISK 19 IF(K.EQ.1) GOTO11 CALL DOUC(19,X,NDD) GOTO13 11 CALL DOUTAT(19,X,NDD,1) 13 PRINT103,K,(SNAME(I,K),I=1,12) PRINT104 C PRINT X DO18I=1,M 18 PRINT105,I,(X(J,I),J=1,NS) C COMPUTE (XT)X AND STORE IN Y DO19I=1,M DO19L=1,M Y(I,L)=0. DO19J=1,NS 19 Y(I,L)=Y(I,L)&X(J,I)*X(J,L) C INVERT (XT)X AND STORE IN Y MM=1 CALL MTINV(Y,ND1,ND2,C,NDH,NDH,M,MM) IF(MM.EQ.1)GOTO20 PRINT106,MM,K STOP C COMPUTE X*((XT)X) INVERSE 20 DO21I=1,NS DO21J=1,M C(I,J)=0. DO21L=1,M 21 C(I,J)=C(I,J)&X(I,L)*Y(L,J) C PUT C ON DISK 19 CALL DOUC(19,C,NDD) C COMPUTE Q MATRIX DO22I=1,NS DO22J=1,NS SUM=0. DO23L=1,M 23 SUM=SUM&C(I,L)*X(J,L) 22 Q(I,J)=Q(I,J)&SUM*WF 10 CONTINUE CALL DROUT C FACTOR Q MATRIX MM=1 CALL MTEVV(Q,ND1,ND1,Z,ND1,ND1,NS,MM) PRINT107 PRINT108,(Q(I,I),I=1,NS) PRINT109 DO24I=1,NC C SET UP FOR CORRELATIONS BETWEEN Z AND EACH SET. SUMZ(I)=0. SSQZ(I)=0. DO25J=1,NS SUMZ(I)=SUMZ(I)&Z(J,I) 25 SSQZ(I)=SSQZ(I)&Z(J,I)**2 24 PRINT105,I,(Z(J,I),J=1,NS) C COMPUTE A(I), I FOR I(TH) SET, STORE A IN Q PRINT114 DO30K=1,NSET M=NCS(K) C READ X AND C FROM DISK. PUT C IN Y IF(K.GT.1)GOTO31 CALL DINAT(19,X,NDD,1) GOTO32 31 CALL DINC(19,X,NDD) 32 CALL DINC(19,Y,NDD) C COMPUTE A AND STORE A IN Q DO35I=1,NC SUM=0. DO37J=1,M Q(I,J)=0. DO36L=1,NS 36 Q(I,J)=Q(I,J)&Z(L,I)*Y(L,J) 37 SUM=SUM&Q(I,J)**2 SUM=SQRT(SUM) DO38LL=1,M 38 Q(I,LL)=Q(I,LL)/SUM 35 CONTINUE C PRINT CANONICAL DIRECTIONS PRINT103,K,(SNAME(I,K),I=1,12) PRINT110 DO39I=1,NC 39 PRINT105,I,(Q(I,J),J=1,M) C COMPUTE ROTATION MATRIX AND STORE IN Y DO40I=1,NC SXY=0. SUMY=0. SSQY=0. DO47J=1,NS Y(J,I)=0. DO45L=1,M 45 Y(J,I)=Y(J,I)&Q(I,L)*X(J,L) C COMPUTE CORRELATIONS BETWEEN Z AND EACH SET OF CANONICAL VARIATE SUMY=SUMY&Y(J,I) SSQY=SSQY&Y(J,I)**2 47 SXY=SXY&Y(J,I)*Z(J,I) PP=XNS*SXY-SUMY*SUMZ(I) QQ=(XNS*SSQY-SUMY**2)*(XNS*SSQZ(I)-SUMZ(I)**2) QQ=SQRT(QQ) CORL(I,K)=PP/QQ 40 CONTINUE IF(ICOR.EQ.0) GOTO43 C PUT Y ON DISK 20 IF(K.EQ.1) GOTO42 CALL DOUC(20,Y,NDD) GOTO43 42 CALL DOUTAT(20,Y,NDD,1) C PRINT Y 43 PRINT111 DO41I=1,NC 41 PRINT105,I,(Y(J,I),J=1,NS) 30 CONTINUE PRINT115 PRINT116 DO46I=1,NSET 46 PRINT117,I,(CORL(J,I),J=1,NC) IF(ICOR)44,500,44 44 CALL DROUT C COMPUTE CORRELATIONS BETWEEN SETS FOR EACH COMPONENTS CALL CCOR(Y,NC,NS,NDD,X,Q,ND1,NSET,XNS) GOTO500 100 FORMAT(18I4) 101 FORMAT(12A6) 102 FORMAT(1H112A6) 103 FORMAT(4H0SETI3,3X,12A6) 104 FORMAT(5X,26HDATA (SUBTRACTED ROW MEAN)) 105 FORMAT(1H04X,I2,2X,10E12.4/(9X,10E12.4)) 106 FORMAT(4H0MM=I1,32HMATRIX CANNOT BE INVERTED IN SETI3) 107 FORMAT(41H0II. CANONICAL VALUES (EIGEN ROOTS OF Q)) 108 FORMAT(5X,10E12.4) 109 FORMAT(47H0 CANONICAL COMPONENTS (EIGEN VECTORS OF Q)) 110 FORMAT(23H0 CANONICAL DIRECTIONS) 111 FORMAT(28H0 CANONICAL ROTATION MATRIX) 112 FORMAT(2I2,F10.5) 114 FORMAT(49H0III. CANONICAL DIRECTIONS AND ROTATION MATRICES) 115 FORMAT(89H0IV. CORRELATIONS BETWEEN EACH SET OF CANONICAL VARIATE 1 AND THE CANONICAL COMPONENTS (Z)) 116 FORMAT(5H0 SET,14X,10HCOMPONENTS/6X,I9,9I12) 117 FORMAT(I4,5X,10E12.4) 999 STOP END $ FORTRAN SUBROUTINE CCOR(Y,NC,NS,NDD,X,Q,ND1,NSET,XNS) DIMENSION Y(ND1,NC),X(ND1,NC),Q(ND1,ND1) DIMENSION SUMY(10),SSQY(10) N=NSET-1 IC=0 NY=1 DO10I=1,N CALL DINAT(20,Y,NDD,NY) NY=NY&NDD DO15J=1,NC SUMY(J)=0. SSQY(J)=0. DO15L=1,NS SUMY(J)=SUMY(J)&Y(L,J) 15 SSQY(J)=SSQY(J)&Y(L,J)**2 II=I&1 DO20JJ=II,NSET IC=IC&1 CALL DINC(20,X,NDD) DO20J=1,NC SUMX=0. SSQX=0. SXY=0. DO21L=1,NS SUMX=SUMX&X(L,J) SSQX=SSQX&X(L,J)**2 21 SXY=SXY&X(L,J)*Y(L,J) PP=XNS*SXY-SUMY(J)*SUMX QQ=(XNS*SSQY(J)-SUMY(J)**2)*(XNS*SSQX-SUMX**2) QQ=SQRT(QQ) 20 Q(J,IC)=PP/QQ 10 CONTINUE PRINT100 100 FORMAT(50H0V. CORRELATIONS BETWEEN SETS FOR EACH COMPONENT) PRINT101,(I,I=1,NC) 101 FORMAT(6H0 SETS13X,10HCOMPONENTS/6H I J,I9,9I12) IC=0 DO30I=1,N II=I&1 DO30J=II,NSET IC=IC&1 30 PRINT102,I,J,(Q(L,IC),L=1,NC) 102 FORMAT(2I3,3X,10E12.4) RETURN END $ FORTRAN READX(X,N,K,IRX,NL) SUBROUTINE READX(X,N,K,IRX,NL) DIMENSION X(NL,K),FMT(12) READ100,(FMT(I),I=1,12) 100 FORMAT(12A6) IF(IRX.EQ.0) GOTO10 DO6I=1,K 6 READFMT,(X(J,I),J=1,N) GOTO20 10 DO8I=1,N 8 READFMT,(X(I,J),J=1,K) 20 RETURN END $ FORTRAN EIGJ00a0 $ INCODE IBMF EIGJ00b0 *EIGENJ EIGENVALUE AND EIGENVECTOR SUBROUTINE EIGJ0020 * CD600D4.009 DATE 05/05/65 EIGJ0030 SUBROUTINE EIGENJ (A,V,NN,IV,E,I3,I4) EIGJ00e0 DIMENSION A(5050),V( I3, I4),E(NN) EIGJ00f0 LOGICAL IV EIGJ0060 * EIGENJ FINDS THE EIGENVALUES AND EIGENVECTORS OF EIGJ0070 * A SYMMETRIC MATRIX (A) USING A MODIFIED THRESHOLD JACOBI METHOD EIGJ0080 * ONLY THE UPPER TRIANGLE IS STORED EIGJ0090 * DIMENSION A(K=M(M&1)/2),V(M,M),E(M) WHERE M=MAXIMUM ORDER EIGJ0100 * DIMENSION A(K=M(M&1)/2),V(M,M),E(M) WHERE M=MAXIMUM ORDER EIGJ0110 N=NN EIGJ0120 IND=1 EIGJ0130 NM2=N-2 EIGJ0140 AMAX=0.0 EIGJ0150 NM=N-1 EIGJ0160 IF(.NOT.IV)GO TO 5 EIGJ0170 * SET UP VECTOR MATRIX EIGJ0180 1 DO 3 I=1,N EIGJ0190 DO 2 J=1,N EIGJ0200 2 V(I,J)=0.0 EIGJ0210 3 V(I,I)=1.0 EIGJ0220 5 K=2 EIGJ0230 DO 28 I=1,NM EIGJ0240 IP=I&1 EIGJ0250 DO 27 J=IP,N EIGJ0260 Y=A(K) EIGJ0270 X=ABS(Y) EIGJ0280 IF(X-AMAX)27,27,6 EIGJ0290 6 GO TO (7,8),IND EIGJ0300 7 AMAX=X EIGJ0310 ITEST=1 EIGJ0360 X=A(II)-A(JJ) EIGJ0370 IF(X)10,11,10 EIGJ0380 10 X=Y/X EIGJ0390 Y=.5*X EIGJ0400 IF(ABS(Y)-.414213562)14,11,11 EIGJ0410 11 C=.707106781 EIGJ0420 IF(Y)12,12,13 EIGJ0430 12 S=-C EIGJ0440 GO TO 15 EIGJ0450 13 S=C EIGJ0460 GO TO 15 EIGJ0470 14 X=Y*Y EIGJ0480 C=1.&X EIGJ0490 S=2.*Y/C EIGJ0500 C=(1.-X)/C EIGJ0510 15 X=S*S EIGJ0520 Y=C*C EIGJ0530 XY=S*C EIGJ0540 AXY=2.*A(K)*XY EIGJ0550 Z=A(II)*Y&AXY&A(JJ)*X EIGJ0560 W=A(II)*X-AXY&A(JJ)*Y EIGJ0570 A(K)=A(K)*(Y-X)&XY*(A(JJ)-A(II)) EIGJ0580 A(II)=Z EIGJ0590 A(JJ)=W EIGJ0600 IF(NM2)24,24,16 EIGJ0610 16 IT=I EIGJ0620 JT=J EIGJ0630 * TRANSFORM A EIGJ0640 DO 23 M=1,NM2 EIGJ0650 NP=N-M EIGJ0660 IF(JT-K)20,17,18 EIGJ0670 17 IT=IT&1 EIGJ0680 JT=JT&NP EIGJ0690 ITEST=2 EIGJ0700 18 IF(IT-K)20,19,19 EIGJ0710 19 IT=IT&1 EIGJ0720 JT=JT&1 EIGJ0730 ITEST=3 EIGJ0740 20 X=A(IT)*C&A(JT)*S EIGJ0750 A(JT)=A(JT)*C-A(IT)*S EIGJ0760 A(IT)=X EIGJ0770 GO TO (21,22,23),ITEST EIGJ0780 21 IT=IT&NP EIGJ0790 JT=JT&NP EIGJ0800 GO TO 23 EIGJ0810 22 IT=IT&1 EIGJ0820 JT=JT&NP-1 EIGJ0830 23 CONTINUE EIGJ0840 24 IF(.NOT.IV)GO TO 27 EIGJ0850 * TRANSFORM V EIGJ0860 25 DO 26 M=1,N EIGJ0870 X=V(M,I)*C&V(M,J)*S EIGJ0880 V(M,J)=V(M,J)*C-V(M,I)*S EIGJ0890 26 V(M,I)=X EIGJ0900 27 K=K&1 EIGJ0910 28 K=K&1 EIGJ0920 * SEQUENCE TO ADJUST THRESHOLD FOR EACH SWEEP EIGJ0930 GO TO (29,31),IND EIGJ0940 29 IF(AMAX)30,35,30 EIGJ0950 30 AT=AMAX EIGJ0960 IND=2 EIGJ0970 F=.25 EIGJ0980 AMAX=F*AT EIGJ0990 MAX=6 EIGJ1000 GO TO 5 EIGJ1010 31 MAX=MAX-1 EIGJ1020 IF(MAX)34,33,32 EIGJ1030 32 F=2.0*F*F EIGJ1040 AMAX=AT*F EIGJ1050 33 C=0.0 EIGJ1060 GO TO 5 EIGJ1070 34 IF(C)33,35,33 EIGJ1080 * COPY EIGENVALUES INTO E EIGJ1090 35 MM=1 EIGJ1100 NM=N&1 EIGJ1110 DO 36 M=1,N EIGJ1120 E(M)=A(MM) EIGJ1130 36 MM=MM&NM-M EIGJ1140 RETURN EIGJ1150 END EIGJ11g0 $ FORTRAN $ INCODE IBMF SUBROUTINE MTEVV(X,I1,I2,D,I3,I4,N,M) DIMENSION A(5050),X(I1,I2),D(I3,I4),E(100),NS(100),MS(100) K=1 DO1J=1,N DO1I=J,N A(K)=X(I,J) 1 K=K&1 CALL EIGENJ(A,D,N,M,E,I3,I4) C THE FOLLOWING SYSTEM SORT PACKAGE ROUTINES SORT THE FLOATING POINT C ARRAY E OF N ELEMENTS IN DESCENDING ORDER. CALL SORTFL(-1) IF(M)4,3,4 3 CALL SORT(E(1),E(2),N) DO 50 I=1,N 50 X(I,I)=E(I) RETURN 4 DO 5 I=1,N MS(I)=I 5 NS(I)=I C THE FOLLOWING ROUTINE SORTS THE FLOATING POINT ARRAY E OF N C ELEMENTS IN DESCENDING ORDER AND ARRAY MS IS THE PASSIVE LIST C WHICH WILL BE REARRANGED IN THE SAME ORDER AS E. CALL SORT(E(1),E(2),N,MS(1)) C THE FOLLOWING ROUTINES SORT THE FIXED POINT ARRAY MS OF N ELEMENTS C IN ASCENDING ORDER WITH ARRAY NS AS THE PASSIVE LIST. CALL SORTFX(1) CALL SORT(MS(1),MS(2),N,NS(1)) CALL ARREV(N,NS,D,I3,I4) DO 6 I=1,N 6 X(I,I)=E(I) RETURN END $ FORTRAN $ INCODE IBMF SUBROUTINE ARREV(N,NS,D,I3,I4) DIMENSION NS(N),D(I3,I4) DO 10 I=1,N 1 IF(NS(I)-I)2,10,2 2 CALL XCH(N,I,NS(I),D,NS,I3,I4) GOTO1 10 CONTINUE RETURN END $ FORTRAN $ INCODE IBMF SUBROUTINE XCH(N,L,M,D,NS,I3,I4) DIMENSION NS(N),D(I3,I4) I=L J=M DO 1 K=1,N T=D(K,I) D(K,I)=D(K,J) 1 D(K,J)=T NS(I)=NS(J) NS(J)=J RETURN END $ FORTRAN $ INCODE IBMF CMTINV MTINV SUBROUTINE MTINV (A,I1,I2,B,I3,I4,N,M) MTPK1115 DIMENSION A(I1,I2), B(I3,I4) MTPK1116 C MTPK1117 C THIS SUBROUTINE CALCULATES THE INVERSE OF MATRIX A BY THE GAUSS- MTPK1118 C JORDAN ELIMINATION SCHEME. ALL CALCULATIONS ARE DONE IN DOUBLE- MTPK1119 C PRECISION ARITHMETIC MTPK1120 C MTPK1121 DOUBLE PRECISION B,BMAX,BMULT MTPK1122 EQUIVALENCE (BMAX,BMULT) MTPK1123 C MTPK1124 C SINGLE TO DOUBLE PRECISION TRANSFER MTPK1125 C MTPK1126 101 DO 103 I = 1,N MTPK1127 102 DO 103 J = 1,N MTPK1128 103 B(I,J) = A(I,J) MTPK1129 C MTPK1130 C FIND THE LARGEST ELEMENT IN THE N-K BY N-K LOWER RIGHT SUBMATRIX. MTPK1133 C MTPK1134 200 DO 430 K = 1,N MTPK1131 A(K,1) = FLOAT(K) MTPK1135 A(K,2) = A(K,1) 203 BMAX = DABS(B(K,K)) MTPK1137 210 DO 219 I = K,N MTPK1138 211 DO 219 J = K,N MTPK1139 IF (BMAX - DABS(B(I,J)) ) 213,219,219 213 BMAX = DABS ( B(I,J)) MTPK1141 214 A(K,1) = FLOAT(I) MTPK1142 215 A(K,2) = FLOAT(J) MTPK1143 219 CONTINUE MTPK1145 216 IF (BMAX - 1.D-38) 801,301,301 C MTPK1146 C EXCHANGE ROWS AND COLUMNS TO PUT B(I,J) ON DIAGONAL. MTPK1147 C MTPK1148 301 I = IFIX(A(K,1)) MTPK1149 302 J = IFIX(A(K,2)) MTPK1150 313 DO 314 M = 1,N MTPK1151 BMULT = B(I,M) MTPK1152 B(I,M) = B(K,M) MTPK1153 314 B(K,M) = BMULT MTPK1155 303 DO 304 M = 1,N MTPK1156 BMULT = B(M,J) MTPK1157 B(M,J) = B(M,K) MTPK1158 304 B(M,K) = BMULT MTPK1159 C MTPK1160 C ACTUAL INVERSION STEP. MTPK1161 C BEGIN ROW ITERATION MTPK1162 C MTPK1163 401 DO 425 I = 1,N MTPK1164 C IGNORE ROW K. MTPK1165 402 IF (I - K) 405,425,405 C SET ROW MULTIPLIER. MTPK1167 405 BMULT = B(I,K) / B(K,K) C MODIFY ROW ELEMENTS. MTPK1169 410 DO 420 J = 1,N MTPK1170 C IGNORE COLUMN K MTPK1171 411 IF (J - K) 414,412,414 412 B(I,J) = -BMULT MTPK1173 413 GO TO 420 MTPK1174 414 B(I,J) = B(I,J) - B(K,J) * BMULT MTPK1175 420 CONTINUE MTPK1176 425 CONTINUE MTPK1177 C MTPK1178 C DIVIDE PIVOT ROW BY PIVOT ELEMENT MTPK1179 C MTPK1180 426 BMULT = 1.D0/ B(K,K) 427 B(K,K) = 1.D0 C ACTUAL DIVISION MTPK1183 428 DO 429 J = 1,N MTPK1184 429 B(K,J) = B(K,J) * BMULT MTPK1185 430 CONTINUE MTPK1186 C MTPK1187 C REARRANGE AND REASSEMBLE INVERSE MATRIX. MTPK1188 C MTPK1189 501 DO 512 IJK = 1,N MTPK1190 502 I = N-IJK&1 MTPK1191 503 L = IFIX(A(I,2)) MTPK1192 504 J = IFIX(A(I,1)) MTPK1193 505 DO 508 M = 1,N MTPK1194 506 BMULT = B(I,M) MTPK1195 507 B(I,M) = B(L,M) MTPK1196 508 B(L,M) = BMULT MTPK1197 509 DO 512 M = 1,N MTPK1198 510 BMULT = B(M,I) MTPK1199 511 B(M,I) = B(M,J) MTPK1200 512 B(M,J) = BMULT MTPK1201 C MTPK1202 C MOVE A INVERSE BACK INTO A MTPK1203 C MTPK1204 601 DO 603 I = 1,N MTPK1205 602 DO 603 J = 1,N MTPK1206 603 A(I,J) = B(I,J) 701 M = 1 MTPK1208 702 RETURN MTPK1209 C MTPK1210 C ERROR RETURN AT SOME TIME BMAX WAS .LT. 0.1D-38 C MTPK1212 801 M = 2 MTPK1213 802 RETURN MTPK1214 999 END MTPK1215