MINICERN/mathlib/gen/v/funpre.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1996/04/01 15:02:57  mclareni
   6 * Mathlib gen
   7 *
   8 *
   9 #include "gen/pilot.h"
  10       SUBROUTINE FUNPRE (FUNC,XFCUM,X2LOW,X2HIGH)
  11 C         F. JAMES,    MAY, 1976
  12 C         MODIFIED OCT, 1980 TO ADD PRINTOUT OF INTEGRAL
  13 C         MODIFIED DEC., 1980 TO DELETE LEADING AND TRAILING
  14 C            RANGES OF X WHERE FUNCTION IS ZERO.
  15 C         MODIFIED JUNE,1982 TO FIX POSSIBLE INFINITE LOOP.
  16 C
  17 C         PREPARES THE USER FUNCTION "FUNC" FOR FUNRAN
  18 C         BY FINDING THE PERCENTILES
  19 C         (IN EFFECT, INVERTING THE CUMULATIVE DISTRIBUTION)
  20       EXTERNAL FUNC
  21       COMMON/FUNINT/TFTOT
  22       DIMENSION XFCUM(100)
  23       DATA NBINS / 99/
  24       DATA NZ / 10/
  25       DATA MAXZ / 20/
  26       DATA NITMAX / 6 /
  27       DATA IFUNC/0/
  28       IFUNC = IFUNC + 1
  29 C         FIND MACHINE ACCURACY
  30       COMP1 = 1.0
  31       DO 200 I= 1, 100
  32       COMP1 = COMP1*0.5
  33       COMP2 = 1.0 - COMP1
  34       IF(COMP2 .EQ. 1.0) GOTO 210
  35   200 CONTINUE
  36       COMP1 = 1.0E-10
  37   210 PRECIS = COMP1
  38 C         FIND RANGE WHERE FUNCTION IS NON-ZERO.
  39       CALL FUNZER(FUNC,X2LOW,X2HIGH,XLOW,XHIGH)
  40       XRANGE = XHIGH-XLOW
  41       IF(XRANGE .LE. 0.) GOTO 900
  42       RTEPS = MAX(0.0001,PRECIS*10.)
  43       TFTOT = GAUSS(FUNC,XLOW,XHIGH,RTEPS)
  44 C         PRINT OUT VALUE OF NORMALIZATION INTEGRAL
  45       WRITE(6,1003) IFUNC,XLOW,XHIGH,TFTOT
  46       RTEPS = 0.001
  47       IF(TFTOT .LE. 0.) GOTO 900
  48       TPCTIL = TFTOT/NBINS
  49       TZ = TPCTIL/NZ
  50       TZMAX = TZ * 2.
  51        XFCUM(1) = XLOW
  52       XFCUM(NBINS+1) = XHIGH
  53       X = XLOW
  54       F = FUNC(X)
  55       IF(F .LT. 0.) GOTO 900
  56       NBINM1 = NBINS - 1
  57 C         LOOP OVER BINS (HUNDREDTH PERCENTILES)
  58       DO 600 IBIN = 1, NBINM1
  59       TCUM = 0.
  60       X1 = X
  61       F1 = F
  62       DXMAX = (XHIGH -X) / NZ
  63       FMIN = TZ/DXMAX
  64       FMINZ = FMIN
  65 C         LOOP OVER TRAPEZOIDS WITHIN A SUPPOSED PERCENTILE
  66       DO 500 IZ= 1, MAXZ
  67       XINCR = TZ/MAX(F1,FMIN,FMINZ)
  68   350 X = X1 + XINCR
  69       F = FUNC(X)
  70       IF(F .LT. 0.) GOTO 900
  71       TINCR = (X-X1) * 0.5 * (F+F1)
  72       IF(TINCR .LT. TZMAX) GOTO 370
  73       XINCR = XINCR * 0.5
  74       GOTO 350
  75   370 CONTINUE
  76       TCUM = TCUM + TINCR
  77       IF(TCUM .GE. TPCTIL*0.99) GOTO 520
  78       FMINZ = TZ*F/ (TPCTIL-TCUM)
  79       F1 = F
  80       X1 = X
  81   500 CONTINUE
  82       WRITE(6,2000)
  83  2000 FORMAT('0      FAILURE TO FIND TRAPEZOID   HELP')
  84 C         END OF TRAPEZOID LOOP
  85 C         ADJUST TRAPEZOID INTEGRAL BY GAUSS WITH NEWTON CORRECTION
  86   520 CONTINUE
  87       X1 = XFCUM(IBIN)
  88       XBEST = X
  89       DTBEST = TPCTIL
  90       TPART = TPCTIL
  91 C         ALLOW FOR MAXIMUM NITMAX MORE ITERATIONS ON GAUSS
  92       DO 550 IHOME= 1, NITMAX
  93       XINCR = (TPCTIL-TPART) / MAX(F,FMIN)
  94   535 X = XBEST + XINCR
  95       X2 = X
  96       TPART2 = GAUSS(FUNC,X1,X2,RTEPS)
  97       DTPAR2 = TPART2-TPCTIL
  98       DTABS = ABS(DTPAR2)
  99       IF(ABS(XINCR) .LT. PRECIS) GOTO 545
 100       IF(DTABS .LT. DTBEST) GOTO 545
 101       XINCR = XINCR * 0.5
 102       GOTO 535
 103   545 DTBEST = DTABS
 104       XBEST = X
 105       IF(DTABS .LT. RTEPS*TPCTIL) GOTO 580
 106       TPART = TPART2
 107       F = FUNC(X)
 108   550 CONTINUE
 109       IHOME = NITMAX
 110 C
 111   580 CONTINUE
 112       XFCUM(IBIN+1) = X
 113       F = FUNC(X)
 114       IF(F .LT. 0.) GOTO 900
 115   600 CONTINUE
 116 C         END OF LOOP OVER BINS
 117       X1 = XFCUM(NBINS)
 118       X2 = XHIGH
 119       TPART = GAUSS(FUNC,X1,X2,RTEPS)
 120       ABERR = ABS(TPART-TPCTIL)/TFTOT
 121       WRITE(6,1001) IFUNC,ABERR
 122       IF(ABERR .GT. RTEPS)  WRITE(6,1002)
 123       RETURN
 124   900 WRITE(6,1000) X,F,XLOW,XHIGH
 125       RETURN
 126  1000 FORMAT('0FUNPRE FINDS NEGATIVE FUNCTION VALUE OR RANGE OF X'/
 127      + ,' X=',E15.6,', F=',E15.6,20X,'XLOW=',E15.6,' XHIGH=',E15.6/)
 128  1001 FORMAT(' SUBROUTINE FUNPRE HAS PREPARED USER FUNCTION NUMBER',I4,
 129      + ' FOR FUNRAN'/' MAXIMUM RELATIVE ERROR IN CUMULATIVE',
 130      + 'DISTRIBUTION WILL BE',E15.5)
 131  1002 FORMAT('+',80X,'WARNING,THIS MAY BE TOO BIG'//)
 132  1003 FORMAT('0SUBROUTINE FUNPRE FINDS THE INTEGRAL OF USER FUNCTION',
 133      +  I2,' FROM ',E12.5,' TO ',E12.5,' IS ',E14.6)
 134       END