EXTERNAL P80SPA DIMENSION TITLE(20) COMMON/FMATOR/ICAR,ISF,ISL,NCID,IVF,IVL,IV,IVKW,IFULL,IPLDIM, 1 IPLOBS,IDATA,MEDSCO,MEDVAR,IMED,MEDLAB,NCAT(16) WRITE(6,5000) READ(5,1100) TITLE WRITE(6,5100) TITLE READ(5,1500) N1,NVAR,NCID N=N1 READ(5,1500) ISF,ISL,ICONF,IVF,IVL READ(5,1500) ICAR,IPLDIM,IPLOBS,IDATA,MEDSCO,MEDVAR,IMED,MEDLAB IREST=NVAR NCTOT=NVAR*NCID IFULL=0 IF (N.EQ.0) IFULL=1 N1=NCID**NVAR IF (N1.EQ.0) N1=1 IF (NCID.GT.0) IREST=0 20 IF (IREST.LE.0) GOTO 40 READ(5,1500) (NCAT(J),J=1,16) WRITE(1,1500) (NCAT(J),J=1,16) DO 30 J=1,16 IF (NCAT(J).EQ.0) GOTO 40 N1=N1*NCAT(J) 30 NCTOT=NCTOT+NCAT(J) IREST=IREST-16 GOTO 20 40 MIN=NCTOT IF (N.EQ.0) N=N1 IF (N.LT.NCTOT) MIN=N IF (ISL.LE.0) ISL=2 IF (ISL.GE.MIN) ISL=MIN-1 IF (ISF.LE.0) ISF=1 IF (ISF.GT.ISL) RETURN IF (IVL.LT.ISF.OR.IVL.GT.ISL) IVL=ISL IF (IVF.LT.ISF) IVF=ISF IF (IVF.GT.ISL) ICONF=0 IF (IMED.GT.ISL) IMED=ISL IF (IMED.EQ.0) IMED=ISL WRITE(6,5200) N,NVAR,ISF,ISL IF (ICONF.NE.1) ICONF=0 IF (ICONF.EQ.0) WRITE(6,5300) IF (ICONF.EQ.1) WRITE(6,5400)IVF,IVL IV=IVL-IVF+1 IVKW=IV*(IV+1)/2 ISUM=1+N*(NVAR+MIN+3)+NCTOT*(MIN+2)+3*NVAR+ 1 ICONF*((N+NCTOT)*(IV+IVKW)+IVKW+NCTOT) CALL P80DEC(P80SPA,A,ISUM,N,NVAR,NCTOT,MIN,ICONF) C CALL DYNAM(P80SPA,A,ISUM,N,NVAR,NCTOT,MIN,ICONF) 1100 FORMAT(20A4) 1500 FORMAT(16I5) 5000 FORMAT(1H1///1H0,5X,17HA N A P R O F - S,76X, 1 29HBERT BETTONVIL JAN DE LEEUW/ 2 1H0,5X,12HOCTOBER 1980,81X,24HDEPARTMENT OF DATATHEORY/ 3 1H0,98X,13HBREESTRAAT 70/1H0,5X, 4 53HCORRESPONDENCE ANALYSIS OF PROFILE-FREQUENCY-MATRICES, 5 40X,24HLEIDEN - THE NETHERLANDS///) 5100 FORMAT(12H0 JOB TITLE,24X,20A4) 5200 FORMAT(21H0 NUMBER OF PROFILES,14X,I5/ 1 22H0 NUMBER OF VARIABLES,13X,I5/ 2 30H0 FIRST DIMENSION OF ANALYSIS,5X,I5/ 3 29H0 LAST DIMENSION OF ANALYSIS,6X,I5) 5300 FORMAT(28H0 NO VARIANCES ARE COMPUTED) 5400 FORMAT(48H0 THE VARIANCES ARE COMPUTED FOR THE DIMENSIONS, 1 I5,8H UP TO,I5) STOP END SUBROUTINE P80DEC(P80SUB,B,ISUM,N,NVAR,NCTOT,MIN,ICONF) DIMENSION A(20000) IDIM=20000 IF (ISUM.LE.IDIM) GO TO 10 WRITE (6,1000) ISUM RETURN 10 CONTINUE CALL P80SUB(A,ISUM,N,NVAR,NCTOT,MIN,ICONF) 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 P80SPA(A,ISUM,N,NVAR,NCTOT,MIN,ICONF) DIMENSION A(ISUM) LOGICAL CONF INDNC=1 INDIOB=INDNC+NVAR INDD=INDIOB+N*NVAR INDE=INDD+NCTOT INDPI=INDE+N INDX=INDPI+N INDY=INDX+N*MIN INDXL=INDY+NCTOT*MIN INDVH=INDXL+NCTOT INDPLV=INDVH+N INDIHN=INDY-N INDIHM=INDXL-NCTOT INDHM1=INDXL INDVAR=INDPLV+2*NVAR IREST=ISUM-INDVAR+1 CONF=(ICONF.EQ.1) CALL P80ORG(N,NCTOT,MIN,NVAR,A(INDNC),A(INDIOB),A(INDD),A(INDE), 1 A(INDE),A(INDPI),A(INDX),A(INDY),A(INDXL),A(INDVH), 2 A(INDPLV),A(INDIHN),A(INDIHM),A(INDHM1), 3 CONF,IREST,A(INDVAR)) RETURN END SUBROUTINE P80ORG(N,M,MIN,NVAR,NCAT,IOBS,D,E,IE,PI,X,Y,XLABDA,VH, 1 IPLOT,IHN,IHM,IHM1,CONF,IREST,VARIAN) DIMENSION NCAT(NVAR),IOBS(N,NVAR),D(M),E(N),IE(N),PI(N),X(N,MIN), 1 Y(M,MIN),XLABDA(M),VH(N),IPLOT(2,NVAR),IHN(N),IHM(M), 2 IHM1(M),VARIAN(IREST),FMT(20) LOGICAL CONF COMMON/FMATOR/ICAR,ISF,ISL,NCID,IVF,IVL,IV,IVKW,IFULL,IPLDIM, 1 IPLOBS,IDATA,MEDSCO,MEDVAR,IMED,MEDLAB,IPLT(16) IF (NCID.GT.0) GOTO 10 REWIND 1 READ(1,1000) (NCAT(J),J=1,NVAR) GOTO 30 10 DO 40 J=1,NVAR 40 NCAT(J)=NCID 30 READ(5,1100) FMT K=1 L=2 IF (IPLDIM.GE.0) K=2 IF (IPLOBS.GE.0) L=1 IF (K.LE.L) READ(5,1000) ((IPLOT(I,J),I=K,L),J=1,NVAR) IT=0 IF (IFULL.EQ.0) GOTO 5 DO 15 J=1,NVAR IHM(J)=1 15 IOBS(1,J)=1 DO 20 I=2,N J=NVAR 35 IHM(J)=IHM(J)+1 IF (IHM(J).LE.NCAT(J)) GOTO 45 J=J-1 GOTO 35 45 IF (J.EQ.NVAR) GOTO 25 J=J+1 DO 55 K=J,NVAR 55 IHM(K)=1 25 DO 20 J=1,NVAR 20 IOBS(I,J)=IHM(J) READ(ICAR,FMT) (IE(I),I=1,N) DO 65 I=1,N PI(I)=IE(I) 65 IT=IT+IE(I) GOTO 75 5 DO 50 I=1,N READ(ICAR,FMT) (IOBS(I,J),J=1,NVAR),K PI(I)=K 50 IT=IT+K 75 IF (IT.GT.0) GOTO 60 DO 70 I=1,N 70 PI(I)=1. IT=N 60 IF (IDATA.EQ.0) IDATA=10 IF (IDATA.GT.N) IDATA=N WRITE (6,1900) IDATA DO 100 I=1,IDATA K=PI(I)+.1 100 WRITE (6,2000) I,K,(IOBS(I,K1),K1=1,NVAR) WRITE(6,1400) IT CALL P80SVD(N,M,MIN,NDIM,NVAR,NCAT,IOBS,PI,E,D,X,Y,XLABDA,VH,T,IT) IF (NDIM.LT.0) RETURN IF (NDIM.EQ.0) GOTO 1020 IF (ISL.GT.NDIM) ISL=NDIM IF (ISF.GT.ISL) RETURN IF (IVL.GT.ISL) IVL=ISL IF (IVF.GT.IVL) CONF=.FALSE. WRITE (6,1500) XLABDA(1) IF (NDIM.EQ.1) GOTO 85 WRITE (6,1300) (XLABDA(J),J=2,NDIM) 85 WRITE(6,1600) DO 88 L=1,N K=PI(L) 88 WRITE(6,1700) L,K,(X(L,J),J=ISF,ISL) J=0 DO 90 K=1,NVAR WRITE(6,1800) K K1=NCAT(K) DO 90 L=1,K1 J=J+1 I=D(J)+.1 90 WRITE (6,1700) L,I,(Y(J,K2),K2=ISF,ISL) IF (MEDSCO.LE.0) GOTO 110 WRITE (MEDSCO,5000) (XLABDA(J),J=ISF,IMED) DO 120 I=1,M 120 WRITE (MEDSCO,5000) (Y(I,J),J=ISF,IMED) DO 125 I=1,N 125 WRITE (MEDSCO,5000) (X(I,J),J=ISF,IMED) 110 IF (MEDLAB.LE.0) GOTO 140 H1=0. DO 150 K=1,NVAR H1=H1+1. H2=H1 DO 150 I=1,K1 H2=H2+.001 150 WRITE (MEDLAB,5100) H2 H1=0. DO 160 I=1,N H1=H1+.001 160 WRITE (MEDLAB,5100) H1 140 H1=IT IF (.NOT.CONF) GOTO 190 ICVX=1 ICVY=ICVX+N*IVKW ICVL=ICVY+M*IVKW IIC=ICVL+IVKW IVX=IIC+NVAR+1 IVY=IVX+N*IV CALL P80VAR(N,M,NDIM,D,X,Y,XLABDA,VARIAN(ICVX),VARIAN(ICVY), 1 VARIAN(ICVL),IVF,IVL,IV,IVKW,H1,E,PI,VARIAN(IIC), 2 NVAR+1,NVAR,IOBS,NCAT,VARIAN(IVX),VARIAN(IVY), 3 VH,T,MEDVAR,IMED) 190 DO 205 I=1,N 205 IE(I)=1.+(9.*PI(I)-.1)/H1 IF (ISF.NE.ISL) GOTO 200 IPLDIM=1 IPLOBS=2 GOTO 280 200 I=0 DO 210 K=1,NVAR K1=NCAT(K) DO 210 L=1,K1 I=I+1 210 IHM(I)=K K=ISL-1 DO 220 I=ISF,K J=I+1 DO 220 L=J,ISL WRITE (6,2100) I,L CALL PLOT21(Y(1,I),Y(1,L),M,IHM,IHM1,-3) 220 WRITE (6,2200) K=ISL-1 DO 250 I=ISF,K J=I+1 DO 250 L=J,ISL WRITE (6,2400) I,L CALL PLOT21(X(1,I),X(1,L),N,IE,IHN,1) 250 WRITE (6,2500) 280 IF (IPLDIM.EQ.0.AND.IPLOBS.EQ.0) RETURN IDK=IPLDIM IOK=IPLOBS J1=1 DO 290 K=1,NVAR K1=NCAT(K) IF (IPLDIM.GT.0) GOTO 300 IDK=IPLOT(1,K) 300 IF (IPLOBS.GE.0) GOTO 310 IOK=IPLOT(2,K) 310 CONTINUE IF (IDK.LE.0) GOTO 395 IF (IDK.EQ.2) GOTO 360 L=J1-1 H2=0. DO 325 J=1,K1 L=L+1 325 IF (D(L).GT.H2) H2=D(L) L=J1-1 DO 330 J=1,K1 L=L+1 330 IHM(J)=1.+(9.*D(L)-.1)/H2 DO 340 J=ISF,ISL WRITE (6,2800) K,J CALL PLOT11(K1,Y(J1,J),K1,IHM,IHM1,1) 340 WRITE (6,2900) IF (IDK.EQ.1) GOTO 395 360 DO 380 L=1,K1 380 IHM(L)=L J=ISL-1 DO 390 L=ISF,J L1=L+1 DO 390 L2=L1,ISL WRITE (6,3600) K,L,L2 CALL PLOT21(Y(J1,L),Y(J1,L2),K1,IHM,IHM1,-3) 390 WRITE (6,3700) 395 IF (IOK.LE.0) GOTO 290 320 L=J1-1 DO 350 J=ISF,ISL DO 370 L1=1,N L2=IOBS(L1,K) IF (L2.EQ.0) VH(L1)=X(L1,J) 370 IF (L2.NE.0) VH(L1)=Y(L+L2,J) WRITE (6,3200) K,J CALL PLOT21(VH,X(1,J),N,IE,IHN,2) 350 WRITE (6,3300) IF (IOK.EQ.1) GOTO 290 400 J=ISL-1 DO 410 L=ISF,J L1=L+1 DO 410 L2=L1,ISL WRITE (6,4000) K,L,L2 CALL PLOT22(X(1,L),X(1,L2),N,Y(J1,L),Y(J1,L2),K1, 1 IOBS(1,K),IHM,IHN,IHM1,-1) 410 WRITE (6,4100) 290 J1=J1+K1 1000 FORMAT(16I5) 1100 FORMAT(20A4) RETURN 1020 WRITE (6,1200) 1200 FORMAT(29H0ALL VARIABLES ARE INDEPENENT) 1500 FORMAT(25H0NON-ZERO SINGULAR VALUES,5X,F10.4) 1300 FORMAT(30X,F10.4) 1400 FORMAT(//23H0NUMBER OF OBSERVATIONS,7X,I10) 1600 FORMAT(1H1,14X,14HPROFILE SCORES/ 1 19H0 NUMBER FREQUENCY,12X,6HSCORES/) 1700 FORMAT(I8,I10,4X,(11F10.4)) 1800 FORMAT(/1H0,5X,27HCATEGORY SCORES OF VARIABLE,I5/ 1 19H0CATEGORY FREQUENCY,12X,6HSCORES/) 1900 FORMAT(8H1LIST OF,I5,25H ROWS OF THE DATAMATRIX:/ 1 19H0 NUMBER FREQUENCY,8X,7HPROFILE/) 2000 FORMAT(I8,I10,4X,(22I5)) 2100 FORMAT(1H1,45X,29HCATEGORY-SCORES IN DIMENSIONS,I5,5H AND,I5) 2200 FORMAT(1H0,41X,46HCATEGORIES OF THE SAME VARIABLE HAVE THE SAME , 1 6HLETTER) 2400 FORMAT(1H1,46X,28HPROFILE-SCORES IN DIMENSIONS,I5,5H AND,I5) 2500 FORMAT(1H0,49X,37HPROFILES ARE LABELED BY THEIR WEIGHTS) 2800 FORMAT(1H1,44X,22HCATEGORIES OF VARIABLE,I5,14H IN DIMENSION,I5) 2900 FORMAT(1H0,48X,39HCATEGORIES ARE LABELED BY THEIR WEIGHTS) 3200 FORMAT(1H1,24X,27HCATEGORY-SCORES OF VARIABLE,I5, 1 50H VERSUS CORRESPONDING PROFILE-SCORES IN DIMENSION,I5) 3300 FORMAT(1H0,49X,37HPROFILES ARE LABELED BY THEIR WEIGHTS) 3600 FORMAT(1H1,39X,22HCATEGORIES OF VARIABLE,I5,15H IN DIMENSIONS,I5, 1 5H AND,I5) 3700 FORMAT(1H0,48X,38HCATEGORIES ARE LABELED BY THEIR NUMBER) 4000 FORMAT(1H1,32X,22HCATEGORIES OF VARIABLE,I5, 1 28H AND PROFILES IN DIMENSIONS,I5,5H AND,I5) 4100 FORMAT(1H0,26X,35HCATEGORIES ARE LABELED BY LETTERS, , 1 47HPROFILES BY DIGITS, ACCORDING TO THEIR CATEGORY) 5000 FORMAT(8F10.6) 5100 FORMAT (2X,F6.3) RETURN END SUBROUTINE P80SVD(N,M,MIN,NDIM,NVAR,NCAT,IOBS,PI,E,D,X,Y,XLABDA, 1 VH,T,IT) DIMENSION NCAT(NVAR),IOBS(N,NVAR),PI(N),E(N),D(M),X(N,MIN), 1 Y(M,MIN),XLABDA(MIN),VH(MIN) T=0 DO 15 J=1,M 15 D(J)=0. IF (N.LT.M) GOTO 10 DO 20 I=1,N H1=PI(I) H2=0. J=0 DO 30 K=1,NVAR K1=NCAT(K) L=IOBS(I,K) IF (L.LT.0.OR.L.GT.K1) IOBS(I,K)=0 DO 30 K2=1,K1 J=J+1 H=0. IF (K2.EQ.L) H=H1 D(J)=D(J)+H H2=H2+H 30 X(I,J)=H T=T+H2 20 E(I)=SQRT(H2) DO 40 J=1,M H=SQRT(D(J)) D(J)=H DO 40 I=1,N H1=H*E(I) IF (H1.LE.0.) GOTO 40 X(I,J)=X(I,J)/H1-H1/T 40 CONTINUE COMMENT<<>> >>> CALL SIVAD (N,M,XLABDA,X,M*N,Y,M*M,VH,ICV) GOTO 50 10 DO 60 I=1,N H1=PI(I) H2=0. J=0 DO 70 K=1,NVAR K1=NCAT(K) L=IOBS(I,K) IF (L.LT.0.OR.L.GT.K1) IOBS(I,K)=0 DO 70 K2=1,K1 J=J+1 H=0. IF (K2.EQ.L) H=H1 D(J)=D(J)+H H2=H2+H 70 Y(J,I)=H T=T+H2 60 E(I)=SQRT(H2) DO 80 J=1,M H=SQRT(D(J)) D(J)=H DO 80 I=1,N H1=H*E(I) IF (H1.EQ.0.) GOTO 80 Y(J,I)=Y(J,I)/H1-H1/T 80 CONTINUE C<<< SINGULAR VALUE DECOMPOSITION >>> CALL SIVAD (M,N,XLABDA,Y,M*N,X,N*N,VH,ICV) 50 IF (ICV.EQ.0) GOTO 1010 EPS=SQRT(FLOAT(M*N))*1.E-7 DO 75 J=1,MIN 75 IF (XLABDA(J).LT.EPS) GOTO 85 85 NDIM=J-1 DO 90 I=J,MIN 90 XLABDA(J)=0. T=SQRT(T) DO 120 I=1,M H1=D(I)/T IF (H1.LE.0.) GOTO 120 DO 130 J=1,NDIM 130 Y(I,J)=Y(I,J)*XLABDA(J)/H1 H1=H1*T D(I)=H1*H1 120 CONTINUE T=T*T DO 140 J=1,NDIM 140 XLABDA(J)=XLABDA(J)*XLABDA(J) H1=IT DO 100 I=1,N DO 110 J=1,NDIM 110 X(I,J)=0. J=0 H2=0. DO 150 K=1,NVAR K2=IOBS(I,K) IF (K2.EQ.0) GOTO 150 K2=J+K2 H2=H2+1. DO 170 L1=1,NDIM 170 X(I,L1)=X(I,L1)+Y(K2,L1) 150 J=J+NCAT(K) IF (H2.EQ.0.) GOTO 100 DO 160 L1=1,NDIM 160 X(I,L1)=X(I,L1)/XLABDA(L1)/H2 100 E(I)=H2 RETURN 1010 WRITE(6,1020) 1020 FORMAT(36H0SINGULAR VALUE DECOMPOSITION FAILED) NDIM=-1 RETURN END SUBROUTINE P80VAR(N,M,NDIM,D,X,Y,XLABDA,VARX,VARY,VARLAB, 1 IVF,IVL,IV,IVKW,TOTOBS,E,PI,IC,NVP1,NVAR,IOBS, 2 NCAT,VX,VY,VL,T,MEDVAR,IMED) DIMENSION D(M),X(N,NDIM),Y(M,NDIM),XLABDA(NDIM), 1 VARX(N,IVKW),VARY(M,IVKW),VARLAB(IVKW),E(N),PI(N), 2 IC(NVP1),IOBS(N,NVAR),NCAT(NVAR),VX(N,IV),VY(M,IV), 3 VL(IV) DO 10 J=1,IVKW VARLAB(J)=0. DO 20 I=1,N 20 VARX(I,J)=0. DO 10 I=1,M 10 VARY(I,J)=0. DO 30 K=1,N PK=PI(K)/TOTOBS IF (PK.EQ.0.) GOTO 30 H1=E(K) I=0 L=0 DO 40 J=1,NVAR J1=IOBS(K,J) IF (J1.EQ.0) GOTO 40 I=I+1 IC(I)=L+J1 40 L=L+NCAT(J) IC(I+1)=0 C C C<<< COMPUTATION OF THE DERIVATIVES OF XLABDA(IVF),.....,XLABDA(IVL) C C<<< AND THE MATCHING COLUMNS OF Y WITH RESPECT TO PI(K) C C C J1=0 DO 50 J=IVF,IVL J1=J1+1 XKJ=X(K,J) DO 60 I=1,M 60 VY(I,J1)=((1.-XKJ*XKJ)*Y(I,J)/2.-XKJ)*H1 DO 70 L=1,NDIM H2=H1*XLABDA(L)*XKJ*X(K,L) I=1 80 L1=IC(I) IF (L1.EQ.0) GOTO 90 H2=H2-Y(L1,J)*Y(L1,L) I=I+1 GOTO 80 90 IF (L.EQ.J) GOTO 100 H3=XLABDA(J)-XLABDA(L) IF (ABS(H3).LT.1.E-7) WRITE(6,1600) H2=H2/H3 DO 110 I=1,M 110 VY(I,J1)=VY(I,J1)+H2*Y(I,L) GOTO 70 100 VL(J1)=H2*TOTOBS/T 70 CONTINUE DO 120 I=1,M 120 VY(I,J1)=VY(I,J1)*TOTOBS/T I=1 130 L1=IC(I) IF (L1.EQ.0) GOTO 50 VY(L1,J1)=VY(L1,J1)+(XKJ-Y(L1,J))*TOTOBS/D(L1) I=I+1 GOTO 130 50 CONTINUE C C C<<< COMPUTATION OF THE DERIVATIVES OF X(IVF),....,X(IVL) C C C DO 140 I=1,N H2=E(I) J1=0 DO 150 J=IVF,IVL J1=J1+1 150 VX(I,J1)=-H2*VL(J1)*X(I,J) IF (H2.EQ.0.) GOTO 140 L=0 DO 145 L1=1,NVAR K1=IOBS(I,L1) IF (K1.EQ.0) GOTO 145 K1=L+K1 DO 155 J=1,IV 155 VX(I,J)=VX(I,J)+VY(K1,J) 145 L=L+NCAT(L1) J1=0 DO 135 J=IVF,IVL J1=J1+1 135 VX(I,J1)=VX(I,J1)/XLABDA(J)/H2 140 CONTINUE L=1 J1=1 DO 160 I=1,IVKW VARLAB(I)=VARLAB(I)+PK*VL(L)*VL(J1) DO 170 L1=1,N 170 VARX(L1,I)=VARX(L1,I)+PK*VX(L1,L)*VX(L1,J1) DO 180 L1=1,M 180 VARY(L1,I)=VARY(L1,I)+PK*VY(L1,L)*VY(L1,J1) J1=J1+1 IF (J1.LE.L) GOTO 160 L=L+1 J1=1 160 CONTINUE 30 CONTINUE DO 190 I=1,IVKW VARLAB(I)=VARLAB(I)/TOTOBS DO 200 L1=1,N 200 VARX(L1,I)=VARX(L1,I)/TOTOBS DO 190 L1=1,M 190 VARY(L1,I)=VARY(L1,I)/TOTOBS WRITE(6,1100) J=0 I=1 DO 210 L=1,IV J=J+L WRITE (6,1000) (VARLAB(K),K=I,J) 210 I=J+1 WRITE(6,1200) DO 220 K=1,N WRITE (6,1300) K J=0 I=1 DO 220 L=1,IV J=J+L WRITE (6,1000) (VARX(K,J1),J1=I,J) 220 I=J+1 J=0 DO 230 K=1,NVAR WRITE(6,1400) K K1=NCAT(K) DO 230 K2=1,K1 J=J+1 WRITE(6,1500) K2 I1=0 I2=1 DO 230 L=1,IV I1=I1+L WRITE (6,1000) (VARY(J,J1),J1=I2,I1) 230 I2=I1+1 IF (IMED.LT.IVF) RETURN IF (MEDVAR.LE.0) RETURN IF (IMED.GT.IVL) IMED=IVL L=(IMED-IVF+1)*(IMED-IVF+2)/2 WRITE (MEDVAR,2000) (VARLAB(I),I=1,L) DO 300 I=1,M 300 WRITE(MEDVAR,2000) (VARY(I,J),J=1,L) DO 310 I=1,N 310 WRITE(MEDVAR,2000) (VARX(I,J),J=1,L) 1000 FORMAT(2X,13E10.2) 1100 FORMAT(37H1COVARIANCE-MATRIX OF THE EIGENVALUES/) 1200 FORMAT(/50H0COVARIANCE-MATRICES OF THE SCORES OF THE PROFILES) 1300 FORMAT(15H0PROFILE NUMBER,I5) 1400 FORMAT(/47H0COVARIANCE-MATRICES OF THE CATEGORY-SCORES OF , 1 8HVARIABLE,I5) 1500 FORMAT(16H0CATEGORY NUMBER,I5) 1600 FORMAT(50H0WARNING: EQUAL EIGENVALUES CAUSE DIVISION BY ZERO) 2000 FORMAT(8E10.3) 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=SPANX/66. 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 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 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