C*********************************************************************** C*** * C*** MAIN PROGRAM CANALS * C*** CANONICAL CORRELATION ANALYSIS * C*** MADE BY EEKE VAN DER BURG AND JAN DE LEEUW * C*** FOR DISCRETE NOMINAL, ORDINAL AND NUMERICAL DATA * C*** FIRST VERSION JULY 1978 * C*** SECOND VERSION JULY 1980 * C*** * C*********************************************************************** C*** * C*** LIST OF ARRAYS * C*** XY RAW DATA * C*** X DATA FIRST SET * C*** Y DATA SECOND SET * C*** XA PRODUCT OF X AND A * C*** YB PRODUCT OF Y AND B * C*** AB WEIGHTS * C*** A WEIGHTS OF THE FIRST SET * C*** B WEIGHTS OF THE SECOND SET * C*** IN INDICES OF RAW DATA VECTOR * C*** LB TIE BLOCK LENGTHS OF DATA VECTOR * C*** NK NUMBER OF CATEGORIES OR UPPER BOUND FOR NON MISSING VALUES * C*** IT VARIABLE TYPE * C*** NM VECTOR FOR PLOTS OF OPTIMAL SCORES * C*** H1,H2,H3,H4,H5,H6,H7,H8,H9 AUXILIAIRY ARRAYS * C*** XY AND X AND Y SHARE THE SAME STORAGE,XA AND H5,YB AND H6 IDEM. * C*** AB AND A AND B SHARE THE SAME STORAGE * C*** * C*** LIST OF PARAMETERS * C*** N NUMBER OF OBJECTS * C*** M NUMBER OF VARIABLES * C*** M1 NUMBER OF VARIABLES OF THE FIRST SET * C*** M2 NUMBER OF VARIABLES OF THE SECOND SET * C*** MC MAXIMUM NUMBER OF CATEGORIES FOR ONE VARIABLE * C*** P NUMBER OF DIMENSIONS * C*** NI MAXIMUM NUMBER OF ITERATIONS * C*** ES CRITERIUM STRESS DIFFERENCE * C*** TH CRITERIUM THETA FOR APPROXIMATION WEIGHTS * C*** II LOGICAL UNIT NUMBER FOR INPUT * C*** IO PRINT OUTPUT RAW DATA (1=YES,0=NO) * C*** IU UNIT NUMBER FOR OUTPUT OF OPTIMALLY SCALED VARIABLES * C*** (DATA VECTORS) (0,1=NO OUTPUT) * C*** IX PARAMETER FOR OUTPUT OF EVERY STRESS (1=YES,0=NO) * C*** ID UNIT NUMBER FOR OUTPUT OPTIMAL SCORES,WEIGHTS & CORRELATIONS * C*** (0,1=NO OUTPUT) * C*** IZ PARAMETER FOR OUTPUT CANONICAL SCORES (0=NO,1=PRINT,GT 1=UNIT)* C*** IY PARAMETER FOR PLOTS OF OPTIMAL SCORES (0=N0,1=YES) * C*** LZ PARAMETER FOR PLOT CANONICAL SCORES (0=N0 ,1=YES) * C*** LR PARAMETER FOR OUTPUT PROJECTIONS(0=N0,1=PRINT,GT 1=UNIT) * C*** IK PARAMETER FOR PRINT OF MISSING SCORES (0=N0,1=YES) * C*** IS LOGICAL UNIT FOR STORAGE OF RAW DATA INFORMATION * C*** IC LOGOCAL UNIT FOR CARD READING * C*** IP LOGICAL UNIT FOR PRINTING * C*** * C*** * C*********************************************************************** C*** * C*** THE CANALS PROGRAM SEARCHES FOR * C*** OPTIMALLY SCALED VARIABLES X AND Y AND WEIGHTS A AND B * C*** IN SUCH A WAY THAT THE SUM OF THE SQUARED CORRELATIONS * C*** (WHICH ARE THE CANONICAL CORRELATIONS) * C*** BETWEEN THE WEIGHTED OPTIMALLY SCALED VARIABLES ARE MAXIMIZED * C*** WHILE THE OPTIMALLY SCALED VARIABLES SATISFY THE MEASUREMENT * C*** RESTRICTIONS AND ARE STANDARDIZED AND THE WEIGHTED OPTIMALLY * C*** SCALED VARIABLES OR CANONICAL SCORES ARE STANDARDIZED, SO: * C*** FIND X,Y AND A,B FROM TRACE(XA-YB)T*(XA-YB) IS MINIMAL (1) * C*** WHILE X SATISFIES THE MEASUREMENT RESTRICTIONS AND XT*X=N*I (2) * C*** AND Y SATISFIES THE MEASUREMENT RESTRICTIONS AND YT*Y=N*I (3) * C*** AND EITHER (XA)T*XA=N*I (4) OR (YB)T*YB=N*I (5) * C*** ALTERNATINGLY (1),(2),(3),(5) AND (1),(2),(3),(4) IS SOLVED * C*** OR (1) AND (2) IS SOLVED WHILE (3) AND (5) ARE SATISFIED * C*** AND (1) AND (3) IS SOLVED WHILE (2) AND (4) ARE SATISFIED * C*** TO GET THE PROPER STANDARDIZATION OF XA OR YB * C*** GRAM-SCHMIDT ORTHOGONALIZATION IS USED * C*** CANALS FOLLOWS AN ALTERNATING LEAST SQUARES PRINCIPLE BECAUSE * C*** X (OR Y) AND A (OR B) ARE BOTH LEAST SQUARES SOLUTIONS * C*** AND X (OR Y) AND A (OR B) ARE SOLVED ALTERNATINGLY. * C*** * C*********************************************************************** C*** * C*** MAIN PROGRAM * C*** ASSIGNMENT OF NUMBERS TO AUXILIAIRY DATA SET(IS),CARD READER(IC) * C*** AND PRINTER(IP); IS=1,IC=5,IP=6 * C*** COMPUTATION OF TOTAL STORAGE FROM INPUT PARAMETERS AND * C*** DYNAMIC ARRAY ALLOCATION IN SUBROUTINE DYNAM. C*** * C*********************************************************************** EXTERNAL BINDEX COMMON/PAR/IS,IC,IP INTEGER N,M1,M2,P,P2,MC,II,IO,IU,IX,ID,NI,IS,IC,IP,NJ,J,IQ,IR,M, * IH,IZ,IY,LZ,LR,IK REAL TX,A DOUBLE PRECISION ES,TH DIMENSION TX(20) C C AUXILIAIRY UNIT IS,CARD READER IC, PRINTER IP IS=1 IC=5 IP=6 WRITE(IP,7) WRITE(IP,1) READ(IC,2)NJ WRITE(IP,5)NJ DO 500 J=1,NJ READ(IC,6)TX READ(IC,2)N,M1,M2,MC,P,NI,ES,TH,II,IO,IU,IX,ID,IZ,IY,LZ,LR,IK C*** DEFAULT VALUES OF THE MAXIMUM NUMBER OF ITERATIONS, C*** STRESS DIFFERENCE AND WEIGHT INCREASE IF(NI.EQ.0)NI=20 IF(ES.EQ.0.0)ES=1.0D-3 IF(TH.EQ.0.0)TH=1.0D-4 P=MIN0(P,M1,M2) P2=P*2 WRITE(IP,3)J,TX,N,M1,M2,MC,P,NI,ES,TH WRITE(IP,8)II,IO,IU,IX,ID,IZ,IY,LZ,LR,IK M=M1+M2 IQ=2*N*M+2*N*P*2+4*M*P+N+MC+4*M+3*MC*2+P*2+2*P**2+P*2 IR=((IQ*4)/1024+1)+70 C*** 70 K IS THE TOTAL LENGTH OF THE PROGRAM, INCLUSIVE ERRORMODULES C*** (70K = 17920 HALF WORDS) IH=IQ+17920 WRITE(IP,4) IR,IH C CALL DYNAM(BINDEX,A,IQ,N,M,M1,M2,MC,P,P2,NI,ES,TH,II,IO,IU,IX,ID, 1 IZ,IY,LZ,LR,IK,IS,IC,IP) IF(NJ.EQ.1) GOTO 490 IF(II.NE.IC) REWIND II REWIND IS 490 CONTINUE 500 CONTINUE C*** 1 FORMAT(1H1,///13H C A N A L S/1H ,103X,19H EEKE VAN DER BURG ,/ * 13H VERSION - 3,30X, * 37H CANONICAL CORRELATION ANALYSIS ,33X,1H\/1H ,106X, * 12HJAN DE LEEUW/13H AUGUST 1982, * 28X,47HFOR DISCRETE NOMINAL,ORDINAL AND NUMERICAL DATA/ * 1H ,103X,20HDEPT. OF DATA THEORY// * 1H ,101X,24HTHE UNIVERSITY OF LEIDEN////) 2 FORMAT(4I5/2I5,2F15.10/10I5) 3 FORMAT(16H PROBLEM NUMBER,I5//13H0 TITLE ,35X,20A4/// * 54H0 DATA SPECIFICATIONS: ,/// * 54H NUMBER OF OBJECTS OR INDIVIDUALS ,I5// * 54H NUMBER OF VARIABLES IN THE FIRST SET ,I5// * 54H NUMBER OF VARIABLES IN THE SECOND SET ,I5// * 54H HIGHEST NUMBER OF CATEGORIES (MAXIMUM OF CARD 6) ,I5//// * 54H0 ANALYSIS SPECIFICATIONS: ,/// * 54H NUMBER OF DIMENSIONS ,I5// * 54H MAXIMUM NUMBER OF ITERATIONS (DEFAULT=20) ,I5// * 51H MINIMUM STRESS DIFFERENCE (DEFAULT=0.001) ,F16.10// * 51H MAXIMUM WEIGHT INCREASE (DEFAULT=0.0001) ,F16.10) 8 FORMAT(1H1,/// * 54H I/O AND PLOT SPECIFICATIONS: ,/// * 54H UNIT NUMBER FOR THE INPUT MEDIUM OF THE DATA ,I5// * 54H PRINT OF RAW DATA (1=ALL OBJECTS, 0=10 OBJECTS) ,I5// * 54H UNIT NUMBER OF THE STORE MEDIUM FOR OPTIMALLY ,/ * 54H SCALED VARIABLES ,I5// * 54H PRINT OUTPUT OF THE STRESS (0=NO, 1=YES) ,I5// * 54H UNIT NUMBER OF STORE MEDIUM FOR CATEGORY ,/ * 54H QUANTIFICATIONS, WEIGHTS AND CORRELATIONS ,I5// * 54H OUTPUT CANONICAL SCORES ( 0=NO OUTPUT, 1=PRINT, ,/ * 54H GT 1=UNIT NUMBER OF STORE MEDIUM) ,I5// * 54H PLOT CATEGORY QUANTIFICATIONS ,/ * 54H (0=NO PLOTS,1=ONE OR MORE PLOTS) ,I5// * 54H PLOT CANONICAL SCORES (0=NO PLOT, 1=PLOT) ,I5// * 54H OUTPUT PROJECTIONS OF THE CATEGORY QUANTIFICATIONS ,/ * 54H (0=NO OUTPUT,1=PRINT,GT 1=UNIT NUMBER OF STORE ,/ * 54H MEDIUM) ,I5// * 54H OUTPUT QUANTIFICATIONS OF MISSING DATA ,/ * 54H (0=NO OUTPUT, 1=PRINT) ,I5///) 4 FORMAT(1H0,43H REQUIRED SPACE FOR ARRAYS AND PROGRAM IS :,I5, *10HK BYTES (,I8,13H HALF WORDS)) 5 FORMAT(22H0 NUMBER OF PROBLEMS ,I5//) 6 FORMAT(20A4) 7 FORMAT(1H1,/////////////////// +30X,58H +++++++++++++++++++++++++++++++++++++++++++++++++++++++++/ +30X,58H + +/ +30X,58H + AUGUST 1982 C A N A L S VERSION 3.0 +/ +30X,58H + +/ +30X,58H + +/ +30X,58H + +/ +30X,58H + EEKE VAN DER BURG TEL 148333 EXT 2257 +/ +30X,58H + +/ +30X,58H +++++++++++++++++++++++++++++++++++++++++++++++++++++++++) STOP END SUBROUTINE BINDEX(A,IQ,N,M,M1,M2,MC,P,P2,NI,ES,TH,II,IO,IU,IX,ID, 1 IZ,IY,LZ,LR,IK,IS,IC,IP) C*************************************** C*** * C*** COMPUTATION OF THE INDEX NUMBERS * C*** TO DEVIDE THE ALLOCATED STORAGE * C*** OF ARRAY A INTO SMALLER ARRAYS. * C*** SUBROUTINES CALLED: CANALS * C*** * C*************************************** INTEGER I,IQ,N,M,M1,M2,MC,P,P2,NI,II,IO,IU,IX,ID,IS,IC,IP,IZ,IY, +LZ,LR,IK REAL A DIMENSION A(IQ) DOUBLE PRECISION ES,TH C I1=1+N*M1*2 I2=I1+N*M2*2 I3=I2+N*P*2 I4=I3+N*P*2 I5=I4+M1*P*2 I6=I5+M2*P*2 I7=I6+N I8=I7+MC I9=I8+M I10=I9+M I11=I10+M I12=I11+MC*2 I13=I12+MC*2 I14=I13+MC*2 I15=I14+P*2 I16=I15+2*P**2 I17=I16+2*P I18=I17+M C*** I19=I18+M*P2 C CALL CANALS(A(1),A(1),A(I1),A(I2),A(I3), 1 A(I4),A(I5),A(I6),A(I7),A(I8), 2 A(I9),A(I10),A(I11),A(I12),A(I13), 3 A(I14),A(I2),A(I3),A(I15),A(I16), 4 A(I4),A(I17),A(I18),N,M,M1,M2,MC,P,P2,NI, 5 ES,TH,II,IO,IU,IX,ID,IZ,IY,LZ,LR,IK,IS,IC,IP) RETURN END C*********************************************************************** C*** * C*** ALGORITHM FLOW OF CANALS SUBROUTINE C*** READ X AND Y & STANDARDIZE ON MEAN=0 AND SUM OF SQUARES=N INDATA* C*** INITIALIZE WEIGHTS A WITH RANDOM VALUES RANDMA* C*** COMPUTE X*A PRODCT* C*** RESCALE XA SUCH THAT (XA)T*XA=N*I GRAMS1* C*** INITIALIZE WEIGHTS B WITH RANDOM NUMBERS RANDMA* C*** COMPUTE Y*B PRODCT* C*** START LOOP NUMERIC SOLUTION * C*** COMPUTE WEIGHTS B AND YB ITERWE* C*** COMPUTE STRESS((XA-YB))**2/N*P STRESS* C*** IF STRESS DIFFERENCE LT EPS GOTO ITERATION CYCLE C*** RESCALE YB AND XA SUCH THAT (YB)T*YB=NI GRAMS2* C*** COMPUTE WEIGHTS A AND XA ITERWE* C*** RESCALE XA AND YB SUCH THAT (XA)T*XA=N*I GRAMS2* C*** END LOOP NUMERIC SOLUTION * C*** START ITERATION CYCLE * C*** RESCALE YB SUCH THAT (YB)T*YB=N*I GRAMS2* C*** COMPUTE STRESS STRESS* C*** MODIFY X ,UPDATE A COLUMN OF X MODIFY* C*** REGRESS A COLUMN OF X REGRSX* C*** STANDARDIZE A COLUMN OF X SQDEV * C*** IF STRESS DIFFERENCE SMALL ENOUGH OUTPUT OF * C*** CATEGORY NUMBERS AND CATEGORY QUANTIFICATIONS OUTDAT* C*** COMPUTE X*A PRODCT* C*** COMPUTE WEIGHTS A AND XA ITERWE* C*** RESCALE XA SUCH THAT (XA)T*XA=N*I GRAMS2* C*** COMPUTE STRESS STRESS* C*** MODIFY Y ,UPDATE A COLUMN OF Y MODIFY* C*** REGRESS A COLUMN OF Y REGRSX* C*** STANDARDIZE A COLUMN OF Y SQDEV * C*** IF STRESS DIFFERENCE SMALL ENOUGH OUTPUT OF * C*** CATEGORY NUMBERS AND CATEGORY QUANTIFICATIONS OUTDAT* C*** COMPUTE Y*B PRODCT* C*** COMPUTE WEIGHTS B AND YB ITERWE* C*** IF STRESS DIFFERENCE SMALL ENOUGH GOTO OUT * C*** END ITERATION CYCLE,GOTO START * C*** OUT: PLOT THE CATEGORY QUANTIFICATIONS VERSUS THE CATEGORY * C*** NUMBERS PLOTOP* C*** COMPUTE CANONICAL CORRELATIONS AND ROTATE A AND B CANOCO* C*** SO THAT (XA)T*XA=NP*I AND (YB)Y*YB=N*I * C*** OUTPUT OF WEIGHTS A AND B OUTWEI* C*** COMPUTE X*A PRODCT* C*** COMPUTE Y*B PRODCT* C*** COMPUTE THE CORRELATIONS (X,XA),(X,YB),(Y,XA),(Y,YB) CORRAX* C*** PRINT THE PROJECTIONS OF THE CATEGORY QUANTIFICATIONS * C*** ON THE VARIATES SINGCO* C*** PRINT AND PLOT THE CANONICAL SCORES PLOTCA* C*** * C*********************************************************************** SUBROUTINE CANALS(XY,X,Y,XA,YB, 1 A,B,IN,LB,NM, 2 NK,IT,H1,H2,H3, 3 H4,H5,H6,H7,H8, 4 AB,H9,H10,N,M,M1,M2,MC,P,P2,NI, 5 ES,TH,II,IO,IU,IX,ID,IZ,IY,LZ,LR,IK,IS,IC,IP) LOGICAL LC,LX DOUBLE PRECISION FL,S1,S2,ES,TH INTEGER I,J,L,K,LB,IN,NK,NM,IT,N,P,P2,M,M1,M2,MC,NI,II,IO,IU,IX,ID 1 ,IS,IC,IP,IW,IZ,IY,H9,LZ,LR,IK REAL FM,H10 DOUBLE PRECISION XY,X,Y,XA,YB,A,B,H1,H2,H3,H4,H5,H6,H7,H8,AB DIMENSION XY(N,M),X(N,M1),Y(N,M2),XA(N,P),YB(N,P), 1 A(M1,P),B(M2,P),IN(N),LB(MC),NM(M), 2 NK(M),IT(M),H1(MC),H2(MC),H3(MC), 3 H4(P),H5(N),H6(N),H7(P,P),H8(P), 4 AB(M,P),H9(M),H10(M,P2),FM(60) READ(IC,2) NK READ(IC,2)IT CALL CHECK1(NK,IT,M,MC,IP) DO 1 J=1,M 1 NM(J)=0 IF(IY.EQ.1)READ(IC,2)NM READ(IC,4) FM 2 FORMAT(16I5) 4 FORMAT(20A4) 5 FORMAT(8H0 FORMAT//3(2H ,20A4/)) 6 FORMAT(75H1 ITERATION STRESS AFTER CHANGING NO * ITERATIONS WEIGHTS/) 7 FORMAT(1H ,2X,I5,4X,F20.10,5X,9HWEIGHTS 2,16X,I5) 8 FORMAT(1H ,11X,F20.10,5X,9HWEIGHTS 1,16X,I5) 9 FORMAT(//1H0,15HVARIABLE NUMBER,5X,16HNO OF CATEGORIES,5X, *17HMEASUREMENT LEVEL,5X,25HNUMBER OF MISSING ENTRIES, *5X,29HPLOT CATEGORY QUANTIFICATIONS/) 10 FORMAT(1H ,3X,I5,13X,I5,13X,I5,12H = NOMINAL ,15X,I5,25X,I1) 11 FORMAT(1H ,3X,I5,13X,I5,13X,I5,12H = ORDINAL ,15X,I5,25X,I1) 12 FORMAT(1H ,3X,I5,13X,I5,13X,I5,12H = NUMERICAL,15X,I5,25X,I1) 14 FORMAT(37H0MAXIMUM NUMBER OF ITERATIONS REACHED) 15 FORMAT(1H1,10HFIRST SET://) 16 FORMAT(1H1,11HSECOND SET://) 17 FORMAT(76H0THE FIRST STRESS IS THE STRESS OF A COMPLETE NUMERICAL *SOLUTION OF THE DATA/ *51H0(WITH A FIXED MISSING CATEGORY FOR EVERY VARIABLE)) C*** C*** READ XY, STANDARDIZE AND WRITE NECESSARY INFO FOR SCALING ON UNIT CALL INDATA(XY,H5,IN,LB,H9,NK,N,M,MC,II,IO,IS,IP,FM) WRITE(IP,9) DO 100 J=1,M K=N-H9(J) IF(IT(J).EQ.0)WRITE(IP,10)J,NK(J),IT(J),K,NM(J) IF(IT(J).EQ.1)WRITE(IP,11)J,NK(J),IT(J),K,NM(J) IF(IT(J).EQ.2)WRITE(IP,12)J,NK(J),IT(J),K,NM(J) 100 CONTINUE WRITE(IP,5)FM C*** INITIALIZE MATRIX A,B , COMPUTE XA,YB AND RESCALE XA L=M1*P CALL RANDMA(A,L) CALL PRODCT(X,A,XA,N,M1,P) CALL GRAMS1(XA,A,N,P,M1) L=M2*P CALL RANDMA(B,L) CALL PRODCT(Y,B,YB,N,M2,P) S1=100.0 C*** COMPUTATION OF THE WEIGHTS OF A COMPLETE NUMERICAL SOLUTION 110 CALL ITERWE(XA,YB,Y,B,N,P,M2,TH,IP,IW) CALL STRESS(XA,YB,N,P,S2) C*** WRITE(IP,8)S2,IW IF(DABS(S2-S1).LT.ES) GOTO 120 CALL GRAMS2(YB,B,XA,A,N,P,M2,M1) CALL ITERWE(YB,XA,X,A,N,P,M1,TH,IP,IW) CALL GRAMS2(XA,A,YB,B,N,P,M1,M2) S1=S2 GOTO 110 120 S1=1.5D0 FL=1.0D0/DFLOAT(N) LX=.FALSE. C*** ITERATION CYCLE WRITE(IP,6) DO 400 K=1,NI C*** RESCALE YB AND COMPUTE STRESS CALL GRAMS2(YB,B,XA,A,N,P,M2,M1) CALL STRESS(XA,YB,N,P,S2) IF(DABS(S1-S2).LT.ES) LX=.TRUE. IF(K.EQ.NI)LX=.TRUE. IF(IX.EQ.1)WRITE(IP,7)K,S2,IW IF(IX.EQ.0.AND.(K.EQ.1.OR.LX))WRITE(IP,7)K,S2,IW IF(K.EQ.NI)WRITE(IP,14) IF(LX) WRITE(IP,17) IF(LX)WRITE(IP,15) C*** MODIFY X DO 200 L=1,M1 C*** COMPUTE NEW UPDATE FOR COLUMN OF X CALL MODIFY(X(1,L),XA,YB,A,H5,L,N,P,M1) DO 190 I=1,P 190 H4(I)=A(L,I) C*** REGRESSION ON COLUMN OF X AND STANDARDIZATION CALL REGRSX(IN,LB,X(1,L),H1,H2,H3,H5,H4,N,P,NK(L),IT(L),IS) DO 195 I=1,P 195 A(L,I)=H4(I) CALL SQDEV(X(1,L),N,FL,N,LC) IF(LX)CALL OUTDAT(IN,LB,X(1,L),NK(L),IT(L),N,L,IP,ID,IK) CALL PRODCT(X,A,XA,N,M1,P) 200 CONTINUE C*** COMPUTE WEIGHTS A CALL ITERWE(YB,XA,X,A,N,P,M1,TH,IP,IW) C*** RESCALE XA AND COMPUTE STRESS CALL GRAMS2(XA,A,YB,B,N,P,M1,M2) CALL STRESS(XA,YB,N,P,S1) IF(IX.EQ.1.AND..NOT.LX)WRITE(IP,8)S1,IW IF(LX)WRITE(IP,16) C*** MODIFY Y DO 300 L=1,M2 J=L+M1 C*** COMPUTE NEW UPDATE FOR A COLUMN OF Y CALL MODIFY(Y(1,L),YB,XA,B,H6,L,N,P,M2) DO 290 I=1,P 290 H4(I)=B(L,I) C*** REGRESSION ON COLUMN OF Y AND STANDARDIZATION CALL REGRSX(IN,LB,Y(1,L),H1,H2,H3,H6,H4,N,P,NK(J),IT(J),IS) DO 295 I=1,P 295 B(L,I)=H4(I) CALL SQDEV(Y(1,L),N,FL,N,LC) IF(LX)CALL OUTDAT(IN,LB,Y(1,L),NK(J),IT(J),N,J,IP,ID,IK) CALL PRODCT(Y,B,YB,N,M2,P) 300 CONTINUE C*** COMPUTE WEIGHTS B CALL ITERWE(XA,YB,Y,B,N,P,M2,TH,IP,IW) REWIND IS IF(LX)GOTO 500 400 CONTINUE 500 CALL OUTVAR(XY,N,M,IU,IP) DO 550 J=1,M 550 CALL PLOTOP(XY(1,J),IN,LB,H1,H2,H3,IS,IP,N,NK(J),IT(J),NM(J),J) CALL CANOCO(YB,B,A,H7,H4,H8,N,P,M2,M1,IP) CALL OUTWEI(A,B,M1,M2,P,IP,ID) CALL PRODCT(X,A,XA,N,M1,P) CALL PRODCT(Y,B,YB,N,M2,P) CALL CORRAX(XY,XA,YB,IN,LB,NK,H4,H7,H9,H10,N,M,M1,MC,P, 1 IP,ID,P2,LR,IS) CALL PLOTCA(XA,YB,XY(1,1),N,P,IZ,IP,LZ) CALL OUTPAR(ID,IU,IZ,IP,LR) RETURN END SUBROUTINE CANOCO(YB,B,A,H7,H4,H8,N,P,M2,M1,IP) C*********************************************************** C*** * C*** COMPUTATION OF THE EIGENVECTORS Q AND * C*** THE EIGENVALUES L**2 OF MATRIX(YB)T*YB/N * C*** THE WEIGHTS A AND B ARE ROTATED : * C*** A=A*Q AND B=B*Q*L**-1 * C*** SO THAT XA AND YB ARE BOTH STANDARDIZED ON * C*** SUM OF SQUARES=N * C*** L ARE THE CANONICAL CORRELATIONS * C*** SUBROUTINES CALLED: TRED2,IMTQL2 * C*** * C*********************************************************** C*** INTEGER K,L,I,IR,IP,N,P,M2,M1 DOUBLE PRECISION YB,B,A,H4,H7,H8 DOUBLE PRECISION DFLOAT,DSQRT DIMENSION YB(N,P),B(M2,P),A(M1,P),H4(P),H7(P,P),H8(P) C*** 1 FORMAT(36H1EIGENVALUE COMPUTATION NOT OKAY,IR=,I5) 2 FORMAT(42H1CANONICAL CORRELATIONS FOR EACH DIMENSION/ *1H0,9X,15F8.3) DO 10 K=1,P DO 10 L=K,P H7(K,L)=0.0D0 DO 5 I=1,N 5 H7(K,L)=H7(K,L)+YB(I,K)*YB(I,L) H7(K,L)=H7(K,L)/DFLOAT(N) 10 H7(L,K)=H7(K,L) C*** IF(P.NE.1)GOTO 11 H4(1)=H7(1,1) H7(1,1)=1.0D0 GOTO 12 11 CALL TRED2(P,P,H4,H8,H7) CALL IMTQL2(P,P,H4,H8,H7,IR) IF(IR.GT.0)WRITE(IP,1)IR IF(IR.GT.0)STOP 12 DO 20 J=1,M1 DO 16 K=1,P H8(K)=0.0D0 DO 16 L=1,P 16 H8(K)=H8(K)+A(J,L)*H7(L,K) DO 17 K=1,P 17 A(J,K)=H8(K) 20 CONTINUE DO 35 K=1,P 35 H4(K)=DSQRT(H4(K)) DO 40 J=1,M2 DO 36 K=1,P H8(K)=0.0D0 DO 36 L=1,P 36 H8(K)=H8(K)+B(J,L)*H7(L,K) DO 37 K=1,P 37 B(J,K)=H8(K)/H4(K) 40 CONTINUE WRITE(IP,2) H4 RETURN END SUBROUTINE CATLGT(NX,LB,N,NK,LV,IP) C**************************************** C*** * C*** COMPUTATION OF TIE BLOCK LENGTHS * C*** * C**************************************** INTEGER LB,I,K,NX,LV,IP DIMENSION NX(N),LB(NK) DO 100 K=1,NK 100 LB(K)=0 K=1 DO 300 I=1,N 150 IF(NX(I).GT.K)GOTO 200 LB(K)=LB(K)+1 GOTO 300 200 IF(K.EQ.NK) GOTO 400 K=K+1 GOTO 150 300 CONTINUE 400 WRITE(IP,401)LV,(LB(K),K=1,NK) 401 FORMAT(1H ,5X,I5,10X,15I6/(1H ,20X,15I6)) RETURN END SUBROUTINE CHECK1(NK,IT,M,MC,IP) C**************************************** C*** * C*** CHECK ON THE MAXIMUM NUMBER OF * C*** CATEGORIES AND THE VARIABLE TYPES * C*** * C**************************************** INTEGER M,MC,IP,NK,IT DIMENSION NK(M),IT(M) DO 10 J=1,M IF(NK(J).LE.MC) GOTO 10 WRITE(IP,1) 1 FORMAT(/////47H0 THE MAXIMUM NUMBER OF CATEGORIES IS INCORRECT/ * 47H *********************************************) STOP 10 CONTINUE DO 20 J=1,M IF(IT(J).GE.0.AND.IT(J).LE.3) GOTO 20 WRITE(IP,11) 11 FORMAT(/////30H0 A VARIABLE TYPE IS INCORRECT/ + 30H ****************************) STOP 20 CONTINUE RETURN END SUBROUTINE CORRAX(XY,XA,YB,IN,LB,NK,H4,H7,H9,H10,N,M,M1,MC,P, 1 IP,ID,P2,LR,IS) C*********************************************************** C*** * C*** COMPUTATION OF THE CORRELATIONS BETWEEN THE * C*** OPTIMALLY SCALED VARIABLES AND THE CANONICAL VARIATES* C*** AND OF THE PROJECTIONS ON THE CANONICAL VARIATES * C*** * C*** SUBROUTINES CALLED: PLOT,SINGCO * C*** * C*********************************************************** C*** INTEGER N,M,M1,MC,P,IP,I,J,K,ID,P2,LR REAL H10,H9 DOUBLE PRECISION XY,XA,YB,H4,FN DIMENSION XY(N,M),XA(N,P),YB(N,P),H4(P),H9(M),H10(M,P2),IN(N), 1 LB(MC),NK(M),H7(P,P) C*** 1 FORMAT(69H1CORRELATIONS BETWEEN THE OPTIMALLY SCALED VARIABLES OF *THE FIRST SET/) 2 FORMAT(70H0CORRELATIONS BETWEEN THE OPTIMALLY SCALED VARIABLES OF *THE SECOND SET/) 3 FORMAT(1H ,62HAND THE CANONICAL VARIATES OF THE FIRST SET FOR EACH * DIMENSION//) 4 FORMAT(1H ,63HAND THE CANONICAL VARIATES OF THE SECOND SET FOR EAC *H DIMENSION//) 5 FORMAT(1H ,I5,4X,15F8.3) 6 FORMAT(I5,9F8.3/(5X,9F8.3)) 7 FORMAT(87H1 CORRELATIONS BETWEEN OPTIMALLY SCALED VARIABLES AND CA *NONICAL VARIATES (HORIZONTAL NO,I2,16H AND VERTICAL NO,I2, *19H) OF THE FIRST SET ) 8 FORMAT(87H1 CORRELATIONS BETWEEN OPTIMALLY SCALED VARIABLES AND CA *NONICAL VARIATES (HORIZONTAL NO,I2,16H AND VERTICAL NO,I2, *19H) OF THE SECOND SET) C*** FN=1.0/DFLOAT(N) WRITE(IP,1) WRITE(IP,3) DO 130 J=1,M DO 120 K=1,P H4(K)=0.0D0 DO 110 I=1,N 110 H4(K)=XY(I,J)*XA(I,K)+H4(K) 120 H10(J,K)=H4(K)*FN WRITE(IP,5)J,(H10(J,K),K=1,P) IF(ID.GT.1)WRITE(ID,6)J,(H10(J,K),K=1,P) IF(J.EQ.M1) WRITE(IP,2) IF(J.EQ.M1) WRITE(IP,3) 130 CONTINUE IF(P.EQ.1) GOTO 245 MP=P-1 DO 240 L=1,MP LP=L+1 DO 240 K=LP,P WRITE(IP,7)L,K 240 CALL PLOT(H10(1,L),H10(1,K),H9,M,M,1,IP) 245 CONTINUE WRITE(IP,1) WRITE(IP,4) DO 330 J=1,M DO 320 K=1,P H4(K)=0.0D0 DO 310 I=1,N 310 H4(K)=XY(I,J)*YB(I,K)+H4(K) 320 H10(J,K+P)=H4(K)*FN WRITE(IP,5)J,(H10(J,K+P),K=1,P) IF(ID.GT.1)WRITE(ID,6)J,(H10(J,K+P),K=1,P) IF(J.EQ.M1) WRITE(IP,2) IF(J.EQ.M1) WRITE(IP,4) 330 CONTINUE IF(P.EQ.1) GOTO 450 MP=P-1 DO 440 L=1,MP LP=L+1 DO 440 K=LP,P WRITE(IP,8)L,K 440 CALL PLOT(H10(1,L+P),H10(1,K+P),H9,M,M,1,IP) 450 CALL SINGCO(H7(1,1),H10,XY,IN,LB,NK,N,P,M,MC,ID,IP,P2,LR,IS) RETURN END SUBROUTINE DYNAM(BINDEX,B,IQ,N,M,M1,M2,MC,P,P2,NI,ES,TH,II,IO,IU, 1 IX,ID,IZ,IY,LZ,LR,IK,IS,IC,IP) C************************************************* C*** * C*** ALLOCATION OF FIXED STORAGE IN STEAD OF * C*** THE DYNAMIC ALLOCATION BY THE ASSEMBLER * C*** ROUTINE DYNAM * C*** * C************************************************* C*** INTEGER IQ,N,M,M1,M2,MC,P,P2,NI,II,IO,IU,IX,ID,IZ,IY,IS,IC,IP,LZ, +LR,IK DOUBLE PRECISION ES,TH REAL A,B C-----ALLOCATION OF STORAGE IF DYNAM IS NOT LINKED DIMENSION A(10000) C*** ATTENTION, IQ HAS TO BE COMPARED WITH THE SEIZE OF ARRAY A IF(IQ.LE.10000)GOTO 100 WRITE(IP,10) 10 FORMAT(36H0NOT ENOUGH STORAGE FOR THIS PROBLEM) STOP 100 CALL BINDEX(A,IQ,N,M,M1,M2,MC,P,P2,NI,ES,TH,II,IO,IU,IX,ID,IZ, + IY,LZ,LR,IK,IS,IC,IP) RETURN END SUBROUTINE GRAMS1(XA,A,N,P,M1) C******************************************** C*** * C*** GRAMS1(XA,A,N,P,M1) * C*** GRAM SCHMIDT ORTHOGONALIZATION OF XA * C*** COMPUTATION OF XA*S ANS A*S SUCH THAT * C*** (XA*S)T*(XA*S)=N*I * C*** * C******************************************** C*** DOUBLE PRECISION XA,A,S1,FN INTEGER I,J,K,L,N,P,M1,K1 DIMENSION XA(N,P),A(M1,P) C*** FN=1.0D0/DFLOAT(N) DO 100 K=1,P S1=0.0D0 DO 30 I=1,N 30 S1=S1+XA(I,K)**2 S1=1.0D0/DSQRT(S1*FN) DO 40 I=1,N 40 XA(I,K)=XA(I,K)*S1 DO 50 J=1,M1 50 A(J,K)=A(J,K)*S1 IF(K.EQ.P) GOTO 110 K1=K+1 DO 90 L=K1,P S1=0.0D0 DO 60 I=1,N 60 S1=S1+XA(I,K)*XA(I,L) S1=S1*FN DO 70 I=1,N 70 XA(I,L)=XA(I,L)-XA(I,K)*S1 DO 80 J=1,M1 80 A(J,L)=A(J,L)-A(J,K)*S1 90 CONTINUE 100 CONTINUE 110 CONTINUE RETURN END SUBROUTINE GRAMS2(XA,A,YB,B,N,P,M1,M2) C*********************************************************** C*** * C*** GRAMS2(XA,A,YB,B,N,P,M1,M2) * C*** GRAM SCHMIDT ORTHOGONALIZATION OF XA * C*** COMPUTATION OF XA*S,A*S,YB*ST**-1 AND B*ST**-1 * C*** SUCH THAT (XA*S)T*(XA*S)=N*I * C*** * C*********************************************************** C*** DOUBLE PRECISION XA,A,YB,B,S1,S2,FN INTEGER I,J,K,L,N,P,M1,M2,K1 DIMENSION XA(N,P),A(M1,P),YB(N,P),B(M2,P) C*** FN=1.0D0/DFLOAT(N) DO 100 K=1,P S1=0.0D0 DO 30 I=1,N 30 S1=S1+XA(I,K)**2 S1=DSQRT(S1*FN) S2=1.0D0/S1 DO 40 I=1,N XA(I,K)=XA(I,K)*S2 40 YB(I,K)=YB(I,K)*S1 DO 50 J=1,M1 50 A(J,K)=A(J,K)*S2 DO 55 J=1,M2 55 B(J,K)=B(J,K)*S1 IF(K.EQ.P) GOTO 110 K1=K+1 DO 90 L=K1,P S1=0.0D0 DO 60 I=1,N 60 S1=S1+XA(I,K)*XA(I,L) S1=S1*FN DO 70 I=1,N XA(I,L)=XA(I,L)-XA(I,K)*S1 70 YB(I,K)=YB(I,K)+YB(I,L)*S1 DO 80 J=1,M1 80 A(J,L)=A(J,L)-A(J,K)*S1 DO 85 J=1,M2 85 B(J,K)=B(J,K)+B(J,L)*S1 90 CONTINUE 100 CONTINUE 110 RETURN END SUBROUTINE IMTQL2(NM,N,D,E,Z,IERR) C ***************************************************************** C * * C * I M T Q L 2 * C * * C * PURPOSE: DETERMINES THE EIGENVALUES AND EIGENVECTORS OF A * C * SYMMETRIC TRIDIAGONAL MATRIX USING THE IMPLICIT * C * QL METHOD. * C * * C * SUBROUTINES CALLED: NONE * C * * C * ADAPTED FROM EISPACK-ORIGINAL VERSION NOT IN DOUBLE PRECISION* C * * C ***************************************************************** DIMENSION D(N),E(N),Z(NM,N) DOUBLE PRECISION D,E,Z,B,C,F,G,P,R,S,MACHEP DOUBLE PRECISION DSQRT,DABS,DSIGN C C ********** MACHEP IS A MACHINE DEPENDENT PARAMETER SPCIFYING C THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC. C C ********** C MACHEP = 2.0**(-20) C IERR = 0 IF (N .EQ. 1) GO TO 1001 C DO 100 I = 2, N 100 E(I-1) = E(I) C E(N) = 0.0 C DO 240 L = 1, N J = 0 C ********** LOOK FOR SMALL SUB-DIAGONAL ELEMENT ********** 105 DO 110 M = L, N IF (M .EQ. N) GO TO 120 IF (DABS(E(M)) .LE. MACHEP* (DABS(D(M)) + DABS(D(M+1)))) X GO TO 120 110 CONTINUE C 120 P = D(L) IF (M .EQ. L) GO TO 240 IF (J .EQ.30) GO TO 1000 J = J + 1 C ********** FORM SHIFT ********** G = (D(L+1) - P) / (2.0 * E(L)) R = DSQRT(G*G+1.0) G = D(M) - P + E(L) / (G + DSIGN(R,G)) S = 1.0 C = 1.0 P = 0.0 NML = M - L C ********** FOR I=M-1 STEP -1 UNTIL L DO -- ********** DO 200 II = 1, NML I = M - II F = S * E(I) B = C * E(I) IF (DABS(F) .LT. DABS(G)) GO TO 150 C = G / F R = DSQRT(C*C+1.0) E(I+1) = F * R S = 1.0 / R C = C * S GO TO 160 150 S = F / G R = DSQRT(S*S+1.0) E(I+1) = G* R C = 1.0 / R S = S * C 160 G = D(I+1) - P R = (D(I) - G) * S + 2.0 * C * B P = S * R D(I+1) = G + P G = C * R - B C ********** FORM VECTOR ********** DO 180 K = 1, N F = Z(K,I+1) Z(K,I+1) = S * Z(K,I) + C * F Z(K,I) = C * Z(K,I) - S * F 180 CONTINUE C 200 CONTINUE C D(L) = D(L) - P E(L) = G E(M) = 0.0 GO TO 105 240 CONTINUE C ********** ORDER EIGENVALUES AND EIGENVECTORS ********** DO 300 II = 2, N I = II - 1 K = I P = D(I) C DO 260 J = II, N IF (D(J) .LE. P ) GO TO 260 K = J P = D(J) 260 CONTINUE C IF (K .EQ. I) GO TO 300 D(K) = D(I) D(I) = P C DO 280 J = 1, N P = Z(J,I) Z(J,I) = Z(J,K) Z(J,K) = P 280 CONTINUE C 300 CONTINUE C GO TO 1001 C ********** SET ERROR -- NO CONVERGENCE TO AN C EIGENVALUE AFTER 30 ITERATIONS ********** 1000 IERR = L 1001 RETURN END SUBROUTINE INDATA(XY,NX,IN,LB,NM,NK,N,M,MC,II,IO,IS,IP,FM) C*********************************************************************** C*** * C*** INPUT OF RAW DATA FROM UNIT II * C*** THE DATA ARE READ ROW WISE (ONE RECORD FOR EVERY OBJECT) * C*** ACCORDING TO FORMAT FM(FORTRAN I-FORMAT) * C*** THE INDEX NUMBERS OF THE RAW DATA AND THE TIE BLOCK LENGHTS ARE * C*** STORED ON UNIT NUMBER IS. * C*** AT RETURN THE MATRIX XY CONTAINS THE STANDARDIZED RAW DATA * C*** (MEAN=0 AND SQUARED SUM=N). * C*** SUBROUTINES CALLED: SORT,CATLGT,SQDEV * C*** * C*********************************************************************** C*** INTEGER N,M,MC,II,IO,IS,IP,I,J,K,NM,NK,IN,LB,NX,NJ REAL FM LOGICAL LC DOUBLE PRECISION XY,AR,FL DIMENSION XY(N,M),NX(N),IN(N),LB(MC),NM(M),NK(M),FM(60) C*** 1 FORMAT(45H1RAW DATA , ROWS=OBJECTS , COLUMNS=VARIABLES ,16X,I6, *3X,19HOBJECTS ARE PRINTED/) 2 FORMAT(1H ,I6,3X,25I4/(1H ,9X,25I4)) 3 FORMAT(1H0,15HVARIABLE NUMBER,5X, *49HMARGINAL FREQUENCIES OF THE SUCCEEDING CATEGORIES/) IF(IO.EQ.0)IO=10 IF(IO.EQ.1.OR.N.LT.10)IO=N WRITE(IP,1)IO DO 10 I=1,IO READ(II,FM)NM WRITE(IP,2)I,NM DO 9 J=1,M IF(NM(J).LE.0.OR.NM(J).GT.NK(J))NM(J)=NK(J)+1 9 XY(I,J)=NM(J) 10 CONTINUE IF(IO.EQ.N) GOTO 21 DO 20 I=11,N READ(II,FM)NM DO 19 J=1,M IF(NM(J).LE.0.OR.NM(J).GT.NK(J)) NM(J)=NK(J)+1 19 XY(I,J)=NM(J) 20 CONTINUE 21 WRITE(IP,3) DO 55 J=1,M DO 25 I=1,N 25 NX(I)=XY(I,J) CALL SORT(NX,1,N,IN,N) NJ=NK(J) CALL CATLGT(NX,LB,N,NJ,J,IP) NM(J)=0 DO 30 K=1,NJ 30 NM(J)=NM(J)+LB(K) C*** NM(J) NUMBER OF NON MISSING ENTRIES FOR VARIABLE J WRITE(IS)IN,(LB(K),K=1,NJ) FL=1.0D0/DFLOAT(N) CALL SQDEV(XY(1,J),N,FL,N,LC) 55 CONTINUE REWIND IS RETURN END SUBROUTINE ITERWE(XA,YB,Y,B,N,P,M2,TH,IP,L) C************************************************************************ C*** * C*** COMPUTE MINIMUM TRACE(XA-YB)T(XA-YB) WHILE XAT*XA=N*I * C*** T MEANS TRANSPOSED MATRIX * C*** SOLUTION FOR B: * C*** B=(YT*Y)**-1*Y*XA * C*** TO AVOID COMPUTATION OF THE SSCP MATRIX B IS APPROXIMATED. * C*** B=B+THETA*E(J,K) * C*** E(J,K)=1 FOR ELEMENT (J,K) AND =0 FOR ALL OTHER ELEMENTS. * C*** REFORMULATION OF THE PROBLEM: * C*** COMPUTE MINIMUM TR(XA-YB-Y*THETA*E(J,K))T(XA-YB-Y*THETA*E(J,K)) * C*** WHILE XAT*XA=N*I * C*** BY DIFFERENTIATION WE FIND THE SMALLEST SQUARE SOLUTION FOR THETA* C*** THETA=(INNER PRODUCT OF COLUMN J OF Y AND COLUMN K OF (XA-YB))/N * C*** A MEASURE FOR GOODNESS OF THE APPROXIMATION IS THE MAXIMUM THETA * C*** OVER J AND K. THE MATRIX B IS APPROXIMATED ANOTHER TIME IF * C*** THIS MAXIMUM IS LESS THAN TH (CONVERGENCE CRITERIUM). * C*** * C************************************************************************ C*** DOUBLE PRECISION XA,YB,Y,B,TH,AT,TM,FN INTEGER N,P,M2,I,J,K,L,IP DIMENSION XA(N,P),YB(N,P),Y(N,M2),B(M2,P) C*** L=0 FN=1.0D0/DFLOAT(N) 50 TM=0.0D0 L=L+1 DO 100 J=1,M2 DO 100 K=1,P AT=0.0D0 DO 80 I=1,N 80 AT=AT+(XA(I,K)-YB(I,K))*Y(I,J) AT=AT*FN IF(DABS(AT).LT.1.0D-20) GOTO 100 TM=DMAX1(TM,DABS(AT)) B(J,K)=B(J,K)+AT DO 90 I=1,N 90 YB(I,K)=YB(I,K)+AT*Y(I,J) 100 CONTINUE IF(L.GE.20)GOTO 200 IF(TM.GT.TH) GOTO 50 200 RETURN END SUBROUTINE LIDICA(MARFRE,YBLANK,EPLUS,DPLUS,KJ) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C*** C C*** THIS ROUTINE COMPUTES A WEIGTHED LINEAR REGRESSION OF THE C C*** MODEL VECTOR YBLANK ON CATEGORICAL DISCRETE DATA C C*** C C*** NO SUBROUTINES CALLED J.VAN RIJCKEVORSEL MAY 78 C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DIMENSION MARFRE(KJ),YBLANK(KJ),EPLUS(KJ),DPLUS(KJ) DOUBLE PRECISION CONST,DELTA,GAMMA,HULP1,HULP2,ALPHA,BETA,DPLUS, +EPLUS,YBLANK INTEGER MARFRE C*** CONST=0. DELTA=0. GAMMA=0. ALPHA=0. BETA =0. C*** DO 1 K=1,KJ CONST=CONST+MARFRE(K) DELTA=DELTA+MARFRE(K)*K GAMMA=GAMMA+MARFRE(K)*YBLANK(K) 1 CONTINUE C*** IF (CONST.GT.1.E-20)CONST=1./CONST DELTA=DELTA*CONST GAMMA=GAMMA*CONST HULP1=0. HULP2=0. C*** DO 2 K=1,KJ EPLUS(K)=YBLANK(K)-GAMMA DPLUS(K)=K-DELTA HULP1=HULP1+EPLUS(K)*MARFRE(K)*DPLUS(K) HULP2=HULP2+DPLUS(K)*MARFRE(K)*DPLUS(K) 2 CONTINUE C*** IF (HULP2.GT.1.E-20) ALPHA=HULP1/HULP2 IF(ALPHA.LT.0.) ALPHA=-ALPHA BETA=GAMMA-ALPHA*DELTA C*** DO 3 K=1,KJ IF(MARFRE(K).EQ.0) GOTO 3 YBLANK(K)=ALPHA*K+BETA 3 CONTINUE RETURN END SUBROUTINE MODIFY(XL,XA,YB,A,H5,L,N,P,M1) C*********************************************************** C*** * C*** COMPUTAION OF AN UPDATE OF VECTOR XL * C*** XL=(YB-XA)*A(L)+A(L)T*A(L)*XL * C*** XL THE L'TH COLUMN OF MATRIX X(N,M1) * C*** A(L) THE L'TH ROW OF MATRIX A(M1,P) * C*** H5 CONTAINS THE PREVIOUS VALUES OF XL AT RETURN * C*** IF H5 AND XA START AT THE SAME STORAGE LOCATION * C*** XA IS DESTROYED. * C*** * C*********************************************************** C*** INTEGER L,N,P,K,I DOUBLE PRECISION XL,XA,YB,A,H5,S,HI DIMENSION XL(N),XA(N,P),YB(N,P),A(M1,P),H5(N) S=0.0D0 DO 10 K=1,P 10 S=S+A(L,K)**2 DO 30 I=1,N HI=0.0D0 DO 20 K=1,P 20 HI=HI+(YB(I,K)-XA(I,K))*A(L,K) H5(I)=XL(I) XL(I)=HI+S*XL(I) 30 CONTINUE RETURN END C**** SUBROUTINE MTR (DIST,W,N) C**** ****************************************************************** C**** * * C**** * SUBROUTINE MTR * C**** * * C**** * PURPOSE: * C**** * PERFORM A WEIGHTED MONOTONIC REGRESSION, * C**** * IN ASCENDING ORDER. * C**** * * C**** * USAGE: * C**** * CALL MTR (DIST,DISP,W,N) * C**** * * C**** * SUBPROGRAMS CALLED: * C**** * NONE. * C**** * * C**** * * C**** * PARAMETERS: * C**** * * C**** * NAME TYPE DESCRIPTION * C**** * * C**** * DIST DEFAULT INPUT ARRAY ON WHICH A MONOTONIC**** * C**** * REGRESSION IS TO BE PERFORMED; AND * C**** * OUTPUT ARRAY CONTAINING THE MONOTONIC * C**** * LEAST-SQUARES TRANSFORMATION OF DIST. * C**** * W DEFAULT INPUT ARRAY OF WEIGHTS OF CORRESPONDING * C**** * ELEMENTS OF DIST. * C**** * N INTEGER NUMBER OF ELEMENTS IN DIST,DISP & W. * C**** * * C**** * * C**** * LOCAL VARIABLES: * C**** * * C**** * NAME TYPE DESCRIPTION * C**** * * C**** * SDS DEFAULT IN LOOP 3: NUMBER EVENTUALLY CONTAINING * C**** * THE WEIGHTED SUM OF AN * C**** * ACTIVE (ORDER-VIOLATING) * C**** * BLOCK. * C**** * SW DEFAULT IN LOOP 3: NUMBER EVENTUALLY CONTAINING * C**** * THE SUM OF THE WEIGHTS OF * C**** * AN ORDER-VIOLATING BLOCK. * C**** * M INTEGER EQUALS N, BUT HALF ITS SIZE. * C**** * I INTEGER INDEX OF MAIN LOOP. * C**** * IT INTEGER INDEX OF ORDER-VIOLATING BLOCK * C**** * TRANSFORMATION LOOP. * C**** * * C**** * * C**** * REMARK: * C**** * TO AVOID EXECUTION-TIME INEFFICIENCIES OR IMPRECISE * C**** * RESULTS SDS, SW, DIST, DISP AND W SHOULD BE OF THE SAME * C**** * TYPE, DOUBLE OR EXTENDED PRECISION. * C**** * * C**** * (ERNST VAN WANING). * C**** * * C**** ****************************************************************** INTEGER I,M,J,IT,IT1 DOUBLE PRECISION DIST,W,SDS,SW DIMENSION DIST(N),W(N) M=N DO 1 I=2,M IF (DIST(I) .GE. DIST(I-1)) GOTO 1 SDS=DIST(I)*W(I) SW=W(I) IT=I 3 IT=IT-1 SDS=SDS+DIST(IT)*W(IT) SW=SW+W(IT) DIST(IT)=SDS/SW IF (IT .EQ. 1) GOTO 2 IF (DIST(IT-1) .GT. DIST(IT)) GOTO 3 2 IT1=IT+1 DO 4 J=IT1,I 4 DIST(J)=DIST(IT) 1 CONTINUE RETURN END SUBROUTINE OUTDAT(IN,LB,X,NK,IT,N,LV,IP,ID,IK) C************************************************************* C*** * C*** PRINT OUTPUT OF THE CATEGORY NUMBERS AND CATEGORY * C*** QUANTIFICATIONS OF DATA VECTOR NO LV. * C*** OUTPUT OF VARIABLE NUMBER,CATEGORY NUMBERS AND * C*** CATEGORY QUANTIFICATIONS TO UNIT ID WITH FORMAT * C*** (2I5,F8.3) FOR EACH CATEGORY * C*** * C************************************************************* C*** INTEGER N,NK,IT,IP,ID,LV,IL,IN,L,K,I,IK DOUBLE PRECISION X REAL Z DIMENSION IN(N),LB(NK),X(N) C*** 94 FORMAT(1H+,88X,9H(NOMINAL)/) 95 FORMAT(1H+,88X,9H(ORDINAL)/) 96 FORMAT(1H+,88X,11H(NUMERICAL)/) 98 FORMAT(82H0CATEGORY NUMBERS,MARGINAL FREQUENCIES AND CATEGORY QUAN *TIFICATIONS OF VARIABLE NO,I4) 100 FORMAT(1H ,I5,2X,I5,2X,F10.3) 102 FORMAT(2I5,F8.3) 103 FORMAT(8H MISSING,I5,2X,F10.3,5X,3H(NO,I5,1H)) C*** Z=0.0 WRITE(IP,98)LV IF(IT.EQ.0)WRITE(IP,94) IF(IT.EQ.1)WRITE(IP,95) IF(IT.EQ.2)WRITE(IP,96) L=0 DO 105 K=1,NK IF(LB(K).EQ.0) WRITE(IP,100)K,LB(K),Z IF(LB(K).EQ.0.AND.ID.GT.1) WRITE(ID,102) LV,K,Z IF(LB(K).EQ.0) GOTO 105 L=L+LB(K) IL=IN(L) WRITE(IP,100)K,LB(K),X(IL) IF(ID.GT.1) WRITE(ID,102)LV,K,X(IL) 105 CONTINUE IF(L.EQ.N.OR.IK.EQ.0) GOTO 150 K=L+1 DO 120 L=K,N IL=1 120 WRITE(IP,103)IL,X(IN(L)),IN(L) 150 CONTINUE RETURN END SUBROUTINE OUTPAR(ID,IU,IZ,IP,LR) C******************************************* C*** * C*** COMMENT ON OUTPUT TO OTHER UNITS * C*** THAN THE PRINTER * C*** * C******************************************* C*** 1 FORMAT(1H1) 10 FORMAT(/////49H FOR EACH CATEGORY THE VARIABLE NUMBER, THE CATEG, + 52HORY NUMBER AND THE CATEGORY QUANTIFICATION HAVE BEEN/ + 17H WRITTEN TO UNIT,I3,24H WITH FORMAT (2I5,F8.3)., + //54H FOR EACH VARIABLE THE VARIABLE NUMBER AND THE WEIGHTS, + 26H HAVE BEEN WRITTEN TO UNIT,I3, + 36H WITH FORMAT (I5,9F8.3/,(5X,9F8.3))., + //46H FOR EACH VARIABLE THE VARIABLE NUMBER AND THE, + 51H CORRELATIONS (FOR THE FIRST SET AND THE SECOND SET, + 14H RESPECTIVILY)/26H HAVE BEEN WRITTEN TO UNIT,I3, + 36H WITH FORMAT (I5,9F8.3/,(5X,9F8.3)).,/////) 20 FORMAT(49H THE OPTIMALLY SCALED VARIABLES HAVE BEEN WRITTEN, + 8H TO UNIT,I3,22H WITH FORMAT (10F8.3)./ + 52H THE ROWS ARE OBJECTS AND THE COLUMNS ARE VARIABLES., + /////) 30 FORMAT(52H FOR EACH OBJECT THE OBJECT NUMBER AND THE CANONICAL, + 55H SCORES (FOR EACH DIMENSION OF THE FIRST SET AND SECOND, + 5H SET /40H RESPECTIVILY) HAVE BEEN WRITTEN TO UNIT,I3, + 36H WITH FORMAT (I5,9F8.3/,(5X,9F8.3)).,/////) 40 FORMAT(48H FOR EACH VARIABLE THE VARIABLE NUMBER, CATEGORY, + 59H NUMBER AND THE PROJECTIONS OF THE CATEGORY QUANTIFICATIONS/ + 54H ON THE CANONICAL VARIATES OF THE FIRST AND SECOND SET, + 39H RESPECTIVILY HAVE BEEN WRITTEN TO UNIT,I3,5H WITH/ + 33H FORMAT (2I5,8F8.3/,(10X,8F8.3)).,/////) C*** WRITE(IP,1) IF(ID.GT.1)WRITE(IP,10)ID,ID,ID IF(IU.GT.1)WRITE(IP,20)IU IF(IZ.GT.1)WRITE(IP,30)IZ IF(LR.GT.1)WRITE(IP,40)LR RETURN END SUBROUTINE OUTVAR(XY,N,M,IU,IP) C************************************************* C*** * C*** OUTPUT OF THE OPTIMALLY SCALED * C*** VARIABLES AS A (N X M) MATRIX * C*** TO UNIT IU * C*** * C************************************************* C*** DIMENSION XY(N,M) DOUBLE PRECISION XY INTEGER IU,IP C*** IF(IU.LE.1)RETURN DO 10 I=1,N 10 WRITE(IU,20)(XY(I,J),J=1,M) 20 FORMAT(10F8.3) RETURN END SUBROUTINE OUTWEI(A,B,M1,M2,P,IP,ID) C************************************************* C*** * C*** PRINT OUTPUT OF THE VARIABLE WEIGHTS * C*** FOR EACH DIMENSION * C*** * C************************************************* C*** INTEGER J,K,IP,L,P,ID DOUBLE PRECISION A,B DIMENSION A(M1,P),B(M2,P) WRITE(IP,99) DO 200 J=1,M1 WRITE(IP,98)J,(A(J,K),K=1,P) 200 IF(ID.GT.1)WRITE(ID,97)J,(A(J,K),K=1,P) WRITE(IP,100) DO 250 J=1,M2 L=J+M1 WRITE(IP,98)L,(B(J,K),K=1,P) 250 IF(ID.GT.1)WRITE(ID,97)L,(B(J,K),K=1,P) C 97 FORMAT(I5,9F8.3/(5X,9F8.3)) 98 FORMAT(1H ,I5,4X,15F8.3) 99 FORMAT(53H0VARIABLE WEIGHTS OF THE FIRST SET FOR EACH DIMENSION,/) 100 FORMAT(54H0VARIABLE WEIGHTS OF THE SECOND SET FOR EACH DIMENSION/) RETURN END SUBROUTINE PLOT(X,Y,IND,NITEM,NRC,IDENT,IWRITE) C ****************************************************************** C * * C * P L O T * C * * C * PURPOSE: PRINTPLOT OF Y VERSUS X * C * IF IDENT.GT.0, THE FIRST NRC POINTS ARE IDENTIFIED BY* C * INTEGER NUMBERS 1-9,0,1-9... , THE REST OF THE POINTS* C * BY CHARACTERS A-I,J,A-I... * C * IF IDENT.EQ.0, ALL POINTS ARE PRINTED AS STARS * C * IF IDENT.LT.0, THE FIRST NRC POINTS ARE IDENTIFIED BY* C * STARS, THE OTHERS BY A 'D' * C * * C * SUBROUTINES CALLED: NONE * C * * C * AUTHOR: WILLEM HEISER RELEASED: MARCH 1978 * C * ADAPTED : NOVEMBER 1982 * C * * C ****************************************************************** C C TO ADJUST THE PLOT CHANGE THE CONSTANTS UNDER WHICH C ## IS PLACED ACCORDING TO THE COMMENTS C DIMENSION LINE(71),XPR(11),X(NITEM),Y(NITEM),LABEL(20),IND(NITEM) C ##=NN+1 (NN=INT(NC/7)) C ##=NUMBER OF COLUMNS PER LINE=NC INTEGER BLANK,STAR DATA BLANK,STAR,MORE/1H ,1H*,1HM/,NL/56/,NC/71/,NN/10/ C ##=NC C ##=NUMBER OF LINES PER PAGE DATA LABEL/1H0,1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9, X 1HJ,1HA,1HB,1HC,1HD,1HE,1HF,1HG,1HH,1HI/ C 1 FORMAT(16X,3H.--,10(7H+------),5H+---.) C ##=NN 2 FORMAT(8X,F7.3,1X,1H!,2X,71A1,3X,1H!) C ##=NC 3 FORMAT(8X,F7.3,1X,1H!,76X,1H!) C ##=NC+5 4 FORMAT(15X,11F7.3) C ##=NN+1 C DO 10 I=1,NITEM 10 IND(I)=I C C THE INDICES ARE ORDERED SUCH THAT THEY WOULD ORDER Y IN DESCENDING C ORDER C M=NITEM 20 M=M/2 IF (M.LT.1) GO TO 40 K=NITEM-M J=1 41 I=J 49 L=I+M IF (Y(IND(I)).GE.Y(IND(L))) GO TO 60 II=IND(I) IND(I)=IND(L) IND(L)=II I=I-M IF (I.GE.1) GO TO 49 60 J=J+1 IF (J-K) 41,41,20 40 CONTINUE C C THE PLOT IS ADJUSTED TO THE RANGE OF THE 'LONGEST' AXIS C XMAX=X(1) XMIN=X(1) DO 70 I=2,NITEM IF (X(I).GT.XMAX) XMAX=X(I) IF (X(I).LT.XMIN) XMIN=X(I) 70 CONTINUE SPANX=XMAX-XMIN SPANY=Y(IND(1))-Y(IND(NITEM)) IF (SPANY.LE.SPANX) GO TO 75 XMIN=XMIN-(SPANY-SPANX)/2.0 SPANX=SPANY 75 STEPX=SPANX/(7.0*NN) STEPY=(STEPX*NC)/(NL+0.) HSTEP=STEPY/2.0 TOP=(Y(IND(NITEM))+Y(IND(1))+SPANX)/2.0 C C THE POINTS ARE PLOTTED LINE AFTER LINE C WRITE(IWRITE,1) L=1 I=1 YPR=TOP+STEPY 80 YPR=YPR-STEPY IF (I.GT.NITEM) GO TO 110 IF ((YPR-Y(IND(I))).GT.HSTEP) GO TO 110 90 DO 95 J=1,NC 95 LINE(J)=BLANK 100 JP=(X(IND(I))-XMIN)/STEPX+1.0 IF (LINE(JP).EQ.BLANK) GO TO 101 LINE(JP)=MORE GO TO 103 101 IF (IDENT.EQ.0) GO TO 102 IF (IDENT.GT.0) GO TO 104 IF (IND(I).LE.NRC) LINE(JP)=STAR IF (IND(I).GT.NRC) LINE(JP)=LABEL(15) GO TO 103 104 IF (IND(I).LE.NRC) LINE(JP)=LABEL(MOD(IND(I),10)+1) IF (IND(I).GT.NRC) X LINE(JP)=LABEL(MOD((IND(I)-NRC),10)+11) GO TO 103 102 LINE(JP)=STAR 103 I=I+1 IF (I.GT.NITEM) GO TO 105 IF ((YPR-Y(IND(I))).LT.HSTEP) GO TO 100 105 WRITE(IWRITE,2) YPR,LINE GO TO 115 110 WRITE(IWRITE,3) YPR 115 L=L+1 IF (L.LE.NL) GO TO 80 WRITE(IWRITE,1) C C VALUES OF X-AXIS ARE PRINTED C XPR(1)=XMIN DO 130 J=1,10 130 XPR(J+1)=XPR(J)+STEPX*7.0 WRITE(IWRITE,4) XPR RETURN END SUBROUTINE PLOTCA(XA,YB,H,N,P,IZ,IP,LZ) C************************************************* C*** * C*** PRINT/PLOT/OUTPUT OF CANONICAL SCORES * C*** H WORK ARRAY * C*** IZ=0 NO OUTPUT,IZ=1 OR -1 PRINT * C*** IZ.GT.1 OR IZ.LT.-1 OUTPUT TO UNIT ABS(IZ) * C*** LZ=0 NO PLOTS,LZ=1 PLOTS * C*** * C*** SUBROUTINES CALLED: PLOT * C*** FUNCTIONS CALLED: IABS * C*** * C************************************************* C*** INTEGER N,P,IZ,IP,LZ DOUBLE PRECISION XA,YB REAL H DIMENSION XA(N,P),YB(N,P),H(N,3) IZ=IABS(IZ) IF(IZ.EQ.0) GOTO 102 IF(IZ.GT.1) GOTO 101 C IZ=1 OR -1 WRITE(IP,1) DO 100 I=1,N 100 WRITE(IP,2)I,(XA(I,L),L=1,P),(YB(I,L),L=1,P) GOTO 102 C IZ GT 1 OR LT -1 101 CONTINUE DO 200 I=1,N 200 WRITE(IZ,6)I,(XA(I,L),L=1,P),(YB(I,L),L=1,P) 102 IF(LZ.EQ.0)RETURN C LZ=1 IF(P.EQ.1)RETURN MP=P-1 DO 110 L=1,MP LP=L+1 DO 110 K=LP,P WRITE(IP,3) WRITE(IP,5)L,K DO 105 I=1,N H(I,1)=XA(I,L) 105 H(I,2)=XA(I,K) CALL PLOT(H(1,1),H(1,2),H(1,3),N,N,0,IP) WRITE(IP,4) WRITE(IP,5)L,K DO 107 I=1,N H(I,1)=YB(I,L) 107 H(I,2)=YB(I,K) CALL PLOT(H(1,1),H(1,2),H(1,3),N,N,0,IP) 110 CONTINUE 1 FORMAT(79H1 CANONICAL SCORES FOR EACH DIMENSION OF THE FIRST RESP *ECTIVELY THE SECOND SET//) 2 FORMAT(I5,9F8.3/,(5X,9F8.3)) 3 FORMAT(43H1 CANONICAL SCORES OF THE FIRST SET) 4 FORMAT(44H1 CANONICAL SCORES OF THE SECOND SET) 5 FORMAT(1H+,46X,23H (HORIZONTAL VARIATE NO,I2,24H AND VERTICAL VARI *ATE NO,I2,1H)) 6 FORMAT(I5,9F8.3/,(5X,9F8.3)) RETURN END SUBROUTINE PLOTOP(XY,IN,LB,H1,H2,H3,IS,IP,N,NK,IT,NM,J) C*********************************************************** C*** * C*** PLOT OF THE CATEGORY QUANTIFICATIONS AND CATEGORY * C*** NUMBERS * C*** SUBROUTINES CALLED: PLOT * C*** * C*********************************************************** C*** DIMENSION XY(N),IN(N),LB(NK),H1(NK),H2(NK),H3(NK) DOUBLE PRECISION XY REAL H1,H2,H3 INTEGER IN,LB,NM,IS,IP,IT,J,L,K C*** READ(IS)IN,LB IF(NM.EQ.0) RETURN L=0 K1=0 DO 20 K=1,NK IF(LB(K).EQ.0) GOTO 20 L=L+LB(K) K1=K1+1 H2(K1)=XY(IN(L)) H1(K1)=K 20 CONTINUE WRITE(IP,1)J IF(IT.EQ.0)WRITE(IP,2) IF(IT.EQ.1)WRITE(IP,3) IF(IT.EQ.2)WRITE(IP,4) CALL PLOT(H1,H2,H3,K1,K1,0,IP) 1 FORMAT(92H1 CATEGORY NUMBERS (HORIZONTAL) AGAINST CATEGORY QUAN *TIFICATION (VERTICAL) OF VARIABLE NO,I5) 2 FORMAT(1H+,97X,9H(NOMINAL)) 3 FORMAT(1H+,97X,9H(ORDINAL)) 4 FORMAT(1H+,97X,11H(NUMERICAL)) RETURN END SUBROUTINE PRODCT(X,A,XA,N,M1,P) C********************************************** C*** * C*** COMPUTATION OF THE MATRIX PRODUCT X*A * C*** * C********************************************** C*** INTEGER N,M1,P,I,J,K DOUBLE PRECISION X,A,XA DIMENSION X(N,M1),A(M1,P),XA(N,P) C*** DO 50 I=1,N DO 50 K=1,P 50 XA(I,K)=0.0D0 C*** DO 60 J=1,M1 DO 60 K=1,P IF(DABS(A(J,K)).LT.1.0D-15) GOTO 60 DO 55 I=1,N 55 XA(I,K)=XA(I,K)+X(I,J)*A(J,K) 60 CONTINUE C*** RETURN END SUBROUTINE RANDMA(A,N) C********************************************************* C*** * C*** MACHINE INDEPENDENT RANDOM GENERATOR * C*** N RANDOM NUMBERS ARE PLACED IN ARRAY A * C*** WARNING:THE NECESSARY AVAILABLE INTEGER ARITHMATIC * C*** IS 125*2796203=348525375 (<2**29). * C*** * C********************************************************* DIMENSION A(N) DOUBLE PRECISION A,Z IY=315747 M=2796203 Z=2.0D0/M DO 10 I=1,N IY=IY*125 IY=IY-(IY/M)*M 10 A(I)=-1.0D0+IY*Z RETURN END SUBROUTINE REGRSX(IN,LB,X,H1,H2,H3,H5,A,N,P,NK,IT,IS) C*********************************************************************** C*** * C*** REGRESSION ON VECTOR X * C*** AT RETURN X IS ELEMENT OF THE CONE SPANNED BY THE ORIGINAL DATA. * C*** IN THE NOMINAL CASE REGRESSION MEANS COMPUTATION OF THE * C*** TIE BLOCK MEANS AND REPLACING THE OLD ELEMENTS BY THESE MEANS * C*** IN THE OTHER CASES THE TIE BLOCK MEANS ARE ALSO CACULATED * C*** AND A MONOTONE REGRESSION ON THESE MEANS IS DONE IN SUBROUTINE * C*** MTR (ORDINAL) OR A LINEAR REGRESSION IN SUBR.LIDICA (NUMERICAL). * C*** AFTER REGRESSION THE OLD ELEMENTS OF VECTOR X ARE REPLACED BY * C*** THE REGRESSED MEANS. * C*** SUBROUTINES CALLED: MTR,LIDICA * C*** * C*********************************************************************** C*** INTEGER I,K,L,L1,L2,NK,IT,N,P,IN,LB,IS DOUBLE PRECISION X,H1,H2,H3,H5,A DIMENSION IN(N),LB(NK),X(N),H1(NK),H2(NK),H3(NK),H5(N),A(P) C*** READ(IS)IN,(LB(K),K=1,NK) C*** INDICES AND TIE BLOCK LENGTHS L2=0 DO 250 K=1,NK H2(K)=0.0D0 IF(LB(K).EQ.0) GOTO 250 L1=L2+1 L2=L2+LB(K) DO 240 L=L1,L2 IL=IN(L) H2(K)=H2(K)+X(IL) 240 CONTINUE H2(K)=H2(K)/LB(K) 250 CONTINUE C*** H2 CONTAINS THE TIE BLOCK MEANS C*** IF(IT.EQ.0) GOTO 270 IF(IT.EQ.2) GOTO 260 L=0 DO 251 K=1,NK IF(LB(K).EQ.0) GOTO 251 L=L+1 H3(L)=LB(K) H1(L)=H2(K) 251 CONTINUE C*** H3 REGRESSION WEIGHTS,ALL POSITIVE CALL MTR(H1,H3,L) L=0 DO 252 K=1,NK IF(LB(K).EQ.0) GOTO 252 L=L+1 H2(K)=H1(L) 252 CONTINUE C*** C*** IN CASE ALL THE ENTRIES WERE MADE EQUAL IN THE REGRESSION C*** THE VARIABLE WEIGHTS ARE MADE ZERO AND C*** THE ENTRIES GET THE VALUES OF THE PREVIOUS ITERATION IF(DABS(H1(1)-H1(L)).GT.1.D-15) GOTO 270 IF(NK.EQ.1)GOTO 270 DO 255 I=1,N 255 X(I)=H5(I) DO 256 K=1,P 256 A(K)=0.0D0 RETURN C*** 260 IF(IT.EQ.2) CALL LIDICA(LB,H2,H1,H3,NK) C*** 270 CONTINUE L2=0 DO 290 K=1,NK IF(LB(K).EQ.0)GOTO 290 L1=L2+1 L2=L2+LB(K) DO 280 L=L1,L2 IL=IN(L) X(IL)=H2(K) 280 CONTINUE 290 CONTINUE RETURN END SUBROUTINE SINGCO(H11,H10,XY,IN,LB,NK,N,P,M,MC,ID,IP,P2,LR,IS) C*** C********************************************************** C*** * C*** PRINT, WRITE TO UNIT OF THE PROJECTIONS OF THE * C*** CATEGORY QUANTIFICATIONS ON THE CANONICAL * C*** VARIATES LR=0 NO OUTPUT LR=1 PRINT OUTPUT * C*** LR GT 1 OUTPUT TO UNIT * C*** * C*** IF P=1 H11 NEEDS STORAGE OF H7 AND H8 * C*** IF P>1 H11 NEEDS STORAGE OF H7 ONLY * C********************************************************** C INTEGER IN,LB,NK,N,P,M,MC,ID,IP,P2,LR DIMENSION H10(M,P2),XY(N,M),IN(N),LB(MC),NK(M),H11(P2) REAL H10,H11 DOUBLE PRECISION XY IF(LR.EQ.0)RETURN IF(LR.EQ.1)WRITE(IP,900) Z=0.0 REWIND IS DO 500 J=1,M NC=NK(J) READ(IS)IN,(LB(K),K=1,NC) L=0 DO 400 K=1,NC IF(LB(K).EQ.0.AND.LR.EQ.1) WRITE(IP,902) J,K,(Z,I=1,P2) IF(LB(K).EQ.0.AND.LR.GT.1) WRITE(LR,901) J,K,(Z,I=1,P2) IF(LB(K).EQ.0) GOTO 400 L=L+LB(K) IL=IN(L) DO 300 I=1,P2 300 H11(I)=XY(IL,J)*H10(J,I) IF(LR.EQ.1) WRITE(IP,902) J,K,(H11(I),I=1,P2) IF(LR.GT.1) WRITE(LR,901) J,K,(H11(I),I=1,P2) 400 CONTINUE 500 CONTINUE 900 FORMAT(1H1,73HTHE PROJECTIONS OF THE CATEGORY QUANTIFICATIONS ON T +HE CANONICAL VARIATES//) 901 FORMAT(2I5,8F8.3/,(10X,8F8.3)) 902 FORMAT(1H ,2I5,8F8.3/,(1H ,10X,8F8.3)) RETURN END SUBROUTINE SORT(A,II,JJ,IN,N) C*********************************************************************** C*** * C*** A VECTOR TO BE SORTED * C*** IN INDEX VECTOR,CONTAINS THE ORIGINAL INDEX NUMBERS AFTER RETURN * C*** AT THE BEGINNING OF THE SUBROUTINE IN(K) EQUALS K * C*** IF THE INDICES ARE KNOWN FROM OUTSIDE THE SUBR,DO LOOP 3 CAN BE * C*** REMOVED * C*** II FIRST ,JJ LAST ELEMENT TO BE SORTED * C*** N CAN BE REMOVED FROM THE CALL LIST,IF A AND IN ARE DIMENSIONED 1* C*** * C*********************************************************************** C*** C*** INTEGER IU,IL,IN,II,JJ,IH,IG,I,J,M,K,IJ INTEGER A,T,TT C*** DIMENSION IU(16),IL(16),A(N),IN(N) C*** C*** M=1 I=II J=JJ C*** DO 3 K=II,JJ 3 IN(K)=K C*** 5 IF (I.GE.J) GOTO 70 10 K=I IJ=(I+J)/2 T=A(IJ) IH=IN(IJ) IF (A(I).LE.T)GOTO 20 A(IJ)=A(I) A(I)=T T=A(IJ) IN(IJ)=IN(I) IN(I)=IH IH=IN(IJ) 20 L=J IF (A(J).GE.T) GOTO 40 A(IJ)=A(J) A(J)=T T=A(IJ) IN(IJ)=IN(J) IN(J)=IH IH=IN(IJ) IF(A(I).LE.T) GOTO 40 A(IJ)=A(I) A(I)=T T=A(IJ) IN(IJ)=IN(I) IN(I)=IH IH=IN(IJ) GOTO 40 30 A(L)=A(K) A(K)=TT IN(L)=IN(K) IN(K)=IG 40 L=L-1 IF (A(L).GT.T) GOTO 40 TT=A(L) IG=IN(L) 50 K=K+1 IF (A(K).LT.T) GOTO 50 IF (K.LE.L) GOTO 30 IF (L-I.LE.J-K) GOTO 60 IL(M)=I IU(M)=L I=K M=M+1 GOTO 80 60 IL(M)=K IU(M)=J J=L M=M+1 GOTO 80 70 M=M-1 IF (M.EQ.0) RETURN I=IL(M) J=IU(M) 80 IF(J-I.GE.11) GOTO 10 IF(I.EQ.II) GOTO 5 I=I-1 90 I=I+1 IF(I.EQ.J) GOTO 70 T=A(I+1) IH=IN(I+1) IF(A(I).LE.T) GOTO 90 K=I 100 A(K+1)=A(K) IN(K+1)=IN(K) K=K-1 IF (T.LT.A(K)) GOTO 100 A(K+1)=T IN(K+1)=IH GOTO 90 END SUBROUTINE SQDEV(A,N,AN,L,LC) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C SUBROUTINE SQDEV(A,N,AN,L,LC) C C PURPOSE C C SET A VECTOR IN DEVIATION OF THE MEAN AND MAKE THE C C SUM OF SQUARES EQUAL TO N OR 1 DEPENDING ON L C C THE LOGICAL LC IS TRUE IF THE STANDARD DEVIATION C C EQUALS ZERO C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C DOUBLE PRECISION A,SY,BY,Y,FN,BS,RY,AN INTEGER I,N,L LOGICAL LC DIMENSION A(N) C C****** COMPUTE SUM AND SUM OF SQUARES ************************ C SY=0.0D0 BY=0.0D0 DO 1 I=1,N Y=A(I) IF(DABS(Y).LT.1.0E-20) GOTO 1 SY=SY+Y BY=BY+Y*Y 1 CONTINUE C C****** COMPUTE THE SQUARE ROOTOF N DIVIDED BY THE SUM OF THE SQUARED C****** DEVIATIONS OF THE MEAN IN RY C FN=DFLOAT(N) LC=.FALSE. C**** BS=SY*AN RY=2*BS*SY SY=N*BS*BS RY=BY+SY-RY 4 IF(RY.LT.1.0D-20) LC=.TRUE. IF(LC) GOTO 5 RY=1.0D0/DSQRT(RY) IF(L.NE.1) RY=DSQRT(FN)*RY C**** C****** MULTIPLY THE SCORES IN DEVIATION OF THE MEAN WITH RY C 3 DO 2 I=1,N Y=A(I) Y=Y-BS Y=Y*RY A(I)=Y 2 CONTINUE RETURN C*** 5 DO 6 I=1,N 6 A(I)=0.0D0 RETURN END SUBROUTINE STRESS(XA,YB,N,P,S) C********************************** C*** * C*** S=TRACE(XA-YB)T(XA-YB)/N*P * C*** T MEANS TRANSPOSE * C*** * C********************************** C*** DOUBLE PRECISION XA,YB,S INTEGER N,P,K,I DIMENSION XA(N,P),YB(N,P) C*** S=0.0D0 DO 50 K=1,P DO 50 I=1,N 50 S=S+(XA(I,K)-YB(I,K))**2 S=S/DFLOAT(N*P) RETURN END SUBROUTINE TRED2(NM,N,D,E,Z) C ****************************************************************** C * * C * T R E D 2 * C * * C * PURPOSE: REDUCES A REAL SYMMETRIC MATRIX TO A SYMMETRIC TRI- * C * DIAGONAL MATRIX USING ORTHOGONAL SIMILARITY TRANS- * C * FORMATIONS; THE ACCUMULATED ORTHOGONAL TRANSFORMA- * C * TIONS ARE RETAINED IN Z. * C * * C * SUBROUTINES CALLED: NONE * C * * C * ADAPTED FROM EISPACK-ORIGINAL VERSION NOT IN DOUBLE PRECISION * C * * C ****************************************************************** DOUBLE PRECISION D,E,Z,H,SCALE,F,G,HH DOUBLE PRECISION DABS,DSQRT,DSIGN DIMENSION Z(NM,N),D(N),E(N) C IF (N .EQ. 1) GO TO 320 C ********** FOR I=N STEP -1 UNTIL 2 DO -- ********** DO 300 II =2, N I = N + 2 - II L = I - 1 H = 0.0 SCALE = 0.0 IF (L .LT.2) GO TO 130 C ********** SCALE ROW (ALGOL TOL THEN NOT NEEDED) ********** DO 120 K =1, L 120 SCALE = SCALE + DABS(Z(I,K)) C IF (SCALE .NE. 0.0) GO TO 140 130 E(I) = Z(I,L) GO TO 290 C 140 DO 150 K = 1, L Z(I,K) = Z(I,K) / SCALE H = H + Z(I,K) * Z(I,K) 150 CONTINUE C F = Z(I,L) G = -DSIGN(DSQRT(H),F) E(I) = SCALE * G H = H - F * G Z(I,L) = F - G F = 0.0 C DO 240 J = 1, L Z(J,I) = Z(I,J) / (SCALE * H) G = 0.0 C ********** FORM ELEMENT OF A*U ********** DO 180 K = 1, J 180 G = G + Z(J,K) * Z(I,K) C JP1 = J + 1 IF (L .LT. JP1) GO TO 220 C DO 200 K = JP1, L 200 G = G + Z(K,J) * Z(I,K) C ********** FORM ELEMENT OF P ********** 220 E(J) = G / H F = F + E(J) * Z(I,J) 240 CONTINUE C HH = F / (H + H) C ********** FORM REDUCED A ********** DO 260 J =1, L F = Z(I,J) G = E(J) - HH * F E(J) = G C DO 260 K = 1, J Z(J,K) = Z(J,K) - F * E(K) -G * Z(I,K) 260 CONTINUE C DO 280 K = 1, L 280 Z(I,K) = SCALE * Z(I,K) C 290 D(I) = H 300 CONTINUE C 320 D(1) = 0.0 E(1) = 0.0 C ********** ACCUMULATION OF TRANSFORMATION MATRICES ********** DO 500 I = 1, N L = I - 1 IF (D(I) .EQ. 0.0) GO TO 380 C C DO 360 J = 1, L G = 0.0 DO 340 K =1, L 340 G = G + Z(I,K) * Z(K,J) C DO 360 K = 1, L Z(K,J) = Z(K,J) - G * Z(K,I) 360 CONTINUE C 380 D(I) = Z(I,I) Z(I,I) = 1.0 IF (L .LT. 1) GO TO 500 C DO 400 J =1, L Z(I,J) = 0.0 Z(J,I) = 0.0 400 CONTINUE C 500 CONTINUE C RETURN END