Trailing-Edge
-
PDP-10 Archives
-
decuslib10-09
-
43,50466/manov4.f4
There are no other files named manov4.f4 in the archive.
C WMU IMPLEMENTATION OF THE WAYNE STATE VERSION OF
C
C THE MIAMI MULTIVARIATE ANALYSIS OF VARIANCE(MANOVA)
C
C MODIFIED FOR WMU COMPUTER CENTER BY: SAM ANEMA
C
C FILES: MANOVA.F4,MANOV1.F4,MANOV2.F4,MANOV3.F4,MANOV4.F4
C
C ADDITIONAL SUBROUTINES: USAGE - RUSS BARR(IN NGLIB).MAC
C
C LOADING PROCEDURE:
C
C R LOADER
C MANOV=MANOVA,USAGE/<MANOV1/>MANOV2/>MANOV3/>MANOV4/>/G
C SAV MANOVA
C
C MANOV4.F4
SUBROUTINE ANALYS
DIMENSION VECHLD(100),INLET(75),SSRES(26,27)
COMMON ORTHES(100,20),DUMMY(26,27),SSHYP(25,26),SSEAD(40,41),
1ESTIM(50,50), NUMERR(8),NERRS(100),NDFCUM(100),NDFTST(
2100),VARNAM(2,50),HEAD(3,100),LEVEL(8),LEVSUB(8,10),LEVCUM(8,10),
3NCELL(100), NTABLE(27),ITABLE(9,9),OBS(100)
4 ,NVAR,NCOVAR,NERWIT,NERRES,HNUM(100),ERROR,NTESTS,RVARC,
5 FIRST,IORD(50),IPOSV(50),IVPOS(50),NCELLS,NVART,NDFTOT,
6SPECOR,VLIST,PRINTR,NFACT , READK, PRINTK,CONTR ,TESTR, MFIRST,
7TRUTH,BLANK,MAXFAC,MAXCEL,MAXPAR,MAXLEV,MAXVAR,ATITLE,AJOBCD,
8AANALY,AFINIS,WITHIN,SPACE(10)
DATA INP,IOUT,IAUX1,IAUX2/5,30,22,23/
DATA IAUX4/24/
EQUIVALENCE (INLET,SSEAD), (VECHLD,SSHYP),(SSRES,DUMMY),
1(TAG,SPACE(1) ), (EONLY,SPACE(2) ),(NSIG,SPACE(7)),(TAGP,SPACE(3))
2,(PRTCLL,SPACE(4)),(PRTSSH,SPACE(5)),(NPCT,SPACE(6))
LOGICAL ERROR,TRUTH,VLIST,SPECOR,READK,CONTR,TESTR,MFIRST,WITHIN
2,EONLY,PRINTK,PRINTR,PRTCLL,PRTSSH,NLIST
INTEGER BLANK
DOUBLE PRECISION VARNAM,VECHLD
REWIND IAUX2
NVARTP = NVART + 1
IF (VLIST) GO TO 120
GO TO 130
100 IF (WITHIN) GO TO 110
REWIND IAUX2
READ(IAUX2)((SSRES(J,K),J = 1,NVART),K = 1,NVARTP)
110 IF (.NOT. VLIST) GO TO 140
120 CALL REORD
IF (ERROR) GO TO 210
130 IF (NVAR .EQ. 1) GO TO 140
REWIND IAUX2
WRITE(IAUX2)((SSRES(J,K),J = 1,NVART),K = 1,NVARTP)
140 REWIND IAUX2
WITHIN = .TRUE.
CALL SETTST
IF (ERROR) GO TO 210
IF (J .EQ. 1) WRITE (IOUT,150)
150 FORMAT (52H0THERE HAS BEEN A DIVISION BY ZERO IN THIS ANALYSIS )
NSIG = NSIG-1
IF (NSIG .GT. 0) GO TO 210
WRITE (IOUT,157)
157 FORMAT(1H1/' WESTERN MICHIGAN UNIVERSITY'//,
1 ' MULTI-VARIATE ANALYSIS OF VARIANCE'/)
155 READ(INP,160) TAG,TAGP,INLET
160 FORMAT (A4,76A1)
IF (TAG .NE. ATITLE) GO TO 165
WRITE(IOUT,163) INLET
163 FORMAT (1H 80A1)
GO TO 155
165 TRUTH = TAG .EQ. AANALY
IF (.NOT. TRUTH) GO TO 210
NLIST = INLET(1) .NE. BLANK
SPECOR = INLET(2) .NE. BLANK
CONTR = INLET(3) .NE. BLANK
DO 166 J = 2,8
IF (INLET(J) .NE. BLANK) GO TO 167
166 CONTINUE
GO TO 169
167 TESTR = INLET(4) .NE. BLANK
READK = INLET(5) .NE. BLANK
MFIRST = INLET(6) .EQ. BLANK
EONLY = INLET(7) .NE. BLANK
PRTCLL = INLET(8) .NE. BLANK
EONLY = EONLY .OR. PRTCLL
SPECOR = SPECOR .OR. EONLY
MFIRST = MFIRST .AND. .NOT. PRTCLL
PRINTK = PRINTK .AND. .NOT. EONLY
PRINTR = PRINTR .AND. .NOT. EONLY
PRTSSH = PRTSSH .AND. .NOT. EONLY
169 NSIG = 1
DO 170 I = 1,9
170 IF (INLET(7) .EQ. NTABLE(I+9) .OR. INLET(8) .EQ. NTABLE(I+9))
1NSIG = I
IF (NLIST)READ(INP,180) NVAR,NCOVAR,(IORD(J),J = 1,38)
VLIST = VLIST .OR. NLIST
180 FORMAT (40I2)
NVCVAR = NVAR+NCOVAR
IF(NLIST.AND.NVCVAR.GT. 38) READ(INP,180)(IORD(J), J = 39,NVCVAR)
IF (NCOVAR .EQ. 0) NCOVAR = 0
DO 190 J = 1,NVCVAR
J1 = IORD(J)
K = IPOSV(J1)
VECHLD (2*J-1) = VARNAM(1,K)
IF (NVCVAR .GT. NVART) GO TO 220
190 VECHLD(2*J) = VARNAM(2,K)
N2 = 2*NVCVAR
NPCT = NPCT+1
NP = NPCT/1000
NR = NPCT - NP*1000
WRITE(IOUT,200)NP,NR,NVAR,NCOVAR, (VECHLD(J),J = 1,N2)
200 FORMAT (8H0PROBLEM I4, 11H REANALYSIS I4, 19H WITH THE FOLLOWING
1I4,13H CRITERIA AND I4, 11H COVARIATES /1X10(3X A6,A4) )
IF (.NOT.(SPECOR .OR. CONTR) ) GO TO 100
210 RETURN
220 WRITE(IOUT,230)
230 FORMAT (35H TOO MANY VARIATES AND COVARIATES )
ERROR = .TRUE.
GO TO 210
END
SUBROUTINE REORD
DIMENSION SSWITH(26,26),SSRES(26,27),Q(500),GM(40)
COMMON ORTHES(100,20),DUMMY(26,27),SSHYP(25,26),SSEAD(40,41),
1ESTIM(50,50), NUMERR(8),NERRS(100),NDFCUM(100),NDFTST(
2100),VARNAM(2,50),HEAD(3,100),LEVEL(8),LEVSUB(8,10),LEVCUM(8,10),
3NCELL(100), NTABLE(27),ITABLE(9,9),OBS(100)
4 ,NVAR,NCOVAR,NERWIT,NERRES,HNUM(100),ERROR,NTESTS,RVARC,
5 FIRST,IORD(50),IPOSV(50),IVPOS(50),NCELLS,NVART,NDFTOT,
6SPECOR,VLIST,PRINTR,NFACT , READK, PRINTK,CONTR ,TESTR, MFIRST,
7TRUTH,BLANK,MAXFAC,MAXCEL,MAXPAR,MAXLEV,MAXVAR,ATITLE,AJOBCD,
8AANALY,AFINIS,WITHIN,SPACE(10)
EQUIVALENCE (SSWITH,SSRES(1,2) ),(SSRES,DUMMY)
1,(Q,DUMMY(203)), (GM,Q(301))
DATA INP,IOUT,IAUX1,IAUX2/5,30,22,23/
DATA IAUX4/24/
DOUBLE PRECISION VARNAM,XHOLDD
LOGICAL ERROR
NVCVAR = NVAR + NCOVAR
DO 190 I = 1,NVCVAR
NV = IORD(I)
IF (NV .EQ. 0 .OR. NV .GT. NVART) GO TO 210
N1 = IPOSV(NV)
IF (N1 .EQ. I) GO TO 190
K = IVPOS(I)
IVPOS(I) = NV
IVPOS(N1) = K
IPOSV(NV) = I
IPOSV(K) = N1
N2 = I
NMIN = MIN0(N1,N2)
N = NMIN-1
IF (N .EQ. 0) GO TO 110
DO 100 J = 1,N
XHOLD = SSWITH(J,N2)
SSWITH(J,N2) = SSWITH(J,N1)
SSWITH(J,N1) = XHOLD
XHOLD = SSRES(N2,J)
SSRES(N2,J) = SSRES(N1,J)
100 SSRES(N1,J) = XHOLD
110 NMAX = MAX0(N1,N2)
N = NMAX+1
IF (N .GT. NVART) GO TO 130
DO 120 J = N,NVART
XHOLD = SSWITH(N1,J)
SSWITH(N1,J) = SSWITH(N2,J)
SSWITH(N2,J) = XHOLD
XHOLD = SSRES(J,N1)
SSRES(J,N1) = SSRES(J,N2)
120 SSRES(J,N2) = XHOLD
130 IF (NMAX-NMIN .EQ. 1) GO TO 150
M1 = NMIN + 1
M2 = NMAX - 1
DO 140 J = M1,M2
XHOLD = SSWITH(J ,NMAX)
SSWITH(J ,NMAX) = SSWITH(NMIN,J )
SSWITH(NMIN,J ) = XHOLD
XHOLD = SSRES(NMAX,J )
SSRES(NMAX,J ) = SSRES(J ,NMIN)
140 SSRES(J ,NMIN) = XHOLD
150 XHOLD = SSWITH(NMAX,NMAX)
SSWITH(NMAX,NMAX) = SSWITH(NMIN,NMIN)
SSWITH(NMIN,NMIN) = XHOLD
XHOLD = SSRES (NMAX,NMAX)
SSRES (NMAX,NMAX) = SSRES (NMIN,NMIN)
SSRES (NMIN,NMIN) = XHOLD
DO 160 K = 1,2
XHOLDD = VARNAM(K,N2)
VARNAM(K,N2) = VARNAM(K,N1)
160 VARNAM(K,N1) = XHOLDD
DO 170 J = 1,NCELLS
XHOLD = ORTHES(J,N2)
ORTHES(J,N2) = ORTHES(J,N1)
170 ORTHES(J,N1) = XHOLD
DO 180 J = 1,NCELLS
XHOLD = ESTIM(J,N2)
ESTIM(J,N2) = ESTIM(J,N1)
180 ESTIM(J,N1) = XHOLD
XHOLD = GM(N2)
GM(N2) = GM(N1)
GM(N1) = XHOLD
190 CONTINUE
200 RETURN
210 ERROR = .TRUE.
WRITE(IOUT,220) NV
220 FORMAT (9H0VARIABLE I4, 15H DOES NOT EXIST )
GO TO 200
END
SUBROUTINE SETTST
DIMENSION TAGR(3),ERRNAM(11,3),HOLD(18),SD(50),TAG(3),
1SSERR(25,25),ESTHLD(50),SSE(40,40) ,SSWITH(26,26),NERSAV(11),
2SSRES(26,27),Q(500)
DATA INP,IOUT,IAUX1,IAUX2/5,30,22,23/
DATA IAUX4/24/
COMMON ORTHES(100,20),DUMMY(26,27),SSHYP(25,26),SSEAD(40,41),
1ESTIM(50,50), NUMERR(8),NERRS(100),NDFCUM(100),NDFTST(
2100),VARNAM(2,50),HEAD(3,100),LEVEL(8),LEVSUB(8,10),LEVCUM(8,10),
3NCELL(100), NTABLE(27),ITABLE(9,9),OBS(100)
4 ,NVAR,NCOVAR,NERWIT,NERRES,HNUM(100),ERROR,NTESTS,RVARC,
5 FIRST,IORD(50),IPOSV(50),IVPOS(50),NCELLS,NVART,NDFTOT,
6SPECOR,VLIST,PRINTR,NFACT , READK, PRINTK,CONTR ,TESTR, MFIRST,
7TRUTH,BLANK,MAXFAC,MAXCEL,MAXPAR,MAXLEV,MAXVAR,ATITLE,AJOBCD,
8AANALY,AFINIS,WITHIN,SPACE(10)
EQUIVALENCE (SSERR,SSHYP(1,2)), (SSE,SSEAD(1,2) ) ,
1 (SSWITH,SSRES(1,2) ),(SSRES,DUMMY), (EONLY,SPACE(2) )
2,(Q,DUMMY(203)),(Q,SD), (Q(51),ESTHLD),(PRTSSH,SPACE(5))
LOGICAL REGRES,ERROR,RVARC,WITHIN,EONLY,PRTSSH
REAL MS,NTABLE
DOUBLE PRECISION VARNAM
NVARTP = NVART + 1
NERTST = 0
NVCVAR = NVAR + NCOVAR
NVAR1 = NVAR+1
REGRES = .FALSE.
KERR = 0
IF (NVAR .EQ. 1 .AND. .NOT. EONLY) WRITE(IOUT,100)
100 FORMAT (7H1SOURCE 27X 2HSS 7X 2HDF 10X 2HMS 14X 1HF10X 11HP LE
1SS THAN )
NERRH = 0
NHYP= NCOVAR
NERR = NERWIT
CALL ALPHAN(3,TAG,12HWITHIN CELLS)
NESKIP = 1
GO TO 160
110 NERRH = 10
CALL ALPHAN(3,TAG,12HRESIDUAL )
NHYP= NCOVAR
NERR = NERRES
NESKIP = 2
GO TO 160
120 NERRH = 9
CALL ALPHAN(3,TAG,12HWITHIN+RESID )
NHYP = NCOVAR
NERR = NERRES+NERWIT
NESKIP = 3
GO TO 160
130 IF (KERR .EQ. 8) GO TO 770
140 KERR = KERR+1
IF (NUMERR(KERR) .EQ. 0) GO TO 130
NERRH = KERR
NCOD = NUMERR(KERR)
CALL ALPHAN(2,TAG,8HERROR )
TAG(3) = HNUM(KERR)
NHYP= NCOVAR
NERR = NDFTST(NCOD)
NESKIP = 4
150 NSTART = NDFCUM(NCOD)
NFIN = NSTART-1+NDFTST(NCOD)
160 DO 170 J = 1,NTESTS
IF (NERRS(J) .EQ. NERRH) GO TO 180
170 CONTINUE
GO TO 760
180 IF (NVAR .GT. 1) GO TO 190
NERTST = NERTST + 1
ERRNAM(NERTST,1) = TAG(1)
ERRNAM(NERTST,2) = TAG(2)
ERRNAM(NERTST,3)=TAG(3)
NERSAV(NERTST) = NERRH
IF (NERTST .NE. 1 .AND. NCOVAR .NE. 0) WRITE(IAUX2)(SSERR(1,J),
1J = 2,NVCVAR)
190 NERR = NERR - NCOVAR
IF (NERR .LT. NVAR) GO TO 820
IF (WITHIN .OR. NESKIP .EQ. 4) GO TO 200
READ(IAUX2) ((SSRES(J,K),J=1,NVART),K=1,NVARTP)
WITHIN = .TRUE.
REWIND IAUX2
200 DO 270 I = 1,NVCVAR
DO 270 J = I,NVCVAR
GO TO (210,220,230,240), NESKIP
210 SUM = SSWITH(I,J)
GO TO 260
220 SUM = SSRES(J,I)
GO TO 260
230 SUM = SSWITH(I,J) + SSRES(J,I)
GO TO 260
240 SUM = 0.0
DO 250 K = NSTART,NFIN
250 SUM = SUM+ORTHES(K,I)*ORTHES(K,J)
260 SSERR(I,J) = SUM
270 SSE(I,J) = SUM
DFERR = NERR
IF (NCOVAR .EQ. 0) GO TO 330
CALL UPRFCT(NCOVAR,SSERR(NVAR+1,NVAR+1) ,M)
IF (M .EQ. 0) GO TO 300
M = NVAR+M
280 WRITE(IOUT,290) TAG, VARNAM(1,M), VARNAM(2,M)
290 FORMAT(1H0/1H03A4, 53H HAS A POSSIBLE LINEAR DEPENDENCY INVOLVING
1 VARIABLE A6,A4,51H AND PRECEEDING COVARIATES FOR THE CURRENT ORDE
2RING )
300 CALL UTIRT(NCOVAR,NVAR,SSERR(NVAR+1,NVAR+1), SSERR(1,NVAR+1) )
DO 320 I = 1,NVAR
DO 320 J = I,NVAR
SUM = 0.0
DO 310 K = NVAR1,NVCVAR
310 SUM = SUM + SSERR(I,K)*SSERR(J,K)
SSHYP(J,I) = SUM
320 SSERR(I,J) = SSERR(I,J) - SUM
CALL UIRT(NCOVAR,NVAR,SSERR(NVAR+1,NVAR+1), SSERR(1,NVAR+1) )
REGRES = .TRUE.
CALL ALPHAN(3,TAGR, 12HREGRESSION )
C
330 IF (.NOT. EONLY) GO TO 335
IF(NVAR .NE. 1)CALL PRINTE(HOLD,ESTHLD,NERRH,0,1)
GO TO 760
335 FAC = 1.0/SQRT(DFERR)
IF (NVAR .EQ. 1)GO TO 410
DO 340 I = 1,NVAR
SD(I) = SQRT(AMAX1(0.0,SSERR(I,I) ) )
340 SSEAD(I,I) = SD(I)*FAC
NVARM1 = NVAR-1
DO 350 I = 1,NVARM1
K = I+1
DO 350 J = K,NVAR
350 SSEAD(J,I) = SSERR(I,J)/(SD(I)*SD(J) )
WRITE(IOUT,360)TAG,NCOVAR
360 FORMAT (1H1 3A4,75H CORRELATIONS OF CRITERIA WITH STANDARD DEVIATI
1ONS ON DIAGONAL ADJUSTED FOR I3,11H COVARIATES )
KEND = 0
370 KBEG = KEND + 1
KEND = KBEG + 9
KEND = MIN0(KEND,NVAR)
WRITE(IOUT,380) (VARNAM(1,J),VARNAM(2,J), J = KBEG,KEND)
380 FORMAT (1H0/10H0VARIABLE 11X,10(A6,A4,1X) )
DO 390 I = KBEG,NVAR
KENDR = MIN0(I,KEND)
390 WRITE(IOUT,400)VARNAM(1,I), VARNAM(2,I), (SSEAD(I,J), J = KBEG,
1KENDR)
400 FORMAT (1H A6,A4,8X, 10F11.3)
IF (KEND .LT. NVAR) GO TO 370
410 DO 420 I = 1,NVAR
DO 420 J = I,NVAR
420 SSEAD(J,I) = SSERR(I,J)
IF (NVAR .NE. 1) GO TO 440
MS = SSERR(1,1)/DFERR
WRITE(IOUT,430) TAG, SSERR(1,1), NERR, MS
430 FORMAT (1H03A4, 11X F15.3, I6, F15.3)
GO TO 460
440 CALL UPRFCT(NVAR,SSERR,M)
IF (M .NE. 0) WRITE(IOUT,450) TAG, VARNAM(1,M) ,
1VARNAM(2,M)
450 FORMAT(1H0/1H03A4, 53H HAS A POSSIBLE LINEAR DEPENDENCY INVOLVING
1 VARIABLE A6,A4,49H AND PRECEDING VARIATES FOR THE CURRENT ORDERI
2NG )
CALL PRINTE(HOLD,ESTHLD,NERRH,0,1)
460 DO 750 LP = 1,NTESTS
M = 0
L = LP
IF (NVAR .NE. 1) L = NTESTS - LP + 1
470 IF (REGRES) GO TO 590
NE = NERRS(L)
480 IF (NE .NE. NERRH) GO TO 730
NSTART = NDFCUM(L)
NFIN = NSTART-1+NDFTST(L)
490 CALL UNPAK(HEAD(1,L),HOLD,NTABLE)
IF (NSTART .NE. 0 .OR. WITHIN) GO TO 500
READ(IAUX2) ((SSRES(J,K),J=1,NVART),K=1,NVARTP)
WITHIN = .TRUE.
REWIND IAUX2
C COMPUTE HYPOTHESIS OR GET RESIDUAL IF REQUIRED.
500 DO 550 I = 1,NVCVAR
DO 550 J = I,NVCVAR
IF (NSTART .NE. 0) GO TO 510
SUM = SSRES(J,I)
GO TO 530
510 SUM = 0.0
DO 520 K = NSTART,NFIN
520 SUM = SUM + ORTHES(K,I)*ORTHES(K,J)
530 IF (NCOVAR .NE. 0)GO TO 540
SSHYP(J,I) = SUM
GO TO 550
540 SSHYP(J,I) = SUM + SSE(I,J)
550 CONTINUE
IF (NCOVAR .EQ. 0) GO TO 580
CALL LWRFCT(NCOVAR,SSHYP(NVAR+1,NVAR+1),M)
CALL LIXR(NCOVAR,NVAR,SSHYP(NVAR+1,NVAR+1), SSHYP(NVAR+1,1) )
DO 570 I = 1,NVAR
DO 570 J = I,NVAR
SUM = 0.0
DO 560 K = NVAR1,NVCVAR
560 SUM = SUM + SSHYP(K,I)*SSHYP(K,J)
570 SSHYP(J,I) = SSHYP(J,I) - SUM - SSEAD(J,I)
580 NHYP = NDFTST(L)
590 IF (NVAR .EQ. 1) GO TO 680
IF (REGRES) GO TO 610
WRITE(IOUT,600)HOLD
600 FORMAT (9H1TEST OF 50A1)
GO TO 630
610 WRITE(IOUT,620)TAG,TAGR
620 FORMAT (9H1TEST OF 3A4,1X,3A4)
NSTART = 0
630 IF (.NOT. PRTSSH) GO TO 639
WRITE(IOUT,631) NCOVAR
631 FORMAT (46H0SUMS OF PRODUCTS FOR HYPOTHESIS ADJUSTED FOR I3,1X,
111HCOVARIATES )
KEND = 0
362 KBEG = KEND + 1
KEND = KBEG + 9
KEND = MIN0(KEND,NVAR)
WRITE(IOUT,363) (VARNAM(1,J),VARNAM(2,J),J=KBEG,KEND)
363 FORMAT (1H0/10H0VARIABLE 11X,10(A6,A4,1X) )
DO 367 I = KBEG,NVAR
KENDR = MIN0(I,KEND)
367 WRITE(IOUT,364) VARNAM(2,I),VARNAM(1,I),(SSHYP(I,J), J= KBEG,
1KENDR)
364 FORMAT (1H A6,A4,8X, 10E11.3)
IF (KEND .LT. NVAR) GO TO 362
IF (.NOT. REGRES) WRITE(IOUT,600)HOLD
IF (REGRES) WRITE(IOUT,620)TAG,TAGR
639 CALL SIGTST(NHYP,NERR,SD,NSTART)
WITHIN = .FALSE.
IF (.NOT. REGRES) GO TO 730
KEND = 0
640 KBEG = KEND+1
KEND = KBEG + 7
KEND = MIN0(NVAR,KEND)
WRITE(IOUT,650) (VARNAM(1,J),VARNAM(2,J),J = KBEG,KEND)
650 FORMAT (1H0/29H0RAW REGRESSION COEFFICIENTS /30X 8HVARIATES /
111H COVARIATES 9X 8(4X A6,A4) )
DO 660 J = NVAR1,NVCVAR
660 WRITE(IOUT,670) VARNAM(1,J), VARNAM(2,J), (SSERR(K,J),K=KBEG,KEND)
670 FORMAT (1H A6,A4,8X, 8F14.3)
IF (KEND .LT. NVAR) GO TO 640
GO TO 720
680 DFHYP = NHYP
MS = SSHYP(1,1)/DFHYP
F = SSHYP(1,1)* DFERR/ (SSEAD(1,1)*DFHYP)
FPRO=FPROB(DFHYP,DFERR,F)
PROB = AMAX1(.001,FPRO )
IF (REGRES) GO TO 700
WRITE(IOUT,690)HOLD, SSHYP(1,1), NHYP, MS, F,PROB
690 FORMAT(1H 18A1, 5X, F15.3, I6,3F15.3)
GO TO 730
700 WRITE(IOUT,710)TAGR,SSHYP(1,1),NHYP,MS,F,PROB
710 FORMAT (1H 3A4,11X, F15.3, I6,3F15.3)
720 REGRES = .FALSE.
GO TO 470
730 IF (M .EQ. 0) GO TO 750
M = NVAR + M
WRITE(IOUT,740)VARNAM(1,M), VARNAM(2,M)
740 FORMAT (1H0/43H0LOSS OF ACCURACY IN ABOVE DUE TO VARIABLE A6,A4,
151H AND PRECEDING COVARIATES FOR THE CURRENT ORDERING )
750 CONTINUE
760 GO TO (110,120,130,130), NESKIP
770 IF (NVAR .NE. 1) GO TO 810
CALL PRINTE(HOLD,ESTHLD,NERRH,NERSAV,NERTST)
IF (EONLY) GO TO 810
IF (NCOVAR .EQ. 0) GO TO 810
N2 = 0
780 N1 = N2+1
N2 = MIN0(N1+7,NERTST)
WRITE(IOUT,790)(ERRNAM(J,1), ERRNAM(J,2),ERRNAM(J,3), J=N1,N2)
790 FORMAT (1H0/28H0RAW REGRESSION COEFFICIENTS/11H0COVARIATES 8X,
18(2X 3A4) )
DO 800 J = 2,NVCVAR
800 WRITE(IOUT,670) VARNAM(1,J), VARNAM(2,J),( SSERR(K,J),K=N1,N2)
IF (N2 .LT. NERTST) GO TO 780
810 RETURN
820 I=1
IF (NVAR .EQ. 1) I = 0
WRITE(IOUT,830)I,TAG
830 FORMAT (I1,15HTOO FEW DF IN 3A4)
GO TO 760
END
SUBROUTINE UPRFCT(N,A,M)
DIMENSION A(25,25)
M=0
N1=N-1
IF(N1) 160,100,100
100 DO 220 K=1,N
AKK=A(K,K)
IF(K .EQ. 1) GO TO 120
DO 110 J=2,K
110 AKK=AKK-A(J-1,K)**2
120 IF(A(K,K)) 140,140,130
130 IF(AKK/A(K,K) .GE. .001) GO TO 150
140 M=K
AKK=0.
150 AKK=SQRT(AKK)
A(K,K)=AKK
IF(K .NE. N) GO TO 170
160 RETURN
170 DO 220 I=K,N1
AKI=A(K,I+1)
IF(K .EQ. 1) GO TO 190
DO 180 J=2,K
180 AKI=AKI-A(J-1,K)*A(J-1,I+1)
190 IF(AKK) 210,200,210
200 A(K,I+1)= 0.0
GO TO 220
210 A(K,I+1)=AKI/AKK
220 CONTINUE
RETURN
END
SUBROUTINE UTISUI (N,A,B)
DIMENSION A(25,25),B(25,25)
DO 200 I=1,N
I1 = I-1
DO 130 J=I,N
IF(I1) 120,120,100
100 DO 110 K=1,I1
110 A(J,I) = A(J,I) - B(K,I)*A(J,K)
120 A(J,I) = A(J,I)/B(I,I)
130 IF (B(I,I) .EQ. 0.0) A(J,I) = 0.0
DO 200 J=1,I
J1 = J-1
IF(J1) 160,160,140
140 DO 150 K=1,J1
150 A(I,J) = A(I,J) - B(K,I)*A(J,K)
160 IF(I1-J) 190,170,170
170 DO 180 K=J,I1
180 A(I,J) = A(I,J) - B(K,I)*A(K,J)
190 A(I,J) = A(I,J)/B(I,I)
200 IF (B(I,I) .EQ. 0.0) A(I,J) = 0.0
RETURN
END
SUBROUTINE UIRT (M,N,S,B)
DIMENSION S(25,25),B(25,25)
DO 150 J=1,N
100 I=M
110 SUM = B(J,I)
REC = 0.0
IF (S(I,I) .NE. 0.) REC = 1./S(I,I)
IF(M-I) 140,140,120
120 IP1 = I+1
DO 130 K=IP1,M
130 SUM = SUM - S(I,K)*B(J,K)
140 SUM = SUM*REC
B(J,I) = SUM
I=I-1
IF(I) 150,150,110
150 CONTINUE
RETURN
END
SUBROUTINE PRINTE( HOLD,ESTHLD,NERRH,NERSAV,NERTST)
DIMENSION HOLD(18), ESTHLD(50) , SSERR(25,25),NERSAV(11),
1SSRES(26,27),Q(500),GM(50)
COMMON ORTHES(100,20),DUMMY(26,27),SSHYP(25,26),SSEAD(40,41),
1ESTIM(50,50), NUMERR(8),NERRS(100),NDFCUM(100),NDFTST(
2100),VARNAM(2,50),HEAD(3,100),LEVEL(8),LEVSUB(8,10),LEVCUM(8,10),
3NCELL(100), NTABLE(27),ITABLE(9,9),OBS(100)
4 ,NVAR,NCOVAR,NERWIT,NERRES,HNUM(100),ERROR,NTESTS,RVARC,
5 FIRST,IORD(50),IPOSV(50),IVPOS(50),NCELLS,NVART,NDFTOT,
6SPECOR,VLIST,PRINTR,NFACT , READK, PRINTK,CONTR ,TESTR, MFIRST,
7TRUTH,BLANK,MAXFAC,MAXCEL,MAXPAR,MAXLEV,MAXVAR,ATITLE,AJOBCD,
8AANALY,AFINIS,WITHIN,SPACE(10)
DATA INP,IOUT,IAUX1,IAUX2/5,30,22,23/
DATA IAUX4/24/
EQUIVALENCE (SSERR,SSHYP(1,2) ),(SSRES,DUMMY)
1,(Q,DUMMY(203)), (GM,Q(0301)), (PRTCLL,SPACE(4))
DOUBLE PRECISION VARNAM
LOGICAL PRTCLL,SPWITH
DATA W /4HW /
NVAR1 = NVAR + 1
NVCVAR = NVAR + NCOVAR
IF (NERTST .EQ. 1 .OR. NCOVAR .EQ. 0) GO TO 110
REWIND IAUX2
DO 100 J = 2,NVCVAR
100 SSERR(NERTST,J) = SSERR(1,J)
110 KEND = 0
IF (NVAR .LE. 10) WRITE(IOUT,120)NCOVAR
120 FORMAT ( 23H0ESTIMATES ADJUSTED FOR I3,11H COVARIATES )
130 KBEG = KEND+1
KEND = KBEG + 7
KEND = MIN0(NVAR,KEND)
IF ((NDFTOT .GT. 10 .AND. KBEG .GT. 1).OR.(KBEG .EQ. 1 .AND. NVAR
1.GT. 10))WRITE(IOUT,140)NCOVAR
140 FORMAT (23H1ESTIMATES ADJUSTED FOR I3,11H COVARIATES )
WRITE(IOUT,150) (VARNAM(1,J),VARNAM(2,J),J = KBEG,KEND)
150 FORMAT (1H0 45X 8HCRITERIA / 9H CONTRAST 11X 8(4XA6,A4))
DO 210 LT = 1,NERTST
IF (NVAR .NE. 1) GO TO 160
NERRH = NERSAV(LT)
IF (LT .NE. NERTST .AND. NCOVAR .NE. 0) READ(IAUX2)(SSERR(LT,J),
1J = 2,NVCVAR)
160 DO 210 L = 1,NTESTS
IF (NERRS(L) .GT. 10 .OR. NERRS(L) .NE. NERRH) GO TO 210
NSTART = NDFCUM(L)
NFIN = NSTART - 1 + NDFTST(L)
IF (NSTART .EQ. 0) GO TO 210
CALL UNPAK(HEAD(1,L), HOLD, NTABLE)
SPWITH = HOLD(1) .EQ. W
WRITE(IOUT,170)HOLD
170 FORMAT (1H 50A1)
L1=0
DO 205 I = NSTART,NFIN
L1 = L1 + 1
DO 190 J = 1,NVAR
ESTHLD(J) = ESTIM(I,J)
IF (NCOVAR .EQ. 0) GO TO 190
JT = J+LT-1
CORR = 0.0
DO 180 K = NVAR1,NVCVAR
IF (SPWITH .AND. PRTCLL) CORR = GM(K)
180 ESTHLD(J) = ESTHLD(J) - (ESTIM(I,K) - CORR) * SSERR(JT,K)
190 CONTINUE
WRITE(IOUT,200) L1, (ESTHLD(K), K = KBEG,KEND)
200 FORMAT (1H I6,12X 8F14.3)
205 CONTINUE
210 CONTINUE
IF (KEND .LT. NVAR) GO TO 130
RETURN
END
SUBROUTINE UNPAK(IOUT,IN,NTABLE)
DIMENSION IN(21), IOUT(3), NTABLE(27)
I=1
DO 100 K=1,3
ITEMP=IOUT(K)
NSCALE=33554432
DO 100 J=1,6
INI=ITEMP/NSCALE
ITEMP=ITEMP-INI*NSCALE
IN(I)=NTABLE(20)
IF (INI .NE. 0) IN(I) = NTABLE(INI)
NSCALE=NSCALE/32
100 I=I+1
RETURN
END
SUBROUTINE LWRFCT(N,A,M)
DIMENSION A(25,25)
M=0
N1=N-1
IF(N1) 160,100,100
100 DO 220 K=1,N
AKK=A(K,K)
IF(K .EQ. 1) GO TO 120
DO 110 J=2,K
110 AKK=AKK-A(K,J-1)**2
120 IF(A(K,K)) 140,140,130
130 IF(AKK/A(K,K) .GE. .00001) GO TO 150
140 M=K
AKK=0.
150 AKK=SQRT(AKK)
A(K,K)=AKK
IF(K .NE. N) GO TO 170
160 RETURN
170 DO 220 I=K,N1
AIK=A(I+1,K)
IF(K .EQ. 1) GO TO 190
DO 180 J=2,K
180 AIK=AIK-A(K,J-1)*A(I+1,J-1)
190 IF(AKK) 210,200,210
200 A(I+1,K)= 0.0
GO TO 220
210 A(I+1,K)=AIK/AKK
220 CONTINUE
RETURN
END
SUBROUTINE LIXR (M,N,S,B)
DIMENSION S(25,25),B(25,25)
DO 130 J=1,N
DO 130 I=1,M
SUM = B(I,J)
IM1 = I-1
IF(IM1) 120,120,100
100 DO 110 K=1,IM1
110 SUM = SUM-S(I,K)*B(K,J)
120 SUM = SUM/S(I,I)
IF (S(I,I) .EQ. 0.0) SUM = 0.0
130 B(I,J) = SUM
RETURN
END
SUBROUTINE SIGTST(NHYP,NERR,SD,NSTART)
DIMENSION FUVAR(50),SD(50), EIG(50), VEC(26,26), PI(50),VECHLD(50)
1,ESTHLD(50), SSERR(50,50),SAVMS(50),Q(500)
COMMON ORTHES(100,20),DUMMY(26,27),SSHYP(25,26),SSEAD(40,41),
1ESTIM(50,50), NUMERR(8),NERRS(100),NDFCUM(100),NDFTST(
2100),VARNAM(2,50),HEAD(3,100),LEVEL(8),LEVSUB(8,10),LEVCUM(8,10),
3NCELL(100), NTABLE(27),ITABLE(9,9),OBS(100)
4 ,NVAR,NCOVAR,NERWIT,NERRES,HNUM(100),ERROR,NTESTS,RVARC,
5 FIRST,IORD(50),IPOSV(50),IVPOS(50),NCELLS,NVART,NDFTOT,
6SPECOR,VLIST,PRINTR,NFACT , READK, PRINTK,CONTR ,TESTR, MFIRST,
7TRUTH,BLANK,MAXFAC,MAXCEL,MAXPAR,MAXLEV,MAXVAR,ATITLE,AJOBCD,
8AANALY,AFINIS,WITHIN,SPACE(10)
DATA INP,IOUT,IAUX1,IAUX2/5,30,22,23/
DATA IAUX4/24/
EQUIVALENCE (SSERR,SSHYP(1,2)),
1(VEC,DUMMY), (PI,VECHLD), (FUVAR,ESTHLD), (Q,DUMMY(203)),
2(Q(101),EIG), (Q(151),VECHLD), (Q(201),ESTHLD), (Q(251),SAVMS)
DOUBLE PRECISION VARNAM
LOGICAL RVARC
TOL = .15
NVAR1 = NVAR + 1
NVCVAR = NVAR + NCOVAR
NROOTS = MIN0(NVAR,NHYP)
XNVAR = NVAR
XNERR = NERR
FAC = SQRT(XNERR)
XNHYP = NHYP
YNHYP = XNHYP
DFHYP = XNVAR*XNHYP
CON = XNERR/XNHYP
DO 100 I = 1,NVAR
SAVMS(I) = SSHYP(I,I)/XNHYP
100 FUVAR(I) = CON*SSHYP(I,I)/SSEAD(I,I)
CALL UTISUI(NVAR,SSHYP,SSERR)
CALL EIGN(NVAR,SSHYP,EIG,VEC,IND)
IF (IND .NE. 0) WRITE(IOUT,110)
110 FORMAT (40H0LOSS OF ACCURACY IN EIGENVALUE PROBLEM )
CALL UIXR(NVAR,NVAR,SSERR,VEC)
PI(NROOTS) = 1.0 + EIG(NROOTS)
IF (NROOTS .EQ. 1) GO TO 130
NROOT1 = NROOTS - 1
DO 120 I = 1,NROOT1
J = NROOTS - I
IF (EIG(I+1) .LT. 0.0) EIG(I+1) = 0.0
120 PI(J) = PI(J+1)*(1.0+EIG(J) )
130 WRITE(IOUT,140)
140 FORMAT (81H0TESTS OF SIGNIFICANCE USING WILKS LAMBDA CRITERION
1 AND CANONICAL CORRELATIONS /15H TEST OF ROOTS 13X 1HF 9X 5HD
2FHYP 7X 5HDFERR 4X 11HP LESS THAN 5X 1HR)
NSIG = 0
DO 160 I = 1,NROOTS
XK = 1.0
DFH = XNHYP*XNVAR
IF (DFH .NE. 2.0) XK = SQRT( (DFH**2-4.0)/(XNVAR**2+XNHYP**2
1-5.0) )
DFERR = (XK*(2.0*XNERR+XNHYP-XNVAR-1.0) - DFHYP + 2.0)/2.0
CON = DFERR/DFHYP
XLAMB = PI(I)**(-1.0/XK)
F = CON* (1.0-XLAMB)/XLAMB
FPRO=FPROB(DFHYP,DFERR,F)
PROB =AMAX1(.001,FPRO )
CANONR = SQRT(EIG(I)/(1.0+EIG(I) ))
WRITE(IOUT,150) I,NROOTS,F,DFHYP,DFERR,PROB, CANONR
150 FORMAT (I3,8H THROUGH I3,5X,5F12.3)
IF (PROB .LE. TOL) NSIG = NSIG + 1
XNHYP = XNHYP - 1.0
DFSUB = NVAR + NHYP + 1 - 2*I
160 DFHYP = DFHYP - DFSUB
N2 = NSIG
N1 = 1
IF (NSIG .EQ. 0) N2 = 1
IF (NSIG .GT. 8) N2 = 8
170 WRITE(IOUT,180)NHYP,NERR, (I,I = 1,N2)
180 FORMAT (1H0/1H0 27X 19HUNIVARIATE F TESTS 14X 47HSTANDARDIZED DIS
1CRIMINANT FUNCTION COEFFICIENTS /9H VARIABLE 11X 2HF( I2,1H, I5,
21H) 4X 7HMEAN SQ 4X 11HP LESS THAN I7,7I8)
DO 200 I = 1,NVAR
FPRO=FPROB(YNHYP,XNERR,FUVAR(I))
PROBF =AMAX1(.001,FPRO )
DO 190 J = 1,N2
190 VECHLD(J) = SD(I)*VEC(I,J)
200 WRITE(IOUT,210)VARNAM(1,I), VARNAM(2,I), FUVAR(I),SAVMS(I),PROBF,
1(VECHLD(J), J = 1,N2)
210 FORMAT (1H A6,A4, 7X, 3F12.3, 5X, 8F8.3)
IF (NSIG .EQ. 0) GO TO 410
220 IF (NSTART .EQ. 0) GO TO 300
WRITE(IOUT,230) (I,I = N1,N2)
230 FORMAT (1H0/21H0DISCRIMINANT SCORES /9H CONTRAST 47X 8I8)
NFIN = NSTART + NHYP - 1
M=0
DO 280 I = NSTART,NFIN
DO 250 J = 1,NVAR
SUM = ESTIM(I,J)
IF (NCOVAR .EQ. 0) GO TO 250
DO 240 K = NVAR1,NVCVAR
240 SUM = SUM - ESTIM(I,K)*SSERR(J,K)
250 ESTHLD(J) = SUM
M = M+1
DO 270 J = N1,N2
SUM = 0.0
DO 260 K = 1,NVAR
260 SUM = SUM + ESTHLD(K)*VEC(K,J)
270 VECHLD(J) = SUM*FAC
280 WRITE(IOUT,290)M, (VECHLD(J), J =N1,N2)
290 FORMAT (I7, 52X 8F8.3)
300 IF (.NOT. RVARC) GO TO 370
WRITE(IOUT,310)(I,I = N1,N2)
310 FORMAT(1H0/52H0CORRELATIONS BETWEEN VARIABLES AND COMPOSITE SCORES
1/9H VARIABLE 47X 8I8)
DO 350 I = 1,NVAR
I1 = I + 1
DO 340 J = N1,N2
SUM = 0.0
DO 320 K1 = 1,I
320 SUM = SUM + SSEAD(I,K1)*VEC(K1,J)
IF (I1 .GT. NVAR) GO TO 340
DO 330 K1= I1,NVAR
330 SUM = SUM + SSEAD(K1,I)*VEC(K1,J)
340 VECHLD(J) = SUM/SD(I)
350 WRITE(IOUT,360)VARNAM(1,I), VARNAM(2,I), (VECHLD(J), J = N1,N2)
360 FORMAT (1H A6,A4, 48X, 8F8.3)
370 IF (N2 .EQ. NSIG) GO TO 410
N1 = N2 + 1
N2 = MIN0(NSIG,N1+7)
WRITE(IOUT,380) (I,I = N1,N2)
380 FORMAT (1H0/1H0 60X 47HSTANDARDIZED DISCRIMINANT FUNCTION COEFFICI
1ENTS / 9H VARIABLE 47X 8I8)
IF (N2 .GT. NSIG) N2 = NSIG
DO 400 I = 1,NVAR
DO 390 J = N1,N2
390 VECHLD(J) = SD(I)*VEC(I,J)
400 WRITE(IOUT,360)VARNAM(1,I), VARNAM(2,I), (VECHLD(J), J = N1,N2)
GO TO 220
410 RETURN
END
SUBROUTINE EIGN(NN,A,EIG,VEC,IND)
DIMENSION A(25,25),GAMMA(50),BETA(50),BETASQ(50),EIG(50)
DIMENSION W(49),VEC(26,26)
DIMENSION P(49),Q(49)
EQUIVALENCE (P(1),BETA(1)),(Q(1),BETA(1))
DIMENSION IPOSV(50),IVPOS(50),IORD(50)
EQUIVALENCE (IPOSV(1),GAMMA(1)),(IVPOS(1),BETA(1)),
1(IORD(1),BETASQ(1))
N=NN
IND=0
IF(N .EQ. 0) GO TO 560
N1=N-1
N2=N-2
ENORM=0.
TRACE=0.
DO 110 J=1,N
DO 100 I=J,N
100 ENORM=ENORM+A(I,J)**2
TRACE=TRACE+A(J,J)
110 ENORM=ENORM-.5*A(J,J)**2
ENORM=ENORM+ENORM
GAMMA(1)=A(1,1)
IF(N2) 280,270,120
120 DO 260 NR=1,N2
B=A(NR+1,NR)
S=0.
DO 130 I=NR,N2
130 S=S+A(I+2,NR)**2
A(NR+1,NR)=0.
IF(S) 250,250,140
140 S=S+B*B
SGN=+1.
IF(B) 150,160,160
150 SGN=-1.
160 SQRTS=SQRT(S)
D=SGN/(SQRTS+SQRTS)
TEMP=SQRT(.5+B*D)
W(NR)=TEMP
A(NR+1,NR)=TEMP
D=D/TEMP
B=-SGN*SQRTS
DO 170 I=NR,N2
TEMP=D*A(I+2,NR)
W(I+1)=TEMP
170 A(I+2,NR)=TEMP
WTAW=0.
DO 220 I=NR,N1
SUM=0.
DO 180 J=NR,I
180 SUM=SUM+A(I+1,J+1)*W(J)
I1=I+1
IF(N1-I1) 210,190,190
190 DO 200 J=I1,N1
200 SUM=SUM+A(J+1,I+1)*W(J)
210 P(I)=SUM
220 WTAW=WTAW+SUM*W(I)
DO 230 I=NR,N1
230 Q(I)=P(I)-WTAW*W(I)
DO 240 J=NR,N1
QJ=Q(J)
WJ=W(J)
DO 240 I=J,N1
240 A(I+1,J+1)=A(I+1,J+1)-2.*(W(I)*QJ+WJ*Q(I))
250 BETA(NR)=B
BETASQ(NR)=B*B
260 GAMMA(NR+1)=A(NR+1,NR+1)
270 B=A(N,N-1)
BETA(N-1)=B
BETASQ(N-1)=B*B
GAMMA(N)=A(N,N)
280 BETASQ(N)=0.
DO 300 I=1,N
DO 290 J=1,N
290 VEC(I,J)=0.
300 VEC(I,I)=1.
M=N
SUM=0.
NPAS=1
GO TO 400
310 SUM=SUM+SHIFT
COSA=1.
G=GAMMA(1)-SHIFT
PP=G
PPBS=PP*PP+BETASQ(1)
PPBR=SQRT(PPBS)
DO 370 J=1,M
COSAP=COSA
IF(PPBS .NE. 0.) GO TO 320
SINA=0.
SINA2=0.
COSA=1.
GO TO 350
320 SINA=BETA(J)/PPBR
SINA2=BETASQ(J)/PPBS
COSA=PP/PPBR
NT=J+NPAS
IF(NT .GE. N) NT=N
330 DO 340 I=1,NT
TEMP=COSA*VEC(I,J)+SINA*VEC(I,J+1)
VEC(I,J+1)=-SINA*VEC(I,J)+COSA*VEC(I,J+1)
340 VEC(I,J)=TEMP
350 DIA=GAMMA(J+1)-SHIFT
U=SINA2*(G+DIA)
GAMMA(J)=G+U
G=DIA-U
PP=DIA*COSA-SINA*COSAP*BETA(J)
IF(J .NE. M) GO TO 360
BETA(J)=SINA*PP
BETASQ(J)=SINA2*PP*PP
GO TO 380
360 PPBS=PP*PP+BETASQ(J+1)
PPBR=SQRT(PPBS)
BETA(J)=SINA*PPBR
370 BETASQ(J)=SINA2*PPBS
380 GAMMA(M+1)=G
NPAS=NPAS+1
IF(BETASQ(M) .GT. 1.E-21) GO TO 410
390 EIG(M+1)=GAMMA(M+1)+SUM
400 BETA(M)=0.
BETASQ(M)=0.
M=M-1
IF(M .EQ. 0) GO TO 430
IF(BETASQ(M) .LE. 1.E-21) GO TO 390
410 A2=GAMMA(M+1)
R2=.5*A2
R1=.5*GAMMA(M)
R12=R1+R2
DIF=R1-R2
TEMP=SQRT(DIF*DIF+BETASQ(M))
R1=R12+TEMP
R2=R12-TEMP
DIF=ABS(A2-R1)-ABS(A2-R2)
IF(DIF .LT. 0.) GO TO 420
SHIFT=R2
GO TO 310
420 SHIFT=R1
GO TO 310
430 EIG(1)=GAMMA(1)+SUM
DO 440 J=1,N
IPOSV(J)=J
IVPOS(J)=J
440 IORD(J)=J
M=N
GO TO 470
450 DO 460 J=1,M
IF(EIG(J) .GE. EIG(J+1)) GO TO 460
TEMP=EIG(J)
EIG(J)=EIG(J+1)
EIG(J+1)=TEMP
ITEMP=IORD(J)
IORD(J)=IORD(J+1)
IORD(J+1)=ITEMP
460 CONTINUE
470 M=M-1
IF(M .NE. 0) GO TO 450
IF(N1 .EQ. 0) GO TO 500
DO 490 L=1,N1
NV=IORD(L)
NP=IPOSV(NV)
IF(NP .EQ. L) GO TO 490
LV=IVPOS(L)
IVPOS(NP)=LV
IPOSV(LV)=NP
DO 480 I=1,N
TEMP=VEC(I,L)
VEC(I,L)=VEC(I,NP)
480 VEC(I,NP)=TEMP
490 CONTINUE
500 ESUM=0.
ESSQ=0.
DO 550 NRR=1,N
K=N1
510 K=K-1
IF(K .LE. 0) GO TO 540
SUM=0.
DO 520 I=K,N1
520 SUM=SUM+VEC(I+1,NRR)*A(I+1,K)
SUM=SUM+SUM
DO 530 I=K,N1
530 VEC(I+1,NRR)=VEC(I+1,NRR)-SUM*A(I+1,K)
GO TO 510
540 ESUM=ESUM+EIG(NRR)
550 ESSQ=ESSQ+EIG(NRR)**2
TEMP=ABS(512.*TRACE)
IF((ABS(TRACE-ESUM)+TEMP)-TEMP .NE. 0.) IND=IND+1
TEMP=1024.*ENORM
IF((ABS(ENORM-ESSQ)+TEMP)-TEMP .NE. 0.) IND=IND+2
560 RETURN
END
SUBROUTINE UIXR (M,N,S,B)
DIMENSION S(25,25),B(26,26)
DO 150 J=1,N
100 I=M
110 SUM = B(I,J)
REC = 0.0
IF (S(I,I) .NE. 0.) REC = 1./S(I,I)
IF(M-I) 140,140,120
120 IP1 = I+1
DO 130 K=IP1,M
130 SUM = SUM - S(I,K)*B(K,J)
140 SUM = SUM*REC
B(I,J) = SUM
I=I-1
IF(I) 150,150,110
150 CONTINUE
RETURN
END
SUBROUTINE UTIRT(M,N,S,B)
DIMENSION S(25,25),B(25,25)
DO 130 J=1,N
DO 130 I=1,M
SUM = B(J,I)
IM1 = I-1
IF(IM1) 120,120,100
100 DO 110 K=1,IM1
110 SUM = SUM-S(K,I)*B(J,K)
120 SUM = SUM/S(I,I)
IF (S(I,I) .EQ. 0.0) SUM = 0.0
130 B(J,I) = SUM
RETURN
END
SUBROUTINE ALPHAN(N,NTABLE,IHOLV)
DIMENSION NTABLE (1),IHOLV(20)
CALL OFILE(1,'TST')
WRITE(1,2)(IHOLV(I),I=1,N)
2 FORMAT(20A5)
END FILE 1
CALL IFILE(1,'TST')
READ(1,3)(NTABLE(I),I=1,N)
3 FORMAT(20A4)
END FILE 1
RETURN
END
FUNCTION FPROB(UV,VV,FV)
U=UV
V=VV
F=FV
IND=0
IF(F.GT.1.)GO TO 10
IND=1
F=1./F
TEMP=U
U=V
V=TEMP
10 Z1=2./(9.*U)
Z2=2./(9.*V)
Z=ABS((1.-Z2)*F**(.33333333)-1+Z1)
Z=Z/SQRT(Z2*F**(.66666667)+Z1)
IF(V.GT.4.)GO TO 30
Z=Z*(1.+.08*Z**4/V**3)
30 Z=(1.+Z*(.196854+Z*(.115194+Z*(.000344+Z*.019527))))**4
FPROB=.5/Z
IF(IND.EQ.1)FPROB=1.-FPROB
RETURN
END