MINICERN/mathlib/gen/d/dadapt.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1996/04/01 15:02:13  mclareni
   6 * Mathlib gen
   7 *
   8 *
   9 #include "gen/pilot.h"
  10       SUBROUTINE DADAPT(F,A,B,NSEG,RELTOL,ABSTOL,RES,ERR)
  11 #if !defined(CERNLIB_DOUBLE)
  12 #include "gen/imp128.inc"
  13       CHARACTER*6 NAME
  14       NAME = 'DADAPT'
  15       CALL MTLPRT(NAME,'D102',
  16      +'not available on this machine - see documentation')
  17       RETURN
  18       END
  19 #endif
  20 #if defined(CERNLIB_DOUBLE)
  21       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
  22
  23 C     RES = Estimated Integral of F from A to B,
  24 C     ERR = Estimated absolute error on RES.
  25 C     NSEG  specifies how the adaptation is to be done:
  26 C        =0   means use previous binning,
  27 C        =1   means fully automatic, adapt until tolerance attained.
  28 C        =n>1 means first split interval into n equal segments,
  29 C             then adapt as necessary to attain tolerance.
  30 C     The specified tolerances are:
  31 C            relative: RELTOL ;  absolute: ABSTOL.
  32 C        It stops when one OR the other is satisfied, or number of
  33 C        segments exceeds NDIM.  Either TOLA or TOLR (but not both!)
  34 C        can be set to zero, in which case only the other is used.
  35
  36       EXTERNAL F
  37
  38       PARAMETER (NDIM=100)
  39       PARAMETER (R1 = 1, HF = R1/2)
  40
  41       DIMENSION XLO(NDIM),XHI(NDIM),TVAL(NDIM),TERS(NDIM)
  42       SAVE XLO,XHI,TVAL,TERS,NTER
  43       DATA NTER /0/
  44
  45       IF(NSEG .LE. 0)  THEN
  46        IF(NTER .EQ. 0) THEN
  47         NSEGD=1
  48         GO TO 2
  49        ENDIF
  50        TVALS=0
  51        TERSS=0
  52        DO 1 I = 1,NTER
  53        CALL DGS56P(F,XLO(I),XHI(I),TVAL(I),TE)
  54        TERS(I)=TE**2
  55        TVALS=TVALS+TVAL(I)
  56        TERSS=TERSS+TERS(I)
  57     1  CONTINUE
  58        ROOT= SQRT(2*TERSS)
  59        GO TO 9
  60       ENDIF
  61       NSEGD=MIN(NSEG,NDIM)
  62     2 XHIB=A
  63       BIN=(B-A)/NSEGD
  64       DO 3 I = 1,NSEGD
  65       XLO(I)=XHIB
  66       XLOB=XLO(I)
  67       XHI(I)=XHIB+BIN
  68       IF(I .EQ. NSEGD) XHI(I)=B
  69       XHIB=XHI(I)
  70       CALL DGS56P(F,XLOB,XHIB,TVAL(I),TE)
  71       TERS(I)=TE**2
  72     3 CONTINUE
  73       NTER=NSEGD
  74       DO 4 ITER = 1,NDIM
  75       TVALS=TVAL(1)
  76       TERSS=TERS(1)
  77       DO 5 I = 2,NTER
  78       TVALS=TVALS+TVAL(I)
  79       TERSS=TERSS+TERS(I)
  80     5 CONTINUE
  81       ROOT= SQRT(2*TERSS)
  82       IF(ROOT .LE. ABSTOL .OR. ROOT .LE. RELTOL*ABS(TVALS)) GO TO 9
  83       IF(NTER .EQ. NDIM) GO TO 9
  84       BIGE=TERS(1)
  85       IBIG=1
  86       DO 6 I = 2,NTER
  87       IF(TERS(I) .GT. BIGE) THEN
  88        BIGE=TERS(I)
  89        IBIG=I
  90       ENDIF
  91     6 CONTINUE
  92       NTER=NTER+1
  93       XHI(NTER)=XHI(IBIG)
  94       XNEW=HF*(XLO(IBIG)+XHI(IBIG))
  95       XHI(IBIG)=XNEW
  96       XLO(NTER)=XNEW
  97       CALL DGS56P(F,XLO(IBIG),XHI(IBIG),TVAL(IBIG),TE)
  98       TERS(IBIG)=TE**2
  99       CALL DGS56P(F,XLO(NTER),XHI(NTER),TVAL(NTER),TE)
 100       TERS(NTER)=TE**2
 101     4 CONTINUE
 102     9 RES=TVALS
 103       ERR=ROOT
 104       RETURN
 105       END
 106 #endif