C     DRIVER PROGRAM FOR GEN_F
C     D. E. RAMIREZ
C     DER@VIRGINIA.   
C     MAY 13, 1999
C     CREATE STORAGE
      INTEGER COUNT, CSIZE, I, IER, MAXR, NTERMS, R
      PARAMETER (CSIZE = 3000, MAXR = 10)
      DOUBLE PRECISION C
      DIMENSION C(0:CSIZE)
      DOUBLE PRECISION CDFERR, CUMF, DENSTY, ERRDEN, G(MAXR), GM, LB,
     +                 NU, PVAL, Y, UB
      INTRINSIC DBLE
C     IMSL SUBROUTINE DFDF
      EXTERNAL DFDF
      PRINT*, 'PROGRAM COMPUTES THE P-VALUES FOR THE '
      PRINT*, 'GENERALIZED F DISTRIBUTION (T/R)/(V/DF)'
 10   CONTINUE
C     INPUT R RANK OF T
      PRINT 800, 'ENTER RANK OF T: 2 <= R <= ', MAXR,'(0 TO EXIT)'
 800  FORMAT(1X,A,I2,2X,A)
      READ*, R
      IF (R.EQ.0) STOP
C     INPUT POSITIVE WEIGHTS
 20   PRINT*, 'ENTER THE POSITIVE WEIGHTS FROM LARGEST TO SMALLEST '
      READ*, (G(I),I=1,R)
      GM=G(1)
      DO 30, I=2, R
         GM=GM*G(I) 
         IF (G(I-1).GE.G(I)) GOTO 30
         PRINT*, 'WEIGHTS ARE NOT SORTED'
         GOTO 20
 30   CONTINUE    
      IF (G(R).LE.0) THEN
         PRINT*, 'WEIGHTS ARE NOT POSITIVE'
         GOTO 20
      ENDIF
C     COMPUTE GEOMETRIC MEAN FOR BOUNDS
      GM=GM**(1.0D0/DBLE(R))
C     INPUT NU
      PRINT*, 'ENTER NU = DEGREES OF FREEDOM FOR V '
      READ*, NU
C     PRINT *, 'ENTER CDF ERROR (SUGGESTED VALUE IS 1E-4)'
C     READ*, CDFERR
      CDFERR=1.0D-4
      PRINT 810, 'CDFERR = ',CDFERR
 40   CONTINUE
      PRINT*
      PRINT*, 'ENTER Y FOR GENERALIZED F DISTRIBUTION (0 TO STOP) '
      READ*, Y
      IF (Y.LE.0) GOTO 10
C     COMPUTE BOUNDS
      LB = DFDF(Y/G(1),DBLE(R),DBLE(NU))
      UB = DFDF(Y/GM,DBLE(R),DBLE(NU))
      PRINT 810, 'BOUNDS FOR P-VALUE', 1.0D0-UB, 1.0D0-LB
 810  FORMAT(/,1X,A,2X,2(F11.7,2X))
C     COMPUTE P-VALUE FOR GENERALIZED F DISTRIBUTION
      IF (Y.NE.0) THEN
         CALL GEN_F(R, G, NU, Y, CDFERR, CUMF, NTERMS, COUNT,
     +            DENSTY, ERRDEN, IER)
         PVAL=1.0D0-CUMF
         IF (IER.EQ.0) THEN
            PRINT 820, 'Y','P-VALUE','NTERMS','DENSITY', 'PDF ERR',
     +                 'EVALS'
 820        FORMAT(/,1X,2(A11,2X),A7,2X,A11,2X,A13,2X,A7)
            PRINT 830, Y,PVAL,NTERMS,DENSTY,ERRDEN,COUNT
 830        FORMAT(/,1X,2(F11.7,2X),I7,2X,F11.7,2X,E13.7,2X,I7)
            GOTO 40
          ELSE
             PRINT*, 'ERROR CODE IS ', IER
             GOTO 10
          ENDIF
      ENDIF
      GOTO 10
      END
************************************************************************
      SUBROUTINE GEN_F(R, G, NU, Y, CDFERR, CUMF, NTERMS, EVALS,
     +               DENSTY, ERRDEN, IER)
C     PROGRAMMED BY D. E. RAMIREZ
C     DER@VIRGINIA.EDU
C     PROGRAM COMPUTES THE LEFT-HAND PROBABILITIE FOR THE
C     GENERALIZED F DISTRIBUTION GIVEN BY F=(T/P)/(V/NU)
C     VERSION MAY 13, 1999
C     IF CSIZE IS CHANGED, THEN CHANGE ALSO IN FUNCTION NEWDEN.
C     IF SMALL IS CHANGED, THEN CHANGE ALSO IN FUNCTION INTERR
      INTEGER CSIZE, MAXP
      DOUBLE PRECISION NEAR1
      PARAMETER (CSIZE = 3000, MAXP = 10, NEAR1 = 0.95D0)
      DOUBLE PRECISION A, AERR, B, C, CDFERR, CUMF, D, DENERR, DENERR0,
     +       DENSTY, DLNGAM, E, E1, ERRDEN, ERROR, G(MAXP), NEWDEN,
     +       NU, NU0, NUPR, NUPRST, ONE, PSTAR, RERR, SMALL, TWO,
     +       VALGSS, VALINT, Y, Y0, ZERO
      DIMENSION C(0:CSIZE), D(CSIZE)
      PARAMETER (SMALL = 1.0D-08)
      INTEGER COUNT, EVALS, I, IER, IRULE, J, L, NTERMS, NTERM0,
     +        P, R, SAVEJ
      INTRINSIC ABS, DBLE, INT, LOG, MIN, MAX, SQRT
      LOGICAL MAXERR
      EXTERNAL NEWDEN
C     IMSL SUBROUTINES
      EXTERNAL DLNGAM, DQDAG, ERSET
      COMMON B, E, NU0, P, NTERM0, C, VALGSS, VALINT, COUNT,
     +       MAXERR, DENERR
      ZERO=0.0D0
      ONE=1.0D0
      TWO=2.0D0
      MAXERR=.TRUE.
C     CONVERT R TO P NEEDED IN COMMON BLOCK FOR NEWDEN
      P=R
      NU0=NU
      IER=0
C     INITIALIZE ERSET TO DISPLAY ONLY TERMINAL ERRORS
      CALL ERSET(3,0,-1)
C     PSTAR IS EVEN
      PSTAR=DBLE(2*((P+1)/2))
C     CHECK THAT P>=2 AND P<=MAXP
      IF (P.LT.2) THEN
         IER=1
         RETURN
      ENDIF
      IF (P.GT.MAXP) THEN
         IER=2
         RETURN
      ENDIF
C     CHECK THAT THE WEIGHTS ARE POSITIVE AND IN DESCENDING ORDER
      DO 20, I=1,P-1
         IF (G(I).LT.G(I+1)) THEN
            IER=3
            RETURN
         ENDIF
 20   CONTINUE
      IF (G(P).LE.ZERO) THEN
         IER=3
         RETURN
      ENDIF
C     CHECK THAT NU IS AT LEAST ONE
      IF (NU.LT.1) THEN
         IER=4
         RETURN
      ENDIF
C     CHECK THAT Y IS POSITIVE
      IF (Y.LE.ZERO) THEN
         IER=5
         RETURN
      ENDIF
C     CHECK THAT CDFERR IS REASONABLE
      IF (CDFERR.LE.ZERO) THEN
         IER=6
         RETURN
      ENDIF
      IF (CDFERR.GT.0.1D0) THEN
         IER=7
         RETURN
      ENDIF
C     CONVERT Y FOR USE WITH T/V
      Y0=P*Y/NU
C     NU+P 
      NUPR=NU0+DBLE(P)
C     NUPRST IS EVEN
      NUPRST=DBLE(2*(INT(NUPR+1)/2))
      B=NEAR1*G(P)
C     EPSILON
      E=ONE-B/G(1)
C     VALUE OF A
      A=ONE
      DO 30, I=1,P
         A=A*(B/G(I))
 30   CONTINUE
      A=SQRT(A)
C     VALUE OF C AND D
      C(0)=A
      E1=ONE-C(0) 
      DO 60, J=1,CSIZE
         D(J)=ZERO
         DO 40, I=1,P
            D(J)=D(J)+(ONE-B/G(I))**J
 40      CONTINUE 
         D(J)=0.5D0*D(J)
         C(J)=ZERO
         DO 50, L=0,J-1
            C(J)=C(J)+D(J-L)*C(L)
 50      CONTINUE
         C(J)=C(J)/DBLE(J)
         E1=E1-C(J)
C        COMPUTE GLOBAL TRUNCATION ERROR BOUND FOR CDF
         DENERR0=DBLE(NU)/B/(DBLE(P)+TWO*(DBLE(J)+ONE))*E1
         IF (DENERR0.LE.CDFERR/Y0) GOTO 70
 60   CONTINUE
C     GLOBAL BOUND FOR (T/P)/(V/NU) HAS BEEN ACHIEVED
 70   CONTINUE
C     SUGGESTED VALUE FOR N
      SAVEJ=MIN(J,CSIZE)
      NTERM0=SAVEJ
C     COMPUTE ONCE AND STORE IN COMMON BLOCK
      VALGSS=LOG(A)+(NTERM0+ONE)*LOG(E)+
     +         DLNGAM((TWO*NTERM0+NU0+P+TWO)/TWO)-LOG(B)-
     +         DLNGAM(DBLE(P)/TWO)-DLNGAM(NU0/TWO)-DLNGAM(NTERM0+TWO)
      VALINT=LOG(A)-DLNGAM(DBLE(P)/TWO)-LOG(B)-
     +         DLNGAM((NU0)/TWO)+DLNGAM(NUPRST/TWO)
C     INITIALIZE FOR DQDAG
      IRULE=2
      AERR=1.0D-5
      RERR=1.0D-5
C     BEGIN TO COMPUTE CUMULATIVE DISTRIBUTION AT Y0
      COUNT=0
C     COMPUTE LOCAL ERROR BOUND ONLY DURING INTEGRATION SUBROUTINE
      MAXERR=.TRUE.
      DENERR=ZERO
      IRULE=2
      CALL DQDAG(NEWDEN,ZERO,Y0,AERR,RERR,IRULE,CUMF,ERROR)
      EVALS=COUNT
      NTERMS=NTERM0
      ERRDEN=MIN(DENERR,DENERR0)
C     SMALL IS PROGRAM EPSILON   
      ERRDEN=MAX(ERRDEN,SMALL)
      MAXERR=.FALSE.
      DENSTY=NEWDEN(Y0)
C     CONVERT T/V TO GENERALIZED F=(T/P)/(V/NU)
      DENSTY=P*DENSTY/NU
      ERRDEN=P*ERRDEN/NU
      RETURN
      END
************************************************************************
      FUNCTION CHECK2(Y,B,E,NTERM0,NU0,P)
      LOGICAL CHECK2
      INTEGER NTERM0, P
      DOUBLE PRECISION B, E, NU0, ONE, T, Y
      INTRINSIC ABS, DBLE, LOG
      ONE=1.0D0
      T=(Y/B)/(ONE+Y/B)
      CHECK2=(DBLE(NTERM0+1).GT.(NU0+P-2.0)/ABS(LOG(E*T)))
      RETURN
      END
************************************************************************
      FUNCTION GSSERR(Y,B,E,NTERM0,NU0,P,VALGSS)
      INTEGER NTERM0, P
      DOUBLE PRECISION B, E, GSSERR, NU0, ONE, TWO, VALGSS, Y
      INTRINSIC DBLE, EXP, LOG
      ONE=1.0D0
      TWO=2.0D0
      GSSERR=VALGSS+(DBLE(P)/TWO+NTERM0)*LOG(Y/B)-
     +         (TWO*NTERM0+NU0+P+TWO)/TWO*LOG(ONE+(ONE-E)*(Y/B))
      GSSERR=EXP(GSSERR)
      RETURN
      END
************************************************************************
      FUNCTION INTERR(Y,B,E,NTERM0,NU0,P,VALINT)
      INTEGER NTERM0, IER, P
      DOUBLE PRECISION B, CS, DCHIDF, E, INTERR, NU0, NUPR, NUPRST, ONE,
     +       PROB, SMALL, TWO, VALET, VALINT,Y, ZERO
      PARAMETER (SMALL = 1.0D-8)
      INTRINSIC ABS, DBLE, EXP, INT, LOG
C     IMSL FUNCTION
      EXTERNAL DCHIDF
      ZERO=0.0D0
      ONE=1.0D0
      TWO=2.0D0
      IF (Y.LE.ZERO) THEN
         INTERR=ZERO
         RETURN
      ENDIF
      NUPR=NU0+P
C     NUPRST IS EVEN
      NUPRST=DBLE(2*(INT(NUPR+1)/2))
      VALET=E*(Y/B)/(ONE+Y/B)
      CS=(TWO*NTERM0+NUPRST-TWO)*ABS(LOG(VALET))
      PROB=DCHIDF(CS,NUPRST)
      PROB=ONE-PROB
      IF (PROB.GT.SMALL) THEN
         INTERR=VALINT+(DBLE(P)/TWO-ONE)*LOG(Y/B)-
     +            ((NU0+DBLE(P))/TWO)*LOG(ONE+Y/B)+
     +            LOG(PROB)+LOG(VALET)-
     +            (NUPRST/TWO)*LOG(VALET*ABS(LOG(VALET)))
         INTERR=EXP(INTERR)
      ELSE
         INTERR=ZERO
      ENDIF
      RETURN
      END
************************************************************************
      FUNCTION DENOFF(Y,XM,NU0)
      INTEGER NTERM0
      DOUBLE PRECISION B, DENOFF, DLNGAM, E, NU0, ONE, TWO,
     +       XM, XN, Y, ZERO
      INTRINSIC EXP, LOG
      EXTERNAL DLNGAM
      ZERO=0.0D0
      ONE=1.0D0
      TWO=2.0D0
      XN=NU0
      IF (XM.GT.TWO) THEN
         IF (Y.LE.ZERO) THEN
            DENOFF=ZERO
            RETURN
         ENDIF
         DENOFF=DLNGAM((XM+XN)/TWO)-DLNGAM(XM/TWO)-DLNGAM(XN/TWO)+
     +            (XM/TWO)*LOG(XM/XN)+((XM-TWO)/TWO)*LOG(Y)-              
     +            ((XM+XN)/TWO)*LOG(ONE+XM*Y/XN)  
      ELSE
         DENOFF=DLNGAM((XM+XN)/TWO)-ZERO-DLNGAM(XN/TWO)+
     +            (XM/TWO)*LOG(XM/XN)+ZERO-
     +            ((XM+XN)/TWO)*LOG(ONE+XM*Y/XN)
      ENDIF
      DENOFF=EXP(DENOFF)
      RETURN
      END
************************************************************************
      FUNCTION NEWDEN(Y)
      INTEGER COUNT, CSIZE, J, NTERM0, P
      PARAMETER (CSIZE = 3000)
      DOUBLE PRECISION B, C(0:CSIZE), DENERR, DENOFF, E, GSSERR,
     +       INTERR, NEWDEN, NU0, TEMPMIN, TWO, VALGSS,
     +       VALINT, XM, XN, Y, ZERO
      LOGICAL MAXERR, CHECK2
      INTRINSIC MAX, MIN
      EXTERNAL CHECK2, DENOFF, GSSERR, INTERR 
      COMMON B, E, NU0, P, NTERM0, C, VALGSS, VALINT, COUNT,
     +       MAXERR, DENERR
      ZERO=0.0D0
      TWO=2.0D0
      XN=NU0
      NEWDEN=ZERO
      DO 10, J=NTERM0,0,-1
         XM=P+TWO*J
         NEWDEN=NEWDEN+C(J)/B*XN/XM*DENOFF(XN*Y/XM/B,XM,NU0)
 10   CONTINUE
      IF (MAXERR) THEN
         IF (CHECK2(Y,B,E,NTERM0,NU0,P)) THEN
            TEMPMIN=MIN(INTERR(Y,B,E,NTERM0,NU0,P,VALINT),
     +                    GSSERR(Y,B,E,NTERM0,NU0,P,VALGSS))
C     TO SEE LOCAL ERROR CALCULATIONS ACTIVATE PRINT STATEMENTS
C     PRINT 700, 'COUNT Y GSSERR INTERR = ', COUNT,Y,
C    +GSSERR(Y,B,E,NTERM0,NU0,P,VALGSS),
C    +INTERR(Y,B,E,NTERM0,NU0,P,VALINT)
         ELSE
            TEMPMIN=GSSERR(Y,B,E,NTERM0,NU0,P,VALGSS)
C     PRINT 700, 'COUNT Y GSSERR = ', COUNT,Y,
C    +GSSERR(Y,B,E,NTERM0,NU0,P,VALGSS)
         ENDIF
C700  FORMAT(1X,A,I7,2X,3(E13.7),2X))
         DENERR=MAX(DENERR, TEMPMIN)
         COUNT=COUNT+1
      ENDIF
      RETURN
      END