[u/mrichter/AliRoot.git] / MINICERN / mathlib / gen / d / dadapt.F

*
* $Id$
*
* $Log$
* Revision 1.1.1.1  1996/04/01 15:02:13  mclareni
* Mathlib gen
*
*
#include "gen/pilot.h"
      SUBROUTINE DADAPT(F,A,B,NSEG,RELTOL,ABSTOL,RES,ERR)
#if !defined(CERNLIB_DOUBLE)
#include "gen/imp128.inc"
      CHARACTER*6 NAME
      NAME = 'DADAPT'
      CALL MTLPRT(NAME,'D102',
     +'not available on this machine - see documentation')
      RETURN
      END
#endif
#if defined(CERNLIB_DOUBLE)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
 
C     RES = Estimated Integral of F from A to B,
C     ERR = Estimated absolute error on RES.
C     NSEG  specifies how the adaptation is to be done:
C        =0   means use previous binning,
C        =1   means fully automatic, adapt until tolerance attained.
C        =n>1 means first split interval into n equal segments,
C             then adapt as necessary to attain tolerance.
C     The specified tolerances are:
C            relative: RELTOL ;  absolute: ABSTOL.
C        It stops when one OR the other is satisfied, or number of
C        segments exceeds NDIM.  Either TOLA or TOLR (but not both!)
C        can be set to zero, in which case only the other is used.
 
      EXTERNAL F
 
      PARAMETER (NDIM=100)
      PARAMETER (R1 = 1, HF = R1/2)
 
      DIMENSION XLO(NDIM),XHI(NDIM),TVAL(NDIM),TERS(NDIM)
      SAVE XLO,XHI,TVAL,TERS,NTER
      DATA NTER /0/
 
      IF(NSEG .LE. 0)  THEN
       IF(NTER .EQ. 0) THEN
        NSEGD=1
        GO TO 2
       ENDIF
       TVALS=0
       TERSS=0
       DO 1 I = 1,NTER
       CALL DGS56P(F,XLO(I),XHI(I),TVAL(I),TE)
       TERS(I)=TE**2
       TVALS=TVALS+TVAL(I)
       TERSS=TERSS+TERS(I)
    1  CONTINUE
       ROOT= SQRT(2*TERSS)
       GO TO 9
      ENDIF
      NSEGD=MIN(NSEG,NDIM)
    2 XHIB=A
      BIN=(B-A)/NSEGD
      DO 3 I = 1,NSEGD
      XLO(I)=XHIB
      XLOB=XLO(I)
      XHI(I)=XHIB+BIN
      IF(I .EQ. NSEGD) XHI(I)=B
      XHIB=XHI(I)
      CALL DGS56P(F,XLOB,XHIB,TVAL(I),TE)
      TERS(I)=TE**2
    3 CONTINUE
      NTER=NSEGD
      DO 4 ITER = 1,NDIM
      TVALS=TVAL(1)
      TERSS=TERS(1)
      DO 5 I = 2,NTER
      TVALS=TVALS+TVAL(I)
      TERSS=TERSS+TERS(I)
    5 CONTINUE
      ROOT= SQRT(2*TERSS)
      IF(ROOT .LE. ABSTOL .OR. ROOT .LE. RELTOL*ABS(TVALS)) GO TO 9
      IF(NTER .EQ. NDIM) GO TO 9
      BIGE=TERS(1)
      IBIG=1
      DO 6 I = 2,NTER
      IF(TERS(I) .GT. BIGE) THEN
       BIGE=TERS(I)
       IBIG=I
      ENDIF
    6 CONTINUE
      NTER=NTER+1
      XHI(NTER)=XHI(IBIG)
      XNEW=HF*(XLO(IBIG)+XHI(IBIG))
      XHI(IBIG)=XNEW
      XLO(NTER)=XNEW
      CALL DGS56P(F,XLO(IBIG),XHI(IBIG),TVAL(IBIG),TE)
      TERS(IBIG)=TE**2
      CALL DGS56P(F,XLO(NTER),XHI(NTER),TVAL(NTER),TE)
      TERS(NTER)=TE**2
    4 CONTINUE
    9 RES=TVALS
      ERR=ROOT
      RETURN
      END
#endif
Commit	Line	Data
fe4da5cc	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