d5/dde/gadapf_8F_source.html

C **********************************************************************


      FUNCTION gadapf(X,A0,B0,F,EPS)

      IMPLICIT NONE

      common/gadap1/ num,ifu

      INTEGER num,ifu

      SAVE /gadap1/


      INTEGER n,l,i

      REAL gadapf,x,a0,b0,f,eps,a,b,f1,f2,f3,s,red,c,w1,u2,ss,sold,sum

      REAL dsum

      REAL aa,bb,f1f,f2f,f3f


      EXTERNAL f

      dimension a(300),b(300),f1(300),f2(300),f3(300),s(300),n(300)

    1 FORMAT(16h gadap:i too big)

      dsum(f1f,f2f,f3f,aa,bb)=5./18.*(bb-aa)*(f1f+1.6*f2f+f3f)

      IF(eps.LT.1.0e-8) eps=1.0e-8

      red=1.4

      l=1

      i=1

      sum=0.

      c=sqrt(15.)/5.

      a(1)=a0

      b(1)=b0

      f1(1)=f(x,0.5*(1+c)*a0+0.5*(1-c)*b0)

      f2(1)=f(x,0.5*(a0+b0))

      f3(1)=f(x,0.5*(1-c)*a0+0.5*(1+c)*b0)

      ifu=3

      s(1)=  dsum(f1(1),f2(1),f3(1),a0,b0)

  100 CONTINUE

      l=l+1

      n(l)=3

      eps=eps*red

      a(i+1)=a(i)+c*(b(i)-a(i))

      b(i+1)=b(i)

      a(i+2)=a(i)+b(i)-a(i+1)

      b(i+2)=a(i+1)

      a(i+3)=a(i)

      b(i+3)=a(i+2)

      w1=a(i)+(b(i)-a(i))/5.

      u2=2.*w1-(a(i)+a(i+2))/2.

      f1(i+1)=f(x,a(i)+b(i)-w1)

      f2(i+1)=f3(i)

      f3(i+1)=f(x,b(i)-a(i+2)+w1)

      f1(i+2)=f(x,u2)

      f2(i+2)=f2(i)

      f3(i+2)=f(x,b(i+2)+a(i+2)-u2)

      f1(i+3)=f(x,a(i)+a(i+2)-w1)

      f2(i+3)=f1(i)

      f3(i+3)=f(x,w1)

      ifu=ifu+6

      IF(ifu.GT.5000) goto 130

      s(i+1)=  dsum(f1(i+1),f2(i+1),f3(i+1),a(i+1),b(i+1))

      s(i+2)=  dsum(f1(i+2),f2(i+2),f3(i+2),a(i+2),b(i+2))

      s(i+3)=  dsum(f1(i+3),f2(i+3),f3(i+3),a(i+3),b(i+3))

      ss=s(i+1)+s(i+2)+s(i+3)

      i=i+3

      IF(i.GT.300)goto 120

      sold=s(i-3)

      IF(abs(sold-ss).GT.eps*(1.+abs(ss))/2.) goto 100

      sum=sum+ss

      i=i-4

      n(l)=0

      l=l-1

  110 CONTINUE

      IF(l.EQ.1) goto 130

      n(l)=n(l)-1

      eps=eps/red

      IF(n(l).NE.0) goto 100

      i=i-1

      l=l-1

      goto 110

  120 WRITE(6,1)

  130 gadapf=sum

      eps=eps/red

      RETURN

      END