Trailing-Edge
-
PDP-10 Archives
-
decuslib20-05
-
decus/20-0137/manova/manov4.for
There is 1 other file named manov4.for in the archive. Click here to see a list.
00100 C WMU IMPLEMENTATION OF THE WAYNE STATE VERSION OF
00200 C
00300 C THE MIAMI MULTIVARIATE ANALYSIS OF VARIANCE(MANOVA)
00400 C
00500 C MODIFIED FOR WMU COMPUTER CENTER BY: SAM ANEMA
00600 C
00700 C FILES: MANOVA.F4,MANOV1.F4,MANOV2.F4,MANOV3.F4,MANOV4.F4
00800 C
00900 C ADDITIONAL SUBROUTINES: USAGE - RUSS BARR(IN NGLIB).MAC
01000 C
01100 C LOADING PROCEDURE:
01200 C
01300 C R LOADER
01400 C MANOV=MANOVA,USAGE/<MANOV1/>MANOV2/>MANOV3/>MANOV4/>/G
01500 C SAV MANOVA
01600 C
01700 C MANOV4.F4
01800 SUBROUTINE ANALYS
01900 DIMENSION VECHLD(100),INLET(75),SSRES(26,27)
02000 COMMON ORTHES(100,20),DUMMY(26,27),SSHYP(25,26),SSEAD(40,41),
02100 1ESTIM(50,50), NUMERR(8),NERRS(100),NDFCUM(100),NDFTST(
02200 2100),VARNAM(2,50),HEAD(3,100),LEVEL(8),LEVSUB(8,10),LEVCUM(8,10),
02300 3NCELL(100), NTABLE(27),ITABLE(9,9),OBS(100)
02400 4 ,NVAR,NCOVAR,NERWIT,NERRES,HNUM(100),ERROR,NTESTS,RVARC,
02500 5 FIRST,IORD(50),IPOSV(50),IVPOS(50),NCELLS,NVART,NDFTOT,
02600 6SPECOR,VLIST,PRINTR,NFACT , READK, PRINTK,CONTR ,TESTR, MFIRST,
02700 7TRUTH,BLANK,MAXFAC,MAXCEL,MAXPAR,MAXLEV,MAXVAR,ATITLE,AJOBCD,
02800 8AANALY,AFINIS,WITHIN,SPACE(10)
02900 DATA INP,IOUT,IAUX1,IAUX2/5,30,22,23/
03000 DATA IAUX4/24/
03100 EQUIVALENCE (INLET,SSEAD), (VECHLD,SSHYP),(SSRES,DUMMY),
03200 1(TAG,SPACE(1) ), (EONLY,SPACE(2) ),(NSIG,SPACE(7)),(TAGP,SPACE(3))
03300 2,(PRTCLL,SPACE(4)),(PRTSSH,SPACE(5)),(NPCT,SPACE(6))
03400 LOGICAL ERROR,TRUTH,VLIST,SPECOR,READK,CONTR,TESTR,MFIRST,WITHIN
03500 2,EONLY,PRINTK,PRINTR,PRTCLL,PRTSSH,NLIST
03600 INTEGER BLANK
03700 DOUBLE PRECISION VARNAM,VECHLD
03800 REWIND IAUX2
03900 NVARTP = NVART + 1
04000 IF (VLIST) GO TO 120
04100 GO TO 130
04200 100 IF (WITHIN) GO TO 110
04300 REWIND IAUX2
04400 READ(IAUX2)((SSRES(J,K),J = 1,NVART),K = 1,NVARTP)
04500 110 IF (.NOT. VLIST) GO TO 140
04600 120 CALL REORD
04700 IF (ERROR) GO TO 210
04800 130 IF (NVAR .EQ. 1) GO TO 140
04900 REWIND IAUX2
05000 WRITE(IAUX2)((SSRES(J,K),J = 1,NVART),K = 1,NVARTP)
05100 140 REWIND IAUX2
05200 WITHIN = .TRUE.
05300 CALL SETTST
05400 IF (ERROR) GO TO 210
05500 IF (J .EQ. 1) WRITE (IOUT,150)
05600 150 FORMAT (52H0THERE HAS BEEN A DIVISION BY ZERO IN THIS ANALYSIS )
05700 NSIG = NSIG-1
05800 IF (NSIG .GT. 0) GO TO 210
05900 WRITE (IOUT,157)
06000 157 FORMAT(1H1/' WESTERN MICHIGAN UNIVERSITY'//,
06100 1 ' MULTI-VARIATE ANALYSIS OF VARIANCE'/)
06200 155 READ(INP,160) TAG,TAGP,INLET
06300 160 FORMAT (A4,76A1)
06400 IF (TAG .NE. ATITLE) GO TO 165
06500 WRITE(IOUT,163) INLET
06600 163 FORMAT (1H 80A1)
06700 GO TO 155
06800 165 TRUTH = TAG .EQ. AANALY
06900 IF (.NOT. TRUTH) GO TO 210
07000 NLIST = INLET(1) .NE. BLANK
07100 SPECOR = INLET(2) .NE. BLANK
07200 CONTR = INLET(3) .NE. BLANK
07300 DO 166 J = 2,8
07400 IF (INLET(J) .NE. BLANK) GO TO 167
07500 166 CONTINUE
07600 GO TO 169
07700 167 TESTR = INLET(4) .NE. BLANK
07800 READK = INLET(5) .NE. BLANK
07900 MFIRST = INLET(6) .EQ. BLANK
08000 EONLY = INLET(7) .NE. BLANK
08100 PRTCLL = INLET(8) .NE. BLANK
08200 EONLY = EONLY .OR. PRTCLL
08300 SPECOR = SPECOR .OR. EONLY
08400 MFIRST = MFIRST .AND. .NOT. PRTCLL
08500 PRINTK = PRINTK .AND. .NOT. EONLY
08600 PRINTR = PRINTR .AND. .NOT. EONLY
08700 PRTSSH = PRTSSH .AND. .NOT. EONLY
08800 169 NSIG = 1
08900 DO 170 I = 1,9
09000 170 IF (INLET(7) .EQ. NTABLE(I+9) .OR. INLET(8) .EQ. NTABLE(I+9))
09100 1NSIG = I
09200 IF (NLIST)READ(INP,180) NVAR,NCOVAR,(IORD(J),J = 1,38)
09300 VLIST = VLIST .OR. NLIST
09400 180 FORMAT (40I2)
09500 NVCVAR = NVAR+NCOVAR
09600 IF(NLIST.AND.NVCVAR.GT. 38) READ(INP,180)(IORD(J), J = 39,NVCVAR)
09700 IF (NCOVAR .EQ. 0) NCOVAR = 0
09800 DO 190 J = 1,NVCVAR
09900 J1 = IORD(J)
10000 K = IPOSV(J1)
10100 VECHLD (2*J-1) = VARNAM(1,K)
10200 IF (NVCVAR .GT. NVART) GO TO 220
10300 190 VECHLD(2*J) = VARNAM(2,K)
10400 N2 = 2*NVCVAR
10500 NPCT = NPCT+1
10600 NP = NPCT/1000
10700 NR = NPCT - NP*1000
10800 WRITE(IOUT,200)NP,NR,NVAR,NCOVAR, (VECHLD(J),J = 1,N2)
10900 200 FORMAT (8H0PROBLEM I4, 11H REANALYSIS I4, 19H WITH THE FOLLOWING
11000 1I4,13H CRITERIA AND I4, 11H COVARIATES /1X10(3X A6,A4) )
11100 IF (.NOT.(SPECOR .OR. CONTR) ) GO TO 100
11200 210 RETURN
11300 220 WRITE(IOUT,230)
11400 230 FORMAT (35H TOO MANY VARIATES AND COVARIATES )
11500 ERROR = .TRUE.
11600 GO TO 210
11700 END
11800 SUBROUTINE REORD
11900 DIMENSION SSWITH(26,26),SSRES(26,27),Q(500),GM(40)
12000 COMMON ORTHES(100,20),DUMMY(26,27),SSHYP(25,26),SSEAD(40,41),
12100 1ESTIM(50,50), NUMERR(8),NERRS(100),NDFCUM(100),NDFTST(
12200 2100),VARNAM(2,50),HEAD(3,100),LEVEL(8),LEVSUB(8,10),LEVCUM(8,10),
12300 3NCELL(100), NTABLE(27),ITABLE(9,9),OBS(100)
12400 4 ,NVAR,NCOVAR,NERWIT,NERRES,HNUM(100),ERROR,NTESTS,RVARC,
12500 5 FIRST,IORD(50),IPOSV(50),IVPOS(50),NCELLS,NVART,NDFTOT,
12600 6SPECOR,VLIST,PRINTR,NFACT , READK, PRINTK,CONTR ,TESTR, MFIRST,
12700 7TRUTH,BLANK,MAXFAC,MAXCEL,MAXPAR,MAXLEV,MAXVAR,ATITLE,AJOBCD,
12800 8AANALY,AFINIS,WITHIN,SPACE(10)
12900 EQUIVALENCE (SSWITH,SSRES(1,2) ),(SSRES,DUMMY)
13000 1,(Q,DUMMY(203)), (GM,Q(301))
13100 DATA INP,IOUT,IAUX1,IAUX2/5,30,22,23/
13200 DATA IAUX4/24/
13300 DOUBLE PRECISION VARNAM,XHOLDD
13400 LOGICAL ERROR
13500 NVCVAR = NVAR + NCOVAR
13600 DO 190 I = 1,NVCVAR
13700 NV = IORD(I)
13800 IF (NV .EQ. 0 .OR. NV .GT. NVART) GO TO 210
13900 N1 = IPOSV(NV)
14000 IF (N1 .EQ. I) GO TO 190
14100 K = IVPOS(I)
14200 IVPOS(I) = NV
14300 IVPOS(N1) = K
14400 IPOSV(NV) = I
14500 IPOSV(K) = N1
14600 N2 = I
14700 NMIN = MIN0(N1,N2)
14800 N = NMIN-1
14900 IF (N .EQ. 0) GO TO 110
15000 DO 100 J = 1,N
15100 XHOLD = SSWITH(J,N2)
15200 SSWITH(J,N2) = SSWITH(J,N1)
15300 SSWITH(J,N1) = XHOLD
15400 XHOLD = SSRES(N2,J)
15500 SSRES(N2,J) = SSRES(N1,J)
15600 100 SSRES(N1,J) = XHOLD
15700 110 NMAX = MAX0(N1,N2)
15800 N = NMAX+1
15900 IF (N .GT. NVART) GO TO 130
16000 DO 120 J = N,NVART
16100 XHOLD = SSWITH(N1,J)
16200 SSWITH(N1,J) = SSWITH(N2,J)
16300 SSWITH(N2,J) = XHOLD
16400 XHOLD = SSRES(J,N1)
16500 SSRES(J,N1) = SSRES(J,N2)
16600 120 SSRES(J,N2) = XHOLD
16700 130 IF (NMAX-NMIN .EQ. 1) GO TO 150
16800 M1 = NMIN + 1
16900 M2 = NMAX - 1
17000 DO 140 J = M1,M2
17100 XHOLD = SSWITH(J ,NMAX)
17200 SSWITH(J ,NMAX) = SSWITH(NMIN,J )
17300 SSWITH(NMIN,J ) = XHOLD
17400 XHOLD = SSRES(NMAX,J )
17500 SSRES(NMAX,J ) = SSRES(J ,NMIN)
17600 140 SSRES(J ,NMIN) = XHOLD
17700 150 XHOLD = SSWITH(NMAX,NMAX)
17800 SSWITH(NMAX,NMAX) = SSWITH(NMIN,NMIN)
17900 SSWITH(NMIN,NMIN) = XHOLD
18000 XHOLD = SSRES (NMAX,NMAX)
18100 SSRES (NMAX,NMAX) = SSRES (NMIN,NMIN)
18200 SSRES (NMIN,NMIN) = XHOLD
18300 DO 160 K = 1,2
18400 XHOLDD = VARNAM(K,N2)
18500 VARNAM(K,N2) = VARNAM(K,N1)
18600 160 VARNAM(K,N1) = XHOLDD
18700 DO 170 J = 1,NCELLS
18800 XHOLD = ORTHES(J,N2)
18900 ORTHES(J,N2) = ORTHES(J,N1)
19000 170 ORTHES(J,N1) = XHOLD
19100 DO 180 J = 1,NCELLS
19200 XHOLD = ESTIM(J,N2)
19300 ESTIM(J,N2) = ESTIM(J,N1)
19400 180 ESTIM(J,N1) = XHOLD
19500 XHOLD = GM(N2)
19600 GM(N2) = GM(N1)
19700 GM(N1) = XHOLD
19800 190 CONTINUE
19900 200 RETURN
20000 210 ERROR = .TRUE.
20100 WRITE(IOUT,220) NV
20200 220 FORMAT (9H0VARIABLE I4, 15H DOES NOT EXIST )
20300 GO TO 200
20400 END
20500 SUBROUTINE SETTST
20600 DIMENSION TAGR(3),ERRNAM(11,3),HOLD(18),SD(50),TAG(3),
20700 1SSERR(25,25),ESTHLD(50),SSE(40,40) ,SSWITH(26,26),NERSAV(11),
20800 2SSRES(26,27),Q(500)
20900 DATA INP,IOUT,IAUX1,IAUX2/5,30,22,23/
21000 DATA IAUX4/24/
21100 COMMON ORTHES(100,20),DUMMY(26,27),SSHYP(25,26),SSEAD(40,41),
21200 1ESTIM(50,50), NUMERR(8),NERRS(100),NDFCUM(100),NDFTST(
21300 2100),VARNAM(2,50),HEAD(3,100),LEVEL(8),LEVSUB(8,10),LEVCUM(8,10),
21400 3NCELL(100), NTABLE(27),ITABLE(9,9),OBS(100)
21500 4 ,NVAR,NCOVAR,NERWIT,NERRES,HNUM(100),ERROR,NTESTS,RVARC,
21600 5 FIRST,IORD(50),IPOSV(50),IVPOS(50),NCELLS,NVART,NDFTOT,
21700 6SPECOR,VLIST,PRINTR,NFACT , READK, PRINTK,CONTR ,TESTR, MFIRST,
21800 7TRUTH,BLANK,MAXFAC,MAXCEL,MAXPAR,MAXLEV,MAXVAR,ATITLE,AJOBCD,
21900 8AANALY,AFINIS,WITHIN,SPACE(10)
22000 EQUIVALENCE (SSERR,SSHYP(1,2)), (SSE,SSEAD(1,2) ) ,
22100 1 (SSWITH,SSRES(1,2) ),(SSRES,DUMMY), (EONLY,SPACE(2) )
22200 2,(Q,DUMMY(203)),(Q,SD), (Q(51),ESTHLD),(PRTSSH,SPACE(5))
22300 LOGICAL REGRES,ERROR,RVARC,WITHIN,EONLY,PRTSSH
22400 REAL MS,NTABLE
22500 DOUBLE PRECISION VARNAM
22600 NVARTP = NVART + 1
22700 NERTST = 0
22800 NVCVAR = NVAR + NCOVAR
22900 NVAR1 = NVAR+1
23000 REGRES = .FALSE.
23100 KERR = 0
23200 IF (NVAR .EQ. 1 .AND. .NOT. EONLY) WRITE(IOUT,100)
23300 100 FORMAT (7H1SOURCE 27X 2HSS 7X 2HDF 10X 2HMS 14X 1HF10X 11HP LE
23400 1SS THAN )
23500 NERRH = 0
23600 NHYP= NCOVAR
23700 NERR = NERWIT
23800 CALL ALPHAN(3,TAG,12HWITHIN CELLS)
23900 NESKIP = 1
24000 GO TO 160
24100 110 NERRH = 10
24200 CALL ALPHAN(3,TAG,12HRESIDUAL )
24300 NHYP= NCOVAR
24400 NERR = NERRES
24500 NESKIP = 2
24600 GO TO 160
24700 120 NERRH = 9
24800 CALL ALPHAN(3,TAG,12HWITHIN+RESID )
24900 NHYP = NCOVAR
25000 NERR = NERRES+NERWIT
25100 NESKIP = 3
25200 GO TO 160
25300 130 IF (KERR .EQ. 8) GO TO 770
25400 140 KERR = KERR+1
25500 IF (NUMERR(KERR) .EQ. 0) GO TO 130
25600 NERRH = KERR
25700 NCOD = NUMERR(KERR)
25800 CALL ALPHAN(2,TAG,8HERROR )
25900 TAG(3) = HNUM(KERR)
26000 NHYP= NCOVAR
26100 NERR = NDFTST(NCOD)
26200 NESKIP = 4
26300 150 NSTART = NDFCUM(NCOD)
26400 NFIN = NSTART-1+NDFTST(NCOD)
26500 160 DO 170 J = 1,NTESTS
26600 IF (NERRS(J) .EQ. NERRH) GO TO 180
26700 170 CONTINUE
26800 GO TO 760
26900 180 IF (NVAR .GT. 1) GO TO 190
27000 NERTST = NERTST + 1
27100 ERRNAM(NERTST,1) = TAG(1)
27200 ERRNAM(NERTST,2) = TAG(2)
27300 ERRNAM(NERTST,3)=TAG(3)
27400 NERSAV(NERTST) = NERRH
27500 IF (NERTST .NE. 1 .AND. NCOVAR .NE. 0) WRITE(IAUX2)(SSERR(1,J),
27600 1J = 2,NVCVAR)
27700 190 NERR = NERR - NCOVAR
27800 IF (NERR .LT. NVAR) GO TO 820
27900 IF (WITHIN .OR. NESKIP .EQ. 4) GO TO 200
28000 READ(IAUX2) ((SSRES(J,K),J=1,NVART),K=1,NVARTP)
28100 WITHIN = .TRUE.
28200 REWIND IAUX2
28300 200 DO 270 I = 1,NVCVAR
28400 DO 270 J = I,NVCVAR
28500 GO TO (210,220,230,240), NESKIP
28600 210 SUM = SSWITH(I,J)
28700 GO TO 260
28800 220 SUM = SSRES(J,I)
28900 GO TO 260
29000 230 SUM = SSWITH(I,J) + SSRES(J,I)
29100 GO TO 260
29200 240 SUM = 0.0
29300 DO 250 K = NSTART,NFIN
29400 250 SUM = SUM+ORTHES(K,I)*ORTHES(K,J)
29500 260 SSERR(I,J) = SUM
29600 270 SSE(I,J) = SUM
29700 DFERR = NERR
29800 IF (NCOVAR .EQ. 0) GO TO 330
29900 CALL UPRFCT(NCOVAR,SSERR(NVAR+1,NVAR+1) ,M)
30000 IF (M .EQ. 0) GO TO 300
30100 M = NVAR+M
30200 280 WRITE(IOUT,290) TAG, VARNAM(1,M), VARNAM(2,M)
30300 290 FORMAT(1H0/1H03A4, 53H HAS A POSSIBLE LINEAR DEPENDENCY INVOLVING
30400 1 VARIABLE A6,A4,51H AND PRECEEDING COVARIATES FOR THE CURRENT ORDE
30500 2RING )
30600 300 CALL UTIRT(NCOVAR,NVAR,SSERR(NVAR+1,NVAR+1), SSERR(1,NVAR+1) )
30700 DO 320 I = 1,NVAR
30800 DO 320 J = I,NVAR
30900 SUM = 0.0
31000 DO 310 K = NVAR1,NVCVAR
31100 310 SUM = SUM + SSERR(I,K)*SSERR(J,K)
31200 SSHYP(J,I) = SUM
31300 320 SSERR(I,J) = SSERR(I,J) - SUM
31400 CALL UIRT(NCOVAR,NVAR,SSERR(NVAR+1,NVAR+1), SSERR(1,NVAR+1) )
31500 REGRES = .TRUE.
31600 CALL ALPHAN(3,TAGR, 12HREGRESSION )
31700 C
31800 330 IF (.NOT. EONLY) GO TO 335
31900 IF(NVAR .NE. 1)CALL PRINTE(HOLD,ESTHLD,NERRH,0,1)
32000 GO TO 760
32100 335 FAC = 1.0/SQRT(DFERR)
32200 IF (NVAR .EQ. 1)GO TO 410
32300 DO 340 I = 1,NVAR
32400 SD(I) = SQRT(AMAX1(0.0,SSERR(I,I) ) )
32500 340 SSEAD(I,I) = SD(I)*FAC
32600 NVARM1 = NVAR-1
32700 DO 350 I = 1,NVARM1
32800 K = I+1
32900 DO 350 J = K,NVAR
33000 350 SSEAD(J,I) = SSERR(I,J)/(SD(I)*SD(J) )
33100 WRITE(IOUT,360)TAG,NCOVAR
33200 360 FORMAT (1H1 3A4,75H CORRELATIONS OF CRITERIA WITH STANDARD DEVIATI
33300 1ONS ON DIAGONAL ADJUSTED FOR I3,11H COVARIATES )
33400 KEND = 0
33500 370 KBEG = KEND + 1
33600 KEND = KBEG + 9
33700 KEND = MIN0(KEND,NVAR)
33800 WRITE(IOUT,380) (VARNAM(1,J),VARNAM(2,J), J = KBEG,KEND)
33900 380 FORMAT (1H0/10H0VARIABLE 11X,10(A6,A4,1X) )
34000 DO 390 I = KBEG,NVAR
34100 KENDR = MIN0(I,KEND)
34200 390 WRITE(IOUT,400)VARNAM(1,I), VARNAM(2,I), (SSEAD(I,J), J = KBEG,
34300 1KENDR)
34400 400 FORMAT (1H A6,A4,8X, 10F11.3)
34500 IF (KEND .LT. NVAR) GO TO 370
34600 410 DO 420 I = 1,NVAR
34700 DO 420 J = I,NVAR
34800 420 SSEAD(J,I) = SSERR(I,J)
34900 IF (NVAR .NE. 1) GO TO 440
35000 MS = SSERR(1,1)/DFERR
35100 WRITE(IOUT,430) TAG, SSERR(1,1), NERR, MS
35200 430 FORMAT (1H03A4, 11X F15.3, I6, F15.3)
35300 GO TO 460
35400 440 CALL UPRFCT(NVAR,SSERR,M)
35500 IF (M .NE. 0) WRITE(IOUT,450) TAG, VARNAM(1,M) ,
35600 1VARNAM(2,M)
35700 450 FORMAT(1H0/1H03A4, 53H HAS A POSSIBLE LINEAR DEPENDENCY INVOLVING
35800 1 VARIABLE A6,A4,49H AND PRECEDING VARIATES FOR THE CURRENT ORDERI
35900 2NG )
36000 CALL PRINTE(HOLD,ESTHLD,NERRH,0,1)
36100 460 DO 750 LP = 1,NTESTS
36200 M = 0
36300 L = LP
36400 IF (NVAR .NE. 1) L = NTESTS - LP + 1
36500 470 IF (REGRES) GO TO 590
36600 NE = NERRS(L)
36700 480 IF (NE .NE. NERRH) GO TO 730
36800 NSTART = NDFCUM(L)
36900 NFIN = NSTART-1+NDFTST(L)
37000 490 CALL UNPAK(HEAD(1,L),HOLD,NTABLE)
37100 IF (NSTART .NE. 0 .OR. WITHIN) GO TO 500
37200 READ(IAUX2) ((SSRES(J,K),J=1,NVART),K=1,NVARTP)
37300 WITHIN = .TRUE.
37400 REWIND IAUX2
37500 C COMPUTE HYPOTHESIS OR GET RESIDUAL IF REQUIRED.
37600 500 DO 550 I = 1,NVCVAR
37700 DO 550 J = I,NVCVAR
37800 IF (NSTART .NE. 0) GO TO 510
37900 SUM = SSRES(J,I)
38000 GO TO 530
38100 510 SUM = 0.0
38200 DO 520 K = NSTART,NFIN
38300 520 SUM = SUM + ORTHES(K,I)*ORTHES(K,J)
38400 530 IF (NCOVAR .NE. 0)GO TO 540
38500 SSHYP(J,I) = SUM
38600 GO TO 550
38700 540 SSHYP(J,I) = SUM + SSE(I,J)
38800 550 CONTINUE
38900 IF (NCOVAR .EQ. 0) GO TO 580
39000 CALL LWRFCT(NCOVAR,SSHYP(NVAR+1,NVAR+1),M)
39100 CALL LIXR(NCOVAR,NVAR,SSHYP(NVAR+1,NVAR+1), SSHYP(NVAR+1,1) )
39200 DO 570 I = 1,NVAR
39300 DO 570 J = I,NVAR
39400 SUM = 0.0
39500 DO 560 K = NVAR1,NVCVAR
39600 560 SUM = SUM + SSHYP(K,I)*SSHYP(K,J)
39700 570 SSHYP(J,I) = SSHYP(J,I) - SUM - SSEAD(J,I)
39800 580 NHYP = NDFTST(L)
39900 590 IF (NVAR .EQ. 1) GO TO 680
40000 IF (REGRES) GO TO 610
40100 WRITE(IOUT,600)HOLD
40200 600 FORMAT (9H1TEST OF 50A1)
40300 GO TO 630
40400 610 WRITE(IOUT,620)TAG,TAGR
40500 620 FORMAT (9H1TEST OF 3A4,1X,3A4)
40600 NSTART = 0
40700 630 IF (.NOT. PRTSSH) GO TO 639
40800 WRITE(IOUT,631) NCOVAR
40900 631 FORMAT (46H0SUMS OF PRODUCTS FOR HYPOTHESIS ADJUSTED FOR I3,1X,
41000 111HCOVARIATES )
41100 KEND = 0
41200 362 KBEG = KEND + 1
41300 KEND = KBEG + 9
41400 KEND = MIN0(KEND,NVAR)
41500 WRITE(IOUT,363) (VARNAM(1,J),VARNAM(2,J),J=KBEG,KEND)
41600 363 FORMAT (1H0/10H0VARIABLE 11X,10(A6,A4,1X) )
41700 DO 367 I = KBEG,NVAR
41800 KENDR = MIN0(I,KEND)
41900 367 WRITE(IOUT,364) VARNAM(2,I),VARNAM(1,I),(SSHYP(I,J), J= KBEG,
42000 1KENDR)
42100 364 FORMAT (1H A6,A4,8X, 10E11.3)
42200 IF (KEND .LT. NVAR) GO TO 362
42300 IF (.NOT. REGRES) WRITE(IOUT,600)HOLD
42400 IF (REGRES) WRITE(IOUT,620)TAG,TAGR
42500 639 CALL SIGTST(NHYP,NERR,SD,NSTART)
42600 WITHIN = .FALSE.
42700 IF (.NOT. REGRES) GO TO 730
42800 KEND = 0
42900 640 KBEG = KEND+1
43000 KEND = KBEG + 7
43100 KEND = MIN0(NVAR,KEND)
43200 WRITE(IOUT,650) (VARNAM(1,J),VARNAM(2,J),J = KBEG,KEND)
43300 650 FORMAT (1H0/29H0RAW REGRESSION COEFFICIENTS /30X 8HVARIATES /
43400 111H COVARIATES 9X 8(4X A6,A4) )
43500 DO 660 J = NVAR1,NVCVAR
43600 660 WRITE(IOUT,670) VARNAM(1,J), VARNAM(2,J), (SSERR(K,J),K=KBEG,KEND)
43700 670 FORMAT (1H A6,A4,8X, 8F14.3)
43800 IF (KEND .LT. NVAR) GO TO 640
43900 GO TO 720
44000 680 DFHYP = NHYP
44100 MS = SSHYP(1,1)/DFHYP
44200 F = SSHYP(1,1)* DFERR/ (SSEAD(1,1)*DFHYP)
44300 FPRO=FPROB(DFHYP,DFERR,F)
44400 PROB = AMAX1(.001,FPRO )
44500 IF (REGRES) GO TO 700
44600 WRITE(IOUT,690)HOLD, SSHYP(1,1), NHYP, MS, F,PROB
44700 690 FORMAT(1H 18A1, 5X, F15.3, I6,3F15.3)
44800 GO TO 730
44900 700 WRITE(IOUT,710)TAGR,SSHYP(1,1),NHYP,MS,F,PROB
45000 710 FORMAT (1H 3A4,11X, F15.3, I6,3F15.3)
45100 720 REGRES = .FALSE.
45200 GO TO 470
45300 730 IF (M .EQ. 0) GO TO 750
45400 M = NVAR + M
45500 WRITE(IOUT,740)VARNAM(1,M), VARNAM(2,M)
45600 740 FORMAT (1H0/43H0LOSS OF ACCURACY IN ABOVE DUE TO VARIABLE A6,A4,
45700 151H AND PRECEDING COVARIATES FOR THE CURRENT ORDERING )
45800 750 CONTINUE
45900 760 GO TO (110,120,130,130), NESKIP
46000 770 IF (NVAR .NE. 1) GO TO 810
46100 CALL PRINTE(HOLD,ESTHLD,NERRH,NERSAV,NERTST)
46200 IF (EONLY) GO TO 810
46300 IF (NCOVAR .EQ. 0) GO TO 810
46400 N2 = 0
46500 780 N1 = N2+1
46600 N2 = MIN0(N1+7,NERTST)
46700 WRITE(IOUT,790)(ERRNAM(J,1), ERRNAM(J,2),ERRNAM(J,3), J=N1,N2)
46800 790 FORMAT (1H0/28H0RAW REGRESSION COEFFICIENTS/11H0COVARIATES 8X,
46900 18(2X 3A4) )
47000 DO 800 J = 2,NVCVAR
47100 800 WRITE(IOUT,670) VARNAM(1,J), VARNAM(2,J),( SSERR(K,J),K=N1,N2)
47200 IF (N2 .LT. NERTST) GO TO 780
47300 810 RETURN
47400 820 I=1
47500 IF (NVAR .EQ. 1) I = 0
47600 WRITE(IOUT,830)I,TAG
47700 830 FORMAT (I1,15HTOO FEW DF IN 3A4)
47800 GO TO 760
47900 END
48000 SUBROUTINE UPRFCT(N,A,M)
48100 DIMENSION A(25,25)
48200 M=0
48300 N1=N-1
48400 IF(N1) 160,100,100
48500 100 DO 220 K=1,N
48600 AKK=A(K,K)
48700 IF(K .EQ. 1) GO TO 120
48800 DO 110 J=2,K
48900 110 AKK=AKK-A(J-1,K)**2
49000 120 IF(A(K,K)) 140,140,130
49100 130 IF(AKK/A(K,K) .GE. .001) GO TO 150
49200 140 M=K
49300 AKK=0.
49400 150 AKK=SQRT(AKK)
49500 A(K,K)=AKK
49600 IF(K .NE. N) GO TO 170
49700 160 RETURN
49800 170 DO 220 I=K,N1
49900 AKI=A(K,I+1)
50000 IF(K .EQ. 1) GO TO 190
50100 DO 180 J=2,K
50200 180 AKI=AKI-A(J-1,K)*A(J-1,I+1)
50300 190 IF(AKK) 210,200,210
50400 200 A(K,I+1)= 0.0
50500 GO TO 220
50600 210 A(K,I+1)=AKI/AKK
50700 220 CONTINUE
50800 RETURN
50900 END
51000 SUBROUTINE UTISUI (N,A,B)
51100 DIMENSION A(25,25),B(25,25)
51200 DO 200 I=1,N
51300 I1 = I-1
51400 DO 130 J=I,N
51500 IF(I1) 120,120,100
51600 100 DO 110 K=1,I1
51700 110 A(J,I) = A(J,I) - B(K,I)*A(J,K)
51800 120 A(J,I) = A(J,I)/B(I,I)
51900 130 IF (B(I,I) .EQ. 0.0) A(J,I) = 0.0
52000 DO 200 J=1,I
52100 J1 = J-1
52200 IF(J1) 160,160,140
52300 140 DO 150 K=1,J1
52400 150 A(I,J) = A(I,J) - B(K,I)*A(J,K)
52500 160 IF(I1-J) 190,170,170
52600 170 DO 180 K=J,I1
52700 180 A(I,J) = A(I,J) - B(K,I)*A(K,J)
52800 190 A(I,J) = A(I,J)/B(I,I)
52900 200 IF (B(I,I) .EQ. 0.0) A(I,J) = 0.0
53000 RETURN
53100 END
53200 SUBROUTINE UIRT (M,N,S,B)
53300 DIMENSION S(25,25),B(25,25)
53400 DO 150 J=1,N
53500 100 I=M
53600 110 SUM = B(J,I)
53700 REC = 0.0
53800 IF (S(I,I) .NE. 0.) REC = 1./S(I,I)
53900 IF(M-I) 140,140,120
54000 120 IP1 = I+1
54100 DO 130 K=IP1,M
54200 130 SUM = SUM - S(I,K)*B(J,K)
54300 140 SUM = SUM*REC
54400 B(J,I) = SUM
54500 I=I-1
54600 IF(I) 150,150,110
54700 150 CONTINUE
54800 RETURN
54900 END
55000 SUBROUTINE PRINTE( HOLD,ESTHLD,NERRH,NERSAV,NERTST)
55100 DIMENSION HOLD(18), ESTHLD(50) , SSERR(25,25),NERSAV(11),
55200 1SSRES(26,27),Q(500),GM(50)
55300 COMMON ORTHES(100,20),DUMMY(26,27),SSHYP(25,26),SSEAD(40,41),
55400 1ESTIM(50,50), NUMERR(8),NERRS(100),NDFCUM(100),NDFTST(
55500 2100),VARNAM(2,50),HEAD(3,100),LEVEL(8),LEVSUB(8,10),LEVCUM(8,10),
55600 3NCELL(100), NTABLE(27),ITABLE(9,9),OBS(100)
55700 4 ,NVAR,NCOVAR,NERWIT,NERRES,HNUM(100),ERROR,NTESTS,RVARC,
55800 5 FIRST,IORD(50),IPOSV(50),IVPOS(50),NCELLS,NVART,NDFTOT,
55900 6SPECOR,VLIST,PRINTR,NFACT , READK, PRINTK,CONTR ,TESTR, MFIRST,
56000 7TRUTH,BLANK,MAXFAC,MAXCEL,MAXPAR,MAXLEV,MAXVAR,ATITLE,AJOBCD,
56100 8AANALY,AFINIS,WITHIN,SPACE(10)
56200 DATA INP,IOUT,IAUX1,IAUX2/5,30,22,23/
56300 DATA IAUX4/24/
56400 EQUIVALENCE (SSERR,SSHYP(1,2) ),(SSRES,DUMMY)
56500 1,(Q,DUMMY(203)), (GM,Q(0301)), (PRTCLL,SPACE(4))
56600 DOUBLE PRECISION VARNAM
56700 LOGICAL PRTCLL,SPWITH
56800 DATA W /4HW /
56900 NVAR1 = NVAR + 1
57000 NVCVAR = NVAR + NCOVAR
57100 IF (NERTST .EQ. 1 .OR. NCOVAR .EQ. 0) GO TO 110
57200 REWIND IAUX2
57300 DO 100 J = 2,NVCVAR
57400 100 SSERR(NERTST,J) = SSERR(1,J)
57500 110 KEND = 0
57600 IF (NVAR .LE. 10) WRITE(IOUT,120)NCOVAR
57700 120 FORMAT ( 23H0ESTIMATES ADJUSTED FOR I3,11H COVARIATES )
57800 130 KBEG = KEND+1
57900 KEND = KBEG + 7
58000 KEND = MIN0(NVAR,KEND)
58100 IF ((NDFTOT .GT. 10 .AND. KBEG .GT. 1).OR.(KBEG .EQ. 1 .AND. NVAR
58200 1.GT. 10))WRITE(IOUT,140)NCOVAR
58300 140 FORMAT (23H1ESTIMATES ADJUSTED FOR I3,11H COVARIATES )
58400 WRITE(IOUT,150) (VARNAM(1,J),VARNAM(2,J),J = KBEG,KEND)
58500 150 FORMAT (1H0 45X 8HCRITERIA / 9H CONTRAST 11X 8(4XA6,A4))
58600 DO 210 LT = 1,NERTST
58700 IF (NVAR .NE. 1) GO TO 160
58800 NERRH = NERSAV(LT)
58900 IF (LT .NE. NERTST .AND. NCOVAR .NE. 0) READ(IAUX2)(SSERR(LT,J),
59000 1J = 2,NVCVAR)
59100 160 DO 210 L = 1,NTESTS
59200 IF (NERRS(L) .GT. 10 .OR. NERRS(L) .NE. NERRH) GO TO 210
59300 NSTART = NDFCUM(L)
59400 NFIN = NSTART - 1 + NDFTST(L)
59500 IF (NSTART .EQ. 0) GO TO 210
59600 CALL UNPAK(HEAD(1,L), HOLD, NTABLE)
59700 SPWITH = HOLD(1) .EQ. W
59800 WRITE(IOUT,170)HOLD
59900 170 FORMAT (1H 50A1)
60000 L1=0
60100 DO 205 I = NSTART,NFIN
60200 L1 = L1 + 1
60300 DO 190 J = 1,NVAR
60400 ESTHLD(J) = ESTIM(I,J)
60500 IF (NCOVAR .EQ. 0) GO TO 190
60600 JT = J+LT-1
60700 CORR = 0.0
60800 DO 180 K = NVAR1,NVCVAR
60900 IF (SPWITH .AND. PRTCLL) CORR = GM(K)
61000 180 ESTHLD(J) = ESTHLD(J) - (ESTIM(I,K) - CORR) * SSERR(JT,K)
61100 190 CONTINUE
61200 WRITE(IOUT,200) L1, (ESTHLD(K), K = KBEG,KEND)
61300 200 FORMAT (1H I6,12X 8F14.3)
61400 205 CONTINUE
61500 210 CONTINUE
61600 IF (KEND .LT. NVAR) GO TO 130
61700 RETURN
61800 END
61900 SUBROUTINE LWRFCT(N,A,M)
62000 DIMENSION A(25,25)
62100 M=0
62200 N1=N-1
62300 IF(N1) 160,100,100
62400 100 DO 220 K=1,N
62500 AKK=A(K,K)
62600 IF(K .EQ. 1) GO TO 120
62700 DO 110 J=2,K
62800 110 AKK=AKK-A(K,J-1)**2
62900 120 IF(A(K,K)) 140,140,130
63000 130 IF(AKK/A(K,K) .GE. .00001) GO TO 150
63100 140 M=K
63200 AKK=0.
63300 150 AKK=SQRT(AKK)
63400 A(K,K)=AKK
63500 IF(K .NE. N) GO TO 170
63600 160 RETURN
63700 170 DO 220 I=K,N1
63800 AIK=A(I+1,K)
63900 IF(K .EQ. 1) GO TO 190
64000 DO 180 J=2,K
64100 180 AIK=AIK-A(K,J-1)*A(I+1,J-1)
64200 190 IF(AKK) 210,200,210
64300 200 A(I+1,K)= 0.0
64400 GO TO 220
64500 210 A(I+1,K)=AIK/AKK
64600 220 CONTINUE
64700 RETURN
64800 END
64900 SUBROUTINE LIXR (M,N,S,B)
65000 DIMENSION S(25,25),B(25,25)
65100 DO 130 J=1,N
65200 DO 130 I=1,M
65300 SUM = B(I,J)
65400 IM1 = I-1
65500 IF(IM1) 120,120,100
65600 100 DO 110 K=1,IM1
65700 110 SUM = SUM-S(I,K)*B(K,J)
65800 120 SUM = SUM/S(I,I)
65900 IF (S(I,I) .EQ. 0.0) SUM = 0.0
66000 130 B(I,J) = SUM
66100 RETURN
66200 END
66300 SUBROUTINE SIGTST(NHYP,NERR,SD,NSTART)
66400 DIMENSION FUVAR(50),SD(50), EIG(50), VEC(26,26), PI(50),VECHLD(50)
66500 1,ESTHLD(50), SSERR(50,50),SAVMS(50),Q(500)
66600 COMMON ORTHES(100,20),DUMMY(26,27),SSHYP(25,26),SSEAD(40,41),
66700 1ESTIM(50,50), NUMERR(8),NERRS(100),NDFCUM(100),NDFTST(
66800 2100),VARNAM(2,50),HEAD(3,100),LEVEL(8),LEVSUB(8,10),LEVCUM(8,10),
66900 3NCELL(100), NTABLE(27),ITABLE(9,9),OBS(100)
67000 4 ,NVAR,NCOVAR,NERWIT,NERRES,HNUM(100),ERROR,NTESTS,RVARC,
67100 5 FIRST,IORD(50),IPOSV(50),IVPOS(50),NCELLS,NVART,NDFTOT,
67200 6SPECOR,VLIST,PRINTR,NFACT , READK, PRINTK,CONTR ,TESTR, MFIRST,
67300 7TRUTH,BLANK,MAXFAC,MAXCEL,MAXPAR,MAXLEV,MAXVAR,ATITLE,AJOBCD,
67400 8AANALY,AFINIS,WITHIN,SPACE(10)
67500 DATA INP,IOUT,IAUX1,IAUX2/5,30,22,23/
67600 DATA IAUX4/24/
67700 EQUIVALENCE (SSERR,SSHYP(1,2)),
67800 1(VEC,DUMMY), (PI,VECHLD), (FUVAR,ESTHLD), (Q,DUMMY(203)),
67900 2(Q(101),EIG), (Q(151),VECHLD), (Q(201),ESTHLD), (Q(251),SAVMS)
68000 DOUBLE PRECISION VARNAM
68100 LOGICAL RVARC
68200 TOL = .15
68300 NVAR1 = NVAR + 1
68400 NVCVAR = NVAR + NCOVAR
68500 NROOTS = MIN0(NVAR,NHYP)
68600 XNVAR = NVAR
68700 XNERR = NERR
68800 FAC = SQRT(XNERR)
68900 XNHYP = NHYP
69000 YNHYP = XNHYP
69100 DFHYP = XNVAR*XNHYP
69200 CON = XNERR/XNHYP
69300 DO 100 I = 1,NVAR
69400 SAVMS(I) = SSHYP(I,I)/XNHYP
69500 100 FUVAR(I) = CON*SSHYP(I,I)/SSEAD(I,I)
69600 CALL UTISUI(NVAR,SSHYP,SSERR)
69700 CALL EIGN(NVAR,SSHYP,EIG,VEC,IND)
69800 IF (IND .NE. 0) WRITE(IOUT,110)
69900 110 FORMAT (40H0LOSS OF ACCURACY IN EIGENVALUE PROBLEM )
70000 CALL UIXR(NVAR,NVAR,SSERR,VEC)
70100 PI(NROOTS) = 1.0 + EIG(NROOTS)
70200 IF (NROOTS .EQ. 1) GO TO 130
70300 NROOT1 = NROOTS - 1
70400 DO 120 I = 1,NROOT1
70500 J = NROOTS - I
70600 IF (EIG(I+1) .LT. 0.0) EIG(I+1) = 0.0
70700 120 PI(J) = PI(J+1)*(1.0+EIG(J) )
70800 130 WRITE(IOUT,140)
70900 140 FORMAT (81H0TESTS OF SIGNIFICANCE USING WILKS LAMBDA CRITERION
71000 1 AND CANONICAL CORRELATIONS /15H TEST OF ROOTS 13X 1HF 9X 5HD
71100 2FHYP 7X 5HDFERR 4X 11HP LESS THAN 5X 1HR)
71200 NSIG = 0
71300 DO 160 I = 1,NROOTS
71400 XK = 1.0
71500 DFH = XNHYP*XNVAR
71600 IF (DFH .NE. 2.0) XK = SQRT( (DFH**2-4.0)/(XNVAR**2+XNHYP**2
71700 1-5.0) )
71800 DFERR = (XK*(2.0*XNERR+XNHYP-XNVAR-1.0) - DFHYP + 2.0)/2.0
71900 CON = DFERR/DFHYP
72000 XLAMB = PI(I)**(-1.0/XK)
72100 F = CON* (1.0-XLAMB)/XLAMB
72200 FPRO=FPROB(DFHYP,DFERR,F)
72300 PROB =AMAX1(.001,FPRO )
72400 CANONR = SQRT(EIG(I)/(1.0+EIG(I) ))
72500 WRITE(IOUT,150) I,NROOTS,F,DFHYP,DFERR,PROB, CANONR
72600 150 FORMAT (I3,8H THROUGH I3,5X,5F12.3)
72700 IF (PROB .LE. TOL) NSIG = NSIG + 1
72800 XNHYP = XNHYP - 1.0
72900 DFSUB = NVAR + NHYP + 1 - 2*I
73000 160 DFHYP = DFHYP - DFSUB
73100 N2 = NSIG
73200 N1 = 1
73300 IF (NSIG .EQ. 0) N2 = 1
73400 IF (NSIG .GT. 8) N2 = 8
73500 170 WRITE(IOUT,180)NHYP,NERR, (I,I = 1,N2)
73600 180 FORMAT (1H0/1H0 27X 19HUNIVARIATE F TESTS 14X 47HSTANDARDIZED DIS
73700 1CRIMINANT FUNCTION COEFFICIENTS /9H VARIABLE 11X 2HF( I2,1H, I5,
73800 21H) 4X 7HMEAN SQ 4X 11HP LESS THAN I7,7I8)
73900 DO 200 I = 1,NVAR
74000 FPRO=FPROB(YNHYP,XNERR,FUVAR(I))
74100 PROBF =AMAX1(.001,FPRO )
74200 DO 190 J = 1,N2
74300 190 VECHLD(J) = SD(I)*VEC(I,J)
74400 200 WRITE(IOUT,210)VARNAM(1,I), VARNAM(2,I), FUVAR(I),SAVMS(I),PROBF,
74500 1(VECHLD(J), J = 1,N2)
74600 210 FORMAT (1H A6,A4, 7X, 3F12.3, 5X, 8F8.3)
74700 IF (NSIG .EQ. 0) GO TO 410
74800 220 IF (NSTART .EQ. 0) GO TO 300
74900 WRITE(IOUT,230) (I,I = N1,N2)
75000 230 FORMAT (1H0/21H0DISCRIMINANT SCORES /9H CONTRAST 47X 8I8)
75100 NFIN = NSTART + NHYP - 1
75200 M=0
75300 DO 280 I = NSTART,NFIN
75400 DO 250 J = 1,NVAR
75500 SUM = ESTIM(I,J)
75600 IF (NCOVAR .EQ. 0) GO TO 250
75700 DO 240 K = NVAR1,NVCVAR
75800 240 SUM = SUM - ESTIM(I,K)*SSERR(J,K)
75900 250 ESTHLD(J) = SUM
76000 M = M+1
76100 DO 270 J = N1,N2
76200 SUM = 0.0
76300 DO 260 K = 1,NVAR
76400 260 SUM = SUM + ESTHLD(K)*VEC(K,J)
76500 270 VECHLD(J) = SUM*FAC
76600 280 WRITE(IOUT,290)M, (VECHLD(J), J =N1,N2)
76700 290 FORMAT (I7, 52X 8F8.3)
76800 300 IF (.NOT. RVARC) GO TO 370
76900 WRITE(IOUT,310)(I,I = N1,N2)
77000 310 FORMAT(1H0/52H0CORRELATIONS BETWEEN VARIABLES AND COMPOSITE SCORES
77100 1/9H VARIABLE 47X 8I8)
77200 DO 350 I = 1,NVAR
77300 I1 = I + 1
77400 DO 340 J = N1,N2
77500 SUM = 0.0
77600 DO 320 K1 = 1,I
77700 320 SUM = SUM + SSEAD(I,K1)*VEC(K1,J)
77800 IF (I1 .GT. NVAR) GO TO 340
77900 DO 330 K1= I1,NVAR
78000 330 SUM = SUM + SSEAD(K1,I)*VEC(K1,J)
78100 340 VECHLD(J) = SUM/SD(I)
78200 350 WRITE(IOUT,360)VARNAM(1,I), VARNAM(2,I), (VECHLD(J), J = N1,N2)
78300 360 FORMAT (1H A6,A4, 48X, 8F8.3)
78400 370 IF (N2 .EQ. NSIG) GO TO 410
78500 N1 = N2 + 1
78600 N2 = MIN0(NSIG,N1+7)
78700 WRITE(IOUT,380) (I,I = N1,N2)
78800 380 FORMAT (1H0/1H0 60X 47HSTANDARDIZED DISCRIMINANT FUNCTION COEFFICI
78900 1ENTS / 9H VARIABLE 47X 8I8)
79000 IF (N2 .GT. NSIG) N2 = NSIG
79100 DO 400 I = 1,NVAR
79200 DO 390 J = N1,N2
79300 390 VECHLD(J) = SD(I)*VEC(I,J)
79400 400 WRITE(IOUT,360)VARNAM(1,I), VARNAM(2,I), (VECHLD(J), J = N1,N2)
79500 GO TO 220
79600 410 RETURN
79700 END
79800 SUBROUTINE EIGN(NN,A,EIG,VEC,IND)
79900 DIMENSION A(25,25),GAMMA(50),BETA(50),BETASQ(50),EIG(50)
80000 DIMENSION W(49),VEC(26,26)
80100 DIMENSION P(49),Q(49)
80200 EQUIVALENCE (P(1),BETA(1)),(Q(1),BETA(1))
80300 DIMENSION IPOSV(50),IVPOS(50),IORD(50)
80400 EQUIVALENCE (IPOSV(1),GAMMA(1)),(IVPOS(1),BETA(1)),
80500 1(IORD(1),BETASQ(1))
80600 N=NN
80700 IND=0
80800 IF(N .EQ. 0) GO TO 560
80900 N1=N-1
81000 N2=N-2
81100 ENORM=0.
81200 TRACE=0.
81300 DO 110 J=1,N
81400 DO 100 I=J,N
81500 100 ENORM=ENORM+A(I,J)**2
81600 TRACE=TRACE+A(J,J)
81700 110 ENORM=ENORM-.5*A(J,J)**2
81800 ENORM=ENORM+ENORM
81900 GAMMA(1)=A(1,1)
82000 IF(N2) 280,270,120
82100 120 DO 260 NR=1,N2
82200 B=A(NR+1,NR)
82300 S=0.
82400 DO 130 I=NR,N2
82500 130 S=S+A(I+2,NR)**2
82600 A(NR+1,NR)=0.
82700 IF(S) 250,250,140
82800 140 S=S+B*B
82900 SGN=+1.
83000 IF(B) 150,160,160
83100 150 SGN=-1.
83200 160 SQRTS=SQRT(S)
83300 D=SGN/(SQRTS+SQRTS)
83400 TEMP=SQRT(.5+B*D)
83500 W(NR)=TEMP
83600 A(NR+1,NR)=TEMP
83700 D=D/TEMP
83800 B=-SGN*SQRTS
83900 DO 170 I=NR,N2
84000 TEMP=D*A(I+2,NR)
84100 W(I+1)=TEMP
84200 170 A(I+2,NR)=TEMP
84300 WTAW=0.
84400 DO 220 I=NR,N1
84500 SUM=0.
84600 DO 180 J=NR,I
84700 180 SUM=SUM+A(I+1,J+1)*W(J)
84800 I1=I+1
84900 IF(N1-I1) 210,190,190
85000 190 DO 200 J=I1,N1
85100 200 SUM=SUM+A(J+1,I+1)*W(J)
85200 210 P(I)=SUM
85300 220 WTAW=WTAW+SUM*W(I)
85400 DO 230 I=NR,N1
85500 230 Q(I)=P(I)-WTAW*W(I)
85600 DO 240 J=NR,N1
85700 QJ=Q(J)
85800 WJ=W(J)
85900 DO 240 I=J,N1
86000 240 A(I+1,J+1)=A(I+1,J+1)-2.*(W(I)*QJ+WJ*Q(I))
86100 250 BETA(NR)=B
86200 BETASQ(NR)=B*B
86300 260 GAMMA(NR+1)=A(NR+1,NR+1)
86400 270 B=A(N,N-1)
86500 BETA(N-1)=B
86600 BETASQ(N-1)=B*B
86700 GAMMA(N)=A(N,N)
86800 280 BETASQ(N)=0.
86900 DO 300 I=1,N
87000 DO 290 J=1,N
87100 290 VEC(I,J)=0.
87200 300 VEC(I,I)=1.
87300 M=N
87400 SUM=0.
87500 NPAS=1
87600 GO TO 400
87700 310 SUM=SUM+SHIFT
87800 COSA=1.
87900 G=GAMMA(1)-SHIFT
88000 PP=G
88100 PPBS=PP*PP+BETASQ(1)
88200 PPBR=SQRT(PPBS)
88300 DO 370 J=1,M
88400 COSAP=COSA
88500 IF(PPBS .NE. 0.) GO TO 320
88600 SINA=0.
88700 SINA2=0.
88800 COSA=1.
88900 GO TO 350
89000 320 SINA=BETA(J)/PPBR
89100 SINA2=BETASQ(J)/PPBS
89200 COSA=PP/PPBR
89300 NT=J+NPAS
89400 IF(NT .GE. N) NT=N
89500 330 DO 340 I=1,NT
89600 TEMP=COSA*VEC(I,J)+SINA*VEC(I,J+1)
89700 VEC(I,J+1)=-SINA*VEC(I,J)+COSA*VEC(I,J+1)
89800 340 VEC(I,J)=TEMP
89900 350 DIA=GAMMA(J+1)-SHIFT
90000 U=SINA2*(G+DIA)
90100 GAMMA(J)=G+U
90200 G=DIA-U
90300 PP=DIA*COSA-SINA*COSAP*BETA(J)
90400 IF(J .NE. M) GO TO 360
90500 BETA(J)=SINA*PP
90600 BETASQ(J)=SINA2*PP*PP
90700 GO TO 380
90800 360 PPBS=PP*PP+BETASQ(J+1)
90900 PPBR=SQRT(PPBS)
91000 BETA(J)=SINA*PPBR
91100 370 BETASQ(J)=SINA2*PPBS
91200 380 GAMMA(M+1)=G
91300 NPAS=NPAS+1
91400 IF(BETASQ(M) .GT. 1.E-21) GO TO 410
91500 390 EIG(M+1)=GAMMA(M+1)+SUM
91600 400 BETA(M)=0.
91700 BETASQ(M)=0.
91800 M=M-1
91900 IF(M .EQ. 0) GO TO 430
92000 IF(BETASQ(M) .LE. 1.E-21) GO TO 390
92100 410 A2=GAMMA(M+1)
92200 R2=.5*A2
92300 R1=.5*GAMMA(M)
92400 R12=R1+R2
92500 DIF=R1-R2
92600 TEMP=SQRT(DIF*DIF+BETASQ(M))
92700 R1=R12+TEMP
92800 R2=R12-TEMP
92900 DIF=ABS(A2-R1)-ABS(A2-R2)
93000 IF(DIF .LT. 0.) GO TO 420
93100 SHIFT=R2
93200 GO TO 310
93300 420 SHIFT=R1
93400 GO TO 310
93500 430 EIG(1)=GAMMA(1)+SUM
93600 DO 440 J=1,N
93700 IPOSV(J)=J
93800 IVPOS(J)=J
93900 440 IORD(J)=J
94000 M=N
94100 GO TO 470
94200 450 DO 460 J=1,M
94300 IF(EIG(J) .GE. EIG(J+1)) GO TO 460
94400 TEMP=EIG(J)
94500 EIG(J)=EIG(J+1)
94600 EIG(J+1)=TEMP
94700 ITEMP=IORD(J)
94800 IORD(J)=IORD(J+1)
94900 IORD(J+1)=ITEMP
95000 460 CONTINUE
95100 470 M=M-1
95200 IF(M .NE. 0) GO TO 450
95300 IF(N1 .EQ. 0) GO TO 500
95400 DO 490 L=1,N1
95500 NV=IORD(L)
95600 NP=IPOSV(NV)
95700 IF(NP .EQ. L) GO TO 490
95800 LV=IVPOS(L)
95900 IVPOS(NP)=LV
96000 IPOSV(LV)=NP
96100 DO 480 I=1,N
96200 TEMP=VEC(I,L)
96300 VEC(I,L)=VEC(I,NP)
96400 480 VEC(I,NP)=TEMP
96500 490 CONTINUE
96600 500 ESUM=0.
96700 ESSQ=0.
96800 DO 550 NRR=1,N
96900 K=N1
97000 510 K=K-1
97100 IF(K .LE. 0) GO TO 540
97200 SUM=0.
97300 DO 520 I=K,N1
97400 520 SUM=SUM+VEC(I+1,NRR)*A(I+1,K)
97500 SUM=SUM+SUM
97600 DO 530 I=K,N1
97700 530 VEC(I+1,NRR)=VEC(I+1,NRR)-SUM*A(I+1,K)
97800 GO TO 510
97900 540 ESUM=ESUM+EIG(NRR)
98000 550 ESSQ=ESSQ+EIG(NRR)**2
98100 TEMP=ABS(512.*TRACE)
98200 IF((ABS(TRACE-ESUM)+TEMP)-TEMP .NE. 0.) IND=IND+1
98300 TEMP=1024.*ENORM
98400 IF((ABS(ENORM-ESSQ)+TEMP)-TEMP .NE. 0.) IND=IND+2
00100 560 RETURN
00200 END
00300 SUBROUTINE UIXR (M,N,S,B)
00400 DIMENSION S(25,25),B(26,26)
00500 DO 150 J=1,N
00600 100 I=M
00700 110 SUM = B(I,J)
00800 REC = 0.0
00900 IF (S(I,I) .NE. 0.) REC = 1./S(I,I)
01000 IF(M-I) 140,140,120
01100 120 IP1 = I+1
01200 DO 130 K=IP1,M
01300 130 SUM = SUM - S(I,K)*B(K,J)
01400 140 SUM = SUM*REC
01500 B(I,J) = SUM
01600 I=I-1
01700 IF(I) 150,150,110
01800 150 CONTINUE
01900 RETURN
02000 END
02100 SUBROUTINE UTIRT(M,N,S,B)
02200 DIMENSION S(25,25),B(25,25)
02300 DO 130 J=1,N
02400 DO 130 I=1,M
02500 SUM = B(J,I)
02600 IM1 = I-1
02700 IF(IM1) 120,120,100
02800 100 DO 110 K=1,IM1
02900 110 SUM = SUM-S(K,I)*B(J,K)
03000 120 SUM = SUM/S(I,I)
03100 IF (S(I,I) .EQ. 0.0) SUM = 0.0
03200 130 B(J,I) = SUM
03300 RETURN
03400 END
04700 FUNCTION FPROB(UV,VV,FV)
04800 U=UV
04900 V=VV
05000 F=FV
05100 IND=0
05200 IF(F.GT.1.)GO TO 10
05300 IND=1
05400 F=1./F
05500 TEMP=U
05600 U=V
05700 V=TEMP
05800 10 Z1=2./(9.*U)
05900 Z2=2./(9.*V)
06000 Z=ABS((1.-Z2)*F**(.33333333)-1+Z1)
06100 Z=Z/SQRT(Z2*F**(.66666667)+Z1)
06200 IF(V.GT.4.)GO TO 30
06300 Z=Z*(1.+.08*Z**4/V**3)
06400 30 Z=(1.+Z*(.196854+Z*(.115194+Z*(.000344+Z*.019527))))**4
06500 FPROB=.5/Z
06600 IF(IND.EQ.1)FPROB=1.-FPROB
06700 RETURN
06800 END