This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / 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