EXTERNAL C82SPA DIMENSION FMT(20),TITLE(20) LOGICAL CONF COMMON/FMATOR/IE1,IE2,ALPHA,BETA,ICAR,IPLTRA,IPL2DI,IPLEXT, 1 MEDSCA,MEDVAR,FMT,CONF,ILIST,IPERM,MEDPER WRITE(6,5000) READ (5,1200) TITLE READ (5,1000) K1,K2 READ (5,1400) SN1,IE1,IE2,ICONF READ (5,1000) ICAR,MEDSCA,MEDVAR,ILIST,IPERM,MEDPER READ (5,1000) IPLTRA,IPL2DI,IPLEXT READ (5,1200) FMT ALPHA=(1.+SN1)/2. BETA=(1.-SN1)/2. WRITE (6,5100) TITLE,K1,K2,SN1 IF (SN1.EQ.-1.) WRITE(6,5500) IF (SN1.EQ.1.) WRITE(6,5600) IF (SN1.EQ.0.) WRITE(6,5700) MINK=K1 IF (K1.GT.K2) MINK=K2 IF (IE1.GE.MINK) GOTO 20 IF (IE1.LE.0) IE1=1 IF (IE2.LT.IE1) IE2=IE1+1 IF (IE2.GE.MINK) IE2=MINK-1 NALDIM=IE2-IE1+1 CONF=(ICONF.EQ.1) IF (CONF) GOTO 10 ICONF=0 MEDVAR=0 10 IF (IPLTRA.LT.0.OR.IPLTRA.GT.4) IPLTRA=0 IF (IPL2DI.LT.0.OR.IPL2DI.GT.4) IPL2DI=0 IF (IPLEXT.NE.1.OR.K1.NE.K2) IPLEXT=0 IF (ICAR.LE.0) ICAR=5 WRITE (6,5900) IE1,IE2,ICONF IF (ICONF.EQ.1) WRITE (6,6000) IF (ICONF.EQ.0) WRITE (6,6100) WRITE (6,6200) ICAR,MEDSCA IF (MEDSCA.EQ.0) WRITE (6,7000) WRITE (6,7100) MEDVAR IF (MEDVAR.EQ.0) WRITE (6,7000) IF (ILIST.EQ.0) ILIST=10 IF (ILIST.EQ.-1.OR.ILIST.GT.K1) ILIST=K1 IF (IPERM.GT.IE2) IPERM=IE2 WRITE (6,5200) ILIST,IE1,IPERM IF (IPERM.LT.IE1) WRITE (6,5300) IF (IPERM.LT.IE1) MEDPER=0 WRITE (6,5400) MEDPER IF (MEDPER.EQ.0) WRITE (6,7000) WRITE (6,7200) IPLTRA IF (IPLTRA.EQ.0) WRITE (6,6300) IF (IPLTRA.EQ.1) WRITE (6,6400) IF (IPLTRA.EQ.2) WRITE (6,6500) IF (IPLTRA.EQ.3) WRITE (6,6600) IF (IPLTRA.EQ.4) WRITE (6,6700) IF (IE1.GE.IE2) IPL2DI=0 WRITE (6,7300) IPL2DI IF (IPL2DI.EQ.0) WRITE (6,6300) IF (IPL2DI.EQ.1) WRITE (6,6400) IF (IPL2DI.EQ.2) WRITE (6,6500) IF (IPL2DI.EQ.3) WRITE (6,6600) IF (IPL2DI.EQ.4) WRITE (6,6700) WRITE (6,7400) IPLEXT IF (IPLEXT.EQ.0) WRITE (6,6300) IF (IPLEXT.EQ.1) WRITE (6,7500) WRITE (6,6800) FMT ISUM=1+K1*K2+K1*MINK+K2*MINK+3*K1+3*K2+MINK+ 1 ICONF*(K1+K2+1)*NALDIM*(NALDIM+3)/2 CALL C82DEC(C82SPA,A,ISUM,K1,K2,MINK,NALDIM) RETURN 20 WRITE (6,5800) 1000 FORMAT(16I5) 1200 FORMAT(20A4) 1400 FORMAT(F5.2,15I5) 5000 FORMAT(1H1///1H0,5X,11HA N A C O R,82X, 1 14HBERT BETTONVIL/ 2 1H0,5X,14HSEPTEMBER 1981,79X,24HDEPARTMENT OF DATATHEORY/ 3 1H0,98X,13HBREESTRAAT 70/ 4 1H0,5X,45HCORRESPONDENCE ANALYSIS OF CONTINGENCY TABLES, 5 48X,24HLEIDEN - THE NETHERLANDS////) 5100 FORMAT(///1H0,2X,9HJOB TITLE,36X,20A4// 1 1H0,2X,19HDATA SPECIFICATIONS/ 2 1H0,2X,14HNUMBER OF ROWS,31X,I5/ 3 1H ,2X,17HNUMBER OF COLUMNS,28X,I5// 4 1H0,2X,23HANALYSIS SPECIFICATIONS/ 5 1H0,2X,24HINDEX OF STANDARDIZATION,21X,F8.2) 5500 FORMAT(1H+,60X,40HTHE COLUMNS ARE IN THE CENTROIDS OF THE , 1 13HMATCHING ROWS) 5600 FORMAT(1H+,60X,46HTHE ROWS ARE IN THE CENTROIDS OF THE MATCHING , 1 7HCOLUMNS) 5700 FORMAT(1H+,60X,42HROWS ARE COLUMNS ARE TREATED SYMMETRICALLY) 5900 FORMAT(41H FIRST OF THE DIMENSIONS TO BE ANALYZED,7X,I5/ 1 40H LAST OF THE DIMENSIONS TO BE ANALYZED,8X,I5/ 2 27H COMPUTATION OF VARIANCES,21X,I5) 6000 FORMAT(1H+,60X,22HVARIANCES ARE COMPUTED/) 6100 FORMAT(1H+,60X,26HVARIANCES ARE NOT COMPUTED/) 6200 FORMAT(1H0,2X,18HI/0 SPECIFICATIONS/ 1 1H0,2X,26HUNIT NUMBER FOR INPUT DATA,19X,I5/ 2 1H ,2X,29HUNIT NUMBER FOR SCORES OUTPUT,16X,I5) 7000 FORMAT(1H+,60X,14HNO SUCH OUTPUT) 7100 FORMAT(1H ,2X,32HUNIT NUMBER FOR VARIANCES OUTPUT,13X,I5) 5200 FORMAT(1H ,2X,40HNUMBER OF LISTED ROWS OF THE DATA-MATRIX, 1 5X,I5/1H ,2X,35HREORDERED DATA-MATRIX IN DIMENSIONS, 2 I4,6H UP TO,I5) 5300 FORMAT(1H+,60X,31HDATA-MATRICES ARE NOT REORDERED) 5400 FORMAT(1H ,2X,30HUNIT NUMBER FOR REORDERED DATA,15X,I5) 7200 FORMAT(/1H0,2X,19HPLOT SPECIFICATIONS/ 1 1H0,2X,28HPLOTS OF THE TRANSFORMATIONS,17X,I5) 7300 FORMAT(1H ,2X,21HTWO-DIMENSIONAL PLOTS,24X,I5) 7400 FORMAT(1H ,2X,34HPLOTS OF ROW- VERSUS COLUMN-SCORES,11X,I5) 7500 FORMAT(1H+,60X,20HTHESE PLOTS ARE MADE) 6300 FORMAT(1H+,60X,13HNO SUCH PLOTS) 6400 FORMAT(1H+,60X,17HONLY FOR THE ROWS) 6500 FORMAT(1H+,60X,20HONLY FOR THE COLUMNS) 6600 FORMAT(1H+,60X,35HSEPARATE PLOTS FOR ROWS AND COLUMNS) 6700 FORMAT(1H+,60X,35HCOMBINED PLOTS FOR ROWS AND COLUMNS) 6800 FORMAT(//1H0,2X,43HFORMAT OF THE ROWS OF THE INPUT DATA MATRIX, 1 2X,20A4) 5800 FORMAT(34H0 NO SOLUTION IN GIVEN DIMENSIONS) STOP END SUBROUTINE C82DEC(C82SUB,B,ISUM,K1,K2,MINK,NALDIM) DIMENSION A(20000) IDIM=20000 IF (ISUM.LE.IDIM) GO TO 10 WRITE (6,1000) ISUM RETURN 10 CONTINUE CALL C82SUB(A,ISUM,K1,K2,MINK,NALDIM) 1000 FORMAT (50H0ARRAY 'A' IN SUBROUTINE 'DECLAR' IS TOO SMALL FOR, 1 9H THIS JOB/21H 'A' MUST BE AT LEAST,I8,5H LONG// 2 45HEXECUTION OF THE PROGRAM HAS BEEN TERMINATED)) RETURN END SUBROUTINE C82SPA(A,ISUM,K1,K2,MINK,NALDIM) DIMENSION A(ISUM) C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C THIS SUBROUTINE DISTRIBUTES THE ARRAY A OVER THE VARIOUS ARRAYS C C NEEDED IN SUBROUTINE C82ORG C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C INDF=1 INDD1=INDF+K1*K2 INDD2=INDD1+K1 INDZ1=INDD2+K2 INDZ2=INDZ1+K1*MINK INDXL=INDZ2+K2*MINK INDI1=INDXL+MINK INDI2=INDI1+K1 INDH1=INDI2+K2 INDH2=INDH1+K1 INDVAR=INDH2+K2 IREST=ISUM-INDVAR+1 CALL C82ORG(K1,K2,MINK,NALDIM,IREST,A(INDF),A(INDD1),A(INDD2), 1 A(INDZ1),A(INDZ2),A(INDXL),A(INDF),A(INDD1), 2 A(INDD2),A(INDI1),A(INDI2),A(INDH1),A(INDH2),A(INDH1), 3 A(INDH2),A(INDVAR)) RETURN END SUBROUTINE C82ORG(K1,K2,MINK,NALDIM,IREST,F,D1,D2,Z1,Z2,XLABDA, 1 INTF,ID1,ID2,IDEN1,IDEN2,HULP1, 2 HULP2,IHULP1,IHULP2,VARS) DIMENSION F(K1,K2),D1(K1),D2(K2),Z1(K1,MINK),Z2(K2,MINK), 1 XLABDA(MINK),INTF(K1,K2),ID1(K1),ID2(K2),IDEN1(K1), 2 IDEN2(K2),HULP1(K1),HULP2(K2), 3 IHULP1(K1),IHULP2(K2),VARS(IREST),FMT(20) LOGICAL CONF COMMON/FMATOR/IE1,IE2,ALPHA,BETA,ICAR,IPLTRA,IPL2DI,IPLEXT, 1 MEDSCA,MEDVAR,FMT,CONF,ILIST,IPERM,MEDPER C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C IN THIS SUBROUTINE A FREQUENCY-TABLE F (K1 K2) WITH MARGINALS C C D1 AND D2 IS DECOMPOSED INTO Z1, Z2 AND XLABDA SUCH THAT C C F = D1.U.U'.D2 + D1.Z1.XLABDA**(1-ALPHA-BETA).Z1'.D2 C C BY MEANS OF SINGULAR VALUE DECOMPOSITION. C C IF CONF=.TRUE. THE VARIANCES OF THE Z-SCORES ARE ALSO C C COMPUTED C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C C C READ IN F, COMPUTE MARGINAL FREQUENCIES AND APPLY SVD C C C TOTAL=0. DO 10 J=1,K2 10 D2(J)=0. MIN=MINK-1 IF (K1.LT.K2) GOTO 50 DO 30 I=1,K1 H1=0. READ (ICAR,FMT) (INTF(I,J),J=1,K2) DO 20 J=1,K2 H2=INTF(I,J) Z1(I,J)=H2 H1=H1+H2 20 D2(J)=D2(J)+H2 TOTAL=TOTAL+H1 30 D1(I)=H1 DO 40 I=1,K1 H1=D1(I) DO 40 J=1,K2 H2=SQRT(H1*D2(J)) IF (H2.EQ.0.) GOTO 40 Z1(I,J)=Z1(I,J)/H2-H2/TOTAL 40 CONTINUE CALL SIVAD(K1,K2,XLABDA,Z1,K1*K2,Z2,K2*K2,HULP2,ICV) GOTO 90 50 DO 70 I=1,K1 H1=0. READ (ICAR,FMT) (INTF(I,J),J=1,K2) DO 60 J=1,K2 H2=INTF(I,J) Z2(J,I)=H2 H1=H1+H2 60 D2(J)=D2(J)+H2 TOTAL=TOTAL+H1 70 D1(I)=H1 DO 80 I=1,K1 H1=D1(I) DO 80 J=1,K2 H2=SQRT(H1*D2(J)) IF (H2.EQ.0.) GOTO 80 Z2(J,I)=Z2(J,I)/H2-H2/TOTAL 80 CONTINUE CALL SIVAD(K2,K1,XLABDA,Z2,K1*K2,Z1,K1*K1,HULP1,ICV) 90 WRITE (6,5000) ILIST J=1 100 K=J+19 IF (K.GT.K2) K=K2 WRITE (6,6700) DO 110 I=1,ILIST 110 WRITE (6,5300) (INTF(I,J1),J1=J,K) J=K+1 IF (J.LE.K2) GOTO 100 IF (ICV.EQ.0.OR.TOTAL.LE.0.) GOTO 790 C C C DETERMINE THE NUMBER OF NONZERO SINGULAR VALUES C C C EPS=SQRT(FLOAT(K1*K2))*1.E-7 DO 120 J=1,MIN 120 IF (XLABDA(J).LT.EPS) GOTO 130 GOTO 140 130 MIN=J-1 140 IF (MIN.EQ.0) GOTO 780 IF (IE1.GT.MIN) GOTO 770 IF (IE2.GT.MIN) IE2=MIN IF (IE1.EQ.IE2) IPL2DI=0 NALDIM=IE2-IE1+1 C C C STANDARDIZE Z1 AND Z2 C C C DO 160 I=1,K1 H1=D1(I)/TOTAL D1(I)=H1 IF (H1.EQ.0.) GOTO 160 H1=SQRT(H1) DO 150 J=1,MIN 150 Z1(I,J)=Z1(I,J)/H1 160 CONTINUE DO 180 I=1,K2 H1=D2(I)/TOTAL D2(I)=H1 IF (H1.EQ.0.) GOTO 180 H1=SQRT(H1) DO 170 J=1,MIN 170 Z2(I,J)=Z2(I,J)/H1 180 CONTINUE IF (.NOT.CONF) GOTO 220 C C C COMPUTE (CO)VARIANCES C C C MINKW=NALDIM*(NALDIM+1)/2 INDVL=1 INDVZ1=INDVL+MINKW INDVZ2=INDVZ1+K1*MINKW INDDXL=INDVZ2+K2*MINKW INDDZ1=INDDXL+NALDIM INDDZ2=INDDZ1+K1*NALDIM DO 190 I=1,K1 DO 190 J=1,K2 190 F(I,J)=FLOAT(INTF(I,J))/TOTAL CALL C82VAR(K1,K2,MIN,IE1,NALDIM,MINKW,F,Z1,Z2,D1,D2,XLABDA,ALPHA, 1 BETA,VARS(INDVL),VARS(INDVZ1),VARS(INDVZ2), 2 VARS(INDDXL),VARS(INDDZ1),VARS(INDDZ2)) J=INDDXL-1 DO 200 I=1,J 200 VARS(I)=VARS(I)/TOTAL IF (IPERM.LT.IE1) GOTO 220 DO 210 I=1,K1 DO 210 J=1,K2 210 INTF(I,J)=F(I,J)*TOTAL+.1 220 IF (ALPHA.EQ.0.) GOTO 240 DO 230 J=IE1,IE2 H1=XLABDA(J)**ALPHA DO 230 I=1,K1 230 Z1(I,J)=Z1(I,J)*H1 240 IF (BETA.EQ.0.) GOTO 260 DO 250 J=IE1,IE2 H1=XLABDA(J)**BETA DO 250 I=1,K2 250 Z2(I,J)=Z2(I,J)*H1 260 CONTINUE C C C PRINT RESULTS C C C I=TOTAL+.5 WRITE (6,1100) I WRITE (6,1200) XLABDA(1) DO 270 I=2,MIN 270 WRITE (6,1300) XLABDA(I) WRITE (6,5100) WRITE (6,5150) DO 280 I=1,K1 IDEN1(I)=I 280 ID1(I)=TOTAL*D1(I)+.5 J=IE1 290 K=J+9 IF (K.GT.IE2) K=IE2 WRITE (6,6700) DO 300 I=1,K1 300 WRITE (6,1400) I,ID1(I),(Z1(I,J1),J1=J,K) J=K+1 IF (J.LE.IE2) GOTO 290 WRITE (6,5200) WRITE (6,5150) DO 310 I=1,K2 IDEN2(I)=I 310 ID2(I)=TOTAL*D2(I)+.5 J=IE1 320 K=J+9 IF (K.GT.IE2) K=IE2 WRITE (6,6700) DO 330 I=1,K2 330 WRITE (6,1400) I,ID2(I),(Z2(I,J1),J1=J,K) J=K+1 IF (J.LE.IE2) GOTO 320 IF (MEDSCA.LE.0) GOTO 360 WRITE(MEDSCA,1700) (XLABDA(J),J=IE1,IE2) DO 340 K=1,K1 340 WRITE (MEDSCA,1700) (Z1(K,J),J=IE1,IE2) DO 350 K=1,K2 350 WRITE (MEDSCA,1700) (Z2(K,J),J=IE1,IE2) 360 IF (.NOT.CONF) GOTO 550 IF (MEDVAR.LE.0) GOTO 390 K=1 L=MINKW WRITE (MEDVAR,1800) (VARS(J),J=K,L) DO 370 I=1,K1 K=K+MINKW L=L+MINKW 370 WRITE(MEDVAR,1800) (VARS(J1),J1=K,L) DO 380 I=1,K2 K=K+MINKW L=L+MINKW 380 WRITE(MEDVAR,1800) (VARS(J1),J1=K,L) C C C PRINT VARIANCES C C C 390 IF (NALDIM.GE.2) WRITE(6,5500) IF (NALDIM.EQ.1) WRITE(6,5550) IE1 N=0 DO 400 I=1,NALDIM N=N+I 400 HULP1(I)=VARS(N) WRITE(6,5600) (HULP1(I),I=1,NALDIM) IF (NALDIM.EQ.1) GOTO 440 IF (NALDIM.GT.2) WRITE (6,5700) IF (NALDIM.EQ.2) WRITE (6,5750) K=NALDIM*(NALDIM+1)/2 CALL SUBTRI(NALDIM,VARS,HULP1,K) L=1 440 DO 490 I=1,K1 IF (NALDIM.GE.2) WRITE(6,5900) I IF (NALDIM.EQ.1) WRITE(6,5950) I DO 450 J=1,NALDIM N=N+J 450 HULP1(J)=VARS(N) WRITE(6,5600) (HULP1(J),J=1,NALDIM) IF (NALDIM.EQ.1) GOTO 490 IF (NALDIM.GT.2) WRITE (6,6000) I IF (NALDIM.EQ.2) WRITE (6,6050) I L=L+K CALL SUBTRI(NALDIM,VARS(L),HULP1,K) 490 CONTINUE DO 540 I=1,K2 IF (NALDIM.GT.1) WRITE(6,6100) I IF (NALDIM.EQ.1) WRITE(6,6150) I DO 500 J=1,NALDIM N=N+J 500 HULP1(J)=VARS(N) WRITE(6,5600) (HULP1(J),J=1,NALDIM) IF (NALDIM.EQ.1) GOTO 540 IF (NALDIM.GT.2) WRITE (6,6200) I IF (NALDIM.EQ.2) WRITE (6,6250) I L=L+K CALL SUBTRI(NALDIM,VARS(L),HULP1,K) 540 CONTINUE C C C PLOT RESULTS C C C 550 IF (IPLTRA.LE.0.OR.IPLTRA.GT.4) GOTO 620 GOTO (560,580,560,600),IPLTRA 560 DO 570 I=IE1,IE2 WRITE (6,3100) I 570 CALL PLOT11(K1,Z1(1,I),K1,IDEN1,IHULP1,-1) IF (IPLTRA.EQ.1) GOTO 620 580 DO 590 I=IE1,IE2 WRITE (6,3300) I 590 CALL PLOT11(K2,Z2(1,I),K2,IDEN2,IHULP2,-1) GOTO 620 600 DO 610 I=IE1,IE2 WRITE (6,3500) I CALL PLOT12(K1,Z1(1,I),K1,K2,Z2(1,I),K2,IDEN1,IDEN2,IHULP1, 1 IHULP2,-3) 610 WRITE (6,3550) 620 IF (IPL2DI.LE.0.OR.IPL2DI.GT.4) GOTO 690 GOTO (630,650,630,670),IPL2DI 630 I2=IE2-1 DO 640 I=IE1,I2 J1=I+1 DO 640 J=J1,IE2 WRITE (6,3200) I,J 640 CALL PLOT21(Z1(1,I),Z1(1,J),K1,IDEN1,IHULP1,-1) IF (IPL2DI.EQ.1) GOTO 690 650 I1=IE2-1 DO 660 I=IE1,I1 J1=I+1 DO 660 J=J1,IE2 WRITE (6,3400) I,J 660 CALL PLOT21(Z2(1,I),Z2(1,J),K2,IDEN2,IHULP2,-1) GOTO 690 670 I1=IE2-1 DO 680 I=IE1,I1 J1=I+1 DO 680 J=J1,IE2 WRITE (6,3600) I,J CALL PLOT22(Z1(1,I),Z1(1,J),K1,Z2(1,I),Z2(1,J),K2,IDEN1, 1 IDEN2,IHULP1,IHULP2,-1) 680 WRITE (6,3550) 690 IF (IPLEXT.NE.1) GOTO 710 DO 700 I=IE1,IE2 WRITE (6,3700) I 700 CALL PLOT21(Z1(1,I),Z2(1,I),K1,IDEN1,IHULP1,-1) 710 IF (IPERM.LT.IE1) GOTO 760 DO 750 I=IE1,IPERM DO 715 J=1,K1 715 Z1(J,I)=-Z1(J,I) DO 720 J=1,K2 720 Z2(J,I)=-Z2(J,I) CALL ORDER (Z1(1,I),K1,IHULP1) CALL ORDER (Z2(1,I),K2,IHULP2) IF (MEDPER.LE.0) GOTO 728 WRITE (MEDPER,1900) (IHULP1(J1),J1=1,K1) WRITE (MEDPER,1900) (ID1(IHULP1(J1)),J1=1,K1) WRITE (MEDPER,1900) (IHULP2(J1),J1=1,K2) WRITE (MEDPER,1900) (ID2(IHULP2(J1)),J1=1,K2) DO 722 J1=1,K1 I1=IHULP1(J1) 722 WRITE (MEDPER,1900) (INTF(I1,IHULP2(L)),L=1,K2) 728 WRITE (6,6400) I J=1 730 K=J+19 IF (K.GT.K2) K=K2 WRITE (6,6450) (IHULP2(J1),J1=J,K) WRITE (6,6500) (ID2(IHULP2(J1)),J1=J,K) WRITE (6,6700) DO 740 J1=1,K1 I1=IHULP1(J1) 740 WRITE (6,6600) I1,ID1(I1),(INTF(I1,IHULP2(L)),L=J,K) J=K+1 750 IF (J.LE.K2) GOTO 730 760 CONTINUE RETURN 770 WRITE(6,5400) RETURN 780 WRITE (6,1600) RETURN 790 WRITE (6,1500) C C C FORMATS C C C 1100 FORMAT(1H1/1H0,2X,28HTOTAL NUMBER OF OBSERVATIONS,2X,I10) 1200 FORMAT(26H0 NONZERO SINGULAR VALUES,9X,F8.4) 1300 FORMAT(35X,F8.4) 1400 FORMAT(1H0,I5,3X,I10,10F11.4) 1500 FORMAT(51H0PROGRAM STOPPED BECAUSE THE INPUT-MATRIX CANNOT BE, 1 9H ANALYZED) 1600 FORMAT(26H0NO SCORES CAN BE ATTACHED/ 1 33H0ROWS AND COLUMNS ARE INDEPENDANT) 1700 FORMAT(8F10.6) 1800 FORMAT(8E10.3) 1900 FORMAT(16I5) 3100 FORMAT(24H1ROW-SCORES IN DIMENSION,I5) 3200 FORMAT(25H1ROW-SCORES IN DIMENSIONS,I5, 1 13H (X-AXIS) AND,I5,9H (Y-AXIS)) 3300 FORMAT(27H1COLUMN-SCORES IN DIMENSION,I5) 3400 FORMAT(28H1COLUMN-SCORES IN DIMENSIONS,I5, 1 13H (X-AXIS) AND,I5,9H (Y-AXIS)) 3500 FORMAT(36H1ROW- AND COLUMN-SCORES IN DIMENSION,I5) 3550 FORMAT(33H0THE ROWS ARE LABELLED BY DIGITS,, 1 23H THE COLUMNS BY LETTERS) 3600 FORMAT(37H1ROW- AND COLUMN-SCORES IN DIMENSIONS,I5, 1 13H (X-AXIS) AND,I5,9H (Y-AXIS)) 3700 FORMAT(47H1ROWSCORES (X-AXIS) VS COLUMNSCORES (Y-AXIS) IN, 1 10H DIMENSION,I5) 5000 FORMAT(1H1/1H ,2X,9HTHE FIRST,I5,24H ROWS OF THE DATA MATRIX/) 5100 FORMAT(1H1/1H ,2X,10HROW SCORES/1H0,2X,3HROW,5X,8HMARGINAL) 5150 FORMAT(8H NUMBER,3X,9HFREQUENCY,3X,6HSCORES) 5200 FORMAT(1H1/1H ,2X,13HCOLUMN SCORES/8H0 COLUMN,3X,8HMARGINAL) 5300 FORMAT(2X,25I5) 5400 FORMAT(32H0NO SOLUTION IN GIVEN DIMENSIONS) 5500 FORMAT(1H1/35H VARIANCES OF THE SINGULAR VALUES//) 5550 FORMAT(1H1/36H VARIANCE OF SINGULAR VALUE NUMBER,I5//) 5600 FORMAT(2X,10(2X,E8.2)) 5700 FORMAT(1H0/44H CORRELATION MATRIX OF THE SINGULAR VALUES//) 5750 FORMAT(1H0/37H CORRELATION OF THE SINGULAR VALUES//) 5900 FORMAT(1H0/3H , 1 37HVARIANCES OF THE SCORES OF ROW NUMBER,I5//) 5950 FORMAT(1H0/3H , 1 35HVARIANCE OF THE SCORE OF ROW NUMBER,I5//) 6000 FORMAT(1H0/1H ,2X, 1 46HCORRELATION MATRIX OF THE SCORES OF ROW NUMBER,I5//) 6050 FORMAT(1H0/1H ,2X, 1 39HCORRELATION OF THE SCORES OF ROW NUMBER,I5//) 6100 FORMAT(1H0/1H ,2X, 1 40HVARIANCES OF THE SCORES OF COLUMN NUMBER,I5//) 6150 FORMAT(1H0/1H ,2X, 1 38HVARIANCE OF THE SCORE OF COLUMN NUMBER,I5//) 6200 FORMAT(1H0/1H ,2X, 1 49HCORRELATION MATRIX OF THE SCORES OF COLUMN NUMBER,I5//) 6250 FORMAT(1H0/1H ,2X, 1 42HCORRELATION OF THE SCORES OF COLUMN NUMBER,I5//) 6400 FORMAT(1H1/1H ,2X,42HTHE DATA-MATRIX PERMUTED ACCORDING TO THE , 1 19HSCORES IN DIMENSION,I5//) 6450 FORMAT(15H0NUMBER COLUMN,20I5) 6500 FORMAT(15H0 ROW MARGINAL,20I5) 6600 FORMAT(1H ,I5,I7,2X,20I5) 6700 FORMAT(1H ) RETURN END SUBROUTINE SUBTRI(N,A,B,N2) DIMENSION A(N2),B(N) C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C THE ELEMENTS OF AN N N SYMMETRIC MATRIX ARE SUPPOSED TO BE C C REPRESENTED BY A OF SIZE N2=N (N+1)/2 IN THE ORDER: C C (1,1),(1,2),(2,2),(1,3),(2,3),(3,3),(1,4),.....,(N,N). C C THE ELEMENTS OF B ARE SUPPOSED TO BE POSITIVE. THIS SUBROUTINE C C RETURNS THE OFF-DIAGONAL ELEMENTS OF THE MATRIX C C (B**-.5) A (B**-.5) IN THE CORRESPONDING ELEMENTS OF A. C C IF B EQUALS THE DIAGONAL ELEMENTS OF A, AND A IS A COVARIANCE- C C MATRIX, THEN A WILL BE CHANGED INTO A CORRELATION-MATRIX. C C MOREOVER, A IS PRINTED IN THIS SUBROUTINE. C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C SR=SQRT(B(1)) L=2 K=2 I=2 J=2 10 A(I)=A(I)/SR I=I+J J=J+1 IF (J.LE.N) GOTO 10 SR=SQRT(B(L)) J=L-1 DO 20 I=1,J A(K)=A(K)/SR 20 K=K+1 I=K+L K=K+1 L=L+1 J=L IF (J.LE.N) GOTO 10 L=0 K=2 IR=0 I=2 30 L=L+I IF (L-K.GE.10) L=K+9 WRITE (6,1000) (A(J),J=K,L) K=K+I I=I+1 IF (I.LE.N) GOTO 30 WRITE (6,2000) IR=IR+10 IF (I-IR.LT.3) RETURN I=IR+2 K=I*(I+1)/2-1 L=K-I GOTO 30 1000 FORMAT(2X,10F10.4) 2000 FORMAT(1H ) END SUBROUTINE C82VAR(K1,K2,MIN,IE1,NALDIM,MINKW,F,Z1,Z2,D1,D2,XLABDA, 1 ALPHA,BETA,VARXL,VARZ1,VARZ2,DERXL,DERZ1,DERZ2) DIMENSION F(K1,K2),Z1(K1,MIN),Z2(K2,MIN),D1(K1),D2(K2), 1 XLABDA(MIN),VARXL(MINKW),VARZ1(MINKW,K1), 2 VARZ2(MINKW,K2),DERXL(NALDIM),DERZ1(K1,NALDIM), 3 DERZ2(K2,NALDIM) IE2=IE1+NALDIM-1 DO 10 I=1,MINKW 10 VARXL(I)=0. I=IE1 J=IE1 L=1 40 DO 20 K=1,K1 20 VARZ1(L,K)=-Z1(K,I)*Z1(K,J)/4. DO 30 K=1,K2 30 VARZ2(L,K)=-Z2(K,I)*Z2(K,J)/4. L=L+1 J=J+1 IF (J.LE.I) GOTO 40 I=I+1 J=IE1 IF (I.LE.IE2) GOTO 40 DO 50 I=1,K1 DO 50 J=1,K2 IF (F(I,J).LE.0.) GOTO 50 L=0 DO 60 K=IE1,IE2 L=L+1 60 DERXL(L)=Z1(I,K)*Z2(J,K)- 1 (Z1(I,K)*Z1(I,K)+Z2(J,K)*Z2(J,K))*XLABDA(K)/2. CALL C82DER(K1,K2,MIN,IE1,NALDIM,Z1,Z2,ALPHA,DERZ1,D1, 1 XLABDA,DERXL,I,J) CALL C82DER(K2,K1,MIN,IE1,NALDIM,Z2,Z1,BETA,DERZ2,D2,XLABDA, 1 DERXL,J,I) I1=1 I2=1 L=1 70 VARXL(L)=VARXL(L)+F(I,J)*DERXL(I1)*DERXL(I2) DO 80 K=1,K1 80 VARZ1(L,K)=VARZ1(L,K)+F(I,J)*DERZ1(K,I1)*DERZ1(K,I2) DO 90 K=1,K2 90 VARZ2(L,K)=VARZ2(L,K)+F(I,J)*DERZ2(K,I1)*DERZ2(K,I2) L=L+1 I2=I2+1 IF (I2.LE.I1) GOTO 70 I1=I1+1 I2=1 IF (I1.LE.NALDIM) GOTO 70 50 CONTINUE RETURN END SUBROUTINE C82DER(N,M,MIN,IE1,NALDIM,Z1,Z2,ALPHA,DER,D,XLABDA, 1 DERXL,I,J) DIMENSION Z1(N,MIN),Z2(M,MIN),DER(N,NALDIM),D(N),XLABDA(MIN), 1 DERXL(NALDIM) C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C IF Z1 AND Z2 ARE THE ROW- AND COLUMN-QUANTIFICATIONS OF AN C C N M MATRIX F, WITH D THE ROW-MARGINALS AND XLABDA THE SINGULAR C C VALUES, THIS SUBROUTINE RETURNS IN DER THE DERIVATIVES OF Z1 C C TO THE ELEMENT (I,J) OF F. C C C C SUBROUTINES CALLED: NONE. C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C DO 10 K=1,N H1=-1. DO 20 K1=1,MIN 20 H1=H1-Z1(I,K1)*Z1(K,K1) IF (K.EQ.I) H1=H1+1./D(I) L=IE1 DO 10 K1=1,NALDIM DER(K,K1)=H1*Z2(J,L)/XLABDA(L)+ 1 Z1(I,L)*Z1(I,L)*Z1(K,L)/2.+ 2 ALPHA*Z1(K,L)*DERXL(K1)/XLABDA(L) 10 L=L+1 L=IE1 DO 30 K1=1,NALDIM DER(I,K1)=DER(I,K1)-Z1(I,L)/D(I) 30 L=L+1 K1=IE1+NALDIM IF (K1.GT.MIN) K1=MIN K=0 I1=1 40 IF (I1.GE.IE1) K=K+1 H1=XLABDA(I1)*Z2(J,I1) H6=XLABDA(I1)*XLABDA(I1) J1=IE1 IF (I1.GE.IE1) J1=I1+1 L=J1-IE1+1 50 H7=XLABDA(J1)*XLABDA(J1) H2=XLABDA(J1)*Z2(J,J1) H2=H2*Z1(I,I1)+H1*Z1(I,J1)-H1*H2 H3=Z1(I,I1)*Z1(I,J1) H4=H6-H7 IF (I1.LT.IE1) GOTO 60 H5=(H2-H7*H3)/H4 DO 70 K2=1,N 70 DER(K2,K)=DER(K2,K)+H5*Z1(K2,J1) 60 IF (L.GT.NALDIM) GOTO 80 H5=(H2-H6*H3)/H4 DO 90 K2=1,N 90 DER(K2,L)=DER(K2,L)-H5*Z1(K2,I1) 80 L=L+1 J1=J1+1 IF (J1.LE.MIN) GOTO 50 IF (K.EQ.0) GOTO 100 H2=XLABDA(I1)**ALPHA DO 110 K2=1,N 110 DER(K2,K)=DER(K2,K)*H2 100 I1=I1+1 IF (I1.LT.K1) GOTO 40 IF (IE1+NALDIM.LE.MIN) RETURN K=K+1 H2=XLABDA(MIN)**ALPHA DO 120 K2=1,N 120 DER(K2,K)=DER(K2,K)*H2 RETURN END SUBROUTINE PLOT11(N1,YR1,NR1,IDENR1,INDR1,ILABEL) DIMENSION YR1(NR1),IDENR1(NR1),INDR1(NR1) DIMENSION LINE(67),XVA(12),NVA(12) C>>>> DIMENSION LINE(NC),XVA,NVA(NN+1) EQUIVALENCE (XVA(1),NVA(1)) C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C THIS SUBROUTINE GIVES A ONE-DIMENSIONAL PLOT OF THE NR1 POINTS C C YR1 THAT ARE REGARDED AS REPLICATIONS OF N1 POINTS. C C THE POINTS ARE LABELLED BY 1,...,9,0 (ILABEL=+/-1) OR A,...,J C C (ILABEL=+/-2), CORRESPONDING WITH THE CONTENT OF IDENR1. C C COINCIDING POINTS ARE LABELLED BY 'M' (ILABEL.LT.0) OR THE SUM C C OF THE ORIGINAL LABELS (ILABEL.GT.0). C C C C THE ARRAY INDR1 IS LOCAL. C C C C SUBROUTINES CALLED: ORDER,WRILIN,SUBLIN,STEPX1. C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C THE VALUES OF YR1 ARE SORTED IN DESCENDING ORDER C C C CALL ORDER(YR1,NR1,INDR1) YMIN=YR1(INDR1(NR1)) YMAX=YR1(INDR1(1)) C C THE PLOT IS ADJUSTED TO THE LENGTH OF THE AXIS C C SPANY=YMAX-YMIN I=N1-1 SPANX=I STEPX=STEPX1(I) STEPY=SPANY/53. C C THE POINTS ARE PLOTTED LINE AFTER LINE C WRITE (6,1) IR1=1 NEXTL=1 L=1 YPR=YMAX 45 IF (L.EQ.NEXTL) GOTO 50 WRITE (6,3) YPR GOTO 100 50 DO 55 J=1,67 55 LINE(J)=0 60 K1=INDR1(IR1) NC1=(FLOAT(MOD(K1-1,N1)+1)-1.)/STEPX+1.5 CALL SUBLIN(LINE(NC1),IDENR1(K1),ILABEL) IR1=IR1+1 NEXTL=55 IF (IR1.LE.NR1) NEXTL=(YMAX-YR1(INDR1(IR1)))/STEPY+1.5 IF (L.EQ.NEXTL) GOTO 60 IF (IABS(ILABEL).EQ.2) GOTO 70 DO 65 J=1,67 K1=LINE(J) 65 IF (K1.LT.0) LINE(J)=-K1 70 CALL WRILIN(YPR,LINE) 100 L=L+1 YPR=YPR-STEPY IF (L.LE.54) GOTO 45 C C VALUES OF X-AXIS ARE PRINTED C XVA(1)=1. DO 130 J=1,11 130 XVA(J+1)=XVA(J)+STEPX*6.0 WRITE (6,1) IF (SPANX.LE.66.) GOTO 150 WRITE (6,4) XVA RETURN 150 H1=66.1/66. DO 160 J=1,12 160 NVA(J)=XVA(J)*H1 WRITE (6,5) NVA 1 FORMAT(32X,3H.--,11(6H+-----),4H+--.) 3 FORMAT(25X,F6.2,1X,1H!,71X,1H!) 4 FORMAT(32X,12F6.2) 5 FORMAT(30X,12I6) RETURN END SUBROUTINE PLOT12(N1,YR1,NR1,N2,YR2,NR2,IDENR1,IDENR2,INDR1, 1 INDR2,ILABEL) DIMENSION YR1(NR1),IDENR1(NR1),INDR1(NR1),YR2(NR2),IDENR2(NR2), 1 INDR2(NR2) DIMENSION LINE(67),XVA(12),NVA(12) C>>>> DIMENSION LINE(NC),XVA,NVA(NN+1) EQUIVALENCE (XVA(1),NVA(1)) C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C THIS SUBROUTINE GIVES A ONE-DIMENSIONAL PLOT OF THE NR1 POINTS C C Y1 THAT ARE REGARDED AS REPLICATIONS OF N1 POINTS, AND THE NR2 C C POINTS YR2 THAT ARE REGARDED AS REPLICATIONS OF N2 POINTS. C C THE POINTS ARE LABELLED BY 1,...,9,0 (ILABEL=+/-1); A,...,J C C (ILABEL=+/-2); OR THE FIRST SET BY LETTERS AND THE SECOND BY C C DIGITS (ILABEL=+/-3), CORRESPONDING WITH THE CONTENT OF IDENR1 C C AND IDENR2. C C COINCIDING POINTS ARE LABELLED BY 'M' (ILABEL.LT.0 OR THE POINTS C C ARE IN DIFFERENT SETS) OR THE SUM OF THEIR ORIGINAL LABELS C C (ILABEL.GT.0 AND THE POINTS ARE IN THE SAME SET). C C C C THE ARRAYS INDR1 AND INDR2 ARE LOCAL. C C C C SUBROUTINES CALLED: ORDER,WRILIN,ADDLIN,SUBLIN,STEPX1. C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C THE VALUES OF YR1 AND YR2 ARE SORTED IN DESCENDING ORDER C C C CALL ORDER(YR1,NR1,INDR1) YMIN=YR1(INDR1(NR1)) YMAX=YR1(INDR1(1)) CALL ORDER(YR2,NR2,INDR2) YM=YR2(INDR2(NR2)) IF (YM.LT.YMIN) YMIN=YM YM=YR2(INDR2(1)) IF (YM.GT.YMAX) YMAX=YM C C THE PLOT IS ADJUSTED TO THE LENGTH OF THE AXIS C C SPANY=YMAX-YMIN I=N1-1 IF (N2.GT.N1) I=N2-1 SPANX=I STEPX=STEPX1(I) STEPY=SPANY/53. C C THE POINTS ARE PLOTTED LINE AFTER LINE C WRITE (6,1) IR1=1 IR2=1 NLR1=(YMAX-YR1(INDR1(1)))/STEPY+1.5 NLR2=(YMAX-YR2(INDR2(1)))/STEPY+1.5 NEXTL=1 L=1 YPR=YMAX 45 IF (L.EQ.NEXTL) GOTO 50 WRITE (6,3) YPR GOTO 100 50 DO 55 J=1,67 55 LINE(J)=0 60 IF (L.LT.NLR1) GOTO 75 K1=INDR1(IR1) NC1=(FLOAT(MOD(K1-1,N1)+1)-1.)/STEPX+1.5 CALL ADDLIN(LINE(NC1),IDENR1(K1),ILABEL) IR1=IR1+1 NLR1=55 IF (IR1.LE.NR1) NLR1=(YMAX-YR1(INDR1(IR1)))/STEPY+1.5 GOTO 60 75 IF (L.LT.NLR2) GOTO 80 K1=INDR2(IR2) NC1=(FLOAT(MOD(K1-1,N2)+1)-1.)/STEPX+1.5 CALL SUBLIN(LINE(NC1),IDENR2(K1),ILABEL) IR2=IR2+1 NLR2=55 IF (IR2.LE.NR2) NLR2=(YMAX-YR2(INDR2(IR2)))/STEPY+1.5 GOTO 75 80 IF (IABS(ILABEL).NE.1) GOTO 85 DO 65 J=1,67 K1=LINE(J) 65 IF (K.LT.0) LINE(J)=-K 85 IF (IABS(ILABEL).NE.2) GOTO 70 DO 90 J=1,67 K1=LINE(J) 90 IF (K.GT.0.AND.K.NE.11) LINE(J)=-K 70 CALL WRILIN(YPR,LINE) NEXTL=NLR1 IF (NLR1.GT.NLR2) NEXTL=NLR2 100 L=L+1 YPR=YPR-STEPY IF (L.LE.54) GOTO 45 C C VALUES OF X-AXIS ARE PRINTED C XVA(1)=1. DO 130 J=1,11 130 XVA(J+1)=XVA(J)+STEPX*6.0 WRITE (6,1) IF (SPANX.LE.66.) GOTO 150 WRITE (6,4) XVA RETURN 150 H1=66.1/66. DO 160 J=1,12 160 NVA(J)=XVA(J)*H1 WRITE (6,5) NVA 1 FORMAT(32X,3H.--,11(6H+-----),4H+--.) 3 FORMAT(25X,F6.2,1X,1H!,71X,1H!) 4 FORMAT(32X,12F6.2) 5 FORMAT(30X,12I6) RETURN END SUBROUTINE PLOT21(X1,Y1,N1,IDENY1,INDY1,ILABEL) DIMENSION X1(N1),Y1(N1),IDENY1(N1),INDY1(N1) DIMENSION LINE(67),XVA(12) C>>>> DIMENSION LINE(NC),XVA(NN+1) C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C THIS SUBROUTINE GIVES A TWO-DIMENSIONAL PLOT OF THE N1 POINTS C C (X1,Y1), LABELLED ACCORDING TO THE CONTENTS OF IDENY1. C C THE POINTS ARE LABELLED WITH 1,...,9,0 (ABS(ILABEL).LE.2) OR C C WITH A,B,...,J (ABS(ILABEL).GE.3). C C THE AXES ARE SCALED EQUALLY (ILABEL IS ODD) OR INDEPENDENTLY C C (ILABEL IS EVEN). C C COINCIDING POINTS ARE LABELLED BY 'M' (ILABEL IS NEGATIVE) OR C C BY THE SUM OF THE CONCERNING LABELS (ILABEL IS POSITIVE). C C C C THE ARRAY INDY1 IS LOCAL. C C C C SUBROUTINES CALLED: ORDER,WRILIN,SUBLIN. C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C THE VALUES OF Y1 ARE SORTED IN DESCENDING ORDER. C C C CALL ORDER(Y1,N1,INDY1) YMIN=AMIN1(0.,Y1(INDY1(N1))) YMAX=AMAX1(0.,Y1(INDY1(1))) C C THE PLOT IS ADJUSTED TO THE LENGTH OF THE AXIS C C SPANY=YMAX-YMIN XMIN=0. XMAX=0. DO 10 I=1,N1 H1=X1(I) IF (H1.LT.XMIN) XMIN=H1 10 IF (H1.GT.XMAX) XMAX=H1 SPANX=XMAX-XMIN IF (MOD(ILABEL,2).EQ.0) GOTO 20 IF (SPANY.GT.SPANX) SPANX=SPANY SPANY=SPANX XMIN=(XMIN+XMAX-SPANX)/2. YMAX=(YMIN+YMAX+SPANY)/2. 20 STEPX=SPANX/66. STEPY=SPANY/53. C C THE POINTS ARE PLOTTED LINE AFTER LINE C WRITE (6,1) I1=1 NEXTL=(YMAX-Y1(INDY1(1)))/STEPY+1.5 L=1 YPR=YMAX 45 IF (L.EQ.NEXTL) GOTO 50 WRITE (6,3) YPR GOTO 100 50 DO 55 J=1,67 55 LINE(J)=0 60 K1=INDY1(I1) NC1=(X1(K1)-XMIN)/STEPX+1.5 CALL SUBLIN(LINE(NC1),IDENY1(K1),ILABEL) I1=I1+1 NEXTL=55 IF (I1.LE.N1) NEXTL=(YMAX-Y1(INDY1(I1)))/STEPY+1.5 IF (L.EQ.NEXTL) GOTO 60 IF (IABS(ILABEL).GE.3) GOTO 70 DO 65 J=1,67 K1=LINE(J) 65 IF (K1.LT.0) LINE(J)=-K1 70 CALL WRILIN(YPR,LINE) 100 L=L+1 YPR=YPR-STEPY IF (L.LE.54) GOTO 45 C C VALUES OF X-AXIS ARE PRINTED C XVA(1)=XMIN DO 130 J=1,11 130 XVA(J+1)=XVA(J)+STEPX*6.0 WRITE (6,1) WRITE (6,4) XVA 1 FORMAT(32X,3H.--,11(6H+-----),4H+--.) 3 FORMAT(25X,F6.2,1X,1H!,71X,1H!) 4 FORMAT(32X,12F6.2) RETURN END SUBROUTINE PLOT22(X1,Y1,N1,X2,Y2,N2,IDENY1,IDENY2,INDY1,INDY2, 1 ILABEL) DIMENSION X1(N1),Y1(N1),IDENY1(N1),INDY1(N1),X2(N2),Y2(N2), 1 IDENY2(N2),INDY2(N2) DIMENSION LINE(67),XVA(12) C>>>> DIMENSION LINE(NC),XVA(NN+1) C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C THIS SUBROUTINE GIVES A TWO-DIMENSIONAL PLOT OF THE N1 POINTS C C (X1,Y1) AND THE N2 POINTS (X2,Y2). C C THE (X1,Y1)-POINTS ARE LABELLED BY 1,...,9,0; THE (X2,Y2)-POINTS C C BY A,...,J, CORRESPONDING WITH THE CONTENTS OF IDENY1 AND C C IDENY2. C C COINCIDING POINTS ARE LABELLED BY 'M' (ILABEL.LT.0 OR THE C C POINTS HAVE A DIFFERENT TYPE OF LABEL) OR THE SUM OF THEIR C C ORIGINAL LABELS (ILABEL.GT.0 AND THE POINTS HAVE THE SAME KIND C C OF LABEL). C C THE AXES ARE SCALED EQUALLY (ILABEL=+/-1) OR INDEPENDENTLY C C (ILABEL=+/-2). C C C C THE ARRAYS INDY1 AND INDY2 ARE LOCAL. C C C C SUBROUTINES CALLED: ORDER,WRILIN,ADDLIN,SUBLIN. C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C THE VALUES OF Y1 AND Y2 ARE SORTED IN DESCENDING ORDER. C C C CALL ORDER(Y1,N1,INDY1) CALL ORDER(Y2,N2,INDY2) YMIN=AMIN1(0.,Y1(INDY1(N1)),Y2(INDY2(N2))) YMAX=AMAX1(0.,Y1(INDY1(1)),Y2(INDY2(1))) C C THE PLOT IS ADJUSTED TO THE LENGTH OF THE AXIS C C SPANY=YMAX-YMIN XMIN=0. XMAX=0. DO 10 I=1,N1 H1=X1(I) IF (H1.LT.XMIN) XMIN=H1 10 IF (H1.GT.XMAX) XMAX=H1 DO 15 I=1,N2 H1=X2(I) IF (H1.LT.XMIN) XMIN=H1 15 IF (H1.GT.XMAX) XMAX=H1 SPANX=XMAX-XMIN IF (IABS(ILABEL).EQ.2) GOTO 20 IF (SPANY.GT.SPANX) SPANX=SPANY SPANY=SPANX XMIN=(XMIN+XMAX-SPANX)/2. YMAX=(YMIN+YMAX+SPANY)/2. 20 STEPX=SPANX/66. STEPY=SPANY/53. C C THE POINTS ARE PLOTTED LINE AFTER LINE C WRITE (6,1) I1=1 I2=1 NL1=(YMAX-Y1(INDY1(1)))/STEPY+1.5 NL2=(YMAX-Y2(INDY2(1)))/STEPY+1.5 L=1 YPR=YMAX 45 NEXTL=NL1 IF (NL2.LT.NL1) NEXTL=NL2 IF (L.EQ.NEXTL) GOTO 50 WRITE (6,3) YPR GOTO 100 50 DO 55 J=1,67 55 LINE(J)=0 60 IF (NEXTL.LT.NL1) GOTO 80 65 K1=INDY1(I1) NC1=(X1(K1)-XMIN)/STEPX+1.5 CALL ADDLIN(LINE(NC1),IDENY1(K1),ILABEL) I1=I1+1 NL1=55 IF (I1.LE.N1) NL1=(YMAX-Y1(INDY1(I1)))/STEPY+1.5 IF (L.EQ.NL1) GOTO 65 80 IF (NEXTL.LT.NL2) GOTO 90 85 K1=INDY2(I2) NC1=(X2(K1)-XMIN)/STEPX+1.5 CALL SUBLIN(LINE(NC1),IDENY2(K1),ILABEL) I2=I2+1 NL2=55 IF (I2.LE.N2) NL2=(YMAX-Y2(INDY2(I2)))/STEPY+1.5 IF (L.EQ.NL2) GOTO 85 90 CALL WRILIN(YPR,LINE) 100 L=L+1 YPR=YPR-STEPY IF (L.LE.54) GOTO 45 C C VALUES OF X-AXIS ARE PRINTED C XVA(1)=XMIN DO 130 J=1,11 130 XVA(J+1)=XVA(J)+STEPX*6.0 WRITE (6,1) WRITE (6,4) XVA 1 FORMAT(32X,3H.--,11(6H+-----),4H+--.) 3 FORMAT(25X,F6.2,1X,1H!,71X,1H!) 4 FORMAT(32X,12F6.2) RETURN END SUBROUTINE ADDLIN(LINE,ID,ICO) IF (LINE.EQ.11.OR.LINE.LT.0) GOTO 10 IF (LINE.GT.0.AND.ICO.LT.0) GOTO 10 LINE=MOD(LINE+ID+9,10)+1 RETURN 10 LINE=11 RETURN END SUBROUTINE SUBLIN(LINE,ID,ICO) IF (LINE.GT.0) GOTO 10 IF (LINE.LT.0.AND.ICO.LT.0) GOTO 10 LINE=-MOD(ID-LINE+9,10)-1 RETURN 10 LINE=11 RETURN END SUBROUTINE WRILIN(YPR,LINE) DIMENSION LINE(67),LOGLIN(67),LABEL(22) LOGICAL LOGLIN,LABEL DATA LABEL(1),LABEL(2),LABEL(3),LABEL(4),LABEL(5),LABEL(6), 1 LABEL(7),LABEL(8),LABEL(9),LABEL(10),LABEL(11),LABEL(12), 2 LABEL(13),LABEL(14),LABEL(15),LABEL(16),LABEL(17),LABEL(18), 3 LABEL(19),LABEL(20),LABEL(21),LABEL(22)/ 4 1HJ,1HI,1HH,1HG,1HF,1HE,1HD,1HC,1HB,1HA,1H ,1H1,1H2,1H3,1H4, 5 1H5,1H6,1H7,1H8,1H9,1H0,1HM/ DO 10 I=1,67 10 LOGLIN(I)=LABEL(LINE(I)+11) WRITE (6,1) YPR,(LOGLIN(I),I=1,67) 1 FORMAT(25X,F6.2,1X,1H!,2X,67A1,2X,1H!) RETURN END FUNCTION STEPX1(I) STEPX=11/I*6 IF (STEPX.EQ.0.) STEPX=6/(1+(I-1)/11) IF (STEPX.EQ.0.) STEPX=66./FLOAT(I) STEPX1=1./STEPX RETURN END SUBROUTINE ORDER(Y,N,IND) DIMENSION Y(N),IND(N) DO 5 I=1,N 5 IND(I)=I M=N 10 M=M/2 IF (M.LT.1) GOTO 50 K=N-M J=1 20 I=J 30 L=I+M IF (Y(IND(I)).GE.Y(IND(L))) GOTO 40 II=IND(I) IND(I)=IND(L) IND(L)=II I=I-M IF (I.GE.1) GOTO 30 40 J=J+1 IF (J-K) 20,20,10 50 CONTINUE RETURN END SUBROUTINE SIVAD(M,N,Q,U,IU,V,IV,E,ICV) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C**** C C**** NAME: SIVAD(M,N,Q,U,IU,V,IV,E,ICV) C C**** C C**** C C**** PURPOSE : COMPUTES AND ORDERS THE SINGULAR VALUES OF A REAL C C**** MATRIX A AND GIVES A COMPLETE ORTHOGONAL DECOMPO- C C**** SITION: C C**** A = U * Q * V(T) C C**** C C**** C C**** PARAMETERS : M NUMBER OF ROWS OF THE MATRIX A C C**** C C**** N NUMBER OF COLUMNS OF THE MATRIX A C C**** C C**** IU LENGTH OF THE INPUT- AND OUTPUTVECTOR C C**** U,SO IU EQUALS M*N C C**** C C**** IV LENGTH OF THE OUTPUTVECTOR V,SO IV C C**** EQUALS N*N C C**** C C**** ICV OUTPUT PARAMETER,CONTAINS THE MAXIMUM C C**** NUMBER OF ITERATIONS USED,IF ICV = 0 C C**** THERE WAS NO CONVERGENCE FOR ONE OF THE C C**** SINGULAR VALUES AND THE SUBROUTINE C C**** HAS STOPPED C C**** C C**** C C**** DIMENSIONS: U(IU) CONTAINS ON INPUT THE M*N MATRIX TO BE C C**** COMPOSED,ON OUTPUT IT CONTAINS THE M*N C C**** ORTHONORMALIZED MATRIX U C C**** C C**** V(IV) REPRESENTS THE ORTHOGONAL MATRIX V C C**** C C**** Q(N) A VECTOR HOLDING THE NONNEGATIVE C C**** SINGULAR VALUES OF A,IN DECREASING C C**** SEQUENCE C C**** C C**** E(N) A VECTOR USED IN THE SUBROUTINE,IT C C**** SHOULD BE DECLARED IN THE MAIN PROGRAM C C**** C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C**** C**** C**** * * * * * * START OF PROGRAM * * * * * * * * * * * * C**** C**** DIMENSION U(IU),V(IV),Q(N),E(N) DATA ZERO,ONE,TWO/0.0,1.0,2.0/ EPS=2.0**(-24) TOL=8.0**(-54) C**** C**** C**** * * * * * * HOUSEHOLDER'S REDUCTION TO BIDIAGONAL FORM C**** C**** G=ZERO X=ZERO DO 130 I=1,N E(I)=G S=ZERO L=I+1 DO 20 J=I,M II=J+(I-1)*M S=S+U(II)*U(II) 20 CONTINUE IF(S.GE.TOL) GOTO 30 G=ZERO GOTO 60 30 CONTINUE JJ=I+(I-1)*M F=U(JJ) G=-SQRT(S) IF(F.LT.ZERO) G=-G H=F*G-S U(JJ)=F-G IF(L.GT.N) GOTO 60 DO 50 J=L,N S=ZERO DO 40 K=I,M S=S+U(K+(I-1)*M)*U(K+(J-1)*M) 40 CONTINUE F=S/H DO 50 K=I,M KJ=K+(J-1)*M U(KJ)=U(KJ)+F*U(K+(I-1)*M) 50 CONTINUE 60 CONTINUE Q(I)=G S=ZERO IF(L.GT.N) GOTO 75 DO 70 J=L,N IJ=I+(J-1)*M S=S+U(IJ)*U(IJ) 70 CONTINUE 75 CONTINUE IF(S.GE.TOL) GOTO 80 G=ZERO GOTO 120 80 CONTINUE I1=(M+1)*I F=U(I1) G=-SQRT(S) IF(F.LT.ZERO) G=-G H=F*G-S U(I1)=F-G DO 90 J=L,N E(J)=U(I+(J-1)*M)/H 90 CONTINUE DO 110 J=L,M S=ZERO DO 100 K=L,N S=S+U(J+(K-1)*M)*U(I+(K-1)*M) 100 CONTINUE DO 110 K=L,N JK=J+(K-1)*M U(JK)=U(JK)+S*E(K) 110 CONTINUE 120 CONTINUE Y=ABS(Q(I))+ABS(E(I)) IF(Y.GT.X) X=Y 130 CONTINUE C**** C**** C**** * * * * * * ACCUMULATION OF RIGHT HAND TRANSFORMATIONS C**** C**** DO 210 II=1,N I=(N+1)-II IF(G.NE.ZERO) GOTO 150 GOTO 190 150 CONTINUE IPLUS=(M+1)*I H=U(IPLUS)*G DO 160 J=L,N V(J+(I-1)*N)=U(I+(J-1)*M)/H 160 CONTINUE DO 180 J=L,N S=ZERO DO 170 K=L,N S=S+U(I+(K-1)*M)*V(K+(J-1)*N) 170 CONTINUE DO 180 K=L,N KJ=K+(J-1)*N V(KJ)=V(KJ)+S*V(K+(I-1)*N) 180 CONTINUE 190 CONTINUE IF(L.GT.N) GOTO 205 DO 200 J=L,N V(J+(I-1)*N)=ZERO V(I+(J-1)*N)=ZERO 200 CONTINUE 205 CONTINUE V(I+(I-1)*N)=ONE G=E(I) L=I 210 CONTINUE C**** C**** C**** * * * * * * ACCUMULATION OF LEFT HAND TRANSFORMATIONS C**** C**** DO 320 II=1,N I=(N+1)-II JJ=I+(I-1)*M L=I+1 G=Q(I) IF(L.GT.N) GOTO 245 DO 240 J=L,N U(I+(J-1)*M)=ZERO 240 CONTINUE 245 CONTINUE IF(G.NE.ZERO) GOTO 250 GOTO 290 250 CONTINUE H=U(JJ)*G IF(L.GT.N) GOTO 275 DO 270 J=L,N S=ZERO DO 260 K=L,M S=S+U(K+(I-1)*M)*U(K+(J-1)*M) 260 CONTINUE F=S/H DO 270 K=I,M KJ=K+(J-1)*M U(KJ)=U(KJ)+F*U(K+(I-1)*M) 270 CONTINUE 275 CONTINUE DO 280 J=I,M JI=J+(I-1)*M U(JI)=U(JI)/G 280 CONTINUE GOTO 310 290 CONTINUE DO 300 J=I,M U(J+(I-1)*M)=ZERO 300 CONTINUE 310 CONTINUE U(JJ)=U(JJ)+ONE 320 CONTINUE C**** C**** C**** * * * * * * DIAGONALIZATION OF THE BIDIAGONAL FORM * * C**** C**** EPS=EPS*X ICV=0 ICONV=1 DO 470 KK=1,N K=(N+1)-KK IF(ICV.LT.ICONV) ICV=ICONV ICONV=0 340 CONTINUE ICONV=ICONV+1 IF(ICONV.GT.50) GOTO 550 DO 350 LL=1,K L=(K+1)-LL IF(ABS(E(L)).LE.EPS) GOTO 390 LMIN=L-1 IF(ABS(Q(LMIN)).LE.EPS) GOTO 360 350 CONTINUE 360 CONTINUE C**** C**** C**** * * * * * * CANCELLATION OF E(L) IF L>1 * * * * * * * C**** C**** C=ZERO S=ONE L1=L-1 DO 380 I=L,K F=S*E(I) E(I)=C*E(I) IF(ABS(F).LE.EPS) GOTO 390 G=Q(I) IF(ABS(F).LT.ABS(G)) GOTO 362 IF (F.EQ.ZERO) GOTO 365 H=G/F H=ABS(F)*SQRT(ONE+H*H) GOTO 368 362 CONTINUE H=F/G H=ABS(G)*SQRT(ONE+H*H) GOTO 368 365 CONTINUE H=ZERO GOTO 369 368 C=G/H S=-F/H 369 Q(I)=H DO 370 J=1,M JL=J+(L1-1)*M JI=J+(I-1)*M Y=U(JL) Z=U(JI) U(JL)=Y*C+Z*S U(JI)=-Y*S+Z*C 370 CONTINUE 380 CONTINUE 390 CONTINUE C**** C**** C**** * * * * * * TEST CONVERGENCE * * * * * * * * * * * * C**** C**** Z=Q(K) IF(L.EQ.K) GOTO 450 C**** C**** C**** * * * * * * SHIFT FROM BOTTOM 2 * 2 MINOR * * * * * * C**** C**** X=Q(L) KMIN=K-1 Y=Q(KMIN) G=E(KMIN) H=E(K) F=((Y-Z)*(Y+Z)+(G-H)*(G+H))/(TWO*H*Y) G=SQRT(F*F+ONE) R=F+G IF(F.LT.ZERO) R=F-G F=((X-Z)*(X+Z)+H*(Y/R-H))/X C**** C**** C**** * * * * * * NEXT QR TRANSFORMATION * * * * * * * * * C**** C**** C=ONE S=ONE L2=L+1 DO 430 I=L2,K G=E(I) Y=Q(I) H=S*G G=C*G IMIN=I-1 IF(ABS(F).LT.ABS(H)) GOTO 392 IF(F.EQ.ZERO) GOTO 395 Z=H/F Z=ABS(F)*SQRT(ONE+Z*Z) GOTO 398 392 CONTINUE Z=F/H Z=ABS(H)*SQRT(ONE+Z*Z) GOTO 398 395 CONTINUE Z=ZERO GOTO 399 398 C=F/Z S=H/Z 399 E(IMIN)=Z F=X*C+G*S G=-X*S+G*C H=Y*S Y=Y*C DO 400 J=1,N JM=J+(IMIN-1)*N JI=J+(I-1)*N X=V(JM) Z=V(JI) V(JM)=X*C+Z*S V(JI)=-X*S+Z*C 400 CONTINUE IF(ABS(F).LT.ABS(H)) GOTO 412 IF(F.EQ.ZERO) GOTO 415 Z=H/F Z=ABS(F)*SQRT(ONE+Z*Z) GOTO 418 412 CONTINUE Z=F/H Z=ABS(H)*SQRT(ONE+Z*Z) GOTO 418 415 CONTINUE Z=ZERO GOTO 419 418 C=F/Z S=H/Z 419 Q(IMIN)=Z F=C*G+S*Y X=-S*G+C*Y DO 420 J=1,M JM=J+(IMIN-1)*M JI=J+(I-1)*M Y=U(JM) Z=U(JI) U(JM)=Y*C+Z*S U(JI)=-Y*S+Z*C 420 CONTINUE 430 CONTINUE E(L)=ZERO E(K)=F Q(K)=X GOTO 340 C**** C**** C**** * * * * * * CONVERGENCE * * * * * * * * * * * * * * * C**** C**** 450 CONTINUE IF(Z.GE.ZERO) GOTO 470 C**** C**** C**** * * * * * * Q(K) IS MADE NON-NEGATIVE * * * * * * * * C**** C**** Q(K)=-Z DO 460 J=1,N JK=J+(K-1)*N V(JK)=-V(JK) 460 CONTINUE 470 CONTINUE C**** C**** C**** * * * * * * ORDERING OF SINGULAR VALUES * * * * * * * C**** C**** NMIN=N-1 DO 540 I=1,NMIN IPLUS=I+1 480 CONTINUE DO 530 J=IPLUS,N IF(Q(I).GE.Q(J)) GOTO 530 R=Q(I) Q(I)=Q(J) Q(J)=R DO 490 K=1,N KI=K+(I-1)*N KJ=K+(J-1)*N R=V(KI) V(KI)=V(KJ) V(KJ)=R 490 CONTINUE DO 510 K=1,M KI=K+(I-1)*M KJ=K+(J-1)*M R=U(KI) U(KI)=U(KJ) U(KJ)=R 510 CONTINUE JPLUS=J+1 IF(JPLUS.GT.N) GOTO 540 GOTO 480 530 CONTINUE 540 CONTINUE GOTO 570 550 CONTINUE ICV=0 570 CONTINUE RETURN END