PROGRAM ANACOR EXTERNAL C82SPA CHARACTER*80 FOMA,TITLE CHARACTER*60 INDA,UTSC,UTVR,UTPR LOGICAL CONF COMMON/LINOUT/ICAR,IPRI COMMON/FMATOR/IE1,IE2,ALPHA,BETA,INPU,IPLTRA,IPL2DI,IPLEXT, 1 MEDSCA,MEDVAR,FOMA,CONF,ILIST,IPERM,MEDPER DIMENSION ITEX(20),ITEM(16),IFMT(60) IPRI = 6 ICAR = 11 WRITE (IPRI,4050) WRITE(IPRI,9109) 9109 FORMAT('WE WILL FIRST BUILD THE PARAMETER FILE INTERACTIVELY',/) WRITE(IPRI,9110) 9110 FORMAT('A TITLE FOR THE JOB') READ (5,*) TITLE WRITE(IPRI,9111) 9111 FORMAT('THE NUMBER OF ROWS') READ (5,*) K1 WRITE(IPRI,9112) 9112 FORMAT('THE NUMBER OF COLUMNS') READ (5,*) K2 WRITE(IPRI,9113) 9113 FORMAT('INDEX OF STANDARIZATION',/ + ,'-1.0 >>> COLUMNS ARE IN THE CENTROIDS OF THE MATCHING ROWS',/ 1 ,'+1.0 >>> ROWS ARE IN THE CENTROIDS OF THE MATCHING COLUMNS',/ 2 ,' 0.0 >>> ROWS AND COLUMNS ARE TREATED SYMMETRICALLY') READ (5,*) SN1 WRITE(IPRI,9114) 9114 FORMAT('FIRST DIMENSION TO BE ANALYZED') READ (5,*) IE1 WRITE(IPRI,9115) 9115 FORMAT('LAST DIMENSION TO BE ANALYZED') READ (5,*) IE2 WRITE(IPRI,9116) 9116 FORMAT('COMPUTATION OF VARIANCES, ONE IS YES, ZERO IS NO') READ (5,*) ICONF WRITE(IPRI,9117) 9117 FORMAT('NUMBER OF LISTED ROWS, ZERO MEANS TEN ROWS') READ (5,*) ILIST WRITE(IPRI,9118) 9118 FORMAT('PERMUTATION UP TO DIMENSION',/ 1 ,'0 MEANS DATA MATRICES ARE NOT REORDERED') READ (5,*) IPERM WRITE(IPRI,9119) 9119 FORMAT('PLOTS OF TRANSFORMATIONS',/ + ,'0 >>> NO SUCH PLOTS',/ 1 ,'1 >>> ONLY FOR THE ROWS',/ 2 ,'2 >>> ONLY FOR THE COLUMNS',/ 3 ,'3 >>> SEPARATE PLOTS FOR ROWS AND COLUMNS',/ 4 ,'4 >>> COMBINED PLOTS FOR ROWS AND COLUMNS') READ (5,*) IPLTRA WRITE(IPRI,9120) 9120 FORMAT('TWO DIMENSIONAL PLOTS',/ + ,'0 >>> NO SUCH PLOTS',/ 1 ,'1 >>> ONLY FOR THE ROWS',/ 2 ,'2 >>> ONLY FOR THE COLUMNS',/ 3 ,'3 >>> SEPARATE PLOTS FOR ROWS AND COLUMNS',/ 4 ,'4 >>> COMBINED PLOTS FOR ROWS AND COLUMNS') READ (5,*) IPL2DI WRITE(IPRI,9121) 9121 FORMAT('PLOTS OF ROW- VERSUS COLUMN SCORES',/ + '0 >>> NO SUCH PLOTS',/ 1 '1 >>> THESE PLOTS ARE MADE') READ (5,*) IPLEXT WRITE(IPRI,9124) 9124 FORMAT('GIVE THE FORMAT FOR THE DATA') READ (5,*) FOMA WRITE(IPRI,9126) 9126 FORMAT('NAME DATA FILE',/ 1 ,'THIS MUST BE AN EXISTING ASCII FILE') READ (5,*) INDA WRITE(IPRI,9127) 9127 FORMAT('NAME SCORES OUTPUT FILE',/ 1 ,'0 MEANS NO SUCH OUTPUT') READ (5,*) UTSC WRITE(IPRI,9128) 9128 FORMAT('NAME VARIANCE OUTPUT FILE',/ 1 ,'0 MEANS NO SUCH OUTPUT') READ (5,*) UTVR WRITE(IPRI,9129) 9129 FORMAT('NAME PERMUTATION OUTPUT FILE',/ 1 ,'0 MEANS NO SUCH OUTPUT') READ (5,*) UTPR OPEN(UNIT = 12, FILE = INDA, STATUS = 'OLD') INPU = 12 IF (UTSC .EQ. '0') THEN MEDSCA = 0 ELSE OPEN(UNIT = 13, FILE = UTSC, STATUS = 'NEW') MEDSCA = 13 END IF IF (UTVR .EQ. '0') THEN MEDVAR = 0 ELSE OPEN(UNIT = 14, FILE = UTVR, STATUS = 'NEW') MEDVAR = 14 END IF IF (UTPR .EQ. '0') THEN MEDPER = 0 ELSE OPEN(UNIT = 15, FILE = UTPR, STATUS = 'NEW') MEDPER = 15 END IF ALPHA=(1.+SN1)/2. BETA=(1.-SN1)/2. 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 INPU = 12 IF (MEDSCA.NE.0) MEDSCA = 13 IF (MEDVAR.NE.0) MEDVAR = 14 IF (ILIST.EQ.0) ILIST=10 IF (ILIST.EQ.-1.OR.ILIST.GT.K1) ILIST=K1 IF (IPERM.GT.IE2) IPERM=IE2 IF (IPERM.LT.IE1) MEDPER=0 IF (MEDPER.NE.0) MEDPER = 15 IF (IE1.GE.IE2) IPL2DI=0 ISUM=1+K1*K2+K1*MINK+K2*MINK+3*K1+3*K2+MINK+ 1 ICONF*(K1+K2+1)*NALDIM*(NALDIM+3)/2+(2*(K1+K2)) C CALL C82DEC(C82SPA,A,ISUM,K1,K2,MINK,NALDIM) CLOSE(INPU, STATUS = 'KEEP') IF (MEDSCA .EQ. 13) CLOSE(MEDSCA, STATUS = 'KEEP') IF (MEDVAR .EQ. 14) CLOSE(MEDVAR, STATUS = 'KEEP') IF (MEDPER .EQ. 15) CLOSE(MEDPER, STATUS = 'KEEP') STOP 20 WRITE (IPRI,5800) CLOSE(INPU, STATUS = 'KEEP') IF (MEDSCA .EQ. 13) CLOSE(MEDSCA, STATUS = 'DELETE') IF (MEDVAR .EQ. 14) CLOSE(MEDVAR, STATUS = 'DELETE') IF (MEDPER .EQ. 15) CLOSE(MEDPER, STATUS = 'DELETE') 1000 FORMAT(16I5) 1200 FORMAT(A80) 1400 FORMAT(F5.2,15I5) 4050 FORMAT (1H ,/////, +10X,58H +++++++++++++++++++++++++++++++++++++++++++++++++++++++++/ +10X,58H + +/ +10X,58H + August 1988 ANACOR Version 1.0 +/ +10X,58H + +/ +10X,58H + Most of the code for ANACOR was written by Bert +/ +10X,58H + Bettonvil in 1980, while he was working at the +/ +10X,58H + Department of Data Theory, University of Leiden. +/ +10X,58H + During a visit to UCLA Stef van Buuren turned the +/ +10X,58H + mainframe version into a PC version by making all +/ +10X,58H + output fit into 80 columns. Jan de Leeuw, Departments +/ +10X,58H + of Psychology and Mathematics, UCLA, added the +/ +10X,58H + interactive front end in August 1988. Binaries of +/ +10X,58H + this program for VMS, UNIX, and MSDOS are also +/ +10X,58H + available. The Macintosh version was compiled with +/ +10X,58H + Language Systems FORTRAN under MPW. The MSDOS version +/ +10X,58H + with MS FORTRAN, the UNIX version with f77, and the +/ +10X,58H + VMS version with VAX FORTRAN. +/ +10X,58H + +/ +10X,58H +++++++++++++++++++++++++++++++++++++++++++++++++++++++++/ + /////) 5800 FORMAT(' No solution in given dimensions') STOP END !!S C82 SUBROUTINE C82DEC(C82SUB,B,ISUM,K1,K2,MINK,NALDIM) COMMON/LINOUT/ICAR,IPRI DIMENSION A(20000) IDIM=20000 IF (ISUM.LE.IDIM) GO TO 10 WRITE (IPRI,1000) ISUM RETURN 10 CONTINUE CALL C82SUB(A,ISUM,K1,K2,MINK,NALDIM) 1000 FORMAT (//' FATAL ERROR: Memory allocation in C82DEC failed'/ + ' This job needs at least IDIM=',I8// + ' EXECUTION OF THE PROGRAM HAS BEEN TERMINATED'//) RETURN END !!S C82SPA 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 IAUX1=INDH2+K2 IAUX2=IAUX1+K1+K2 INDVAR=IAUX2+K1+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(IAUX1),A(IAUX2),A(INDVAR)) RETURN END !!S C82ORG SUBROUTINE C82ORG(K1,K2,MINK,NALDIM,IREST,F,D1,D2,Z1,Z2,XLABDA, 1 INTF,ID1,ID2,IDEN1,IDEN2,HULP1, 2 HULP2,IHULP1,IHULP2,AUX1,AUX2,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),AUX1(K1+K2),AUX2(K1+K2), 4 VARS(IREST) LOGICAL CONF CHARACTER*80 FOMA,TITLE COMMON/LINOUT/ICAR,IPRI COMMON/FMATOR/IE1,IE2,ALPHA,BETA,INPU,IPLTRA,IPL2DI,IPLEXT, 1 MEDSCA,MEDVAR,FOMA,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 (INPU,FOMA) (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 (INPU,FOMA) (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 (IPRI,5000) ILIST J=1 100 K=J+19 IF (K.GT.K2) K=K2 WRITE (IPRI,6700) DO 110 I=1,ILIST 110 WRITE (IPRI,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 (IPRI,1100) I WRITE (IPRI,1200) XLABDA(1) DO 270 I=2,MIN 270 WRITE (IPRI,1300) XLABDA(I) WRITE (IPRI,5100) WRITE (IPRI,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 (IPRI,6700) DO 300 I=1,K1 300 WRITE (IPRI,1400) I,ID1(I),(Z1(I,J1),J1=J,K) J=K+1 IF (J.LE.IE2) GOTO 290 WRITE (IPRI,5200) WRITE (IPRI,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 (IPRI,6700) DO 330 I=1,K2 330 WRITE (IPRI,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 (IPRI,5500) IF (NALDIM.EQ.1) WRITE (IPRI,5550) IE1 N=0 DO 400 I=1,NALDIM N=N+I 400 HULP1(I)=VARS(N) WRITE (IPRI,5600) (HULP1(I),I=1,NALDIM) IF (NALDIM.EQ.1) GOTO 440 IF (NALDIM.GT.2) WRITE (IPRI,5700) IF (NALDIM.EQ.2) WRITE (IPRI,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 (IPRI,5900) I IF (NALDIM.EQ.1) WRITE (IPRI,5950) I DO 450 J=1,NALDIM N=N+J 450 HULP1(J)=VARS(N) WRITE (IPRI,5600) (HULP1(J),J=1,NALDIM) IF (NALDIM.EQ.1) GOTO 490 IF (NALDIM.GT.2) WRITE (IPRI,6000) I IF (NALDIM.EQ.2) WRITE (IPRI,6050) I L=L+K CALL SUBTRI(NALDIM,VARS(L),HULP1,K) 490 CONTINUE DO 540 I=1,K2 IF (NALDIM.GT.1) WRITE (IPRI,6100) I IF (NALDIM.EQ.1) WRITE (IPRI,6150) I DO 500 J=1,NALDIM N=N+J 500 HULP1(J)=VARS(N) WRITE (IPRI,5600) (HULP1(J),J=1,NALDIM) IF (NALDIM.EQ.1) GOTO 540 IF (NALDIM.GT.2) WRITE (IPRI,6200) I IF (NALDIM.EQ.2) WRITE (IPRI,6250) I L=L+K CALL SUBTRI(NALDIM,VARS(L),HULP1,K) 540 CONTINUE C C C PLOT RESULTS C C C C Patch for 80 column plot output SvB 550 DO 541 I=1,K1 541 HULP1(I)=FLOAT(I) DO 542 I=1,K2 542 HULP2(I)=FLOAT(I) C End patch C IF (IPLTRA.LE.0.OR.IPLTRA.GT.4) GOTO 620 GOTO (560,580,560,600),IPLTRA 560 DO 570 I=IE1,IE2 WRITE (IPRI,3100) I 570 CALL PLOT(HULP1,Z1(1,I),IDEN1,K1,K1,1,IPRI) IF (IPLTRA.EQ.1) GOTO 620 580 DO 590 I=IE1,IE2 WRITE (IPRI,3300) I 590 CALL PLOT(HULP2,Z2(1,I),IDEN2,K2,K2,1,IPRI) GOTO 620 600 DO 610 I=IE1,IE2 WRITE (IPRI,3500) I DO 611 J=1,K1 611 AUX1(J)=Z1(J,I) DO 612 J=1,K2 612 AUX1(K1+J)=Z2(J,I) CALL PLOT(HULP1,AUX1,IDEN1,K1+K2,K1,1,IPRI) C Bad programming practices. Relies on contageous arrays.SvB. C 610 WRITE (IPRI,3550) 610 CONTINUE 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 (IPRI,3200) I,J 640 CALL PLOT(Z1(1,I),Z1(1,J),IDEN1,K1,K1,1,IPRI) IF (IPL2DI.EQ.1) GOTO 690 650 I1=IE2-1 DO 660 I=IE1,I1 J1=I+1 DO 660 J=J1,IE2 WRITE (IPRI,3400) I,J 660 CALL PLOT(Z2(1,I),Z2(1,J),IDEN2,K2,K2,1,IPRI) GOTO 690 670 I1=IE2-1 DO 680 I=IE1,I1 J1=I+1 DO 680 J=J1,IE2 WRITE (IPRI,3600) I,J DO 671 JJ=1,K1 AUX1(JJ)=Z1(JJ,I) 671 AUX2(JJ)=Z1(JJ,J) DO 672 JJ=1,K2 AUX1(K1+JJ)=Z2(JJ,I) 672 AUX2(K1+JJ)=Z2(JJ,J) CALL PLOT(AUX1,AUX2,IDEN1,K1+K2,K1,1,IPRI) 680 CONTINUE 690 IF (IPLEXT.NE.1) GOTO 710 DO 700 I=IE1,IE2 WRITE (IPRI,3700) I 700 CALL PLOT(Z1(1,I),Z2(1,I),IDEN1,MINK,MINK,1,IPRI) 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 (IPRI,6400) I J=1 730 K=J+19 IF (K.GT.K2) K=K2 WRITE (IPRI,6450) (IHULP2(J1),J1=J,K) WRITE (IPRI,6500) (ID2(IHULP2(J1)),J1=J,K) WRITE (IPRI,6700) DO 740 J1=1,K1 I1=IHULP1(J1) 740 WRITE (IPRI,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 (IPRI,5400) RETURN 780 WRITE (IPRI,1600) RETURN 790 WRITE (IPRI,1500) C C C FORMATS C C C 1100 FORMAT(/' Total number of observations',4X,I10) 1200 FORMAT(' Nonzero singular values',11X,F8.3) 1300 FORMAT(35X,F8.3) 1400 FORMAT(1H ,I5,1X,I7,3X,9F7.2) 1500 FORMAT(' Program stopped because the input matrix cannot be', 1 'analyzed') 1600 FORMAT(' No scores can be attached'/ 1 ' Rows and columns are independent') 1700 FORMAT(8F10.6) 1800 FORMAT(8E10.3) 1900 FORMAT(16I5) 3100 FORMAT(1H ,32X,'Row scores in dimension',I2) 3200 FORMAT(1H ,18X,'Row scores in dimensions',I2,' (X axis) and',I2, + ' (Y axis)') 3300 FORMAT(1H ,32X,'Column scores in dimension',I2) 3400 FORMAT(1H ,18X,'Column scores in dimensions',I2, 1 ' (X axis) and',I2,' (Y axis)') 3500 FORMAT(1H ,22X,'Row (N) and column (C) scores in dimension',I2) C 3550 FORMAT(33H THE ROWS ARE LABELLED BY DIGITS,, C 1 23H THE COLUMNS BY LETTERS) 3600 FORMAT(1H ,10X,'Row (N) and column (C) scores in dimensions', 1 I2,' (X axis) and',I2,' (Y axis)') 3700 FORMAT(1H ,13X,'Row scores (X axis) vs column scores (Y axis) ', 1 'in dimension',I2) 5000 FORMAT(1H /1H ,9HThe first,I4,24H rows of the data matrix/) 5100 FORMAT(1H /1H ,10HROW SCORES/1H ,3HRow,5X,8HMarginal) 5150 FORMAT(7H number,2X,9Hfrequency,3X,6HScores) 5200 FORMAT(1H /1H ,2X,13HCOLUMN SCORES/' Column',2X,8HMarginal) 5300 FORMAT(1H ,16I5/(1H ,16I5)) 5400 FORMAT(32H No solution in given dimensions) 5500 FORMAT(1H / + ' Variances of the singular values '/) 5550 FORMAT(1H / + ' Variance of singular value number ',I5/) 5600 FORMAT(10F8.3/(10F8.3)) 5700 FORMAT(1H / + ' Correlation matrix of the singular values '/) 5750 FORMAT(1H / + ' Correlation of the singular values '/) 5900 FORMAT(1H / + ' Variances of the scores of row number ',I5/) 5950 FORMAT(1H / + ' Variance of the score of row number ',I5/) 6000 FORMAT(1H / + ' Correlation matrix of the scores of row number ',I5/) 6050 FORMAT(1H / + ' Correlation of the scores of row number ',I5/) 6100 FORMAT(1H / + ' Variances of the scores of column number ',I5/) 6150 FORMAT(1H / + ' Variance of the score of column number ',I5/) 6200 FORMAT(1H / + ' Correlation matrix of the scores of column number ',I5/) 6250 FORMAT(1H / + ' Correlation of the scores of column number ',I5/) 6400 FORMAT(1H / + ' The data matrix permuted according to the scores in ', + 'dimension',I2/) 6450 FORMAT(15H Column ,12I5/(15X,12I5)) 6500 FORMAT(15H Row marginal,12I5/(15X,12I5)) 6600 FORMAT(1H ,I5,I7,2X,12I5/(15X,12I5)) 6700 FORMAT(1H ) RETURN END !!S SUBTRI SUBROUTINE SUBTRI(N,A,B,N2) DIMENSION A(N2),B(N) COMMON/LINOUT/ICAR,IPRI 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 (IPRI,1000) (A(J),J=K,L) K=K+I I=I+1 IF (I.LE.N) GOTO 30 WRITE (IPRI,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(10F8.2) 2000 FORMAT(1H ) END !!S DERIVS 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 !!S ORDER 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 !!S PLOT 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(64),XPR(10),LABEL(20),X(NITEM),Y(NITEM),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/20/,NC/64/,NN/9/ C ##=NC C ##=NUMBER OF LINES PER PAGE DATA LABEL/1H0,1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9, + 1HJ,1HA,1HB,1HC,1HD,1HE,1HF,1HG,1HH,1HI/ C 1 FORMAT(9X, 3H.--,9 (7H+------),5H+---.) C ##=NN 2 FORMAT(1X,F7.2,1X,1H|,2X,64A1,3X,1H|) C ##=NC 3 FORMAT(1X,F7.2,1X,1H|,69X,1H|) C ##=NC+5 4 FORMAT(8X, 10F7.2) 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-1.0) STEPY=(STEPX*NC)/(NL-0.5) 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,NN 130 XPR(J+1)=XPR(J)+STEPX*7.0 WRITE(IWRITE,4) XPR RETURN END !!S SIVAD 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=0.30E-38 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