Google
 

Trailing-Edge - PDP-10 Archives - decus_20tap5_198111 - 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