MINICERN/mathlib/gen/v/funlux.F

   1       SUBROUTINE FUNLXP (FUNC,XFCUM,X2LOW,X2HIGH)
   2 C         F. JAMES,   Sept, 1994
   3 C
   4 C         Prepares the user function FUNC for FUNLUX
   5 C         Inspired by and mostly copied from FUNPRE and FUNRAN
   6 C         except that
   7 C    1. FUNLUX uses RANLUX underneath,
   8 C    2. FUNLXP expands the first and last bins to cater for
   9 C              functions with long tails on left and/or right,
  10 C    3. FUNLXP calls FUNPCT to do the actual finding of percentiles.
  11 C    4. both FUNLXP and FUNPCT use RADAPT for Gaussian integration.
  12 C
  13       EXTERNAL FUNC
  14       COMMON/FUNINT/TFTOT
  15       DIMENSION XFCUM(200)
  16       PARAMETER (RTEPS=0.0002)
  17       SAVE IFUNC
  18       DATA IFUNC/0/
  19       IFUNC = IFUNC + 1
  20 C         FIND RANGE WHERE FUNCTION IS NON-ZERO.
  21       CALL FUNLZ(FUNC,X2LOW,X2HIGH,XLOW,XHIGH)
  22       XRANGE = XHIGH-XLOW
  23       IF(XRANGE .LE. 0.)  THEN
  24         WRITE (6,'(A,2G15.5)') ' FUNLXP finds function range .LE.0',
  25      +  XLOW,XHIGH
  26         GO TO 900
  27       ENDIF
  28       CALL RADAPT(FUNC,XLOW,XHIGH,1,RTEPS,0.,TFTOT ,UNCERT)
  29       WRITE(6,1003) IFUNC,XLOW,XHIGH,TFTOT
  30  1003 FORMAT(' FUNLXP: integral of USER FUNCTION',
  31      +  I3,' from ',E12.5,' to ',E12.5,' is ',E14.6)
  32 C
  33 C      WRITE (6,'(A,A)') ' FUNLXP preparing ',
  34 C     + 'first the whole range, then left tail, then right tail.'
  35       CALL FUNPCT(FUNC,IFUNC,XLOW,XHIGH,XFCUM,1,99,TFTOT,IERR)
  36       IF (IERR .GT. 0)  GO TO 900
  37       X2 = XFCUM(3)
  38       CALL RADAPT(FUNC,XLOW,X2,1,RTEPS,0.,TFTOT1 ,UNCERT)
  39       CALL FUNPCT(FUNC,IFUNC,XLOW,X2 ,XFCUM,101,49,TFTOT1,IERR)
  40       IF (IERR .GT. 0)  GO TO 900
  41       X3 = XFCUM(98)
  42       CALL RADAPT(FUNC,X3,XHIGH,1,RTEPS,0.,TFTOT2 ,UNCERT)
  43       CALL FUNPCT(FUNC,IFUNC,X3,XHIGH,XFCUM,151,49,TFTOT2,IERR)
  44       IF (IERR .GT. 0)  GO TO 900
  45       WRITE(6,1001) IFUNC,XLOW,XHIGH
  46  1001 FORMAT(' FUNLXP has prepared USER FUNCTION',I3,
  47      + ' between',G12.3,' and',G12.3,' for FUNLUX')
  48       RETURN
  49   900 CONTINUE
  50       WRITE (6,'(A)') ' Fatal error in FUNLXP. FUNLUX will not work.'
  51       END
  52 C
  53       SUBROUTINE FUNPCT(FUNC,IFUNC,XLOW,XHIGH,XFCUM,NLO,NBINS,TFTOT,
  54      +                  IERR)
  55 C        Array XFCUM is filled from NLO to NLO+NBINS, which makes
  56 C        the number of values NBINS+1, or the number of bins NBINS
  57       EXTERNAL FUNC
  58       DIMENSION XFCUM(*)
  59       PARAMETER (RTEPS=0.005, NZ=10, MAXZ=20, NITMAX=6,PRECIS=1.E-6)
  60 CC      DOUBLE PRECISION TPCTIL, TZ, TCUM, XINCR, DTABS,
  61 CC     +  TINCR, TZMAX, XBEST, DTBEST, DTPAR2
  62 C
  63       IERR = 0
  64       IF (TFTOT .LE. 0.) GO TO 900
  65       TPCTIL = TFTOT/NBINS
  66       TZ = TPCTIL/NZ
  67       TZMAX = TZ * 2.
  68       XFCUM(NLO) = XLOW
  69       XFCUM(NLO+NBINS) = XHIGH
  70       X = XLOW
  71       F = FUNC(X)
  72       IF (F .LT. 0.) GO TO 900
  73 C         Loop over percentile bins
  74       DO 600 IBIN = NLO, NLO+NBINS-2
  75       TCUM = 0.
  76       X1 = X
  77       F1 = F
  78       DXMAX = (XHIGH -X) / NZ
  79       FMIN = TZ/DXMAX
  80       FMINZ = FMIN
  81 C         Loop over trapezoids within a supposed percentile
  82       DO 500 IZ= 1, MAXZ
  83       XINCR = TZ/MAX(F1,FMIN,FMINZ)
  84   350 X = X1 + XINCR
  85       F = FUNC(X)
  86       IF (F .LT. 0.) GO TO 900
  87       TINCR = (X-X1) * 0.5 * (F+F1)
  88       IF (TINCR .LT. TZMAX) GO TO 370
  89       XINCR = XINCR * 0.5
  90       GO TO 350
  91   370 CONTINUE
  92       TCUM = TCUM + TINCR
  93       IF (TCUM .GE. TPCTIL*0.99) GO TO 520
  94       FMINZ = TZ*F/ (TPCTIL-TCUM)
  95       F1 = F
  96       X1 = X
  97   500 CONTINUE
  98       WRITE(6,'(A)') ' FUNLUX:  WARNING. FUNPCT fails trapezoid.'
  99 C         END OF TRAPEZOID LOOP
 100 C         Adjust interval using Gaussian integration with
 101 C             Newton corrections since F is the derivative
 102   520 CONTINUE
 103       X1 = XFCUM(IBIN)
 104       XBEST = X
 105       DTBEST = TPCTIL
 106       TPART = TPCTIL
 107 C         Allow for maximum NITMAX more iterations on RADAPT
 108       DO 550 IHOME= 1, NITMAX
 109   535 XINCR = (TPCTIL-TPART) / MAX(F,FMIN)
 110       X = XBEST + XINCR
 111       X2 = X
 112         IF (IHOME .GT. 1 .AND. X2 .EQ. XBEST) THEN
 113         WRITE (6,'(A,G12.3)')
 114      +  ' FUNLUX: WARNING from FUNPCT: insufficient precision at X=',X
 115         GO TO 580
 116         ENDIF
 117       REFX = ABS(X)+PRECIS
 118       CALL RADAPT(FUNC,X1,X2,1,RTEPS,0.,TPART2,UNCERT)
 119       DTPAR2 = TPART2-TPCTIL
 120       DTABS = ABS(DTPAR2)
 121       IF(ABS(XINCR)/REFX .LT. PRECIS) GOTO 545
 122       IF(DTABS .LT. DTBEST) GOTO 545
 123       XINCR = XINCR * 0.5
 124       GOTO 535
 125   545 DTBEST = DTABS
 126       XBEST = X
 127       TPART = TPART2
 128       F = FUNC(X)
 129       IF(F .LT. 0.) GOTO 900
 130       IF(DTABS .LT. RTEPS*TPCTIL) GOTO 580
 131   550 CONTINUE
 132       WRITE (6,'(A,I4)')
 133      +   ' FUNLUX: WARNING from FUNPCT: cannot converge, bin',IBIN
 134 C
 135   580 CONTINUE
 136       XINCR = (TPCTIL-TPART) / MAX(F,FMIN)
 137       X = XBEST + XINCR
 138       XFCUM(IBIN+1) = X
 139       F = FUNC(X)
 140       IF(F .LT. 0.) GOTO 900
 141   600 CONTINUE
 142 C         END OF LOOP OVER BINS
 143       X1 = XFCUM(NLO+NBINS-1)
 144       X2 = XHIGH
 145       CALL RADAPT(FUNC,X1,X2,1,RTEPS,0.,TPART ,UNCERT)
 146       ABERR = ABS(TPART-TPCTIL)/TFTOT
 147 C      WRITE(6,1001) IFUNC,XLOW,XHIGH
 148       IF(ABERR .GT. RTEPS)  WRITE(6,1002) ABERR
 149       RETURN
 150   900 WRITE(6,1000) X,F
 151       IERR = 1
 152       RETURN
 153  1000 FORMAT(/' FUNLUX fatal error in FUNPCT: function negative:'/
 154      + ,' at X=',E15.6,', F=',E15.6/)
 155 C 1001 FORMAT(' FUNPCT has prepared USER FUNCTION',I3,
 156 C     + ' between',G12.3,' and',G12.3,' for FUNLUX.')
 157  1002 FORMAT(' WARNING: Relative error in cumulative distribution',
 158      + ' may be as big as',F10.7)
 159       END
 160       SUBROUTINE FUNLUX(ARRAY,XRAN,LEN)
 161 C         Generation of LEN random numbers in any given distribution,
 162 C         by 4-point interpolation in the inverse cumulative distr.
 163 C         which was previously generated by FUNLXP
 164       COMMON/FUNINT/XUNI
 165       DIMENSION ARRAY(200)
 166       DIMENSION XRAN(LEN)
 167 C  Bin width for main sequence, and its inverse
 168       PARAMETER (GAP= 1./99.,  GAPINV=99.)
 169 C  Top of left tail, bottom of right tail (each tail replaces 2 bins)
 170       PARAMETER (TLEFT= 2./99.,BRIGHT=97./99.)
 171 C  Bin width for minor sequences (tails), and its inverse
 172       PARAMETER (GAPS=TLEFT/49.,  GAPINS=1./GAPS)
 173 C
 174 C   The array ARRAY is assumed to have the following structure:
 175 C        ARRAY(1-100) contains the 99 bins of the inverse cumulative
 176 C                     distribution of the entire function.
 177 C        ARRAY(101-150) contains the 49-bin blowup of main bins
 178 C                       1 and 2 (left tail of distribution)
 179 C        ARRAY(151-200) contains the 49-bin blowup of main bins
 180 C                       98 and 99 (right tail of distribution)
 181 C
 182       CALL RANLUX(XRAN,LEN)
 183
 184       DO 500 IBUF= 1, LEN
 185       X = XRAN(IBUF)
 186       J = INT(  X    *GAPINV) + 1
 187       IF (J .LT. 3)  THEN
 188          J1 = INT( X *GAPINS)
 189              J = J1 + 101
 190              J = MAX(J,102)
 191              J = MIN(J,148)
 192          P = (   X -GAPS*(J1-1)) * GAPINS
 193          A = (P+1.0) * ARRAY(J+2) - (P-2.0)*ARRAY(J-1)
 194          B = (P-1.0) * ARRAY(J) - P * ARRAY(J+1)
 195          XRAN(IBUF) = A*P*(P-1.0)*0.16666667  + B*(P+1.0)*(P-2.0)*0.5
 196       ELSE IF (J .GT. 97)  THEN
 197          J1 = INT((X-BRIGHT)*GAPINS)
 198              J = J1 + 151
 199              J = MAX(J,152)
 200              J = MIN(J,198)
 201          P = (X -BRIGHT -GAPS*(J1-1)) * GAPINS
 202          A = (P+1.0) * ARRAY(J+2) - (P-2.0)*ARRAY(J-1)
 203          B = (P-1.0) * ARRAY(J) - P * ARRAY(J+1)
 204          XRAN(IBUF) = A*P*(P-1.0)*0.16666667  + B*(P+1.0)*(P-2.0)*0.5
 205       ELSE
 206 CC      J = MAX(J,2)
 207 CC      J = MIN(J,98)
 208          P = (   X -GAP*(J-1)) * GAPINV
 209          A = (P+1.0) * ARRAY(J+2) - (P-2.0)*ARRAY(J-1)
 210          B = (P-1.0) * ARRAY(J) - P * ARRAY(J+1)
 211          XRAN(IBUF) = A*P*(P-1.0)*0.16666667  + B*(P+1.0)*(P-2.0)*0.5
 212       ENDIF
 213   500 CONTINUE
 214       XUNI = X
 215       RETURN
 216       END
 217       SUBROUTINE FUNLZ(FUNC,X2LOW,X2HIGH,XLOW,XHIGH)
 218 C         FIND RANGE WHERE FUNC IS NON-ZERO.
 219 C         WRITTEN 1980, F. JAMES
 220 C         MODIFIED, NOV. 1985, TO FIX BUG AND GENERALIZE
 221 C         TO FIND SIMPLY-CONNECTED NON-ZERO REGION (XLOW,XHIGH)
 222 C         ANYWHERE WITHIN THE GIVEN REGION (X2LOW,H2HIGH).
 223 C            WHERE 'ANYWHERE' MEANS EITHER AT THE LOWER OR UPPER
 224 C            EDGE OF THE GIVEN REGION, OR, IF IN THE MIDDLE,
 225 C            COVERING AT LEAST 1% OF THE GIVEN REGION.
 226 C         OTHERWISE IT IS NOT GUARANTEED TO FIND THE NON-ZERO REGION.
 227 C         IF FUNCTION EVERYWHERE ZERO, FUNLZ SETS XLOW=XHIGH=0.
 228       EXTERNAL FUNC
 229       XLOW = X2LOW
 230       XHIGH = X2HIGH
 231 C         FIND OUT IF FUNCTION IS ZERO AT ONE END OR BOTH
 232       XMID = XLOW
 233       IF (FUNC(XLOW) .GT. 0.) GO TO 120
 234       XMID = XHIGH
 235       IF (FUNC(XHIGH) .GT. 0.)  GO TO 50
 236 C         FUNCTION IS ZERO AT BOTH ENDS,
 237 C         LOOK FOR PLACE WHERE IT IS NON-ZERO.
 238       DO 30 LOGN= 1, 7
 239       NSLICE = 2**LOGN
 240       DO 20 I= 1, NSLICE, 2
 241       XMID = XLOW + I * (XHIGH-XLOW) / NSLICE
 242       IF (FUNC(XMID) .GT. 0.)  GO TO 50
 243    20 CONTINUE
 244    30 CONTINUE
 245 C         FALLING THROUGH LOOP MEANS CANNOT FIND NON-ZERO VALUE
 246       WRITE(6,554)
 247       WRITE(6,555) XLOW, XHIGH
 248       XLOW = 0.
 249       XHIGH = 0.
 250       GO TO 220
 251 C
 252    50 CONTINUE
 253 C         DELETE 'LEADING' ZERO RANGE
 254       XH = XMID
 255       XL = XLOW
 256       DO 70 K= 1, 20
 257       XNEW = 0.5*(XH+XL)
 258       IF (FUNC(XNEW) .EQ. 0.) GO TO 68
 259       XH = XNEW
 260       GO TO 70
 261    68 XL = XNEW
 262    70 CONTINUE
 263       XLOW = XL
 264       WRITE(6,555) X2LOW,XLOW
 265   120 CONTINUE
 266       IF (FUNC(XHIGH) .GT. 0.) GO TO 220
 267 C         DELETE 'TRAILING' RANGE OF ZEROES
 268       XL = XMID
 269       XH = XHIGH
 270       DO 170 K= 1, 20
 271       XNEW = 0.5*(XH+XL)
 272       IF (FUNC(XNEW) .EQ. 0.) GO TO 168
 273       XL = XNEW
 274       GO TO 170
 275   168 XH = XNEW
 276   170 CONTINUE
 277       XHIGH = XH
 278       WRITE(6,555) XHIGH, X2HIGH
 279 C
 280   220 CONTINUE
 281       RETURN
 282   554 FORMAT('0CANNOT FIND NON-ZERO FUNCTION VALUE')
 283   555 FORMAT(' FUNCTION IS ZERO FROM X=',E12.5,' TO ',E12.5)
 284       END