Trailing-Edge
-
PDP-10 Archives
-
decuslib20-05
-
decus/20-0137/bmd/bmd06m.for
There is 1 other file named bmd06m.for in the archive. Click here to see a list.
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