Trailing-Edge
-
PDP-10 Archives
-
decuslib20-05
-
decus/20-0137/bmd/bmd02r.for
There is 1 other file named bmd02r.for in the archive. Click here to see a list.
00100 C STEPWISE REGRESSION - MAIN PROGRAM MAY 2, 1966
00200 C
00300 C THIS PROGRAM IS A SIFTED VERSION OF THE ORIGINAL FORTRAN II
00400 C PROGRAM, BMD02R WITH SOME MODIFICATIONS TO MAKE IT OPERABLE.
00500 C IT WAS THEN CONVERTED TO 360 FORTRAN IV (H-LEVEL)
00600 C
00700 C PROBLEM CARD FORMAT
00800 C COL NAME
00900 C 1- 6 XMAN PROBLM
01000 C 10-15 CODE ALPHANUMERIC PROBLEM NAME
01100 C 17-20 N NUMBER OF CASES
01200 C 24-25 NOV NUMBER OF ORIGINAL VARIABLES
01300 C 29-30 NTGC NUMBER OF TRANSGENERATION CARDS
01400 C 34-35 NVA NUMBER OF VARIABLES ADDED BY TRANSGENERATION
01500 C 39-40 NAIT ALTERNATE INPUT TAPE NUMBER
01600 C 44-45 NSPC NUMBER OF SUBPROBLEM CARDS
01700 C 48-49 NLV NUMBER OF LABELED VARIABLES
01800 C 51-53 SDAM YES IF ST. DEV. AND MEANS TO BE PRINTED
01900 C 55-57 COVP YES IF COVARIANCE MATRIX TO BE PRINTED
02000 C 59-61 CORP YES IF CORRELATION MATRIX TO BE PRINTED
02100 C 63-65 ZEROI YES IF ZERO REGRESSION INTERCEPT DESIRED
02200 C 67-69 WIND NO IF ALT. INPUT TAPE NOT TO BE REWOUND
02300 C 71-72 NVFC NUMBER OF VARIABLE FORMAT CARDS
02400 C
02500 C SUB-PROBLEM CARD FORMAT
02600 C
02700 C COL NAME
02800 C 1- 6 WMAN SUBPRO
02900 C 9-10 KDEP DEPENDENT VARIABLE NUMBER
03000 C 13-15 MAXSTP MAXIMUM NUMBER OF STEPS
03100 C 20-25 FINC F FOR INCLUSION
03200 C 30-35 FOUT F FOR DELETION
03300 C 40-45 TOL TOLERANCE
03400 C 49-50 NVIP NUMBER OF VARIABLES TO BE PLOTTED
03500 C 53-55 CDF YES IF CONTROL DELETE CARDS ARE INCLUDED
03600 C 58-60 RESID YES IF RESIDUALS ARE TO BE PRINTED
03700 C 63-65 SUMTAB YES IF SUMMARY TABLE DESIRED
03800 C
03900 COMMON NTGC,TRANS(100),KTRANS(3,100),X(81),A(50,50),IP,FINC,RESDF
04000 1,FOUT,KAY,FLAG,N,NINCS,KDEP,IR,P(50,50),Q(50,50),XN,TOL,DF,LLL(10)
04100 2,B(80),COL(11),PINT(11),QINT(11),RES(275),NVIP,ALPHA,RMAX,RMIN,
04200 3RESID,IVPT(33),NVI,KEEP(5),IS
04300 C
04400 DIMENSION BNAME(80),ALBEL(80),ANAME(80),R(80),FE(80),NIEN(80)
04500 DIMENSION KP(50,50),IQ(50,50),STDEV(80),C(80)
04600 DIMENSION XMIN(80),XMAX(80),BES(5),NUSE(20)
04700 DIMENSION XMEAN(80),D(80),F(80),TOLEV(80),INEN(80)
04800 INTEGER KQ(80),JQ(80)
04900 EQUIVALENCE (P,KP),(Q,IQ),(KTRANS(1),XMIN),(KTRANS(81),XMAX),
05000 1(KTRANS(161),STDEV),(TRANS,C)
05100 C
05200 DOUBLE PRECISION XMAN,CODE,YMAN,ALBEL,WMAN,VMAN,Q009HL,Q010HL,ANAM
05300 1E,BNAME,ENTER,PROBLM,SUBPRO,TRNGEN,CONDEL,DXPLTS,FINISH,DUM
05400 DATA PROBLM,SUBPRO,TRNGEN,CONDEL,YES,DXPLTS,FINISH,XNO/'PROBLM',
05500 1'SUBPRO','TRNGEN','CONDEL','YES','IDXPLT','FINISH',' NO'/
05600 DATA Q009HL/6HREMOVE/
05700 DATA Q010HL/6HENTERE/
05800 DOUBLE PRECISION FIN,FOU,TO
05900 8003 FORMAT('1 BMD02R - STEPWISE REGRESSION - REVISED ',
06000 1'JUNE 26, 1969'/
06100 22X,40HHEALTH SCIENCES COMPUTING FACILITY, UCLA//2X,12HPROBLEM CODE
06200 317XA6/2X15HNUMBER OF CASES16XI4/2X28HNUMBER OF ORIGINAL VARIABLES,
06300 45X,I2/2X,25HNUMBER OF VARIABLES ADDED,8 X,I2/2X,25HTOTAL NUMBER OF
06400 5 VARIABLES, 8X,I2/2X,22HNUMBER OF SUB-PROBLEMS,11X,I2)
06500 MAIT=5
06600 CALL USAGEB('BMD02R')
06700 C
06800 C READ PROBLEM CARD
06900 C
07000 1 READ (5,8001)XMAN,CODE,N,NOV,NTGC,NVA,NAIT,NSPC,NLV,SDAM,COVP,CORP
07100 1,ZEROI,WIND,NVFC
07200 8001 FORMAT(A6,3X,A6,1X,I4,5(3X,I2),2X,I2,5(1X,A3),1X,I2)
07300 C
07400 C CHECK PROBLEM CARD FOR VALID PARAMETERS
07500 C
07600 IF(XMAN .EQ. FINISH) GO TO 9001
07700 2 IF(XMAN .NE. PROBLM) GO TO 9002
07800 3 IP = NOV+NVA
07900 WRITE (6,8003)CODE,N,NOV,NVA,IP,NSPC
08000 IF(-NSPC)355,9005,9005
08100 355 IF((NOV-1)*(NOV-81))4,9003,9003
08200 4 IF((IP-1)*(IP-81))5,9003,9003
08300 5 IF(NVFC.GT.0.AND.NVFC.LE.10)GO TO 106
08400 WRITE(6,105)
08500 NVFC=1
08600 105 FORMAT(1H023X71HNUMBER OF VARIABLE FORMAT CARDS INCORRECTLY SPECIF
08700 1IED, ASSUMED TO BE 1.)
08800 106 IF(NTGC) 108,108,7
08900 C
09000 C READ TRANSGENERATION CARDS
09100 7 DO 1000 I=1,NTGC
09200 READ (5,8002)YMAN,(KTRANS(J,I),J=1,3),TRANS(I)
09300 8002 FORMAT(A6,I3,I2,I3,F6.0)
09400 C
09500 C CHECK TRANSGENERATION PARAMETERS
09600 C
09700 IF(YMAN .NE. TRNGEN) GO TO 9006
09800 7777 KIP=IP
09900 IF(NVA)7780,7778,7778
10000 7780 KIP=NOV
10100 7778 IF(KTRANS(1,I)-KIP)7776,7776,9777
10200 7776 IF(KTRANS(3,I)-KIP)7775,7775,9777
10300 7775 IF(KTRANS(2,I)-24)7779,1000,9775
10400 7779 IF((KTRANS(2,I)-19)*(KTRANS(2,I)-18))1000,9775,1000
10500 9777 WRITE (6,8777)I
10600 8777 FORMAT(///51H VARIABLE NUMBER SPECIFIED ON TRANSGENERATION CARD ,
10700 1 I4,12H EXCEEDS P+Q )
10800 NAIT=-24
10900 GO TO 1000
11000 9775 WRITE (6,8775)KTRANS(2,I)
11100 8775 FORMAT(///29H ILLEGAL TRANSGENERATION CODE I3,10H SPECIFIED)
11200 NAIT=-24
11300 1000 CONTINUE
11400 108 IF(NAIT)9999,107,8
11500 107 NAIT=5
11600 GO TO 9
11700 8 IF(MAIT.NE.5.AND.MAIT.NE.NAIT)REWIND MAIT
11800 MAIT=NAIT
11900 IF(WIND.NE.XNO)REWIND NAIT
12000 C
12100 C READ LABEL CARDS
12200 C
12300 9 CALL RDLBL2(NLV,IP,ALBEL)
12400 C
12500 C READ VARIABLE FORMAT CARDS
12600 C
12700 10 NVFC=NVFC*18
12800 READ (5,8004)(RES(I),I=1,NVFC)
12900 8004 FORMAT(18A4)
13000 WRITE(6,8055) (RES(I), I=1,NVFC)
13100 8055 FORMAT(' THE VARIABLE FORMAT IS ',18A4/(1H ,25X,18A4))
13200 C
13300 C INITIALIZE ACCUMULATORS AND MATRIX A
13400 C
13500 XN=N
13600 DO 3000 I=1,IP
13700 XMEAN(I)=0.0
13800 DO 3000 J=1,IP
13900 A(I,J)=0.0
14000 3000 CONTINUE
14100 C
14200 DO 4000 K=1,N
14300 NINCS=K
14400 READ (NAIT,RES)(X(L),L=1,NOV)
14500 IF(NTGC) 11,11,12
14600 12 CALL TRANGN
14700 11 IF(ZEROI .NE. YES) GO TO 121
14800 122 DO 4100 I=1,IP
14900 XMEAN(I)=XMEAN(I)+X(I)
15000 DO 4100 J=I,IP
15100 A(I,J)=A(I,J)+X(I)*X(J)
15200 4100 CONTINUE
15300 XM=N
15400 GO TO 4300
15500 121 XK=K
15600 DO 4201 I=1,IP
15700 XMEAN(I)=XMEAN(I)+X(I)
15800 4201 CONTINUE
15900 DO 4200 I=1,IP
16000 DO 4200 J=I,IP
16100 IF((XK-1.0) .EQ. 0.0) GO TO 4200
16200 A(I,J)=A(I,J)+(XK*X(J)-XMEAN(J))*(XK*X(I)-XMEAN(I))/(XK*(XK-1.0))
16300 C
16400 4200 CONTINUE
16500 XM=N-1
16600 C
16700 C
16800 4300 WRITE(2)(X(I),I=1,IP)
16900 67890 FORMAT(20A4)
17000 C
17100 C
17200 C
17300 C
17400 4000 CONTINUE
17500 REWIND 2
17600 C
17700 C REPLACE XMEAN WITH MEAN VECTOR, A WITH COVARIANCE MATRIX,AND
17800 C COMPUTE STANDARD DEVIATIONS
17900 C
18000 DO 5000 I=1,IP
18100 XMEAN(I) = XMEAN(I) / XN
18200 DO 5000 J=I,IP
18300 A(I,J) = A(I,J) / XM
18400 5000 CONTINUE
18500 DO 5100 I=1,IP
18600 STDEV(I)=SQRT(A(I,I))
18700 5100 CONTINUE
18800 C
18900 C IF ZERO REGRESSION INTERCEPT IS DESIRED, PRINT WARNING.
19000 C
19100 IF(ZEROI .NE. YES) GO TO 5110
19200 5105 WRITE (6,5120)
19300 5120 FORMAT(1H0,6X, 92HWARNING...WHEN THE ZERO REGRESSION INTERCEPT IS
19400 1CHOSEN, ALL VARIANCES, COVARIANCES, STANDARD/17X84HDEVIATIONS AND
19500 2CORRELATIONS ARE COMPUTED ABOUT THE ORIGIN RATHER THAN ABOUT THE M
19600 3EAN/17X,34H(SEE PROGRAM WRITEUP - SECTION 4).)
19700 C
19800 C WRITE OUT MEANS AND STANDARD DEVIATIONS IF REQUESTED
19900 C
20000 5110 IF(SDAM .NE. YES) GO TO 15
20100 14 WRITE (6,8005)
20200 8005 FORMAT(///// 4X,8HVARIABLE,8X,4HMEAN,7X,18HSTANDARD DEVIATION)
20300 C
20400 16 DO 9000 I=1,IP
20500 WRITE (6,8006)ALBEL(I),I,XMEAN(I),STDEV(I)
20600 8006 FORMAT(3X,A6,I3,2X,F14.5,4X,F14.5)
20700 9000 CONTINUE
20800 C
20900 C
21000 C
21100 C
21200 C
21300 C
21400 C PRINT COVARIANCE MATRIX IF REQUESTED
21500 C
21600 15 IF(COVP .NE. YES) GO TO 19
21700 18 WRITE (6,8008)
21800 8008 FORMAT(19H1 COVARIANCE MATRIX)
21900 CALL AOUT
22000 C
22100 C REPLACE UPPER DIAGONAL SECTION OF MATRIX WITH CORRELATION MATRIX
22200 C
22300 19 DO 11000 I=1,IP
22400 DO 11000 J=I,IP
22500 AAA = STDEV(I) * STDEV(J)
22600 IF(AAA .EQ. 0.0) A(I,J) = 0.0
22700 IF(AAA .NE. 0.0)
22800 1A(I,J) = A(I,J) / AAA
22900 A(J,I)=A(I,J)
23000 11000 CONTINUE
23100 C
23200 C PRINT CORRELATION MATRIX IF REQUESTED
23300 C
23400 IF(CORP .NE. YES) GO TO 824
23500 23 WRITE (6,8009)
23600 8009 FORMAT(20H1 CORRELATION MATRIX)
23700 CALL AOUT
23800 824 IF(ZEROI .EQ. YES) GO TO 242
23900 241 XN=XN-1.0
24000 242 DO 13000 M=1,NSPC
24100 C
24200 C RESTORE THE CORRELATION MATRIX
24300 C
24400 DO 23000 I=1,IP
24500 A(I,I)=1.0
24600 K=I+1
24700 DO 23000 J=K,IP
24800 A(I,J)=A(J,I)
24900 23000 CONTINUE
25000 C
25100 C READ SUB PROBLEM CARD
25200 C
25300 READ(5,8010)WMAN,KDEP,MAXSTP,FIN,FOU,TO,NVIP,CDF,RESID,SUMTAB
25400 8010 FORMAT(A6,2X,I2,2X,I3,3(4X,A6),3X,I2,3(2X,A3))
25500 IF(WMAN .NE. SUBPRO) GO TO 9009
25600 FINC=0.01
25700 CALL ATOF(FIN,6,FINC)
25800 FOUT=0.005
25900 CALL ATOF(FOU,6,FOUT)
26000 TOL=0.001
26100 CALL ATOF(TO,6,TOL)
26200 331 IF(MAXSTP) 332,332,34
26300 332 MAXSTP=IP *2
26400 34 DO 14000 I=1,IP
26500 C(I)=0.0
26600 14000 CONTINUE
26700 WRITE (6,8034)M,KDEP,MAXSTP,FINC,FOUT,TOL
26800 8034 FORMAT(2H1 , 10HSUB-PROBLM,I5/3X,18HDEPENDENT VARIABLE,11X,I2/3X,
26900 123HMAXIMUM NUMBER OF STEPS,5X,I3/3X,21HF-LEVEL FOR INCLUSION,4X,
27000 2F8.6/3X,20HF-LEVEL FOR DELETION,5X,F8.6/3X,15HTOLERANCE LEVEL,10X,
27100 3F8.6)
27200 IF(CDF .NE. YES) GO TO 36
27300 35 IF(IP -66) 351,351,352
27400 351 READ (5,8011)VMAN,(C(I),I=1,IP)
27500 IF(VMAN .NE. CONDEL) GO TO 9010
27600 GO TO 36
27700 352 READ (5,8011)VMAN,(C(I),I=1,66),DUM,(C(J),J=67,IP)
27800 8011 FORMAT(A6,66F1.0)
27900 361 IF(VMAN.NE.CONDEL .OR. DUM.NE.CONDEL) GO TO 9010
28000 36 DO 15000 I=1,IP
28100 IF( C(I)) 37,37,15000
28200 37 C(I)=2.0
28300 15000 CONTINUE
28400 C(KDEP)=1.0
28500 DF=0.0
28600 LKL=0
28700 L=0
28800 39 L=L+1
28900 C CALL SUBROUTINE TO ENTER VARIABLE,CALCULATE VALUES TO BE PRINTED
29000 IF ((DF-XN).EQ.0.0) GO TO 1117
29100 CALL STEPRG
29200 IF(FLAG) 391,1117,392
29300 1117 WRITE (6,9118)
29400 9118 FORMAT(//58H F-LEVEL OR TOLERANCE INSUFFICIENT FOR FURTHER COMPUTA
29500 1TION)
29600 GO TO 117
29700 391 ENTER=(+Q009HL)
29800 GO TO 393
29900 392 ENTER=(+Q010HL)
30000 393 RESDF=XN-DF
30100 RESSS =XN*(STDEV(KDEP)**2)*A(KDEP,KDEP)
30200 RESMS = 0.0
30300 IF(RESDF .NE. 0.0)
30400 1RESMS = RESSS / RESDF
30500 REGDF=DF
30600 REGSS=XN*(STDEV(KDEP)**2)-RESSS
30700 REGMS = 0.0
30800 IF(REGDF .NE. 0.0)
30900 1REGMS = REGSS / REGDF
31000 FRATIO = 0.0
31100 STERR = 0.0
31200 IF(RESMS .GT. 0.0) FRATIO = REGMS / RESMS
31300 IF(RESMS .LT. (ABS(XMEAN(KDEP))*1.E-6)) RESMS = 0.0
31400 IF(RESMS .GT. 0.0) STERR = SQRT(RESMS)
31500 HOLD = 1.0 - A(KDEP,KDEP)
31600 XMULTR = 0.0
31700 IF(HOLD .GT. 0.0)
31800 1XMULTR = SQRT(HOLD)
31900 IDF=DF
32000 IRDF=RESDF
32100 WRITE (6,8012)L,ENTER,KAY,XMULTR,STERR,IDF,REGSS,REGMS,FRATIO,IRDF
32200 1,RESSS,RESMS
32300 8012 FORMAT(////4X,11HSTEP NUMBER,2X,I3/4X,9HVARIABLE ,A6,2HD ,I4//4X,
32400 110HMULTIPLE R,12X,F 9.4/4X,18HSTD. ERROR OF EST.,F13.4 // 4X,
32500 220HANALYSIS OF VARIANCE/27X,2HDF,4X,14HSUM OF SQUARES,4X,11HMEAN S
32600 3QUARE,4X,7HF RATIO/12X,10HREGRESSION,3X,I4,F16.3,F14.3,F12.3
32700 4/12X,8HRESIDUAL,5X,I4,F16.3,F14.3)
32800 C
32900 C A VARIABLE IS IN THE EQUATION IF C(I) IS LESS THAN OR EQUAL TO 0.0
33000 C
33100 NVI=0
33200 NVO=0
33300 ALPHA=XMEAN(KDEP)
33400 DO 16000 I=1,IP
33500 IF(I-KDEP) 441,16000,441
33600 441 IF(C(I)) 41,41,43
33700 C
33800 C COMPUTE MULTIPLE REGRESSION EQUATION COEFFICIENTS,STD.ERROR,
33900 C AND F TO REMOVE, FOR VARIABLES IN THE REGRESSION
34000 C
34100 41 NVI=NVI+1
34200 B(NVI) = 0.0
34300 IF(STDEV(I) .NE. 0.0)
34400 1B(NVI)=STDEV(KDEP)*AF(I,KDEP)/STDEV(I)
34500 D(NVI) = 0.0
34600 IF(STDEV(I).NE.0.0 .AND. A(I,I).LT.0.0)
34700 1D(NVI)=(STERR/STDEV(I))*SQRT(-A(I,I)/XN)
34800 F(NVI) = 0.0
34900 IF(D(NVI) .NE. 0.0)
35000 1F(NVI)=(B(NVI)/D(NVI))**2
35100 ALPHA=ALPHA-B(NVI)*XMEAN(I)
35200 ANAME(NVI)=ALBEL(I)
35300 KQ(NVI)=C(I)+9.0
35400 INEN(NVI)=I
35500 GO TO 16001
35600 C
35700 C A VARIABLE IS OUT OF THE REGRESSION IF C(I) IS GREATER THAN OR
35800 C EQUAL TO 1
35900 C
36000 C
36100 C
36200 C COMPUTE PARTIAL CORRELATION COEFFICIENTS, TOLERANCE, AND
36300 C F TO ENTER FOR VARIABLES OUT OF THE REGRESSION
36400 C
36500 43 NVO=NVO+1
36600 BNAME(NVO)=ALBEL(I)
36700 JQ(NVO)=C(I)
36800 NIEN(NVO)=I
36900 TOLEV(NVO)= A(I,I)
37000 R(NVO) = 0.0
37100 FE(NVO) = 0.0
37200 STORE = A(I,I) * A(KDEP,KDEP)
37300 IF(STORE .LE. 0.0) GO TO 16001
37400 R(NVO) = AF(I,KDEP) / SQRT(STORE)
37500 STORE = A(I,I)*A(KDEP,KDEP)-(AF(I,KDEP)**2)
37600 IF(STORE.GT.0.0 .AND. (RESDF-1.0).GT.0.0)
37700 1FE(NVO)=((AF(I,KDEP)**2)*(RESDF-1.0))/ STORE
37800 16001 IF(I-KAY) 16000,16002,16000
37900 16002 IF(C(I)) 16003,16003,16004
38000 16003 FKAY=F(NVI)
38100 GO TO 16000
38200 16004 FKAY=FE(NVO)
38300 16000 CONTINUE
38400 IF(ZEROI .NE. YES) GO TO 443
38500 442 ALPHA=0.0
38600 C
38700 C WRITE HEADING FOR COEFFICIENTS
38800 C
38900 443 WRITE(6,8013)ALPHA
39000 8013 FORMAT(/60X,'.'/21X,'VARIABLES IN EQUATION',18X,'.',19X,'VARIABLES
39100 1 NOT IN EQUATION'/60X,'.'/6X,'VARIABLE',6X,'COEFFICIENT STD. ERRO
39200 2R F TO REMOVE . VARIABLE PARTIAL CORR. TOLERANCE
39300 3 F TO ENTER'/2(60X,'.'/),6X,'(CONSTANT',5X,F11.5,' )',27X,'.')
39400 C
39500 WRITE(1)L,KAY,FLAG,XMULTR,FKAY,NVI
39600 LKL=LKL+1
39700 C
39800 C PRINT THE REGRESSION ANALYSIS TABLE
39900 C
40000 C
40100 C
40200 C
40300 C
40400 C
40500 IF(FLAG.GT.0.0)CALL WHICHX(KAY,NVI,IS,KEEP)
40600 473 NGO=0
40700 44 IF(NVO)46,46,45
40800 45 IF(NVI.LE.0)GO TO 56
40900 47 LNV=MIN0(NVI,NVO)
41000 C
41100 C NVO AND NVI BOTH POSITIVE,PRINT BOTH SIDES OF TABLE
41200 C
41300 C
41400 49 DO 17000 I=1,LNV
41500 WRITE(6,8014)ANAME(I),INEN(I),B(I),D(I),F(I),KQ(I),BNAME(I),NIEN(I
41600 1),R(I),TOLEV(I),FE(I),JQ(I)
41700 8014 FORMAT(5X,A6,1X,I2,1X,F16.5,1X,F11.5,1X,F11.4,' (',I1,') . ',A6
41800 1,1X,I2,1X,F15.5,1X,F13.4,1X,F12.4,' (',I1,')')
41900 17000 CONTINUE
42000 C
42100 C
42200 C
42300 C
42400 C
42500 C
42600 52 NVI=NVI-LNV
42700 NVO=NVO-LNV
42800 NGO=LNV
42900 GO TO 44
43000 C
43100 C NVO ZERO, PRINT LEFT SIDE ONLY
43200 C
43300 46 IF(NVI.LE.0)GO TO 55
43400 C
43500 53 DO 19000 I=1,NVI
43600 II=I+NGO
43700 WRITE(6,8016)ANAME(II),INEN(II),B(II),D(II),F(II),KQ(II)
43800 8016 FORMAT(5X,A6,1X,I2,1X,F16.5,1X,F11.5,1X,F11.4,' (',I1,') .')
43900 19000 CONTINUE
44000 GO TO 55
44100 C
44200 C
44300 C
44400 C
44500 C
44600 C
44700 C
44800 C NVI ZERO,PRINT RIGHT SIDE ONLY
44900 C
45000 C
45100 56 DO 21000 I=1,NVO
45200 II= I+NGO
45300 WRITE(6,8018)BNAME(II),NIEN(II),R(II),TOLEV(II),FE(II),JQ(II)
45400 8018 FORMAT(60X,'. ',A6,1X,I2,1X,F15.5,1X,F13.4,1X,F12.4,' (',I1,')'
45500 1)
45600 21000 CONTINUE
45700 C
45800 C
45900 C
46000 C
46100 C
46200 C
46300 55 IF(L-MAXSTP) 39 ,552,552
46400 552 WRITE (6,8036)
46500 8036 FORMAT(23H SPECIFIED STEP REACHED )
46600 117 ENDFILE 1
46700 REWIND 1
46800 IF(SUMTAB .NE. YES) GO TO 9605
46900 9606 IF(LKL) 9621,9621,9622
47000 9621 WRITE (6,9632)
47100 9632 FORMAT(////49H0SUMMARY TABLE OMITTED DUE TO LACK OF INFORMATION )
47200 GO TO 9605
47300 9622 WRITE (6,9602)
47400 9602 FORMAT(15H1 SUMMARY TABLE// 5X,4HSTEP,16X,8HVARIABLE,15X,8HMULTIPL
47500 1E,18X,8HINCREASE,10X,10HF VALUE TO,5X,21HNUMBER OF INDEPENDENT/4X,
47600 26HNUMBER,10X,7HENTERED,2X,7HREMOVED,9X,1HR,11X,3HRSQ,15X,6HIN RSQ,
47700 39X,15HENTER OR REMOVE,4X,18HVARIABLES INCLUDED//)
47800 R1SQ=0.0
47900 9609 DO 23310 I=1,LKL
48000 READ(1)LMN,KAY,FLAG,XMULTR,FKAY,NVI
48100 RSQ=XMULTR**2
48200 RSQI=RSQ-R1SQ
48300 R1SQ=RSQ
48400 C
48500 C
48600 C
48700 C
48800 C
48900 C
49000 C
49100 C
49200 9611 IF(FLAG) 23314,23310,23313
49300 23313 WRITE (6,9631)LMN,ALBEL(KAY),KAY,XMULTR,RSQ,RSQI,FKAY, NVI
49400 9631 FORMAT(5X,I3,10X,A6,1X,I2,14X,F9.4,F13.4,6X,F11.4,1X,F19.4,15X,I2)
49500 GO TO 23310
49600 23314 WRITE (6,9612)LMN,ALBEL(KAY),KAY,XMULTR,RSQ,RSQI,FKAY, NVI
49700 9612 FORMAT(5X,I3,20X,A6,1X,I2,4X,F9.4,F13.4,6X,F11.4,1X,F19.4,15X,I2)
49800 23310 CONTINUE
49900 9605 REWIND 1
50000 IF(NVIP)8888,8888,75
50100 8888 IF(RESID.EQ.YES)GO TO 58
50200 GO TO 13000
50300 C
50400 C READ INXPLT CARD
50500 C
50600 75 IF(NVIP .GT. 30) GO TO 9012
50700 READ (5,8032)WMAN,(IVPT(J), J=1,NVIP)
50800 8032 FORMAT(A6,33I2)
50900 IF(WMAN .NE. DXPLTS) GO TO 9011
51000 C
51100 755 DO 31000 K=1,NVIP
51200 IF(IVPT(K)-IP) 31000,31000,9011
51300 31000 CONTINUE
51400 4405 DO 20100 J=1,NVIP
51500 XMIN(J)=+999999.9
51600 XMAX(J)=-999999.9
51700 20100 CONTINUE
51800 RMIN=999999.9
51900 RMAX=-999999.9
52000 58 CALL RESIDS
52100 13000 CONTINUE
52200 GO TO 1
52300 9001 WRITE (6,8020)
52400 8020 FORMAT(///24H FINISH CARD ENCOUNTERED)
52500 9999 WRITE (6,8021)
52600 8021 FORMAT(19H PROGRAM TERMINATED)
52700 2210 STOP
52800 9002 WRITE (6,8022)
52900 8022 FORMAT(43H NEITHER PROBLM NOR FINISH CARD ENCOUNTERED)
53000 GO TO 9999
53100 9003 WRITE (6,8023)
53200 8023 FORMAT(50H0NUMBER OF VARIABLES, P OR P+Q, OUTSIDE OF LIMITS.)
53300 GO TO 9999
53400 9004 WRITE (6,8024)
53500 8024 FORMAT(37H CARD INCORRECTLY PUNCHED OR MISSING.)
53600 GO TO 9999
53700 9005 WRITE (6,8025)
53800 8025 FORMAT(31H NO SUB-PROBLEM CARD SPECIFIED.)
53900 GO TO 9999
54000 9006 WRITE (6,8026)
54100 8026 FORMAT(16H0TRANSGENERATION)
54200 GO TO 9004
54300 9009 WRITE (6,8029)
54400 8029 FORMAT(12H0SUB-PROBLEM)
54500 GO TO 9004
54600 9010 WRITE (6,8030)
54700 8030 FORMAT(15H0CONTROL-DELETE)
54800 GO TO 9004
54900 9011 WRITE (6,8033)
55000 8033 FORMAT(11H0INDEX-PLOT)
55100 GO TO 9004
55200 9012 WRITE(6,8035) NVIP
55300 8035 FORMAT('0HTE NUMBER OF VARIABLES SPECIFIED FOR THE INDEX-PLOT CARD
55400 1NUST NOT EXCEED 30,',I11,' IS TOO LARGE.')
55500 GO TO 9999
55600 END
55700 FUNCTION AF(I,J)
55800 C FUNCTION AF FOR BMD02R MAY 2, 1966
55900 DIMENSION KP(50,50),IQ(50,50),STDEV(80),C(80)
56000 COMMON NTGC,TRANS(100),KTRANS(3,100),X(81),A(50,50),IP,FINC,RESDF
56100 1,FOUT,KAY,FLAG,N,NINCS,KDEP,IR,P(50,50),Q(50,50),XN,TOL,DF,LLL(10)
56200 2,B(80),COL(11),PINT(11),QINT(11),RES(275),NVIP,ALPHA,RMAX,RMIN,
56300 3RESID,IVPT(33),NVI,KEEP(5),IS
56400 C
56500 EQUIVALENCE (P,KP),(Q,IQ),(KTRANS(81),STDEV),(TRANS,C)
56600 KKG=MIN0(I,J)
56700 LLG=MAX0(I,J)
56800 AF=A(KKG,LLG)
56900 RETURN
57000 END
57100 SUBROUTINE AOUT
57200 C SUBROUTINE AOUT FOR BMD02R MAY 2, 1966
57300 COMMON NTGC,TRANS(100),KTRANS(3,100),X(81),A(50,50),IP,FINC,RESDF
57400 1,FOUT,KAY,FLAG,N,NINCS,KDEP,IR,P(50,50),Q(50,50),XN,TOL,DF,LLL(10)
57500 2,B(80),COL(11),PINT(11),QINT(11),RES(275),NVIP,ALPHA,RMAX,RMIN,
57600 3RESID,IVPT(33),NVI,KEEP(5),IS
57700 C
57800 DIMENSION FMT1(4),FMT2(4),AI(9)
57900 DIMENSION KP(50,50),IQ(50,50),STDEV(80),C(80)
58000 EQUIVALENCE (P,KP),(Q,IQ),(KTRANS(81),STDEV),(TRANS,C)
58100 DOUBLE PRECISION FMT1,FMT2,AI,FIRST,TOP,SKIP
58200 C THE FOLLOWING DATA STATEMENTS ARE FORMATS FOR THE CORRELATION
58300 C COEFFICIENTS MATRIX.
58400 C
58500 DATA FMT1/' ','H VARIABLE',',(I7,9I11)',') '/,TOP
58600 1,SKIP/' (1H1,9',' (1H0,9'/
58700 DATA FMT2/'(5X,I2, ',' ','F11.3)) ',' '/,
58800 1FIRST/' 2X,(10'/
58900 DATA AI/' 13X,(9',' 24X,(8',' 35X,(7',' 46X,(6',' 57X,(5',
59000 1' 68X,(4',' 79X,(3',' 90X,(2',' 101X,('/
59100 C
59200 MN=0
59300 KK=0
59400 DO 1000 I=1,IP,10
59500 KK=KK+1
59600 NRTEN=I
59700 M=KK*10-IP
59800 IF(M)1,1,2
59900 1 M=9
60000 MN=(M+1)*KK
60100 GO TO 3
60200 2 M=9-M
60300 MN=MN+M+1
60400 3 MM=M+1
60500 DO 2000 LL=1,MM
60600 LLL(LL)=LL+(KK-1)*10
60700 2000 CONTINUE
60800 IF(KK-1) 4,4,5
60900 4 KTOP=0
61000 FMT1(1)=SKIP
61100 GO TO 6
61200 5 KTOP=1
61300 FMT1(1)=TOP
61400 6 WRITE (6,FMT1)(LLL(LK),LK=1,MM)
61500 WRITE(6,9901)
61600 9901 FORMAT(3X,6HNUMBER/)
61700 FMT2(2)=FIRST
61800 K=NRTEN
61900 DO 3000 J=1,NRTEN
62000 WRITE (6,FMT2)J,(A(J,L),L=K,MN)
62100 3000 CONTINUE
62200 NN=K
62300 ID=NRTEN
62400 IF(M)1000,1000,5000
62500 5000 DO 4000 JK=1,M
62600 FMT2(2)=AI(JK)
62700 NN=NN+1
62800 ID=ID+1
62900 WRITE (6,FMT2)ID,(A(ID,L),L=NN,MN)
63000 4000 CONTINUE
63100 1000 CONTINUE
63200 RETURN
63300 END
63400 SUBROUTINE WHICHX(INDEX,NVI,IS,KEEP)
63500 DIMENSION KEEP(5)
63600 IF(NVI.GT.1)GO TO 10
63700 IS=0
63800 DO 5 I=1,5
63900 5 KEEP(I)=0
64000 10 DO 20 I=1,5
64100 IF(KEEP(I).NE.0)GO TO 20
64200 KEEP(I)=INDEX
64300 IS=IS+1
64400 GO TO 30
64500 20 CONTINUE
64600 30 RETURN
64700 END
64800 SUBROUTINE RDLBL2(NLBVAR,NVAR,ARRAY)
64900 C SUBROUTINE RDLBL2 FOR BMD02R MAY 2, 1966
65000 C SUBROUTINE TO READ IN LABELS CARDS, STORE THEM IN ARRAY,
65100 C AND SUBSTITUTE BLANKS FOR UNLABELED VARIABLES
65200 C NVAR IS TOTAL NUMBER OF VARIABLES
65300 C NLBVAR IS NUMBER OF LABELED VARIABLES EXPECTED
65400 C
65500 DIMENSION ARRAY(1),DUMY(7),IDUM(7)
65600 DOUBLE PRECISION ARRAY,BLANKS,DUMY,TEST,ALABEL
65700 DATA ALABEL/'LABELS'/,BLANKS/' '/
65800 C BLANK VARIABLES
65900 DO 1 I=1,NVAR
66000 1 ARRAY(I)=BLANKS
66100 C IF NO LABELS, RETURN
66200 IF(NLBVAR) 9,9,2
66300 2 N=0
66400 C READ 1 LABELS CARD
66500 20 READ (5,3) TEST,(IDUM(J),DUMY(J),J=1,7)
66600 3 FORMAT(A6,7(I4,A6))
66700 C TEST FOR 'LABELS' IN FIRST 6 COLS.
66800 IF(TEST.EQ.ALABEL)GO TO 6
66900 C ERROR--PRINT MESSAGE AND QUIT
67000 4 WRITE (6,5)
67100 5 FORMAT(36H0LABELS CARD NOT FOUND WHEN EXPECTED)
67200 CALL EXIT
67300 STOP
67400 C EXAMINE 7 FIELDS
67500 6 DO 8 J=1,7
67600 K=IDUM(J)
67700 C TEST INDEX. IF 0, IGNORE. IF ILLEGAL, PRINT MESSAGE AND
67800 C IGNORE EXCEPT TO COUNT
67900 IF(K) 11,8,10
68000 10 IF(K-NVAR) 7,7,11
68100 11 WRITE (6,12)K,DUMY(J)
68200 12 FORMAT(18H0LABELS CARD INDEX,I7,18H INCORRECT. LABEL ,A6,9H IGNORE
68300 1D.)
68400 GO TO 13
68500 C MOVE LABEL TO ARRAY
68600 7 ARRAY(K)=DUMY(J)
68700 C STEP NUMBER OF VARIABLES
68800 13 N=N+1
68900 C TEST FOR END. IF END, RETURN. IF NOT, SCAN OTHER FIELDS.
69000 IF(N-NLBVAR) 8,9,9
69100 8 CONTINUE
69200 GO TO 20
69300 9 RETURN
69400 END
69500 SUBROUTINE RESIDS
69600 C SUBROUTINE RESIDS FOR BMD02R MAY 2, 1966
69700 C THIS SUBROUTINE COMPUTES THE RESIDUALS FOR THE REGRESSION.
69800 COMMON NTGC,TRANS(100),KTRANS(3,100),X(81),A(50,50),IP,FINC,RESDF
69900 1,FOUT,KAY,FLAG,N,NINCS,KDEP,IR,P(50,50),Q(50,50),XN,TOL,DF,LLL(10)
70000 2,B(80),COL(11),PINT(11),QINT(11),PES(275),NVIP,ALPHA,RMAX,RMIN,
70100 3RESID,IVPT(33),NVI,KEEP(5),IS
70200 C
70300 DIMENSION FMT1(6),FMT2(6),FMT3(8),FMT4(4),FMT5(5),FMT6(4)
70400 DIMENSION KP(50,50),IQ(50,50),STDEV(80),C(80),XMIN(80),XMAX(80)
70500 EQUIVALENCE (P,KP),(Q,IQ),(KTRANS(1),XMIN),(KTRANS(81),XMAX),
70600 1(KTRANS(161),STDEV),(TRANS,C)
70700 DOUBLE PRECISION FMT1,FMT2,FMT3,FMT4,FMT5,FMT6,SKPONE,SKPTWO,PREON
70800 1 E,PRETWO,ONEPRE,TWOPRE,SECONE,SECTWO,BEFONE,BEFTWO,Q011HL,Q012HL
70900 DATA DOT,YES/'.','YES'/
71000 C THE FOLLOWING DATA STATEMENTS ARE FORMAT STATEMENTS FOR THE
71100 C PLOT OF THE RESIDUALS.
71200 C
71300 DATA FMT1/' ','(26HPLOT O','F RESIDUAL','S (Y-AXIS)',
71400 1',34X,2H..,','3X )) '/,PREONE,PRETWO/'(2H1 1','(2H1 2'/
71500 DATA FMT2/' ','(14H VS. V','ARIABLE ,I','3,9H (X-AX',
71600 1'IS),34X,2H','.. ,3X)) '/,SKPONE,SKPTWO/'(2X, 1','(2X, 2'/
71700 DATA FMT3/' ','6F10.3,2','H..,1X)/',' ',' F7.3,4',
71800 1'F10.3,5X',',2H..,1X',')) '/,ONEPRE,TWOPRE,SECONE,SECTWO/
71900 2'(1(2X,','(2(2X,','1(10X,','2(10X,'/
72000 DATA FMT4/' ','6X,51A1,','3X,2H..,','1X)) '/
72100 DATA FMT5/' ','1X,F6.2,','2H .,50A','1, 3X2H..',' ,1X)) '/
72200 1,BEFONE,BEFTWO/'(1( ','(2( '/
72300 DATA FMT6/' ','8X,1H.,5','0A1,3X,2','H..,1X))'/
72400 DATA Q011HL/6H /
72500 DATA Q012HL/6H* /
72600 YMIN=999999.9
72700 YMAX=-999999.9
72800 IF(RESID.EQ.YES)WRITE(6,406)
72900 406 FORMAT('1 LIST OF RESIDUALS')
73000 C
73100 DO 24000 I=1,N
73200 READ(2)(X(L),L=1,IP)
73300 67890 FORMAT(20A4)
73400 C
73500 C
73600 SUMB=0.0
73700 NVI=0
73800 DO 25000 J=1,IP
73900 IF(C(J)) 581,581,25000
74000 581 NVI=NVI+1
74100 582 SUMB=SUMB+B(NVI)*X(J)
74200 25000 CONTINUE
74300 YHAT=ALPHA+SUMB
74400 RES=X(KDEP)-YHAT
74500 IF(NVIP.LE.0)GO TO 3030
74600 DO 20200 L=1,NVIP
74700 MM=IVPT(L)
74800 IF(XMIN(L)-X(MM)) 81,81,80
74900 80 XMIN(L)=X(MM)
75000 81 IF(XMAX(L)-X(MM)) 97,20200,20200
75100 97 XMAX(L)=X(MM)
75200 20200 CONTINUE
75300 RMAX=AMAX1(RMAX,RES)
75400 RMIN=AMIN1(RMIN,RES)
75500 YMAX=AMAX1(YMAX,YHAT)
75600 YMIN=AMIN1(YMIN,YHAT)
75700 IF(RESID.NE.YES)GO TO 24000
75800 3030 IF(MOD(I,55).NE.1)GO TO 402
75900 IF(I.NE.1)WRITE(6,410)
76000 410 FORMAT('1')
76100 WRITE(6,407)KDEP,(Q011HL,KEEP(L),L=1,IS)
76200 407 FORMAT(/' CASE',9X,'Y',14X,'Y',11X/' NUMBER X(',I2,')',9X,'C
76300 1OMPUTED',7X,'RESIDUAL',7X,5(A1,'X(',I2,')',9X))
76400 WRITE(6,409)
76500 409 FORMAT(' ')
76600 402 WRITE(6,408)I,X(KDEP),YHAT,RES,(X(KEEP(L)),L=1,IS)
76700 408 FORMAT(1X,I4,8(1X,F14.4))
76800 24000 CONTINUE
76900 REWIND 2
77000 IF(NVIP.LE.0)RETURN
77100 SRS=(RMAX-RMIN)/49.0
77200 DO 20300 I=1,11
77300 AI=1+5*(I-1)
77400 COL(I)=RMIN+SRS*(AI-1.0)
77500 20300 CONTINUE
77600 IIX=0
77700 NVIP=NVIP+1
77800 IVPT(NVIP)=81
77900 XMAX(NVIP)=YMAX
78000 XMIN(NVIP)=YMIN
78100 4409 DO 20400 I=1,NVIP,2
78200 II=IVPT(I)
78300 IF(I+1-NVIP) 98,98,118
78400 118 KKK=1
78500 C THE FOLLOWING FMT- STATEMENTS ARE FOR A SINGLE PLOT ON THE PAGE
78600 FMT1(1)=PREONE
78700 FMT2(1)=SKPONE
78800 FMT3(1)=ONEPRE
78900 FMT3(4)=SECONE
79000 FMT4(1)=ONEPRE
79100 FMT5(1)=BEFONE
79200 FMT6(1)=BEFONE
79300 GO TO 99
79400 98 KKK=2
79500 C THE FOLLOWING FMT- STATEMENTS ARE FOR TWO PLOTS ON THE PAGE
79600 FMT1(1)=PRETWO
79700 FMT2(1)=SKPTWO
79800 FMT3(1)=TWOPRE
79900 FMT3(4)=SECTWO
80000 FMT4(1)=TWOPRE
80100 FMT5(1)=BEFTWO
80200 FMT6(1)=BEFTWO
80300 JJ=IVPT(I+1)
80400 J1=I+1
80500 99 SPI=((XMAX(I )-XMIN(I ))/49.0)
80600 IF(KKK-1) 101,101,102
80700 102 SQI=((XMAX(J1)-XMIN(J1))/49.0)
80800 101 DO 20500 J=1,11
80900 AJ=1+5*(J-1)
81000 PINT(J)=XMIN(I )+SPI*(AJ-1.0)
81100 IF(KKK-1) 20500,20500,104
81200 104 QINT(J)=XMIN(J1)+SQI*(AJ-1.0)
81300 20500 CONTINUE
81400 DO 20600 K=1,50
81500 DO 20600 J=1,50
81600 KP(K,J)=IIX
81700 IQ(K,J)=IIX
81800 20600 CONTINUE
81900 C
82000 C
82100 C
82200 DO 20700 J=1,N
82300 READ(2)(X(L),L=1,IP)
82400 C
82500 C
82600 IR = 1
82700 C
82800 C
82900 NVI=0
83000 RESS=0.0
83100 DO 450JPJ=1,IP
83200 IF(C(JPJ))451,451,450
83300 451 NVI=NVI+1
83400 RESS=RESS+B(NVI)*X(JPJ)
83500 450 CONTINUE
83600 YHAT=RESS+ALPHA
83700 RESS=X(KDEP)-YHAT
83800 IF(I.GE.NVIP-1)X(81)=YHAT
83900 IF(SPI.NE.0.0)IR=(X(II)-XMIN(I))/SPI+1.5
84000 IRES = 1
84100 IF(SRS .NE. 0.0)
84200 1IRES=((RESS-RMIN)/SRS)+1.5
84300 KP(IRES,IR)= KP(IRES,IR)+1
84400 IF(KKK-1) 20700,20700,109
84500 109 JQ = 1
84600 IF(SQI .NE. 0.0)
84700 1JQ=((X(JJ)-XMIN(J1))/SQI)+1.5
84800 IQ(IRES,JQ)=IQ(IRES,JQ)+1
84900 20700 CONTINUE
85000 REWIND 2
85100 WRITE (6,FMT1)
85200 IF(KKK-1) 116,116,110
85300 116 IF(II.LT.81)WRITE(6,FMT2)II
85400 IF(II.EQ.81)WRITE(6,105)
85500 105 FORMAT(' VS. COMPUTED Y (X-AXIS)',34X,'..')
85600 GO TO 111
85700 110 IF(JJ.LT.81)WRITE(6,FMT2)II,JJ
85800 IF(JJ.EQ.81)WRITE(6,100)II
85900 100 FORMAT(' VS. VARIABLE ',I3,' (X-AXIS)',34X,'.. VS. COMPUTED Y
86000 1 (X-AXIS)',34X,'..')
86100 111 XX=(+Q011HL)
86200 XY=(+Q012HL)
86300 WRITE (6,9961)
86400 9961 FORMAT(1H )
86500 IF(KKK-1) 9613,9613,9614
86600 9613 WRITE (6,FMT3)(PINT(K),K=1,11,2),(PINT(L),L=2,11,2)
86700 GO TO 9615
86800 9614 WRITE (6,FMT3)(PINT(K),K=1,11,2),(QINT(L),L=1,11,2), (PINT(M),M=2,
86900 111,2),(QINT(J),J=2,11,2)
87000 9615 IF(KKK-1) 9401,9401,9402
87100 9401 MNMN=51
87200 GO TO 9403
87300 9402 MNMN=102
87400 9403 WRITE (6,FMT4)(DOT,J=1,MNMN)
87500 KLN=0
87600 KLM=5
87700 DO 20800 K=1,50
87800 DO 20900 J=1,50
87900 IF(KP(K,J)) 82,82,83
88000 82 P(K,J)=XX
88100 GO TO 86
88200 83 IF(KP(K,J)-10) 84,85,85
88300 C LEFT ADJUST THE INTEGER IN KP(K,J) WHICH IS LESS THAN 10
88400 84 KP(K,J) = INUMB(KP(K,J))
88500 GO TO 86
88600 85 P(K,J)=XY
88700 86 IF(KKK-1) 20900,20900,112
88800 112 IF(IQ(K,J)) 87,87,88
88900 87 Q(K,J)=XX
89000 GO TO 20900
89100 88 IF(IQ(K,J)-10) 89,91,91
89200 C LEFT ADJUST THE INTEGER IN IQ(K,J) WHICH IS LESS THAN 10
89300 89 IQ(K,J) = INUMB(IQ(K,J))
89400 GO TO 20900
89500 91 Q(K,J)=XY
89600 20900 CONTINUE
89700 KLM=KLM-1
89800 IF(KLM-4) 93,94,94
89900 93 IF(KLM) 95,95,96
90000 94 KLN=KLN+1
90100 IF(KKK-1) 9551,9551,9552
90200 9551 WRITE (6,FMT5)COL(KLN),(P(K,J),J=1,50)
90300 GO TO 20800
90400 9552 WRITE (6,FMT5)COL(KLN),(P(K,J),J=1,50),COL(KLN),(Q(K,L),L=1,50)
90500 GO TO 20800
90600 95 KLM=5
90700 96 IF(KKK-1) 9661,9661,9662
90800 9661 WRITE (6,FMT6)(P(K,J),J=1,50)
90900 GO TO 20800
91000 9662 WRITE (6,FMT6)(P(K,J),J=1,50),(Q(K,L),L=1,50)
91100 20800 CONTINUE
91200 WRITE (6,FMT4)(DOT,J=1,MNMN)
91300 IF(KKK-1) 113,113,114
91400 113 WRITE (6,FMT3)(PINT(K),K=1,11,2),(PINT(L),L=2,11,2)
91500 GO TO 20401
91600 114 WRITE (6,FMT3)(PINT(K),K=1,11,2),(QINT(L),L=1,11,2),(PINT(M),M=2,1
91700 11,2),(QINT(J),J=2,11,2)
91800 20401 REWIND 2
91900 20400 CONTINUE
92000 RETURN
92100 END
92200 SUBROUTINE STEPRG
92300 C SUBROUTINE STEPRG FOR BMD02R MAY 2, 1966
92400 COMMON NTGC,TRANS(100),KTRANS(3,100),X(81),A(50,50),IP,FINC,RESDF
92500 1,FOUT,KAY,FLAG,N,NINCS,KDEP,IR,P(50,50),Q(50,50),XN,TOL,DF,LLL(10)
92600 2,B(80),COL(11),PINT(11),QINT(11),RES(275),NVIP,ALPHA,RMAX,RMIN,
92700 3RESID,IVPT(33),NVI,KEEP(5),IS
92800 C
92900 DIMENSION KP(50,50),IQ(50,50),STDEV(80),C(80)
93000 EQUIVALENCE (P,KP),(Q,IQ),(KTRANS(81),STDEV),(TRANS,C)
93100 CC=A(KDEP,KDEP)
93200 AA=XN-DF
93300 BB=FINC+AA-1.0
93400 VIN=0.0
93500 IF(BB.NE.0.0)VIN=FINC*CC/BB+2.0
93600 VOUT=0.0
93700 IF(AA.NE.0.0)VOUT=FOUT*CC/AA-7.0
93800 VMIN= 9999.9
93900 VMAX=0.0
94000 KMIN=0
94100 KMAX=0
94200 DO 1000 K=1,IP
94300 IF(C(K)-1.0) 9,1000,10
94400 9 VSUBK = C(K)
94500 IF(A(K,K) .NE. 0.0)
94600 1VSUBK = C(K) - (AF(K,KDEP)**2) / A(K,K)
94700 IF(VMIN-VSUBK) 1000,1000,1
94800 1 VMIN=VSUBK
94900 KMIN=K
95000 GO TO 1000
95100 10 IF(A(K,K)-TOL) 1000,8,8
95200 8 VSUBK = C(K)
95300 IF(A(K,K) .NE. 0.0)
95400 1VSUBK=C(K)+AF(K,KDEP)**2/A(K,K)
95500 IF(VSUBK-VMAX) 1000,1000,4
95600 4 IF(XN - DF - 3.0 + C(K)) 1000,1000,45
95700 45 VMAX=VSUBK
95800 KMAX=K
95900 1000 CONTINUE
96000 IF(VOUT-VMIN) 2,2,3
96100 3 C(KMIN)=C(KMIN)+9.0
96200 KAY=KMIN
96300 FLAG= -1.0
96400 GO TO 7
96500 2 IF(CC)6,6,25
96600 25 IF(VMAX - VIN) 6,5,5
96700 5 IF(KMAX)6,6,11
96800 11 C(KMAX)=C(KMAX)-9.0
96900 KAY=KMAX
97000 FLAG=1.0
97100 7 CALL STEP
97200 DF=DF+FLAG
97300 RETURN
97400 6 FLAG= 0.0
97500 RETURN
97600 END
97700 SUBROUTINE STEP
97800 C SUBROUTINE STEP FOR BMD02R MAY 2, 1966
97900 DIMENSION KP(50,50),IQ(50,50),STDEV(80),C(80)
98000 DIMENSION U(81)
98100 COMMON NTGC,TRANS(100),KTRANS(3,100),X(81),A(50,50),IP,FINC,RESDF
98200 1,FOUT,KAY,FLAG,N,NINCS,KDEP,IR,P(50,50),Q(50,50),XN,TOL,DF,LLL(10)
98300 2,B(80),COL(11),PINT(11),QINT(11),RES(275),NVIP,ALPHA,RMAX,RMIN,
98400 3RESID,IVPT(33),NVI,KEEP(5),IS
98500 C
98600 EQUIVALENCE (P,KP),(Q,IQ),(KTRANS(81),STDEV),(TRANS,C)
98700 KAY1=KAY-1
98800 KAY2=KAY+1
98900 XAY3=A(KAY,KAY)
99000 IF(KAY1) 3,3,4
99100 4 DO 1000 I=1,KAY1
99200 U(I)=A(I,KAY)
99300 A(I,KAY)=0.0
99400 1000 CONTINUE
99500 3 U(KAY)=-FLAG
99600 A(KAY,KAY)=0.0
99700 IF(KAY2-IP ) 1,1,2
99800 1 DO 2000 I=KAY2,IP
99900 U(I)=A(KAY,I)
00100 A(KAY,I)=0.0
00200 2000 CONTINUE
00300 2 DO 3000 I=1,IP
00400 DO 3000 J=I,IP
00500 IF(XAY3 .NE. 0.0)
00600 1A(I,J)=A(I,J)-(U(I)*U(J))/XAY3
00700 3000 CONTINUE
00800 RETURN
00900 END
01000 SUBROUTINE TRANGN
01100 C SUBROUTINE TRANGN FOR BMD02R MAY 2, 1966
01200 COMMON NTGC,TRANS(100),KTRANS(3,100),X(81),A(50,50),IP,FINC,RESDF
01300 1,FOUT,KAY,FLAG,L,NINCS,KDEP,IR,P(50,50),Q(50,50),XN,TOL,DF,LLL(10)
01400 2,B(80),COL(11),PINT(11),QINT(11),RES(275),NVIP,ALPHA,RMAX,RMIN,
01500 3RESID,IVPT(33),NVI,KEEP(5),IS
01600 C
01700 DIMENSION KP(50,50),IQ(50,50),STDEV(80),C(80)
01800 EQUIVALENCE (P,KP),(Q,IQ),(KTRANS(81),STDEV),(TRANS,C)
01900 ASN(XX)=ATAN(XX/SQRT(1.0-XX**2))
02000 DO 100 I=1,NTGC
02100 M=KTRANS(1,I)
02200 N=KTRANS(3,I)
02300 NTRANS=KTRANS(2,I)
02400 IF(M-81) 91,91,99
02500 91 IF(N-81) 92,92,99
02600 92 IF((NTRANS-25)*NTRANS)50,99,99
02700 99 WRITE (6,199)I
02800 199 FORMAT(22H TRANSGENERATION CARD ,I3,27H MISPUNCHED OR OUT OF ORDER
02900 1)
03000 GO TO 100
03100 50 GO TO (1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,99,99,20,21,22,23
03200 1,24),NTRANS
03300 1 IF(X(N)) 198,107,108
03400 107 X(M)=0.0
03500 GO TO 100
03600 108 X(M)= SQRT(X(N))
03700 GO TO 100
03800 2 IF(X(N)) 198,111,112
03900 111 X(M)=1.0
04000 GO TO 100
04100 112 X(M)= SQRT(X(N))+SQRT(X(N)+1.0)
04200 GO TO 100
04300 3 IF(X(N)) 198,198,114
04400 114 X(M) = ALOG10(X(N))
04500 GO TO 100
04600 4 X(M)=EXP(X(N))
04700 GO TO 100
04800 5 IF(X(N)) 198,107,117
04900 117 IF(X(N)-1.0) 118,119,119
05000 118 E=SQRT(X(N))
05100 X(M)=ASN(E)
05200 GO TO 100
05300 119 X(M)=3.14159265/2.0
05400 GO TO 100
05500 6 FN=L
05600 E = 0.0
05700 B(1) = 0.0
05800 IF((FN+1.0) .EQ. 0.0) GO TO 61
05900 E=X(N)/(FN+1.0)
06000 B(1) = E + 1.0 / (FN + 1.0)
06100 61 IF(E) 198,123,124
06200 123 IF(B(1)) 198,107,127
06300 127 X(M)=ASN(SQRT(B(1)))
06400 GO TO 100
06500 124 IF(B(1)) 198,128,129
06600 128 X(M)=ASN(SQRT(E))
06700 GO TO 100
06800 129 E=SQRT(E)
06900 B(1) = SQRT(B(1))
07000 X(M)=ASN(E)+ASN(B(1))
07100 GO TO 100
07200 7 IF(X(N)) 131,198,131
07300 131 X(M)= 1.0/X(N)
07400 GO TO 100
07500 8 X(M)=X(N)+ TRANS( I )
07600 GO TO 100
07700 9 X(M)=X(N)* TRANS( I )
07800 GO TO 100
07900 10 IF(X(N)) 198,107,133
08000 133 X(M)= X(N)**TRANS( I)
08100 GO TO 100
08200 11 NEWB= TRANS( I)
08300 X(M)= X(N)+X(NEWB)
08400 GO TO 100
08500 12 NEWB= TRANS( I)
08600 X(M)=X(N)-X(NEWB)
08700 GO TO 100
08800 13 NEWB= TRANS( I)
08900 X(M)=X(N)*X(NEWB)
09000 GO TO 100
09100 14 NEWB= TRANS( I)
09200 IF(X(NEWB)) 134,197,134
09300 134 X(M)= X(N)/X(NEWB)
09400 GO TO 100
09500 15 IF(X(N)-TRANS( I)) 107,111,111
09600 16 XNEWB=TRANS( I)
09700 IF(X(N)-(XNEWB)) 107,111,111
09800 17 IF(X(N))198,198,163
09900 163 X(M)= ALOG(X(N))
10000 GO TO 100
10100 20 X(M)= SIN(X(N))
10200 GO TO 100
10300 21 X(M)= COS(X(N))
10400 GO TO 100
10500 22 IF(X(N)-1.57079632) 186,186,198
10600 186 IF(X(N)+1.57079632) 198,187,187
10700 187 X(M)=ATAN(X(N))
10800 GO TO 100
10900 23 NEWB= TRANS( I)
11000 IF(X(N)) 198,198,188
11100 188 X(M)=X(N)**X(NEWB)
11200 GO TO 100
11300 24 IF(TRANS(I)) 198,198,189
11400 189 X(M)= TRANS(I)**X(N)
11500 GO TO 100
11600 197 N=NEWB
11700 198 WRITE (6,201)N,NINCS,KTRANS(2,I),M
11800 201 FORMAT(23H THE VALUE OF VARIABLE ,I4, 9H IN CASE ,I5,55H VIOLATED
11900 1THE RESTRICTIONS FOR TRANSGENERATION OF TYPE ,I3,1H./40H THE PROGR
12000 2AM CONTINUED LEAVING VARIABLE ,I4,11H UNCHANGED.)
12100 100 CONTINUE
12200 RETURN
12300 END
12400 FUNCTION INUMB(I)
12500 C
12600 C THE FUNCTION 'INUMB' PLACES A LEFT JUSTIFIED ALPHANUMERIC CHARACTE
12700 C REPRESENTING THE HEXADECIMAL FORM OF 'I' INTO 'INUMB' IF 'I' IS BE
12800 C '0' AND '15' INCLUSIVE. A '*' IS RETURNED IF I FALLS OUTSIDE THES
12900 C
13000 DIMENSION IT(17)
13100 DATA IT/' ','1','2','3','4','5','6','7','8','9','A','B','C','D','E
13200 .','F','*'/
13300 IF (I.LT.0.OR.I.GT.15) GO TO 20
13400 INUMB=IT(I+1)
13500 RETURN
13600 20 INUMB=IT(17)
13700 RETURN
13800 END
13900 SUBROUTINE ATOF(A,N,F)
14000 DIMENSION A(1)
14100 LOGICAL BLANK
14200 BLANK=.TRUE.
14300 S=1.0
14400 NUMB=0
14500 TEN=1.0
14600 DIV=1.0
14700 DO 10 I=1,N
14800 L=INTCHR(A,I)
14900 IF(L.EQ.36) GO TO 10
15000 BLANK=.FALSE.
15100 IF(L.NE.38) GO TO 2
15200 S=-1.0
15300 GO TO 10
15400 2 IF(L.NE.44) GO TO 4
15500 TEN=10.0
15600 GO TO 10
15700 4 IF(L.GT.9) GO TO 9
15800 NUMB=NUMB*10+L
15900 DIV=DIV*TEN
16000 9 CONTINUE
16100 10 CONTINUE
16200 IF(BLANK)RETURN
16300 F=S*FLOAT(NUMB)/DIV
16400 RETURN
16500 END
16600 FUNCTION INTCHR(STRING,N)
16700 DIMENSION SEQ(50),STRING(1),EBCD(5)
16800 DATA SEQ/1H0,1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9,
16900 X 1HA,1HB,1HC,1HD,1HE,1HF,1HG,1HH,1HI,1HJ,
17000 X 1HK,1HL,1HM,1HN,1HO,1HP,1HQ,1HR,1HS,1HT,
17100 X 1HU,1HV,1HW,1HX,1HY,1HZ,1H ,1H+,1H-,1H*,
17200 X 1H/,1H(,1H),1H,,1H.,1H',1H=,1H$,1H ,1H /
17300 DATA EBCD/1H+,1H(,1H),1H',1H=/
17400 CALL GETCHR(STRING,N,CHR)
17500 IF (CHR.NE.SEQ(37)) GO TO 2
17600 INTCHR = 36
17700 GO TO 10
17800 2 DO 1 I=1,48
17900 IF(SEQ(I).EQ.CHR) GO TO 9
18000 1 CONTINUE
18100 I=51
18200 IF(EBCD(1).EQ.CHR) I=38
18300 IF(EBCD(2).EQ.CHR) I=42
18400 IF(EBCD(3).EQ.CHR) I=43
18500 IF(EBCD(4).EQ.CHR) I=46
18600 IF(EBCD(5).EQ.CHR) I=47
18700 9 INTCHR=I-1
18800 10 RETURN
18900 END
19000
19100
19200
19300
19400
19500
19600
19700
19800
19900
20000
20100
20200