C CANONICAL ANALYSIS JUNE 30, 1966 C THIS IS A SIFTED VERSION OF BMD06M ORIGINALLY WRITTEN IN C FORTRAN II. SOME MODIFICATIONS WERE MADE TO MAKE IT OPERABLE C AND SLIGHTLY MORE EFFICIENT THAN THE SIFTED VERSION. DIMENSION RHO(70,70),DATA(100),SX(70),SX2(70),D(35,35), 1A(35,35),E(35,35),G(35,35),B(35),C(35),QQQ(70,35) DIMENSION L1(35),K1(35),FMT(180) DIMENSION NNEWA(99),LLCODE(99),LLVA(99),BBNEW(99) DIMENSION V(70) COMMON QQQ COMMON N , A , B , C , RHO , D EQUIVALENCE (G,QQQ(1226)),(E,QQQ) EQUIVALENCE(DATA,V) C DATA Q003HL/'YES '/ C 405 FORMAT(42H1BMD06M - CANONICAL ANALYSIS - VERSION OF , 116HDECEMBER 8, 1964/ 2/40H HEALTH SCIENCES COMPUTING FACILITY,UCLA///) C C THIS PROGRAM USES TWO SCRATCH TAPES. THEY ARE DESIGNATED HERE C BY IT2 AND IT3. C INFORM=5 CALL USAGEB('BMD06M') IT2=2 IT3=3 REWIND IT2 REWIND IT3 DOUBLE PRECISION A123,B123,C123,XODE,KODE,TODE,PROB DATA A123/'FINISH'/,B123/'PROBLM'/,C123/'TRNGEN'/ 999 READ (5,5001)TODE,PROB,NCASES,NVAR,IP,IQ,CRITER, PRINTD,NVG,NTAPE, 1KVR 5001 FORMAT(2A6,I4,I3,2I2,F5.0,A3,35X,3I2) IF(TODE.EQ.A123) GO TO 401 IF(TODE.EQ.B123) GO TO 403 402 WRITE (6,404) TODE 404 FORMAT (51H PROGRAM EXPECTED FINISH OR PROBLM CARD. IT READ A , 1 A6,5H CARD ) 401 IF(INFORM-5)27,27,26 26 REWIND INFORM 27 STOP 403 IF((NTAPE-IT2)*(NTAPE-IT3))473,440,473 473 CALL TPWD(NTAPE,INFORM) WRITE (6,405) WRITE (6,9002)PROB 9002 FORMAT(10X,23HTHIS IS PROBLEM NUMBER A6//) WRITE (6,9003)NCASES,IP,IQ,CRITER,NVAR,NVG,KVR 9003 FORMAT(10X,18HNUMBER OF CASES = I4//10X,35HNUMBER OF VARIABLES IN 1FIRST SET = I2// 210X,36HNUMBER OF VARIABLES IN SECOND SET = I2// 310X,75HONLY THE CANONICAL COEFFICIENTS WITH CORRESPONDING CORRELAT 4ION LARGER THAN F6.5,13H ARE COMPUTED// 510X,30HNUMBER OF VARIABLES READ IN = I3// 610X,35HNUMBER OF TRANS-GENERATION CARDS = I2// 710X, 39HNUMBER OF VARIABLE FORMAT CARDS USED = I2//) IF(INFORM-5)901,902,901 902 WRITE (6,9005) 9005 FORMAT(10X,22HINPUT DATA IS BY CARDS//) GO TO 531 C 440 WRITE (6,441) NTAPE 441 FORMAT (19H ILLEGAL TAPE UNIT ,I3,22H LISTED ON PROBLM CARD ) GO TO 401 C 901 WRITE (6,9006)INFORM 9006 FORMAT(10X,47HINPUT DATA IS BY TAPE---LOGICAL TAPE NUMBER IS I2//) 531 IPQ=IP+IQ IF(IQ-35) 532,532,910 532 IF(IP-IQ)533,533,920 533 IA=IP+1 N=IP IF(NVAR*(NVAR-101)) 534,535,535 535 WRITE (6,536)NVAR 536 FORMAT(1H0,20X,65HNUMBER OF ENTERING VARIABLES (T) EXCEEDS PROBLEM 1 LIMITS. YOU HAVE,I5,1H.) GO TO 401 534 CASES=NCASES MERRY=0 LCASE=0 IF(KVR.GT.0.AND.KVR.LE.10) GO TO 540 WRITE(6,6000) 6000 FORMAT(1H0,23X,71HNUMBER OF VARIABLE FORMAT CARDS INCORRECTLY SPEC XIFIED, ASSUMED TO BE 1.) KVR = 1 540 CONTINUE KVR=KVR*18 READ (5,5003)(FMT(I),I=1,KVR) 5003 FORMAT(18A4) WRITE (6,5004) (FMT(I),I=1,KVR) 5004 FORMAT (1X,18HVARIABLE FORMAT IS /,(1H ,18A4)) C 25 DO 4 I=1,IPQ SX(I)=0.0 SX2(I)=0.0 DO 4 J=1,IPQ 4 RHO(I,J)=0.0 IF(NVG) 29,29,28 28 DO 24 I=1,NVG READ (5,5002)XODE,NNEWA(I),LLCODE(I),LLVA(I),BBNEW(I) 5002 FORMAT(A6,I3,I2,I3,F6.0) IF(XODE.EQ.C123) GO TO 23 22 NVG=-101 GO TO 24 23 IF(LLCODE(I)-41)232,24,22 232 IF(LLCODE(I)*(LLCODE(I)-18))24,22,22 24 CONTINUE IF(NVG)234,29,29 234 J=I WRITE (6,5073) XODE,NNEWA(J),LLCODE(J),LLVA(J),BBNEW(J) 5073 FORMAT (1X,A6,I4,I3,I4,F7.0) WRITE (6,4003) 4003 FORMAT (49H ABOVE TRNGEN CARD ERROR. PROGRAM CANNOT PROCEED ) GO TO 401 29 IF((+Q003HL).NE.PRINTD) GO TO 21 295 WRITE (6,1035)PROB 1035 FORMAT (1H1,9X,22HINPUT DATA OF PROBLEM A6///2X,8HCASE NO./) 21 DO 150 II=1,NCASES H=II READ (INFORM,FMT)(DATA(I),I=1,NVAR) IF((+Q003HL).NE.PRINTD) GO TO 33 325 WRITE (6,1040)II,(DATA(I),I=1,NVAR) 1040 FORMAT(1H0,I6,2X,10F11.6/(9X,10F11.6)) 33 IF(NVG) 702,702,34 34 CALL TRNGEN (DATA,IPQ,NVG,NCASES,LCASE,NNEWA,LLCODE,LLVA,BBNEW,MER 1RY,II) IF(LCASE) 35,36,36 35 LCASE=0 GO TO 150 C 36 WRITE (IT2) (DATA(I),I=1,IPQ) 702 DO 8 I=1,IPQ SX(I)=SX(I)+DATA(I) SX2(I)=DATA(I)-SX(I)/H IF(H-1.0) 718,728,718 728 TT=0. GO TO 738 718 TT=SX2(I)*H/(H-1.0) IF(H.EQ.1.) TT=0. 738 DO 8 J=1,I 8 RHO(I,J)=RHO(I,J)+TT*SX2(J) 150 CONTINUE 5759 DO 7589 I=1,IPQ DO 7589 J=1,I 7589 RHO(J,I)=RHO(I,J) IF(NVG)706,706,707 707 REWIND IT2 WRITE (6,1045)IP, PROB,IQ 1045 FORMAT(1H1,10X,38HDATA AFTER TRANSGENERATION, THE FIRST I2, 1 42H BELONG TO THE FIRST SET IN THIS PROBLEM (A6,1H),/ 2 10X,10HTHE OTHER I2, 25H BELONG TO THE SECOND SET//) DO 40 II=1,NCASES READ (IT2) (DATA(I),I=1,IPQ) 40 WRITE (6,1040)II,(DATA(I),I=1,IPQ) REWIND IT2 706 WRITE (6,1100)PROB 1100 FORMAT(1H1,10X,32HSTANDARD DEVIATIONS FOR PROBLEM ,A6/1H0,10X, 1 8HVARIABLE,3X,8HST. DEV./1H0) DO 11 I=1,IPQ DATA(I)=SQRT(RHO(I,I)/(H-1.0)) 11 WRITE (6,1101)I,DATA(I) 1101 FORMAT(10X,I6,F15.6) DO 12 I=1,IPQ DO 12 J=1,IPQ 12 RHO(I,J)=RHO(I,J)/((H-1.0)*DATA(I)*DATA(J)) WRITE (6,1018)PROB 1018 FORMAT(1H1,10X,30HCORRELATION MATRIX OF PROBLEM A6///) DO 140 I=1,IPQ WRITE (6,1017)I,(RHO(I,J),J=1,IPQ) 1017 FORMAT(4H0ROWI3,3X,15F8.4/(10X,15F8.4)) 140 WRITE (IT3) (RHO(I,J),J=1,IPQ) REWIND IT3 DO 205 I=1,IP DO 205 J=1,IP 205 D(I,J)=RHO(I,J) DO 210 II=IA,IPQ I=II-IP DO 210 JJ=IA,IPQ J=JJ-IP 210 A(I,J)=RHO(II,JJ) CALL INV(A,A,IQ,1.E-7,0) DO 215 I=1,IP DO 215 JJ=IA,IPQ J=JJ-IP 215 E(I,J)=RHO(I,JJ) DO 230 I=1,IQ DO 230 J=1,IP G(I,J)=0.0 DO 230 K=1,IQ 230 G(I,J)=G(I,J)+A(I,K)*E(J,K) DO 235 I=1,IP DO 235 J=1,IP A(I,J)=0.0 DO 235 K=1,IQ 235 A(I,J)=A(I,J)+E(I,K)*G(K,J) DO 240 I=1,IP DO 240 J=1,IP 240 RHO(I,J)=A(I,J) CALL INV(D,A,IP,1.E-7,1) CALL HDIAG(A,IP,1,G,NNRR) DO 225 I=1,IP 225 C(I)=A(I,I) WRITE (6,1021)PROB 1021 FORMAT(1H1,10X,34HCANONICAL CORRELATIONS OF PROBLEM A6///) IPP=IP-1 DO 250 I=1,IPP K=I+1 DO 250 J=K,IP IF(C(I)-C(J))255,250,250 255 TEMP=C(I) C(I)=C(J) C(J)=TEMP 250 CONTINUE DO 260 I=1,IP IF(C(I))251,252,252 251 WRITE (6,4002)I,C(I) SX(I)=CRITER-1.0 GO TO 260 252 SX(I)=SQRT(C(I)) 260 WRITE (6,1022)I,SX(I) 1022 FORMAT(10X,I2,10X,F10.7) DO 265 I=1,IPQ 265 READ (IT3) (RHO(I,J),J=1,IPQ) REWIND IT3 WRITE (6,1023)CRITER 1023 FORMAT(1H1,10X,77HCANONICAL COEFFICIENTS---ONLY THOSE WITH CORRESP 1ONDING CANONICAL CORRELATIONS/ 2 38X, 13H LARGER THAN F6.4, 13H ARE COMPUTED/10X,80(1H.)) WRITE (6,2023) 2023 FORMAT(10X,80(1H.)) DO 300 I=1,IP DO 300 J=1,IP 300 A(I,J)=RHO(I,J) N=IPQ DO 777 L=1,IP IF(SX(L)-CRITER) 7777,7777,305 305 WRITE (6,1024)SX(L) 1024 FORMAT(1H0,5X,24HCANONICAL CORRELATION = F10.8) WRITE (6,2024) 2024 FORMAT(5X,34(1H*)//) DO 307 I=1,IP DO 307 J=IA,IPQ RHO(I,J)=-RHO(I,J)/SX(L) 307 RHO(J,I)=RHO(I,J) CALL VCTR (RHO,V,N,0.0) ZZZ=0.0 DO 501 I=1,IP DO 501 J=1,IP 501 ZZZ=ZZZ+V(I)*A(I,J)*V(J) ZZZ=1.0/SQRT(ZZZ) DO 502 I=1,IPQ V(I)=V(I)*ZZZ 502 QQQ(I,L)=V(I) WRITE (6,1025) 1025 FORMAT(10X,43HCOEFFICIENTS FOR THE FIRST SET OF VARIABLES//) WRITE (6,1026)(V(I),I=1,IP) 1026 FORMAT(8(5X,F10.6)) WRITE (6,1027) 1027 FORMAT(1H0/1H0,9X,45HCOEFFICIENTS FOR THE SECOND SET OF VARIABLES 1//) WRITE (6,1026)(V(I),I=IA,IPQ) WRITE (6,2000) 2000 FORMAT(3(1H0/)) DO 776 I=1,IPQ 776 READ (IT3) (RHO(I,J),J=1,IPQ) 8000 FORMAT (20A4) REWIND IT3 777 CONTINUE L=IP+1 7777 L2=L-1 CALL ABC(RHO,QQQ,QQQ,D,IP,IP,L2,RHO) CALL ABC(RHO(IP+1,1),QQQ,QQQ(IP+1,1),D,IP,IQ,L2,RHO(L,1)) CALL ABC(RHO(IP+1,IP+1),QQQ(IP+1,1),QQQ(IP+1,1),D,IQ,IQ,L2,RHO(L,L 1)) L1(1)=2*L2 LL1 = L1(1) DO 7778 I=L,LL1 DO 7778 J=1,L2 7778 RHO(J,I)=RHO(I,J) WRITE(6,7779) 7779 FORMAT(78H1CANONICAL CHECK MATRIX - ONLY 0., 1. AND CANONICAL CORR 1ELATIONS SHOULD APPEAR) DO 7780 I=1,LL1 7780 WRITE (6,1019) I,(RHO(I,J),J=1,LL1) 1019 FORMAT (4H0ROW,I3,3X,15F7.3,/(10X,15F7.3)) GO TO 999 910 WRITE (6,4000)IQ 915 WRITE (6,9004) IF(KVR.GT.0.AND.KVR.LE.10) GO TO 950 WRITE(6,6000) KVR = 1 950 CONTINUE KVR=KVR*18 READ (5,5003)(FMT(I),I=1,KVR) WRITE (6,5004) (FMT(I),I=1,KVR) IF(-NVG)916,918,918 916 DO 917 I=1,NVG 917 READ (5,5002)KODE 918 DO 919 I=1,NCASES 919 READ (INFORM,FMT)(DATA(J),J=1,NVAR) GO TO 999 4000 FORMAT(1H0,34X,50HMAXIMUM NUMBER OF VARIABLES (Q) EXCEEDED. YOU HA 1VE,I5,1H.) 920 WRITE (6,4001)IP,IQ GO TO 915 4001 FORMAT(1H0,9X,41HTHE NUMBER OF VARIABLES IN THE FIRST SET,,I4,47H, 1 IS GREATER THAN THE NUMBER IN THE SECOND SET,,I4,1H.) 4002 FORMAT(1H028X28HCANONICAL CORRELATION NUMBERI3,20H SQUARED IS EQUA 1L TO F10.8/1H019X80HTHE PROGRAM WILL CONTINUE AND DETERMINE COEFFI 2CIENTS FOR THE VALID CORRELATIONS.) 9004 FORMAT(1H ,31X,54HPROGRAM WILL GO TO NEXT PROBLEM (IF ANY) OR TERM 1INATE.) END C SUBROUTINE ABC FOR BMD06M JUNE 30, 1966 SUBROUTINE ABC(R,A,C,D,IP,IQ,L,S) DIMENSION R(70,70),A(70,35),C(70,35),D(35,35),S(70 ,70) DO 1 I=1,IQ DO 1 J=1,L D(I,J)=0.0 DO 1 K=1,IP 1 D(I,J)=D(I,J)+R(I,K)*A(K,J) DO 2 I=1,L DO 2 J=1,L S(I,J)=0.0 DO 2 K=1,IQ 2 S(I,J)=S(I,J)+C(K,J)*D(K,I) RETURN END C SUBROUTINE HDIAG FOR BMD06M JUNE 30, 1966 C MIHDI3, FORTRAN II DIAGONALIZATION OF A REAL SYMMETRIC MATRIX BY C THE JACOBI METHOD. C SIFTED TO FORTRAN IV C CALLING SEQUENCE FOR DIAGONALIZATION C CALL HDIAG( H, N, IEGEN, U, NR) C WHERE H IS THE ARRAY TO BE DIAGONALIZED. C N IS THE ORDER OF THE MATRIX, H. C C IEGEN MUST BE SET UNEQUAL TO ZERO IF ONLY EIGENVALUES ARE C TO BE COMPUTED. C IEGEN MUST BE SET EQUAL TO ZERO IF EIGENVALUES AND EIGENVECTORS C ARE TO BE COMPUTED. C C U IS THE UNITARY MATRIX USED FOR FORMATION OF THE EIGENVECTORS. C C NR IS THE NUMBER OF ROTATIONS. C C A DIMENSION STATEMENT MUST BE INSERTED IN THE SUBROUTINE. C DIMENSION H(N,N), U(N,N), X(N), IQ(N) C C SUBROUTINE PLACES COMPUTER IN THE FLOATING TRAP MODE C C THE SUBROUTINE OPERATES ONLY ON THE ELEMENTS OF H THAT ARE TO THE C RIGHT OF THE MAIN DIAGONAL. THUS, ONLY A TRIANGULAR C SECTION NEED BE STORED IN THE ARRAY H. C C SUBROUTINE HDIAG (H,N,IEGEN,U,NR) DIMENSION H(35,35),U(35,35),X(35),IQ(35) C EXTERNAL SIGN C IF (IEGEN) 15,10,15 10 DO 14 I=1,N DO 14 J=1,N IF(I-J)12,11,12 11 U(I,J)=1.0 GO TO 14 12 U(I,J)=0.0 14 CONTINUE C 15 NR = 0 IF (N-1) 1000,1000,17 C C SCAN FOR LARGEST OFF DIAGONAL ELEMENT IN EACH ROW C X(I) CONTAINS LARGEST ELEMENT IN ITH ROW C IQ(I) HOLDS SECOND SUBSCRIPT DEFINING POSITION OF ELEMENT C 17 NMI1=N-1 DO 30 I=1,NMI1 X(I) = 0.0 IPL1=I+1 DO 30 J=IPL1,N IF ( X(I) - ABS( H(I,J))) 20,20,30 20 X(I)=ABS(H(I,J)) IQ(I)=J 30 CONTINUE C C SET INDICATOR FOR SHUT-OFF.RAP=2**-27,NR=NO. OF ROTATIONS RAP=.745058059E-08 HDTEST=1.0E38 C C FIND MAXIMUM OF X(I) S FOR PIVOT ELEMENT AND C TEST FOR END OF PROBLEM C 40 DO 70 I=1,NMI1 IF (I-1) 60,60,45 45 IF ( XMAX- X(I)) 60,70,70 60 XMAX=X(I) IPIV=I JPIV=IQ(I) 70 CONTINUE C C IS MAX. X(I) EQUAL TO ZERO, IF LESS THAN HDTEST, REVISE HDTEST IF ( XMAX) 1000,1000,80 80 IF (HDTEST) 90,90,85 85 IF (XMAX - HDTEST) 90,90,148 90 HDIMIN = ABS( H(1,1) ) DO 110 I= 2,N IF (HDIMIN- ABS( H(I,I))) 110,110,100 100 HDIMIN=ABS(H(I,I)) 110 CONTINUE C HDTEST=HDIMIN*RAP C C RETURN IF MAX.H(I,J)LESS THAN(2**-27)ABSF(H(K,K)-MIN) IF (HDTEST- XMAX) 148,1000,1000 148 NR = NR+1 C C COMPUTE TANGENT, SINE AND COSINE,H(I,I),H(J,J) 150 TANG=SIGN(2.0,(H(IPIV,IPIV)-H(JPIV,JPIV)))*H(IPIV,JPIV)/(ABS(H(IPI 1V,IPIV)-H(JPIV,JPIV))+SQRT((H(IPIV,IPIV)-H(JPIV,JPIV))**2+4.0*H(IP 2IV,JPIV)**2)) COSINE=1.0/SQRT(1.0+TANG**2) SINE=TANG*COSINE HII=H(IPIV,IPIV) H(IPIV,IPIV)=COSINE**2*(HII+TANG*(2.0*H(IPIV,JPIV)+TANG*H(JPIV,JPI 1V))) H(JPIV,JPIV)=COSINE**2*(H(JPIV,JPIV)-TANG*(2.0*H(IPIV,JPIV)-TANG*H 1 II)) H(IPIV,JPIV)=0.0 C C PSEUDO RANK THE EIGENVALUES C ADJUST SINE AND COS FOR COMPUTATION OF H(IK) AND U(IK) IF ( H(IPIV,IPIV) - H(JPIV,JPIV)) 152,153,153 152 HTEMP = H(IPIV,IPIV) H(IPIV,IPIV) = H(JPIV,JPIV) H(JPIV,JPIV) = HTEMP C RECOMPUTE SINE AND COS HTEMP = SIGN(1.0, -SINE) * COSINE COSINE = ABS(SINE) SINE = HTEMP 153 CONTINUE C C INSPECT THE IQS BETWEEN I+1 AND N-1 TO DETERMINE C WHETHER A NEW MAXIMUM VALUE SHOULD BE COMPUTED SINCE C THE PRESENT MAXIMUM IS IN THE I OR J ROW. C DO 350 I=1,NMI1 IF(I-IPIV)210,350,200 200 IF(I-JPIV)210,350,210 210 IF(IQ(I)-IPIV)230,240,230 230 IF(IQ(I)-JPIV)350,240,350 240 K=IQ(I) 250 HTEMP=H(I,K) H(I,K)=0.0 IPL1=I+1 X(I) =0.0 C C SEARCH IN DEPLETED ROW FOR NEW MAXIMUM C DO 320 J=IPL1,N IF ( X(I)- ABS( H(I,J)) ) 300,300,320 300 X(I) = ABS(H(I,J)) IQ(I)=J 320 CONTINUE H(I,K)=HTEMP 350 CONTINUE C X(IPIV) =0.0 X(JPIV) =0.0 C C CHANGE THE OTHER ELEMENTS OF H C DO 530 I=1,N C IF(I-IPIV)370,530,420 370 HTEMP = H(I,IPIV) H(I,IPIV) = COSINE*HTEMP + SINE*H(I,JPIV) IF ( X(I) - ABS( H(I,IPIV)) )380,390,390 380 X(I) = ABS(H(I,IPIV)) IQ(I) = IPIV 390 H(I,JPIV) = -SINE*HTEMP + COSINE*H(I,JPIV) IF ( X(I) - ABS( H(I,JPIV)) ) 400,530,530 400 X(I) = ABS(H(I,JPIV)) IQ(I) = JPIV GO TO 530 C 420 IF(I-JPIV)430,530,480 430 HTEMP = H(IPIV,I) H(IPIV,I) = COSINE*HTEMP + SINE*H(I,JPIV) IF ( X(IPIV) - ABS( H(IPIV,I)) ) 440,450,450 440 X(IPIV) = ABS(H(IPIV,I)) IQ(IPIV) = I 450 H(I,JPIV) = -SINE*HTEMP + COSINE*H(I,JPIV) IF ( X(I) - ABS( H(I,JPIV)) ) 400,530,530 C 480 HTEMP = H(IPIV,I) H(IPIV,I) = COSINE*HTEMP + SINE*H(JPIV,I) IF ( X(IPIV) - ABS( H(IPIV,I)) ) 490,500,500 490 X(IPIV) = ABS(H(IPIV,I)) IQ(IPIV) = I 500 H(JPIV,I) = -SINE*HTEMP + COSINE*H(JPIV,I) IF ( X(JPIV) - ABS( H(JPIV,I)) ) 510,530,530 510 X(JPIV) = ABS(H(JPIV,I)) IQ(JPIV) = I 530 CONTINUE C C TEST FOR COMPUTATION OF EIGENVECTORS C IF(IEGEN)40,540,40 540 DO 550 I=1,N HTEMP=U(I,IPIV) U(I,IPIV)=COSINE*HTEMP+SINE*U(I,JPIV) 550 U(I,JPIV)=-SINE*HTEMP+COSINE*U(I,JPIV) GO TO 40 1000 RETURN END C SUBROUTINE INV FOR BMD06M JUNE 30, 1966 SUBROUTINE INV(A,B,N,T,KD) DIMENSION U(35),A(35,35),B(35,35) DIMENSION IN(35) DO 1 I=1,N 1 IN(I)=0 K=1 2 DO 3 I=K,N U(I)=A(I,K) 3 A(I,K)=0.0 P=U(K) DO 4 I=1,K U(I)=A(K,I) 4 A(K,I)=0.0 U(K)=-1.0 IN(K)=1 H=0.0 DO 5 I=1,N Y=U(I)/P DO 6 J=1,I 6 A(I,J)=A(I,J)-U(J)*Y IF(IN(I))5,7,5 7 IF(KD)12,12,13 13 Z=B(I,I)-Y*(2.0*B(I,K)-Y*B(K,K)) DO 8 J=1,N B(I,J)=B(I,J)-Y*B(K,J) 8 B(J,I)=B(I,J) B(I,I)=Z 12 IF(A(I,I)-H)5,5,9 9 H=A(I,I) KK=I 5 CONTINUE IF(KD)20,20,21 21 Z=B(K,K)/P P=SQRT(P) DO 14 I=1,N B(I,K)=B(I,K)/P 14 B(K,I)=B(I,K) B(K,K)=Z 20 K=KK IF(H-T)10,10,2 10 DO 91 I=1,N IF(IN(I))90,90,81 90 DO 92 J=1,I A(I,J)=0.0 A(J,I)=0.0 IF(KD)92,92,93 93 B(I,J)=0.0 B(J,I)=0.0 92 CONTINUE GO TO 91 81 DO 11 J=1,I A(I,J)=-A(I,J) 11 A(J,I)=A(I,J) 91 CONTINUE RETURN END C SUBROUTINE TPWD FOR BMD06M JUNE 30, 1966 SUBROUTINE TPWD(NT1,NT2) IF(NT1)40,10,12 10 NT1=5 12 IF(NT1-NT2)14,19,14 14 IF(NT2.EQ.5) GO TO 19 REWIND NT2 19 IF(NT1-5)18,24,18 18 IF(NT1-6)22,40,22 22 REWIND NT1 24 NT2=NT1 28 RETURN 40 WRITE (6,49) STOP 49 FORMAT(25H ERROR ON TAPE ASSIGNMENT) END C SUBROUTINE TRNGEN FOR BMD06M JUNE 30, 1966 SUBROUTINE TRNGEN (DATA,IPQ,NVG,NCASES,LCASE,NNEWA,LLCODE,LLVA,BBN 1EW,MERRY,II) DIMENSION DATA(100),NNEWA(99),LLCODE(99),LLVA(99),BBNEW(99) C EXTERNAL SIGN REAL LOG10 ASN(Q000FL)=ATAN(Q000FL/SQRT(1.0-Q000FL**2)) SAMP=NCASES ITEM=II DO 3 J=1,NVG 305 NEWA=NNEWA(J) LCODE=LLCODE(J) LVA=LLVA(J) BNEW=BBNEW(J) IF(LCODE-10) 4,4,5 5 NEWB=BNEW 4 D=DATA(LVA) IF(LCODE-41)500,180,500 500 GO TO (10,20,30,40,50,60,70,80,90,100,110,120,130,140, 1 150,160,170),LCODE 10 IF(D)99,7,8 7 DATA(NEWA)=0.0 GO TO 3 8 DATA(NEWA)=SQRT(D) GO TO 3 20 IF(D)99,11,12 11 DATA(NEWA)=1.0 GO TO 3 12 DATA(NEWA)=SQRT(D)+SQRT(D+1.0) GO TO 3 30 IF(D)99,99,14 14 DATA(NEWA)=ALOG10(D) GO TO 3 40 DATA(NEWA)=EXP(D) GO TO 3 50 IF(D)99,16,17 16 DATA(NEWA)=0.0 GO TO 3 17 IF(D-1.0)18,19,99 19 DATA(NEWA)=3.14159265/2.0 GO TO 3 18 A=SQRT(D) DATA(NEWA)=ASN(A) GO TO 3 60 A=D/(SAMP+1.0) B=A+1.0/(SAMP+1.0) IF(A)99,23,24 23 IF(B)99,26,27 26 DATA(NEWA)=0.0 GO TO 3 27 DATA(NEWA)=ASN(SQRT(B)) GO TO 3 24 IF(B)99,28,29 28 DATA(NEWA)=ASN(SQRT(A)) GO TO 3 29 A=SQRT(A) B=SQRT(B) DATA(NEWA)=ASN(A)+ASN(B) GO TO 3 70 IF(D)31,99,31 31 DATA(NEWA)=1.0/D GO TO 3 80 DATA(NEWA)=D+BNEW GO TO 3 90 DATA(NEWA)=D*BNEW GO TO 3 100 IF(D)33,32,33 32 DATA(NEWA)=0.0 GO TO 3 33 DATA(NEWA)=D**BNEW GO TO 3 110 DATA(NEWA)=D+DATA(NEWB) GO TO 3 120 DATA(NEWA)=D-DATA(NEWB) GO TO 3 130 DATA(NEWA)=D*DATA(NEWB) GO TO 3 140 IF(DATA(NEWB))34,99,34 34 DATA(NEWA)=D/DATA(NEWB) 150 BNEW=NEWB IF(D-BNEW) 151,152,152 152 DATA(NEWA)=1.0 GO TO 3 151 DATA(NEWA)=0.0 GO TO 3 160 IF(D-DATA(NEWB)) 161,162,162 162 DATA(NEWA)=1.0 GO TO 3 161 DATA(NEWA)=0.0 170 IF(-D)175,99,99 175 DATA(NEWA)=ALOG(D) GO TO 3 180 IF(D)3,503,3 503 IF(SIGN(10.0,D)) 504,3,3 504 DATA(NEWA)=BNEW 3 CONTINUE GO TO 42 99 LCASE=-999 IF(MERRY-J) 402,401,402 402 MERRY=J WRITE (6,1404)J 401 WRITE (6,1405)ITEM WRITE (6,1408) NCASES=NCASES-1 1408 FORMAT(45H0THIS CASE WILL BE DELETED FOR ALL VARIABLES ) 1404 FORMAT(30H0THE INSTRUCTIONS INDICATED ON/25H TRANS GENERATOR CARD 1NO.I2,4H RE-/29H SULTED IN THE VIOLATION OF A/31H RESTRICTION FOR 2THIS TRANSFOR-/31H MATION. THE VIOLATION OCCURRED/27H FOR THE CASE 3 LISTED BELOW./) 1405 FORMAT( 9H CASE NO.I5) 42 RETURN END C SUBROUTINE VCTR FOR BMD06M JUNE 30, 1966 SUBROUTINE VCTR(A,V,N,ALPHA) DIMENSION A(70,70),V(70) A(1,1)=A(1,1)-ALPHA 6 DO 15 I=2,N A(I,I)=A(I,I)-ALPHA 70 II=I-1 7 DO 15 J=1,II 8 IF (A(I,J))9,15,9 9 IF (ABS(A(J,J))-ABS(A(I,J)))11,10,10 10 R=A(I,J)/A(J,J) GO TO 130 11 R=A(J,J)/A(I,J) DO 12 K=1,N C=A(J,K) A(J,K)=A(I,K) 12 A(I,K)=C 130 JJ=J+1 13 DO 14 K=JJ,N 14 A(I,K)=A(I,K)-R*A(J,K) 15 CONTINUE C=A(N,N) V(N)=1.0 DO 29 I=2,N JJ=N-I+1 R=0.0 II=N-I+2 DO 25 K=II,N 25 R=R+A(JJ,K)*V(K) IF (ABS(A(JJ,JJ))-1.0E-10)27,27,28 27 V(JJ)=1.0 C=0.0 DO 26 J=II,N 26 V(J)=0.0 GO TO 29 28 V(JJ)=(C-R)/A(JJ,JJ) 29 CONTINUE RETURN END