dqelg.f

Go to the documentation of this file.

      SUBROUTINE dqelg(N,EPSTAB,RESULT,ABSERR,RES3LA,NRES)
C***BEGIN PROLOGUE  DQELG
C***REFER TO  DQAGIE,DQAGOE,DQAGPE,DQAGSE
C***ROUTINES CALLED  D1MACH
C***REVISION DATE  830518   (YYMMDD)
C***KEYWORDS  EPSILON ALGORITHM, CONVERGENCE ACCELERATION,
C             EXTRAPOLATION
C***AUTHOR  PIESSENS,ROBERT,APPL. MATH. & PROGR. DIV. - K.U.LEUVEN
C           DE DONCKER,ELISE,APPL. MATH & PROGR. DIV. - K.U.LEUVEN
C***PURPOSE  THE ROUTINE DETERMINES THE LIMIT OF A GIVEN SEQUENCE OF
C            APPROXIMATIONS, BY MEANS OF THE EPSILON ALGORITHM OF
C            P.WYNN. AN ESTIMATE OF THE ABSOLUTE ERROR IS ALSO GIVEN.
C            THE CONDENSED EPSILON TABLE IS COMPUTED. ONLY THOSE
C            ELEMENTS NEEDED FOR THE COMPUTATION OF THE NEXT DIAGONAL
C            ARE PRESERVED.
C***DESCRIPTION
C
C           EPSILON ALGORITHM
C           STANDARD FORTRAN SUBROUTINE
C           DOUBLE PRECISION VERSION
C
C           PARAMETERS
C              N      - INTEGER
C                       EPSTAB(N) CONTAINS THE NEW ELEMENT IN THE
C                       FIRST COLUMN OF THE EPSILON TABLE.
C
C              EPSTAB - DOUBLE PRECISION
C                       VECTOR OF DIMENSION 52 CONTAINING THE ELEMENTS
C                       OF THE TWO LOWER DIAGONALS OF THE TRIANGULAR
C                       EPSILON TABLE. THE ELEMENTS ARE NUMBERED
C                       STARTING AT THE RIGHT-HAND CORNER OF THE
C                       TRIANGLE.
C
C              RESULT - DOUBLE PRECISION
C                       RESULTING APPROXIMATION TO THE INTEGRAL
C
C              ABSERR - DOUBLE PRECISION
C                       ESTIMATE OF THE ABSOLUTE ERROR COMPUTED FROM
C                       RESULT AND THE 3 PREVIOUS RESULTS
C
C              RES3LA - DOUBLE PRECISION
C                       VECTOR OF DIMENSION 3 CONTAINING THE LAST 3
C                       RESULTS
C
C              NRES   - INTEGER
C                       NUMBER OF CALLS TO THE ROUTINE
C                       (SHOULD BE ZERO AT FIRST CALL)
C
C***END PROLOGUE  DQELG
C
      DOUBLE PRECISION ABSERR,DABS,DELTA1,DELTA2,DELTA3,DMAX1,D1MACH,
     *  EPMACH,EPSINF,EPSTAB,ERROR,ERR1,ERR2,ERR3,E0,E1,E1ABS,E2,E3,
     *  OFLOW,RES,RESULT,RES3LA,SS,TOL1,TOL2,TOL3
      INTEGER I,IB,IB2,IE,INDX,K1,K2,K3,LIMEXP,N,NEWELM,NRES,NUM
      dimension epstab(52),res3la(3)
C
C           LIST OF MAJOR VARIABLES
C           -----------------------
C
C           E0     - THE 4 ELEMENTS ON WHICH THE COMPUTATION OF A NEW
C           E1       ELEMENT IN THE EPSILON TABLE IS BASED
C           E2
C           E3                 E0
C                        E3    E1    NEW
C                              E2
C           NEWELM - NUMBER OF ELEMENTS TO BE COMPUTED IN THE NEW
C                    DIAGONAL
C           ERROR  - ERROR = ABS(E1-E0)+ABS(E2-E1)+ABS(NEW-E2)
C           RESULT - THE ELEMENT IN THE NEW DIAGONAL WITH LEAST VALUE
C                    OF ERROR
C
C           MACHINE DEPENDENT CONSTANTS
C           ---------------------------
C
C           EPMACH IS THE LARGEST RELATIVE SPACING.
C           OFLOW IS THE LARGEST POSITIVE MAGNITUDE.
C           LIMEXP IS THE MAXIMUM NUMBER OF ELEMENTS THE EPSILON
C           TABLE CAN CONTAIN. IF THIS NUMBER IS REACHED, THE UPPER
C           DIAGONAL OF THE EPSILON TABLE IS DELETED.
C
C***FIRST EXECUTABLE STATEMENT  DQELG
      epmach = d1mach(4)
      oflow = d1mach(2)
      nres = nres+1
      abserr = oflow
      result = epstab(n)
      IF(n.LT.3) GO TO 100
      limexp = 50
      epstab(n+2) = epstab(n)
      newelm = (n-1)/2
      epstab(n) = oflow
      num = n
      k1 = n
      DO 40 i = 1,newelm
        k2 = k1-1
        k3 = k1-2
        res = epstab(k1+2)
        e0 = epstab(k3)
        e1 = epstab(k2)
        e2 = res
        e1abs = dabs(e1)
        delta2 = e2-e1
        err2 = dabs(delta2)
        tol2 = dmax1(dabs(e2),e1abs)*epmach
        delta3 = e1-e0
        err3 = dabs(delta3)
        tol3 = dmax1(e1abs,dabs(e0))*epmach
        IF(err2.GT.tol2.OR.err3.GT.tol3) GO TO 10
C
C           IF E0, E1 AND E2 ARE EQUAL TO WITHIN MACHINE
C           ACCURACY, CONVERGENCE IS ASSUMED.
C           RESULT = E2
C           ABSERR = ABS(E1-E0)+ABS(E2-E1)
C
        result = res
        abserr = err2+err3
C ***JUMP OUT OF DO-LOOP
        GO TO 100
   10   e3 = epstab(k1)
        epstab(k1) = e1
        delta1 = e1-e3
        err1 = dabs(delta1)
        tol1 = dmax1(e1abs,dabs(e3))*epmach
C
C           IF TWO ELEMENTS ARE VERY CLOSE TO EACH OTHER, OMIT
C           A PART OF THE TABLE BY ADJUSTING THE VALUE OF N
C
        IF(err1.LE.tol1.OR.err2.LE.tol2.OR.err3.LE.tol3) GO TO 20
        ss = 0.1d+01/delta1+0.1d+01/delta2-0.1d+01/delta3
        epsinf = dabs(ss*e1)
C
C           TEST TO DETECT IRREGULAR BEHAVIOUR IN THE TABLE, AND
C           EVENTUALLY OMIT A PART OF THE TABLE ADJUSTING THE VALUE
C           OF N.
C
        IF(epsinf.GT.0.1d-03) GO TO 30
   20   n = i+i-1
C ***JUMP OUT OF DO-LOOP
        GO TO 50
C
C           COMPUTE A NEW ELEMENT AND EVENTUALLY ADJUST
C           THE VALUE OF RESULT.
C
   30   res = e1+0.1d+01/ss
        epstab(k1) = res
        k1 = k1-2
        error = err2+dabs(res-e2)+err3
        IF(error.GT.abserr) GO TO 40
        abserr = error
        result = res
   40 CONTINUE
C
C           SHIFT THE TABLE.
C
   50 IF(n.EQ.limexp) n = 2*(limexp/2)-1
      ib = 1
      IF((num/2)*2.EQ.num) ib = 2
      ie = newelm+1
      DO 60 i=1,ie
        ib2 = ib+2
        epstab(ib) = epstab(ib2)
        ib = ib2
   60 CONTINUE
      IF(num.EQ.n) GO TO 80
      indx = num-n+1
      DO 70 i = 1,n
        epstab(i)= epstab(indx)
        indx = indx+1
   70 CONTINUE
   80 IF(nres.GE.4) GO TO 90
      res3la(nres) = result
      abserr = oflow
      GO TO 100
C
C           COMPUTE ERROR ESTIMATE
C
   90 abserr = dabs(result-res3la(3))+dabs(result-res3la(2))
     *  +dabs(result-res3la(1))
      res3la(1) = res3la(2)
      res3la(2) = res3la(3)
      res3la(3) = result
  100 abserr = dmax1(abserr,0.5d+01*epmach*dabs(result))
      RETURN
      END