C      ALGORITHM 691, COLLECTED ALGORITHMS FROM ACM.
C      THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
C      VOL. 17, NO. 2, PP. 218-232.  JUNE, 1991.
This disk contains 8 files:
 
README : this information file
 
stest.for : which contains simple test routines (single precision)
 
          TEST    (main)
          F
 
squad.for : which contains quadrature routines (single precision)
 
          QXGS
          QXGSE
          QXG
          QXGE
          QXLQM
          QXRUL
          QXRRD
          QXCPY
 
saux.for : which contains auxiliary routines (single precision)
 
          QPSRT
          QELG
          R1MACH (machine dependent, TO BE MODIFIED)
 
dtest.for : which contains simple test routines (double precision)
 
          DTEST    (main)
          F
 
dquad.for : which contains quadrature routines (double precision)
 
          DQXGS
          DQXGSE
          DQXG
          DQXGE
          DQXLQM
          DQXRUL
          DQXRRD
          DQXCPY
 
daux.for : which contains auxiliary routines (double precision)
 
          DQPSRT
          DQELG
          D1MACH  (machine dependent, TO BE MODIFIED)
 
CC+------------------------------------------------------------------+CC
CC+                                                                  +CC
CC+                LAWRENCE LIVERMORE NATIONAL LABORATORY            +CC
CC+                                                                  +CC
CC+        MATHEMATICS AND STATISTICS DIVISION SOFTWARE LIBRARY      +CC
CC+                                                                  +CC
CC+------------------------------------------------------------------+CC
CC+                                                                  +CC
CC+  THE SLATEC MATHEMATICAL PROGRAM LIBRARY IS ISSUED BY            +CC
CC+  THE FOLLOWING:                                                  +CC
CC+                                                                  +CC
CC+       AIR FORCE WEAPONS LABORATORY, ALBUQUERQUE                  +CC
CC+       LAWRENCE LIVERMORE NATIONAL LABORATORY, LIVERMORE          +CC
CC+       LOS ALAMOS NATIONAL SCIENTIFIC LABORATORY, LOS ALAMOS      +CC
CC+       NATIONAL BUREAU OF STANDARDS, WASHINGTON                   +CC
CC+       OAK RIDGE NATIONAL LABORATORY, OAK RIDGE                   +CC
CC+       SANDIA NATIONAL LABORATORIES, ALBUQUERQUE                  +CC
CC+       SANDIA NATIONAL LABORATORIES, LIVERMORE                    +CC
CC+                                                                  +CC
CC+  PLEASE REFER QUESTIONS REGARDING THIS SOFTWARE TO THE           +CC
CC+  MSD CONSULTING OFFICE, EXTENSION 3-2976.                        +CC
CC+                                                                  +CC
CC+                                                                  +CC
CC+              * * * * * N O T I C E * * * * *                     +CC
CC+                                                                  +CC
CC+  THIS MATERIAL WAS PREPARED AS AN ACCOUNT OF WORK SPONSORED      +CC
CC+  BY THE UNITED STATES GOVERNMENT.  NEITHER THE UNITED            +CC
CC+  STATES GOVERNMENT, NOR THE DEPARTMENT OF DEFENSE, NOR           +CC
CC+  ANY OF THEIR EMPLOYEES, MAKES ANY WARRANTY, EXPRESSED           +CC
CC+  OR IMPLIED, OR ASSUMES ANY LEGAL LIABILITY OR                   +CC
CC+  RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR               +CC
CC+  USEFULNESS OF ANY INFORMATION, APPARATUS, PRODUCT,              +CC
CC+  OR PROCESS DISCLOSED, OR REPRESENTS ITS USE WOULD NOT           +CC
CC+  INFRINGE UPON PRIVATELY OWNED RIGHTS.                           +CC
CC+                                                                  +CC
CC+------------------------------------------------------------------+CC
*DECK DQPSRT
      SUBROUTINE DQPSRT(LIMIT,LAST,MAXERR,ERMAX,ELIST,IORD,NRMAX)
C***BEGIN PROLOGUE  DQPSRT
C***REFER TO  DQAGE,DQAGIE,DQAGPE,DQAWSE
C***ROUTINES CALLED  (NONE)
C***REVISION DATE  810101   (YYMMDD)
C***KEYWORDS  SEQUENTIAL SORTING
C***AUTHOR  PIESSENS, ROBERT, APPLIED MATH. AND PROGR. DIV. -
C             K. U. LEUVEN
C           DE DONCKER, ELISE, APPLIED MATH. AND PROGR. DIV. -
C             K. U. LEUVEN
C***PURPOSE  THIS ROUTINE MAINTAINS THE DESCENDING ORDERING IN THE
C            LIST OF THE LOCAL ERROR ESTIMATED RESULTING FROM THE
C            INTERVAL SUBDIVISION PROCESS. AT EACH CALL TWO ERROR
C            ESTIMATES ARE INSERTED USING THE SEQUENTIAL SEARCH
C            METHOD, TOP-DOWN FOR THE LARGEST ERROR ESTIMATE AND
C            BOTTOM-UP FOR THE SMALLEST ERROR ESTIMATE.
C***DESCRIPTION
C
C           ORDERING ROUTINE
C           STANDARD FORTRAN SUBROUTINE
C           DOUBLE PRECISION VERSION
C
C           PARAMETERS (MEANING AT OUTPUT)
C              LIMIT  - INTEGER
C                       MAXIMUM NUMBER OF ERROR ESTIMATES THE LIST
C                       CAN CONTAIN
C
C              LAST   - INTEGER
C                       NUMBER OF ERROR ESTIMATES CURRENTLY IN THE LIST
C
C              MAXERR - INTEGER
C                       MAXERR POINTS TO THE NRMAX-TH LARGEST ERROR
C                       ESTIMATE CURRENTLY IN THE LIST
C
C              ERMAX  - DOUBLE PRECISION
C                       NRMAX-TH LARGEST ERROR ESTIMATE
C                       ERMAX = ELIST(MAXERR)
C
C              ELIST  - DOUBLE PRECISION
C                       VECTOR OF DIMENSION LAST CONTAINING
C                       THE ERROR ESTIMATES
C
C              IORD   - INTEGER
C                       VECTOR OF DIMENSION LAST, THE FIRST K ELEMENTS
C                       OF WHICH CONTAIN POINTERS TO THE ERROR
C                       ESTIMATES, SUCH THAT
C                       ELIST(IORD(1)),...,  ELIST(IORD(K))
C                       FORM A DECREASING SEQUENCE, WITH
C                       K = LAST IF LAST.LE.(LIMIT/2+2), AND
C                       K = LIMIT+1-LAST OTHERWISE
C
C              NRMAX  - INTEGER
C                       MAXERR = IORD(NRMAX)
C***END PROLOGUE  DQPSRT
C
      DOUBLE PRECISION ELIST,ERMAX,ERRMAX,ERRMIN
      INTEGER I,IBEG,IDO,IORD,ISUCC,J,JBND,JUPBN,K,LAST,LIMIT,MAXERR,
     1  NRMAX
      DIMENSION ELIST(LAST),IORD(LAST)
C
C           CHECK WHETHER THE LIST CONTAINS MORE THAN
C           TWO ERROR ESTIMATES.
C
C***FIRST EXECUTABLE STATEMENT  DQPSRT
      IF(LAST.GT.2) GO TO 10
      IORD(1) = 1
      IORD(2) = 2
      GO TO 90
C
C           THIS PART OF THE ROUTINE IS ONLY EXECUTED IF, DUE TO A
C           DIFFICULT INTEGRAND, SUBDIVISION INCREASED THE ERROR
C           ESTIMATE. IN THE NORMAL CASE THE INSERT PROCEDURE SHOULD
C           START AFTER THE NRMAX-TH LARGEST ERROR ESTIMATE.
C
   10 ERRMAX = ELIST(MAXERR)
      IF(NRMAX.EQ.1) GO TO 30
      IDO = NRMAX-1
      DO 20 I = 1,IDO
        ISUCC = IORD(NRMAX-1)
C ***JUMP OUT OF DO-LOOP
        IF(ERRMAX.LE.ELIST(ISUCC)) GO TO 30
        IORD(NRMAX) = ISUCC
        NRMAX = NRMAX-1
   20    CONTINUE
C
C           COMPUTE THE NUMBER OF ELEMENTS IN THE LIST TO BE MAINTAINED
C           IN DESCENDING ORDER. THIS NUMBER DEPENDS ON THE NUMBER OF
C           SUBDIVISIONS STILL ALLOWED.
C
   30 JUPBN = LAST
      IF(LAST.GT.(LIMIT/2+2)) JUPBN = LIMIT+3-LAST
      ERRMIN = ELIST(LAST)
C
C           INSERT ERRMAX BY TRAVERSING THE LIST TOP-DOWN,
C           STARTING COMPARISON FROM THE ELEMENT ELIST(IORD(NRMAX+1)).
C
      JBND = JUPBN-1
      IBEG = NRMAX+1
      IF(IBEG.GT.JBND) GO TO 50
      DO 40 I=IBEG,JBND
        ISUCC = IORD(I)
C ***JUMP OUT OF DO-LOOP
        IF(ERRMAX.GE.ELIST(ISUCC)) GO TO 60
        IORD(I-1) = ISUCC
   40 CONTINUE
   50 IORD(JBND) = MAXERR
      IORD(JUPBN) = LAST
      GO TO 90
C
C           INSERT ERRMIN BY TRAVERSING THE LIST BOTTOM-UP.
C
   60 IORD(I-1) = MAXERR
      K = JBND
      DO 70 J=I,JBND
        ISUCC = IORD(K)
C ***JUMP OUT OF DO-LOOP
        IF(ERRMIN.LT.ELIST(ISUCC)) GO TO 80
        IORD(K+1) = ISUCC
        K = K-1
   70 CONTINUE
      IORD(I) = LAST
      GO TO 90
   80 IORD(K+1) = LAST
C
C           SET MAXERR AND ERMAX.
C
   90 MAXERR = IORD(NRMAX)
      ERMAX = ELIST(MAXERR)
      RETURN
      END
*DECK DQELG
      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  CONVERGENCE ACCELERATION,EPSILON ALGORITHM,EXTRAPOLATION
C***AUTHOR  PIESSENS, ROBERT, APPLIED MATH. AND PROGR. DIV. -
C             K. U. LEUVEN
C           DE DONCKER, ELISE, APPLIED MATH. AND PROGR. DIV. -
C             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***END PROLOGUE  DQELG
C
      DOUBLE PRECISION ABSERR,DABS,DELTA1,DELTA2,DELTA3,DMAX1,D1MACH,
     1  EPMACH,EPSINF,EPSTAB,ERROR,ERR1,ERR2,ERR3,E0,E1,E1ABS,E2,E3,
     2  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))
     1  +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
*DECK D1MACH
      DOUBLE PRECISION FUNCTION D1MACH(I)
C***BEGIN PROLOGUE  D1MACH
C***DATE WRITTEN   750101   (YYMMDD)
C***REVISION DATE  890213   (YYMMDD)
C***CATEGORY NO.  R1
C***KEYWORDS  LIBRARY=SLATEC,TYPE=DOUBLE PRECISION(R1MACH-S D1MACH-D),
C             MACHINE CONSTANTS
C***AUTHOR  FOX, P. A., (BELL LABS)
C           HALL, A. D., (BELL LABS)
C           SCHRYER, N. L., (BELL LABS)
C***PURPOSE  RETURNS DOUBLE PRECISION MACHINE DEPENDENT CONSTANTS
C***DESCRIPTION
C
C   D1MACH CAN BE USED TO OBTAIN MACHINE-DEPENDENT PARAMETERS
C   FOR THE LOCAL MACHINE ENVIRONMENT.  IT IS A FUNCTION
C   SUBPROGRAM WITH ONE (INPUT) ARGUMENT, AND CAN BE CALLED
C   AS FOLLOWS, FOR EXAMPLE
C
C        D = D1MACH(I)
C
C   WHERE I=1,...,5.  THE (OUTPUT) VALUE OF D ABOVE IS
C   DETERMINED BY THE (INPUT) VALUE OF I.  THE RESULTS FOR
C   VARIOUS VALUES OF I ARE DISCUSSED BELOW.
C
C   D1MACH( 1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE.
C   D1MACH( 2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE.
C   D1MACH( 3) = B**(-T), THE SMALLEST RELATIVE SPACING.
C   D1MACH( 4) = B**(1-T), THE LARGEST RELATIVE SPACING.
C   D1MACH( 5) = LOG10(B)
C
C   ASSUME DOUBLE PRECISION NUMBERS ARE REPRESENTED IN THE T-DIGIT,
C   BASE-B FORM
C
C              SIGN (B**E)*( (X(1)/B) + ... + (X(T)/B**T) )
C
C   WHERE 0 .LE. X(I) .LT. B FOR I=1,...,T, 0 .LT. X(1), AND
C   EMIN .LE. E .LE. EMAX.
C
C   THE VALUES OF B, T, EMIN AND EMAX ARE PROVIDED IN I1MACH AS
C   FOLLOWS:
C   I1MACH(10) = B, THE BASE.
C   I1MACH(14) = T, THE NUMBER OF BASE-B DIGITS.
C   I1MACH(15) = EMIN, THE SMALLEST EXPONENT E.
C   I1MACH(16) = EMAX, THE LARGEST EXPONENT E.
C
C   TO ALTER THIS FUNCTION FOR A PARTICULAR ENVIRONMENT,
C   THE DESIRED SET OF DATA STATEMENTS SHOULD BE ACTIVATED BY
C   REMOVING THE C FROM COLUMN 1.  ALSO, THE VALUES OF
C   D1MACH(1) - D1MACH(4) SHOULD BE CHECKED FOR CONSISTENCY
C   WITH THE LOCAL OPERATING SYSTEM.
C
C***REFERENCES  FOX P.A., HALL A.D., SCHRYER N.L.,*FRAMEWORK FOR A
C                 PORTABLE LIBRARY*, ACM TRANSACTIONS ON MATHEMATICAL
C                 SOFTWARE, VOL. 4, NO. 2, JUNE 1978, PP. 177-188.
C***ROUTINES CALLED   XERROR
C
C THE CALL TO XERROR HAS BEEN COMMENTED TO GET A STAND ALONE ROUTINE
C IT CAN BE ADAPTED TO FIT THE LOCAL ERROR HANDLING PROCEDURES
C        NOTE ADDED BY F. ROMANI  7/11/89
C
C***END PROLOGUE  D1MACH
C
      INTEGER SMALL(4)
      INTEGER LARGE(4)
      INTEGER RIGHT(4)
      INTEGER DIVER(4)
      INTEGER LOG10(4)
C
      DOUBLE PRECISION DMACH(5)
      SAVE DMACH
C
      EQUIVALENCE (DMACH(1),SMALL(1))
      EQUIVALENCE (DMACH(2),LARGE(1))
      EQUIVALENCE (DMACH(3),RIGHT(1))
      EQUIVALENCE (DMACH(4),DIVER(1))
      EQUIVALENCE (DMACH(5),LOG10(1))
C
C     MACHINE CONSTANTS FOR THE AMIGA
C     ABSOFT FORTRAN COMPILER USING THE 68020/68881 COMPILER OPTION
C
C     DATA SMALL(1), SMALL(2) / Z'00100000', Z'00000000' /
C     DATA LARGE(1), LARGE(2) / Z'7FEFFFFF', Z'FFFFFFFF' /
C     DATA RIGHT(1), RIGHT(2) / Z'3CA00000', Z'00000000' /
C     DATA DIVER(1), DIVER(2) / Z'3CB00000', Z'00000000' /
C     DATA LOG10(1), LOG10(2) / Z'3FD34413', Z'509F79FF' /
C
C     MACHINE CONSTANTS FOR THE AMIGA
C     ABSOFT FORTRAN COMPILER USING SOFTWARE FLOATING POINT
C
C     DATA SMALL(1), SMALL(2) / Z'00100000', Z'00000000' /
C     DATA LARGE(1), LARGE(2) / Z'7FDFFFFF', Z'FFFFFFFF' /
C     DATA RIGHT(1), RIGHT(2) / Z'3CA00000', Z'00000000' /
C     DATA DIVER(1), DIVER(2) / Z'3CB00000', Z'00000000' /
C     DATA LOG10(1), LOG10(2) / Z'3FD34413', Z'509F79FF' /
C
C     MACHINE CONSTANTS FOR THE APOLLO
C
C     DATA SMALL(1), SMALL(2) / 16#00100000, 16#00000000 /
C     DATA LARGE(1), LARGE(2) / 16#7FFFFFFF, 16#FFFFFFFF /
C     DATA RIGHT(1), RIGHT(2) / 16#3CA00000, 16#00000000 /
C     DATA DIVER(1), DIVER(2) / 16#3CB00000, 16#00000000 /
C     DATA LOG10(1), LOG10(2) / 16#3FD34413, 16#509F79FF /
C
C     MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM
C
C     DATA SMALL(1) / ZC00800000 /
C     DATA SMALL(2) / Z000000000 /
C     DATA LARGE(1) / ZDFFFFFFFF /
C     DATA LARGE(2) / ZFFFFFFFFF /
C     DATA RIGHT(1) / ZCC5800000 /
C     DATA RIGHT(2) / Z000000000 /
C     DATA DIVER(1) / ZCC6800000 /
C     DATA DIVER(2) / Z000000000 /
C     DATA LOG10(1) / ZD00E730E7 /
C     DATA LOG10(2) / ZC77800DC0 /
C
C     MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM
C
C     DATA SMALL(1) / O1771000000000000 /
C     DATA SMALL(2) / O0000000000000000 /
C     DATA LARGE(1) / O0777777777777777 /
C     DATA LARGE(2) / O0007777777777777 /
C     DATA RIGHT(1) / O1461000000000000 /
C     DATA RIGHT(2) / O0000000000000000 /
C     DATA DIVER(1) / O1451000000000000 /
C     DATA DIVER(2) / O0000000000000000 /
C     DATA LOG10(1) / O1157163034761674 /
C     DATA LOG10(2) / O0006677466732724 /
C
C     MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS
C
C     DATA SMALL(1) / O1771000000000000 /
C     DATA SMALL(2) / O7770000000000000 /
C     DATA LARGE(1) / O0777777777777777 /
C     DATA LARGE(2) / O7777777777777777 /
C     DATA RIGHT(1) / O1461000000000000 /
C     DATA RIGHT(2) / O0000000000000000 /
C     DATA DIVER(1) / O1451000000000000 /
C     DATA DIVER(2) / O0000000000000000 /
C     DATA LOG10(1) / O1157163034761674 /
C     DATA LOG10(2) / O0006677466732724 /
C
C     MACHINE CONSTANTS FOR THE CDC 170/180 SERIES USING NOS/VE
C
C     DATA SMALL(1) / Z"3001800000000000" /
C     DATA SMALL(2) / Z"3001000000000000" /
C     DATA LARGE(1) / Z"4FFEFFFFFFFFFFFE" /
C     DATA LARGE(2) / Z"4FFE000000000000" /
C     DATA RIGHT(1) / Z"3FD2800000000000" /
C     DATA RIGHT(2) / Z"3FD2000000000000" /
C     DATA DIVER(1) / Z"3FD3800000000000" /
C     DATA DIVER(2) / Z"3FD3000000000000" /
C     DATA LOG10(1) / Z"3FFF9A209A84FBCF" /
C     DATA LOG10(2) / Z"3FFFF7988F8959AC" /
C
C     MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES
C
C     DATA SMALL(1) / 00564000000000000000B /
C     DATA SMALL(2) / 00000000000000000000B /
C     DATA LARGE(1) / 37757777777777777777B /
C     DATA LARGE(2) / 37157777777777777777B /
C     DATA RIGHT(1) / 15624000000000000000B /
C     DATA RIGHT(2) / 00000000000000000000B /
C     DATA DIVER(1) / 15634000000000000000B /
C     DATA DIVER(2) / 00000000000000000000B /
C     DATA LOG10(1) / 17164642023241175717B /
C     DATA LOG10(2) / 16367571421742254654B /
C
C     MACHINE CONSTANTS FOR THE CELERITY C1260
C
C     DATA SMALL(1), SMALL(2) / Z'00100000', Z'00000000' /
C     DATA LARGE(1), LARGE(2) / Z'7FEFFFFF', Z'FFFFFFFF' /
C     DATA RIGHT(1), RIGHT(2) / Z'3CA00000', Z'00000000' /
C     DATA DIVER(1), DIVER(2) / Z'3CB00000', Z'00000000' /
C     DATA LOG10(1), LOG10(2) / Z'3FD34413', Z'509F79FF' /
C
C     MACHINE CONSTANTS FOR THE CONVEX C-1
C
C     DATA SMALL(1), SMALL(2) / '00100000'X,'00000000'X /
C     DATA LARGE(1), LARGE(2) / '7FFFFFFF'X,'FFFFFFFF'X /
C     DATA RIGHT(1), RIGHT(2) / '3CC00000'X,'00000000'X /
C     DATA DIVER(1), DIVER(2) / '3CD00000'X,'00000000'X /
C     DATA LOG10(1), LOG10(2) / '3FF34413'X,'509F79FF'X /
C
C     MACHINE CONSTANTS FOR THE CRAY-1
C
C     DATA SMALL(1) / 201354000000000000000B /
C     DATA SMALL(2) / 000000000000000000000B /
C     DATA LARGE(1) / 577767777777777777777B /
C     DATA LARGE(2) / 000007777777777777774B /
C     DATA RIGHT(1) / 376434000000000000000B /
C     DATA RIGHT(2) / 000000000000000000000B /
C     DATA DIVER(1) / 376444000000000000000B /
C     DATA DIVER(2) / 000000000000000000000B /
C     DATA LOG10(1) / 377774642023241175717B /
C     DATA LOG10(2) / 000007571421742254654B /
C
C     MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200
C
C     NOTE - IT MAY BE APPROPRIATE TO INCLUDE THE FOLLOWING CARD -
C     STATIC DMACH(5)
C
C     DATA SMALL /    20K, 3*0 /
C     DATA LARGE / 77777K, 3*177777K /
C     DATA RIGHT / 31420K, 3*0 /
C     DATA DIVER / 32020K, 3*0 /
C     DATA LOG10 / 40423K, 42023K, 50237K, 74776K /
C
C     MACHINE CONSTANTS FOR THE ELXSI 6400
C     (ASSUMING REAL*8 IS THE DEFAULT DOUBLE PRECISION)
C
C     DATA SMALL(1), SMALL(2) / '00100000'X,'00000000'X /
C     DATA LARGE(1), LARGE(2) / '7FEFFFFF'X,'FFFFFFFF'X /
C     DATA RIGHT(1), RIGHT(2) / '3CB00000'X,'00000000'X /
C     DATA DIVER(1), DIVER(2) / '3CC00000'X,'00000000'X /
C     DATA LOG10(1), LOG10(2) / '3FD34413'X,'509F79FF'X /
C
C     MACHINE CONSTANTS FOR THE HARRIS 220
C
C     DATA SMALL(1), SMALL(2) / '20000000, '00000201 /
C     DATA LARGE(1), LARGE(2) / '37777777, '37777577 /
C     DATA RIGHT(1), RIGHT(2) / '20000000, '00000333 /
C     DATA DIVER(1), DIVER(2) / '20000000, '00000334 /
C     DATA LOG10(1), LOG10(2) / '23210115, '10237777 /
C
C     MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 SERIES
C
C     DATA SMALL(1), SMALL(2) / O402400000000, O000000000000 /
C     DATA LARGE(1), LARGE(2) / O376777777777, O777777777777 /
C     DATA RIGHT(1), RIGHT(2) / O604400000000, O000000000000 /
C     DATA DIVER(1), DIVER(2) / O606400000000, O000000000000 /
C     DATA LOG10(1), LOG10(2) / O776464202324, O117571775714 /
C
C     MACHINE CONSTANTS FOR THE HP 2100
C     THREE WORD DOUBLE PRECISION OPTION WITH FTN4
C
C     DATA SMALL(1), SMALL(2), SMALL(3) / 40000B,       0,       1 /
C     DATA LARGE(1), LARGE(2), LARGE(3) / 77777B, 177777B, 177776B /
C     DATA RIGHT(1), RIGHT(2), RIGHT(3) / 40000B,       0,    265B /
C     DATA DIVER(1), DIVER(2), DIVER(3) / 40000B,       0,    276B /
C     DATA LOG10(1), LOG10(2), LOG10(3) / 46420B,  46502B,  77777B /
C
C     MACHINE CONSTANTS FOR THE HP 2100
C     FOUR WORD DOUBLE PRECISION OPTION WITH FTN4
C
C     DATA SMALL(1), SMALL(2) /  40000B,       0 /
C     DATA SMALL(3), SMALL(4) /       0,       1 /
C     DATA LARGE(1), LARGE(2) /  77777B, 177777B /
C     DATA LARGE(3), LARGE(4) / 177777B, 177776B /
C     DATA RIGHT(1), RIGHT(2) /  40000B,       0 /
C     DATA RIGHT(3), RIGHT(4) /       0,    225B /
C     DATA DIVER(1), DIVER(2) /  40000B,       0 /
C     DATA DIVER(3), DIVER(4) /       0,    227B /
C     DATA LOG10(1), LOG10(2) /  46420B,  46502B /
C     DATA LOG10(3), LOG10(4) /  76747B, 176377B /
C
C     MACHINE CONSTANTS FOR THE HP 9000
C
C     DATA SMALL(1), SMALL(2) / 00040000000B, 00000000000B /
C     DATA LARGE(1), LARGE(2) / 17737777777B, 37777777777B /
C     DATA RIGHT(1), RIGHT(2) / 07454000000B, 00000000000B /
C     DATA DIVER(1), DIVER(2) / 07460000000B, 00000000000B /
C     DATA LOG10(1), LOG10(2) / 07764642023B, 12047674777B /
C
C     MACHINE CONSTANTS FOR THE IBM 360/370 SERIES,
C     THE XEROX SIGMA 5/7/9, THE SEL SYSTEMS 85/86, AND
C     THE PERKIN ELMER (INTERDATA) 7/32.
C
C     DATA SMALL(1), SMALL(2) / Z00100000, Z00000000 /
C     DATA LARGE(1), LARGE(2) / Z7FFFFFFF, ZFFFFFFFF /
C     DATA RIGHT(1), RIGHT(2) / Z33100000, Z00000000 /
C     DATA DIVER(1), DIVER(2) / Z34100000, Z00000000 /
C     DATA LOG10(1), LOG10(2) / Z41134413, Z509F79FF /
C
C     MACHINE CONSTANTS FOR THE IBM PC
C     ASSUMES THAT ALL ARITHMETIC IS DONE IN DOUBLE PRECISION
C     ON 8088, I.E., NOT IN 80 BIT FORM FOR THE 8087.
C
C     DATA SMALL(1)           / 2.23D-308              /
C     DATA LARGE(1)           / 1.79D+308              /
C     DATA RIGHT(1)           / 1.11D-16               /
C     DATA DIVER(1)           / 2.22D-16               /
C     DATA LOG10(1)           / 0.301029995663981195D0 /
C
C     MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR)
C
C     DATA SMALL(1), SMALL(2) / "033400000000, "000000000000 /
C     DATA LARGE(1), LARGE(2) / "377777777777, "344777777777 /
C     DATA RIGHT(1), RIGHT(2) / "113400000000, "000000000000 /
C     DATA DIVER(1), DIVER(2) / "114400000000, "000000000000 /
C     DATA LOG10(1), LOG10(2) / "177464202324, "144117571776 /
C
C     MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR)
C
C     DATA SMALL(1), SMALL(2) / "000400000000, "000000000000 /
C     DATA LARGE(1), LARGE(2) / "377777777777, "377777777777 /
C     DATA RIGHT(1), RIGHT(2) / "103400000000, "000000000000 /
C     DATA DIVER(1), DIVER(2) / "104400000000, "000000000000 /
C     DATA LOG10(1), LOG10(2) / "177464202324, "476747767461 /
C
C     MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING
C     32-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL).
C
C     DATA SMALL(1), SMALL(2) /    8388608,           0 /
C     DATA LARGE(1), LARGE(2) / 2147483647,          -1 /
C     DATA RIGHT(1), RIGHT(2) /  612368384,           0 /
C     DATA DIVER(1), DIVER(2) /  620756992,           0 /
C     DATA LOG10(1), LOG10(2) / 1067065498, -2063872008 /
C
C     DATA SMALL(1), SMALL(2) / O00040000000, O00000000000 /
C     DATA LARGE(1), LARGE(2) / O17777777777, O37777777777 /
C     DATA RIGHT(1), RIGHT(2) / O04440000000, O00000000000 /
C     DATA DIVER(1), DIVER(2) / O04500000000, O00000000000 /
C     DATA LOG10(1), LOG10(2) / O07746420232, O20476747770 /
C
C     MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING
C     16-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL).
C
C     DATA SMALL(1), SMALL(2) /    128,      0 /
C     DATA SMALL(3), SMALL(4) /      0,      0 /
C     DATA LARGE(1), LARGE(2) /  32767,     -1 /
C     DATA LARGE(3), LARGE(4) /     -1,     -1 /
C     DATA RIGHT(1), RIGHT(2) /   9344,      0 /
C     DATA RIGHT(3), RIGHT(4) /      0,      0 /
C     DATA DIVER(1), DIVER(2) /   9472,      0 /
C     DATA DIVER(3), DIVER(4) /      0,      0 /
C     DATA LOG10(1), LOG10(2) /  16282,   8346 /
C     DATA LOG10(3), LOG10(4) / -31493, -12296 /
C
C     DATA SMALL(1), SMALL(2) / O000200, O000000 /
C     DATA SMALL(3), SMALL(4) / O000000, O000000 /
C     DATA LARGE(1), LARGE(2) / O077777, O177777 /
C     DATA LARGE(3), LARGE(4) / O177777, O177777 /
C     DATA RIGHT(1), RIGHT(2) / O022200, O000000 /
C     DATA RIGHT(3), RIGHT(4) / O000000, O000000 /
C     DATA DIVER(1), DIVER(2) / O022400, O000000 /
C     DATA DIVER(3), DIVER(4) / O000000, O000000 /
C     DATA LOG10(1), LOG10(2) / O037632, O020232 /
C     DATA LOG10(3), LOG10(4) / O102373, O147770 /
C
C     MACHINE CONSTANTS FOR THE SUN
C
C     DATA SMALL(1), SMALL(2) / Z'00100000', Z'00000000' /
C     DATA LARGE(1), LARGE(2) / Z'7FEFFFFF', Z'FFFFFFFF' /
C     DATA RIGHT(1), RIGHT(2) / Z'3CA00000', Z'00000000' /
C     DATA DIVER(1), DIVER(2) / Z'3CB00000', Z'00000000' /
C     DATA LOG10(1), LOG10(2) / Z'3FD34413', Z'509F79FF' /
C
C     MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES FTN COMPILER
C
C     DATA SMALL(1), SMALL(2) / O000040000000, O000000000000 /
C     DATA LARGE(1), LARGE(2) / O377777777777, O777777777777 /
C     DATA RIGHT(1), RIGHT(2) / O170540000000, O000000000000 /
C     DATA DIVER(1), DIVER(2) / O170640000000, O000000000000 /
C     DATA LOG10(1), LOG10(2) / O177746420232, O411757177572 /
C
C     MACHINE CONSTANTS FOR VAX 11/780
C     (EXPRESSED IN INTEGER AND HEXADECIMAL)
C     THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSYEMS
C     THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS
C
C     DATA SMALL(1), SMALL(2) /        128,           0 /
C     DATA LARGE(1), LARGE(2) /     -32769,          -1 /
C     DATA RIGHT(1), RIGHT(2) /       9344,           0 /
C     DATA DIVER(1), DIVER(2) /       9472,           0 /
C     DATA LOG10(1), LOG10(2) /  546979738,  -805796613 /
C
C     DATA SMALL(1), SMALL(2) / Z00000080, Z00000000 /
C     DATA LARGE(1), LARGE(2) / ZFFFF7FFF, ZFFFFFFFF /
C     DATA RIGHT(1), RIGHT(2) / Z00002480, Z00000000 /
C     DATA DIVER(1), DIVER(2) / Z00002500, Z00000000 /
C     DATA LOG10(1), LOG10(2) / Z209A3F9A, ZCFF884FB /
C
C     MACHINE CONSTANTS FOR VAX 11/780 (G-FLOATING)
C     (EXPRESSED IN INTEGER AND HEXADECIMAL)
C     THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSYEMS
C     THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS
C
C     DATA SMALL(1), SMALL(2) /         16,           0 /
C     DATA LARGE(1), LARGE(2) /     -32769,          -1 /
C     DATA RIGHT(1), RIGHT(2) /      15552,           0 /
C     DATA DIVER(1), DIVER(2) /      15568,           0 /
C     DATA LOG10(1), LOG10(2) /  1142112243, 2046775455 /
C
C     DATA SMALL(1), SMALL(2) / Z00000010, Z00000000 /
C     DATA LARGE(1), LARGE(2) / ZFFFF7FFF, ZFFFFFFFF /
C     DATA RIGHT(1), RIGHT(2) / Z00003CC0, Z00000000 /
C     DATA DIVER(1), DIVER(2) / Z00003CD0, Z00000000 /
C     DATA LOG10(1), LOG10(2) / Z44133FF3, Z79FF509F /
C
C
C***FIRST EXECUTABLE STATEMENT  D1MACH
C
C THE CALL TO XERROR HAS BEEN COMMENTED TO GET A STAND ALONE ROUTINE
C IT CAN BE ADAPTED TO FIT THE LOCAL ERROR HANDLING PROCEDURES
C        NOTE ADDED BY F. ROMANI  7/11/89
C
C     IF (I .LT. 1  .OR.  I .GT. 5)
C    1   CALL XERROR ('D1MACH -- I OUT OF BOUNDS', 25, 1, 2)
C
      D1MACH = DMACH(I)
      RETURN
C
      END
*DECK DQXGS
      SUBROUTINE DQXGS(F,A,B,EPSABS,EPSREL,RESULT,ABSERR,IER,
     1   LIMIT,LENIW,LENW,LAST,IWORK,WORK)
C
C            THE ROUTINE CALCULATES AN APPROXIMATION RESULT TO A GIVEN
C            DEFINITE INTEGRAL  I = INTEGRAL OF F OVER (A,B),
C            HOPEFULLY SATISFYING FOLLOWING CLAIM FOR ACCURACY
C            ABS(I-RESULT).LE.MAX(EPSABS,EPSREL*ABS(I)).
C
C        PARAMETERS
C         ON ENTRY
C            F      - DOUBLE PRECISION
C                     FUNCTION SUBPROGRAM DEFINING THE INTEGRAND
C                     FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE
C                     DECLARED E X T E R N A L IN THE DRIVER PROGRAM.
C
C            A      - DOUBLE PRECISION
C                     LOWER LIMIT OF INTEGRATION
C
C            B      - DOUBLE PRECISION
C                     UPPER LIMIT OF INTEGRATION
C
C            EPSABS - DOUBLE PRECISION
C                     ABSOLUTE ACCURACY REQUESTED
C            EPSREL - DOUBLE PRECISION
C                     RELATIVE ACCURACY REQUESTED
C                     IF  EPSABS.LE.0
C                     AND EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28),
C                     THE ROUTINE WILL END WITH IER = 6.
C
C         ON RETURN
C            RESULT - DOUBLE PRECISION
C                     APPROXIMATION TO THE INTEGRAL
C
C            ABSERR - DOUBLE PRECISION
C                     ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR,
C                     WHICH SHOULD EQUAL OR EXCEED ABS(I-RESULT)
C
C            IER    - INTEGER
C                     IER = 0 NORMAL AND RELIABLE TERMINATION OF THE
C                             ROUTINE. IT IS ASSUMED THAT THE REQUESTED
C                             ACCURACY HAS BEEN ACHIEVED.
C                     IER.GT.0 ABNORMAL TERMINATION OF THE ROUTINE
C                             THE ESTIMATES FOR INTEGRAL AND ERROR ARE
C                             LESS RELIABLE. IT IS ASSUMED THAT THE
C                             REQUESTED ACCURACY HAS NOT BEEN ACHIEVED.
C            ERROR MESSAGES
C                     IER = 1 MAXIMUM NUMBER OF SUBDIVISIONS ALLOWED
C                             HAS BEEN ACHIEVED. ONE CAN ALLOW MORE SUB-
C                             DIVISIONS BY INCREASING THE VALUE OF LIMIT
C                             (AND TAKING THE ACCORDING DIMENSION
C                             ADJUSTMENTS INTO ACCOUNT. HOWEVER, IF
C                             THIS YIELDS NO IMPROVEMENT IT IS ADVISED
C                             TO ANALYZE THE INTEGRAND IN ORDER TO
C                             DETERMINE THE INTEGRATION DIFFICULTIES. IF
C                             THE POSITION OF A LOCAL DIFFICULTY CAN BE
C                             DETERMINED (E.G. SINGULARITY,
C                             DISCONTINUITY WITHIN THE INTERVAL) ONE
C                             WILL PROBABLY GAIN FROM SPLITTING UP THE
C                             INTERVAL AT THIS POINT AND CALLING THE
C                             INTEGRATOR ON THE SUBRANGES. IF POSSIBLE,
C                             AN APPROPRIATE SPECIAL-PURPOSE INTEGRATOR
C                             SHOULD BE USED, WHICH IS DESIGNED FOR
C                             HANDLING THE TYPE OF DIFFICULTY INVOLVED.
C                         = 2 THE OCCURRENCE OF ROUNDOFF ERROR IS DETEC-
C                             TED, WHICH PREVENTS THE REQUESTED
C                             TOLERANCE FROM BEING ACHIEVED.
C                             THE ERROR MAY BE UNDER-ESTIMATED.
C                         = 3 EXTREMELY BAD INTEGRAND BEHAVIOUR
C                             OCCURS AT SOME POINTS OF THE INTEGRATION
C                             INTERVAL.
C                         = 4 THE ALGORITHM DOES NOT CONVERGE.
C                             ROUNDOFF ERROR IS DETECTED IN THE
C                             EXTRAPOLATION TABLE. IT IS PRESUMED THAT
C                             THE REQUESTED TOLERANCE CANNOT BE
C                             ACHIEVED, AND THAT THE RETURNED RESULT IS
C                             THE BEST WHICH CAN BE OBTAINED.
C                         = 5 THE INTEGRAL IS PROBABLY DIVERGENT, OR
C                             SLOWLY CONVERGENT. IT MUST BE NOTED THAT
C                             DIVERGENCE CAN OCCUR WITH ANY OTHER VALUE
C                             OF IER.
C                         = 6 THE INPUT IS INVALID, BECAUSE
C                             (EPSABS.LE.0 AND
C                              EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28)
C                             OR LIMIT.LT.1 OR LENW.LT.LIMIT*46 OR
C                                LENIW .LT. LIMIT*3
C                             RESULT, ABSERR, LAST ARE SET TO
C                             ZERO.EXCEPT WHEN LIMIT OR LENW OR LENIW
C                             IS INVALID, IWORK(1), WORK(LIMIT*2+1) AND
C                             WORK(LIMIT*3+1) ARE SET TO ZERO, WORK(1)
C                             IS SET TO A AND WORK(LIMIT+1) TO B.
C
C         DIMENSIONING PARAMETERS
C            LIMIT - INTEGER
C                    LIMIT DETERMINES THE MAXIMUM NUMBER OF SUBINTERVALS
C                    IN THE PARTITION OF THE GIVEN INTEGRATION INTERVAL
C                    (A,B), LIMIT.GE.1.
C                    IF LIMIT.LT.1, THE ROUTINE WILL END WITH IER = 6.
C
C            LENW  - INTEGER
C                    DIMENSIONING PARAMETER FOR WORK
C                    LENW MUST BE AT LEAST LIMIT*46
C                    IF LENW.LT.LIMIT*46, THE ROUTINE WILL END
C                    WITH IER = 6.
C
C            LENIW - INTEGER
C                    DIMENSIONING PARAMETER FOR IWORK
C                    LENIW MUST BE AT LEAST LIMIT*3
C                    IF LENW.LT.LIMIT*3, THE ROUTINE WILL END
C                    WITH IER = 6.
C
C            LAST  - INTEGER
C                    ON RETURN, LAST EQUALS THE NUMBER OF SUBINTERVALS
C                    PRODUCED IN THE SUBDIVISION PROCESS, DETERMINES THE
C                    NUMBER OF SIGNIFICANT ELEMENTS ACTUALLY IN THE WORK
C                    ARRAYS.
C
C         WORK ARRAYS
C            IWORK - INTEGER
C                    VECTOR OF DIMENSION AT LEAST 3*LIMIT, THE FIRST K
C                    ELEMENTS OF WHICH CONTAIN POINTERS
C                    TO THE ERROR ESTIMATES OVER THE SUBINTERVALS
C                    SUCH THAT WORK(LIMIT*3+IWORK(1)),... ,
C                    WORK(LIMIT*3+IWORK(K)) FORM A DECREASING
C                    SEQUENCE, WITH K = LAST IF LAST.LE.(LIMIT/2+2),
C                    AND K = LIMIT+1-LAST OTHERWISE
C
C            WORK  - DOUBLE PRECISION
C                    VECTOR OF DIMENSION AT LEAST LENW
C                    ON RETURN
C                    WORK(1), ..., WORK(LAST) CONTAIN THE LEFT
C                     END-POINTS OF THE SUBINTERVALS IN THE
C                     PARTITION OF (A,B),
C                    WORK(LIMIT+1), ..., WORK(LIMIT+LAST) CONTAIN
C                     THE RIGHT END-POINTS,
C                    WORK(LIMIT*2+1), ..., WORK(LIMIT*2+LAST) CONTAIN
C                     THE INTEGRAL APPROXIMATIONS OVER THE SUBINTERVALS,
C                    WORK(LIMIT*3+1), ..., WORK(LIMIT*3+LAST)
C                     CONTAIN THE ERROR ESTIMATES.
C                    WORK(LIMIT*4+1), ... IS THE AREA RESERVED TO STORE
C                     FUNCTIONAL VALUES.
C***ROUTINES CALLED  DQXGSE,XERROR
C
C THE CALL TO XERROR HAS BEEN COMMENTED TO GET A STAND ALONE ROUTINE
C IT CAN BE ADAPTED TO FIT THE LOCAL ERROR HANDLING PROCEDURES
C
      DOUBLE PRECISION A,ABSERR,B,EPSABS,EPSREL,F,RESULT,WORK
      INTEGER IER,IWORK,LAST,LENIW,LENW,LIMIT,LVL,L1,L2,L3,L4,L5
C
      DIMENSION IWORK(LENIW),WORK(LENW)
C
      EXTERNAL F
C
C         CHECK VALIDITY OF LIMIT,LENIW AND LENW.
C
C***FIRST EXECUTABLE STATEMENT  DQXGS
      IER = 6
      LAST = 0
      RESULT = 0.0D+00
      ABSERR = 0.0D+00
      IF(LIMIT.LT.1.OR.LENIW.LT.LIMIT*3.OR.LENW.LT.LIMIT*46) GO TO 10
C
C         PREPARE CALL FOR DQXGSE.
C
      L1 = LIMIT+1
      L2 = LIMIT+L1
      L3 = LIMIT+L2
      L4 = LIMIT+L3
      L5 = 21*LIMIT+L3
C
      CALL DQXGSE(F,A,B,EPSABS,EPSREL,LIMIT,RESULT,ABSERR,
     *  IER,WORK(1),WORK(L1),WORK(L2),WORK(L3),IWORK,LAST,
     *  WORK(L4),WORK(L5),IWORK(L1),IWORK(L2))
C
C         CALL ERROR HANDLER IF NECESSARY.
C
      LVL = 0
10    IF(IER.EQ.6) LVL = 1
C
C THE CALL TO XERROR HAS BEEN COMMENTED TO GET A STAND ALONE ROUTINE
C IT CAN BE ADAPTED TO FIT THE LOCAL ERROR HANDLING PROCEDURES
C
C     IF(IER.NE.0)CALL XERROR('ABNORMAL RETURN FROM DQXGS  ',28,IER,LVL)
      RETURN
      END
*DECK DQXGSE
      SUBROUTINE DQXGSE(F,A,B,EPSABS,EPSREL,LIMIT,RESULT,ABSERR,
     *   IER,ALIST,BLIST,RLIST,ELIST,IORD,LAST,VALP,VALN,LP,LN)
C
C            THE ROUTINE CALCULATES AN APPROXIMATION RESULT TO A GIVEN
C            DEFINITE INTEGRAL I = INTEGRAL OF F OVER (A,B),
C            HOPEFULLY SATISFYING FOLLOWING CLAIM FOR ACCURACY
C            ABS(I-RESULT).LE.MAX(EPSABS,EPSREL*ABS(I)).
C
C        PARAMETERS
C         ON ENTRY
C            F      - DOUBLE PRECISION
C                     FUNCTION SUBPROGRAM DEFINING THE INTEGRAND
C                     FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE
C                     DECLARED E X T E R N A L IN THE DRIVER PROGRAM.
C
C            A      - DOUBLE PRECISION
C                     LOWER LIMIT OF INTEGRATION
C
C            B      - DOUBLE PRECISION
C                     UPPER LIMIT OF INTEGRATION
C
C            EPSABS - DOUBLE PRECISION
C                     ABSOLUTE ACCURACY REQUESTED
C            EPSREL - DOUBLE PRECISION
C                     RELATIVE ACCURACY REQUESTED
C                     IF  EPSABS.LE.0
C                     AND EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28),
C                     THE ROUTINE WILL END WITH IER = 6.
C
C            LIMIT  - INTEGER
C                     GIVES AN UPPERBOUND ON THE NUMBER OF SUBINTERVALS
C                     IN THE PARTITION OF (A,B)
C
C         ON RETURN
C            RESULT - DOUBLE PRECISION
C                     APPROXIMATION TO THE INTEGRAL
C
C            ABSERR - DOUBLE PRECISION
C                     ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR,
C                     WHICH SHOULD EQUAL OR EXCEED ABS(I-RESULT)
C
C            IER    - INTEGER
C                     IER = 0 NORMAL AND RELIABLE TERMINATION OF THE
C                             ROUTINE. IT IS ASSUMED THAT THE REQUESTED
C                             ACCURACY HAS BEEN ACHIEVED.
C                     IER.GT.0 ABNORMAL TERMINATION OF THE ROUTINE
C                             THE ESTIMATES FOR INTEGRAL AND ERROR ARE
C                             LESS RELIABLE. IT IS ASSUMED THAT THE
C                             REQUESTED ACCURACY HAS NOT BEEN ACHIEVED.
C            ERROR MESSAGES
C                         = 1 MAXIMUM NUMBER OF SUBDIVISIONS ALLOWED
C                             HAS BEEN ACHIEVED. ONE CAN ALLOW MORE SUB-
C                             DIVISIONS BY INCREASING THE VALUE OF LIMIT
C                             (AND TAKING THE ACCORDING DIMENSION
C                             ADJUSTMENTS INTO ACCOUNT). HOWEVER, IF
C                             THIS YIELDS NO IMPROVEMENT IT IS ADVISED
C                             TO ANALYZE THE INTEGRAND IN ORDER TO
C                             DETERMINE THE INTEGRATION DIFFICULTIES. IF
C                             THE POSITION OF A LOCAL DIFFICULTY CAN BE
C                             DETERMINED (E.G. SINGULARITY,
C                             DISCONTINUITY WITHIN THE INTERVAL) ONE
C                             WILL PROBABLY GAIN FROM SPLITTING UP THE
C                             INTERVAL AT THIS POINT AND CALLING THE
C                             INTEGRATOR ON THE SUBRANGES. IF POSSIBLE,
C                             AN APPROPRIATE SPECIAL-PURPOSE INTEGRATOR
C                             SHOULD BE USED, WHICH IS DESIGNED FOR
C                             HANDLING THE TYPE OF DIFFICULTY INVOLVED.
C                         = 2 THE OCCURRENCE OF ROUNDOFF ERROR IS DETEC-
C                             TED, WHICH PREVENTS THE REQUESTED
C                             TOLERANCE FROM BEING ACHIEVED.
C                             THE ERROR MAY BE UNDER-ESTIMATED.
C                         = 3 EXTREMELY BAD INTEGRAND BEHAVIOUR
C                             OCCURS AT SOME POINTS OF THE INTEGRATION
C                             INTERVAL.
C                         = 4 THE ALGORITHM DOES NOT CONVERGE.
C                             ROUNDOFF ERROR IS DETECTED IN THE
C                             EXTRAPOLATION TABLE.
C                             IT IS PRESUMED THAT THE REQUESTED
C                             TOLERANCE CANNOT BE ACHIEVED, AND THAT THE
C                             RETURNED RESULT IS THE BEST WHICH CAN BE
C                             OBTAINED.
C                         = 5 THE INTEGRAL IS PROBABLY DIVERGENT, OR
C                             SLOWLY CONVERGENT. IT MUST BE NOTED THAT
C                             DIVERGENCE CAN OCCUR WITH ANY OTHER VALUE
C                             OF IER.
C                         = 6 THE INPUT IS INVALID, BECAUSE
C                             EPSABS.LE.0 AND
C                             EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28).
C                             RESULT, ABSERR, LAST, RLIST(1),
C                             IORD(1) AND ELIST(1) ARE SET TO ZERO.
C                             ALIST(1) AND BLIST(1) ARE SET TO A AND B
C                             RESPECTIVELY.
C
C            ALIST  - DOUBLE PRECISION
C                     VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
C                      LAST  ELEMENTS OF WHICH ARE THE LEFT END POINTS
C                     OF THE SUBINTERVALS IN THE PARTITION OF THE
C                     GIVEN INTEGRATION RANGE (A,B)
C
C            BLIST  - DOUBLE PRECISION
C                     VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
C                      LAST  ELEMENTS OF WHICH ARE THE RIGHT END POINTS
C                     OF THE SUBINTERVALS IN THE PARTITION OF THE GIVEN
C                     INTEGRATION RANGE (A,B)
C
C            RLIST  - DOUBLE PRECISION
C                     VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
C                      LAST  ELEMENTS OF WHICH ARE THE INTEGRAL
C                     APPROXIMATIONS ON THE SUBINTERVALS
C
C            ELIST  - DOUBLE PRECISION
C                     VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
C                      LAST  ELEMENTS OF WHICH ARE THE MODULI OF THE
C                     ABSOLUTE ERROR ESTIMATES ON THE SUBINTERVALS
C
C            IORD   - INTEGER
C                     VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST K
C                     ELEMENTS OF WHICH ARE POINTERS TO THE
C                     ERROR ESTIMATES OVER THE SUBINTERVALS,
C                     SUCH THAT ELIST(IORD(1)), ..., ELIST(IORD(K))
C                     FORM A DECREASING SEQUENCE, WITH K = LAST
C                     IF LAST.LE.(LIMIT/2+2), AND K = LIMIT+1-LAST
C                     OTHERWISE
C
C            LAST   - INTEGER
C                     NUMBER OF SUBINTERVALS ACTUALLY PRODUCED IN THE
C                     SUBDIVISION PROCESS
C
C            VALP   - DOUBLE PRECISION
C            VALN     ARRAYS OF DIMENSION AT LEAST (21,LIMIT) USED TO
C                     SAVE THE FUNCTIONAL VALUES
C
C            LP     - INTEGER
C            LN       VECTORS OF DIMENSION AT LEAST LIMIT, USED TO
C                     STORE THE ACTUAL NUMBER OF FUNCTIONAL VALUES
C                     SAVED IN THE CORRESPONDING COLUMN
C                     OF VALP,VALN
C
C***ROUTINES CALLED  F,D1MACH,DQELG,DQXLQM,DQPSRT,DQXRRD,DQXCPY
C
      DOUBLE PRECISION A,ABSEPS,ABSERR,ALIST,AREA,AREA1,AREA12,AREA2,A1,
     *  A2,B,BLIST,B1,B2,CORREC,DABS,DEFABS,DEFAB1,DEFAB2,D1MACH,DMAX1,
     *  DRES,ELIST,EPMACH,EPSABS,EPSREL,ERLARG,ERLAST,ERRBND,ERRMAX,
     *  ERROR1,ERROR2,ERRO12,ERRSUM,ERTEST,F,OFLOW,RESABS,RESEPS,RESULT,
     *  RES3LA,RLIST,RLIST2,SMALL,UFLOW,
     *  VALP,VALN,VP1,VP2,VN1,VN2
      INTEGER ID,IER,IERRO,IORD,IROFF1,IROFF2,IROFF3,JUPBND,K,KSGN,
     *  KTMIN,LAST,LIMIT,MAXERR,NRES,NRMAX,NUMRL2,
     *  LP,LN,LP1,LP2,LN1,LN2
      LOGICAL EXTRAP,NOEXT
      DIMENSION ALIST(LIMIT),BLIST(LIMIT),ELIST(LIMIT),IORD(LIMIT),
     * RES3LA(3),RLIST(LIMIT),RLIST2(52),
     * VALP(21,LIMIT),VALN(21,LIMIT),LP(LIMIT),LN(LIMIT),
     * VP1(21),VP2(21),VN1(21),VN2(21)
      EXTERNAL F
C
C            MACHINE DEPENDENT CONSTANTS
C            ---------------------------
C
C          EPMACH IS THE LARGEST RELATIVE SPACING.
C          UFLOW IS THE SMALLEST POSITIVE MAGNITUDE.
C          OFLOW IS THE LARGEST POSITIVE MAGNITUDE.
C
C***FIRST EXECUTABLE STATEMENT  DQXGSE
      EPMACH = D1MACH(4)
C
C            TEST ON VALIDITY OF PARAMETERS
C            ------------------------------
      IER = 0
      LAST = 0
      RESULT = 0.0D+00
      ABSERR = 0.0D+00
      ALIST(1) = A
      BLIST(1) = B
      RLIST(1) = 0.0D+00
      ELIST(1) = 0.0D+00
      IF(EPSABS.LE.0.0D+00.AND.EPSREL.LT.DMAX1(0.5D+02*EPMACH,0.5D-28))
     *   IER = 6
      IF(IER.EQ.6) GO TO 999
C
C           FIRST APPROXIMATION TO THE INTEGRAL
C           -----------------------------------
C
      UFLOW = D1MACH(1)
      OFLOW = D1MACH(2)
      IERRO = 0
      LP(1)=1
      LN(1)=1
      VALP(1,1)=F((A+B)*0.5D0)
      VALN(1,1)=VALP(1,1)
      CALL DQXLQM(F,A,B,RESULT,ABSERR,DEFABS,RESABS,
     *           VALP(1,1),VALN(1,1),LP(1),LN(1), 2  )
C
C           TEST ON ACCURACY.
C
      DRES = DABS(RESULT)
      ERRBND = DMAX1(EPSABS,EPSREL*DRES)
      LAST = 1
      RLIST(1) = RESULT
      ELIST(1) = ABSERR
      IORD(1) = 1
      IF(ABSERR.LE.1.0D+02*EPMACH*DEFABS.AND.ABSERR.GT.ERRBND) IER = 2
      IF(LIMIT.EQ.1) IER = 1
      IF(IER.NE.0.OR.(ABSERR.LE.ERRBND.AND.ABSERR.NE.RESABS).OR.
     *  ABSERR.EQ.0.0D+00) GO TO 999
C
C           INITIALIZATION
C           --------------
C
      RLIST2(1) = RESULT
      ERRMAX = ABSERR
      MAXERR = 1
      AREA = RESULT
      ERRSUM = ABSERR
      ABSERR = OFLOW
      NRMAX = 1
      NRES = 0
      NUMRL2 = 2
      KTMIN = 0
      EXTRAP = .FALSE.
      NOEXT = .FALSE.
      IROFF1 = 0
      IROFF2 = 0
      IROFF3 = 0
      KSGN = -1
      IF(DRES.GE.(0.1D+01-0.5D+02*EPMACH)*DEFABS) KSGN = 1
C
C           MAIN DO-LOOP
C           ------------
C
      DO 90 LAST = 2,LIMIT
C
C           BISECT THE SUBINTERVAL WITH THE NRMAX-TH LARGEST ERROR
C           ESTIMATE.
C
        A1 = ALIST(MAXERR)
        B1 = 0.5D+00*(ALIST(MAXERR)+BLIST(MAXERR))
        A2 = B1
        B2 = BLIST(MAXERR)
        ERLAST = ERRMAX
        CALL DQXRRD(F,VALN(1,MAXERR),LN(MAXERR),B1,A1,VN1,VP1,LN1,LP1)
        CALL DQXLQM(F,A1,B1,AREA1,ERROR1,RESABS,DEFAB1,VP1,VN1,LP1,LN1,
     *              2)
        CALL DQXRRD(F,VALP(1,MAXERR),LP(MAXERR),A2,B2,VP2,VN2,LP2,LN2)
        CALL DQXLQM(F,A2,B2,AREA2,ERROR2,RESABS,DEFAB2,VP2,VN2,LP2,LN2,
     *              2)
C
C           IMPROVE PREVIOUS APPROXIMATIONS TO INTEGRAL
C           AND ERROR AND TEST FOR ACCURACY.
C
        AREA12 = AREA1+AREA2
        ERRO12 = ERROR1+ERROR2
        ERRSUM = ERRSUM+ERRO12-ERRMAX
        AREA = AREA+AREA12-RLIST(MAXERR)
        IF(DEFAB1.EQ.ERROR1.OR.DEFAB2.EQ.ERROR2) GO TO 15
        IF(DABS(RLIST(MAXERR)-AREA12).GT.0.1D-04*DABS(AREA12)
     *  .OR.ERRO12.LT.0.99D+00*ERRMAX) GO TO 10
        IF(EXTRAP) IROFF2 = IROFF2+1
        IF(.NOT.EXTRAP) IROFF1 = IROFF1+1
   10   IF(LAST.GT.10.AND.ERRO12.GT.ERRMAX) IROFF3 = IROFF3+1
   15   RLIST(MAXERR) = AREA1
        RLIST(LAST) = AREA2
        ERRBND = DMAX1(EPSABS,EPSREL*DABS(AREA))
C
C           TEST FOR ROUNDOFF ERROR AND EVENTUALLY SET ERROR FLAG.
C
        IF(IROFF1+IROFF2.GE.10.OR.IROFF3.GE.20) IER = 2
        IF(IROFF2.GE.5) IERRO = 3
C
C           SET ERROR FLAG IN THE CASE THAT THE NUMBER OF SUBINTERVALS
C           EQUALS LIMIT.
C
        IF(LAST.EQ.LIMIT) IER = 1
C
C           SET ERROR FLAG IN THE CASE OF BAD INTEGRAND BEHAVIOUR
C           AT A POINT OF THE INTEGRATION RANGE.
C
        IF(DMAX1(DABS(A1),DABS(B2)).LE.(0.1D+01+0.1D+03*EPMACH)*
     *  (DABS(A2)+0.1D+04*UFLOW)) IER = 4
C
C           APPEND THE NEWLY-CREATED INTERVALS TO THE LIST.
C
        IF(ERROR2.GT.ERROR1) GO TO 20
        ALIST(LAST) = A2
        BLIST(MAXERR) = B1
        BLIST(LAST) = B2
        ELIST(MAXERR) = ERROR1
        ELIST(LAST) = ERROR2
        CALL DQXCPY(VALP(1,MAXERR),VP1,LP1)
        LP(MAXERR)=LP1
        CALL DQXCPY(VALN(1,MAXERR),VN1,LN1)
        LN(MAXERR)=LN1
        CALL DQXCPY(VALP(1,LAST),VP2,LP2)
        LP(LAST)=LP2
        CALL DQXCPY(VALN(1,LAST),VN2,LN2)
        LN(LAST)=LN2
        GO TO 30
   20   ALIST(MAXERR) = A2
        ALIST(LAST) = A1
        BLIST(LAST) = B1
        RLIST(MAXERR) = AREA2
        RLIST(LAST) = AREA1
        ELIST(MAXERR) = ERROR2
        ELIST(LAST) = ERROR1
        CALL DQXCPY(VALP(1,MAXERR),VP2,LP2)
        LP(MAXERR)=LP2
        CALL DQXCPY(VALN(1,MAXERR),VN2,LN2)
        LN(MAXERR)=LN2
        CALL DQXCPY(VALP(1,LAST),VP1,LP1)
        LP(LAST)=LP1
        CALL DQXCPY(VALN(1,LAST),VN1,LN1)
        LN(LAST)=LN1
C
C           CALL SUBROUTINE DQPSRT TO MAINTAIN THE DESCENDING ORDERING
C           IN THE LIST OF ERROR ESTIMATES AND SELECT THE SUBINTERVAL
C           WITH NRMAX-TH LARGEST ERROR ESTIMATE (TO BE BISECTED NEXT).
C
   30   CALL DQPSRT(LIMIT,LAST,MAXERR,ERRMAX,ELIST,IORD,NRMAX)
C ***JUMP OUT OF DO-LOOP
        IF(ERRSUM.LE.ERRBND) GO TO 115
C ***JUMP OUT OF DO-LOOP
        IF(IER.NE.0) GO TO 100
        IF(LAST.EQ.2) GO TO 80
        IF(NOEXT) GO TO 90
        ERLARG = ERLARG-ERLAST
        IF(DABS(B1-A1).GT.SMALL) ERLARG = ERLARG+ERRO12
        IF(EXTRAP) GO TO 40
C
C           TEST WHETHER THE INTERVAL TO BE BISECTED NEXT IS THE
C           SMALLEST INTERVAL.
C
        IF(DABS(BLIST(MAXERR)-ALIST(MAXERR)).GT.SMALL) GO TO 90
        EXTRAP = .TRUE.
        NRMAX = 2
C
C           THE BOUND 0.3*ERTEST HAS BEEN INTRODUCED TO PERFORM A
C           MORE CAUTIOUS EXTRAPOLATION THAN IN THE ORIGINAL DQAGSE
C           ROUTINE
C
   40   IF(IERRO.EQ.3.OR.ERLARG.LE.0.3D0*ERTEST) GO TO 60
C
C           THE SMALLEST INTERVAL HAS THE LARGEST ERROR.
C           BEFORE BISECTING DECREASE THE SUM OF THE ERRORS OVER THE
C           LARGER INTERVALS (ERLARG) AND PERFORM EXTRAPOLATION.
C
        ID = NRMAX
        JUPBND = LAST
        IF(LAST.GT.(2+LIMIT/2)) JUPBND = LIMIT+3-LAST
        DO 50 K = ID,JUPBND
          MAXERR = IORD(NRMAX)
          ERRMAX = ELIST(MAXERR)
C ***JUMP OUT OF DO-LOOP
          IF(DABS(BLIST(MAXERR)-ALIST(MAXERR)).GT.SMALL) GO TO 90
          NRMAX = NRMAX+1
   50   CONTINUE
C
C           PERFORM EXTRAPOLATION.
C
   60   NUMRL2 = NUMRL2+1
        RLIST2(NUMRL2) = AREA
        CALL DQELG(NUMRL2,RLIST2,RESEPS,ABSEPS,RES3LA,NRES)
        KTMIN = KTMIN+1
        IF(KTMIN.GT.5.AND.ABSERR.LT.0.1D-02*ERRSUM) IER = 5
        IF(ABSEPS.GE.ABSERR) GO TO 70
        KTMIN = 0
        ABSERR = ABSEPS
        RESULT = RESEPS
        CORREC = ERLARG
        ERTEST = DMAX1(EPSABS,EPSREL*DABS(RESEPS))
C ***JUMP OUT OF DO-LOOP
        IF(ABSERR.LE.ERTEST) GO TO 100
C
C           PREPARE BISECTION OF THE SMALLEST INTERVAL.
C
   70   IF(NUMRL2.EQ.1) NOEXT = .TRUE.
        IF(IER.EQ.5) GO TO 100
        MAXERR = IORD(1)
        ERRMAX = ELIST(MAXERR)
        NRMAX = 1
        EXTRAP = .FALSE.
        SMALL = SMALL*0.5D+00
        ERLARG = ERRSUM
        GO TO 90
   80   SMALL = DABS(B-A)*0.375D+00
        ERLARG = ERRSUM
        ERTEST = ERRBND
        RLIST2(2) = AREA
   90 CONTINUE
C
C           SET FINAL RESULT AND ERROR ESTIMATE.
C           ------------------------------------
C
  100 IF(ABSERR.EQ.OFLOW) GO TO 115
      IF(IER+IERRO.EQ.0) GO TO 110
      IF(IERRO.EQ.3) ABSERR = ABSERR+CORREC
      IF(IER.EQ.0) IER = 3
      IF(RESULT.NE.0.0D+00.AND.AREA.NE.0.0D+00) GO TO 105
      IF(ABSERR.GT.ERRSUM) GO TO 115
      IF(AREA.EQ.0.0D+00) GO TO 130
      GO TO 110
  105 IF(ABSERR/DABS(RESULT).GT.ERRSUM/DABS(AREA)) GO TO 115
C
C           TEST ON DIVERGENCE.
C
  110 IF(KSGN.EQ.(-1).AND.DMAX1(DABS(RESULT),DABS(AREA)).LE.
     * DEFABS*0.1D-01) GO TO 130
      IF(0.1D-01.GT.(RESULT/AREA).OR.(RESULT/AREA).GT.0.1D+03
     * .OR.ERRSUM.GT.DABS(AREA)) IER = 6
      GO TO 130
C
C           COMPUTE GLOBAL INTEGRAL SUM.
C
  115 RESULT = 0.0D+00
      DO 120 K = 1,LAST
         RESULT = RESULT+RLIST(K)
  120 CONTINUE
      ABSERR = ERRSUM
  130 IF(IER.GT.2) IER = IER-1
  999 RETURN
      END
*DECK DQXG
      SUBROUTINE DQXG(F,A,B,EPSABS,EPSREL,KEY,RESULT,ABSERR,IER,
     1   LIMIT,LENIW,LENW,LAST,IWORK,WORK)
C
C            THE ROUTINE CALCULATES AN APPROXIMATION RESULT TO A GIVEN
C            DEFINITE INTEGRAL  I = INTEGRAL OF F OVER (A,B),
C            HOPEFULLY SATISFYING FOLLOWING CLAIM FOR ACCURACY
C            ABS(I-RESULT).LE.MAX(EPSABS,EPSREL*ABS(I)).
C
C        PARAMETERS
C         ON ENTRY
C            F      - DOUBLE PRECISION
C                     FUNCTION SUBPROGRAM DEFINING THE INTEGRAND
C                     FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE
C                     DECLARED E X T E R N A L IN THE DRIVER PROGRAM.
C
C            A      - DOUBLE PRECISION
C                     LOWER LIMIT OF INTEGRATION
C
C            B      - DOUBLE PRECISION
C                     UPPER LIMIT OF INTEGRATION
C
C            EPSABS - DOUBLE PRECISION
C                     ABSOLUTE ACCURACY REQUESTED
C            EPSREL - DOUBLE PRECISION
C                     RELATIVE ACCURACY REQUESTED
C                     IF  EPSABS.LE.0
C                     AND EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28),
C                     THE ROUTINE WILL END WITH IER = 6.
C
C            KEY    - INTEGER
C                     KEY FOR CHOICE OF LOCAL INTEGRATION RULE
C                     RMS FORMULAS ARE USED WITH
C                      13 - 19               POINTS IF KEY.LT.1,
C                      13 - 19 - (27)        POINTS IF KEY = 1,
C                      13 - 19 - (27) - (41) POINTS IF KEY = 2,
C                           19 -  27  - (41) POINTS IF KEY = 3,
C                                 27  -  41  POINTS IF KEY.GT.3.
C
C                         (RULES) USED IF THE FUNCTION APPEARS
C                         ENOUGH REGULAR IN THE SUBINTERVAL
C
C         ON RETURN
C            RESULT - DOUBLE PRECISION
C                     APPROXIMATION TO THE INTEGRAL
C
C            ABSERR - DOUBLE PRECISION
C                     ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR,
C                     WHICH SHOULD EQUAL OR EXCEED ABS(I-RESULT)
C
C            IER    - INTEGER
C                     IER = 0 NORMAL AND RELIABLE TERMINATION OF THE
C                             ROUTINE. IT IS ASSUMED THAT THE REQUESTED
C                             ACCURACY HAS BEEN ACHIEVED.
C                     IER.GT.0 ABNORMAL TERMINATION OF THE ROUTINE
C                             THE ESTIMATES FOR RESULT AND ERROR ARE
C                             LESS RELIABLE. IT IS ASSUMED THAT THE
C                             REQUESTED ACCURACY HAS NOT BEEN ACHIEVED.
C            ERROR MESSAGES
C                     IER = 1 MAXIMUM NUMBER OF SUBDIVISIONS ALLOWED
C                             HAS BEEN ACHIEVED. ONE CAN ALLOW MORE SUB-
C                             DIVISIONS BY INCREASING THE VALUE OF LIMIT
C                             (AND TAKING THE ACCORDING DIMENSION
C                             ADJUSTMENTS INTO ACCOUNT. HOWEVER, IF
C                             THIS YIELDS NO IMPROVEMENT IT IS ADVISED
C                             TO ANALYZE THE INTEGRAND IN ORDER TO
C                             DETERMINE THE INTEGRATION DIFFICULTIES. IF
C                             THE POSITION OF A LOCAL DIFFICULTY CAN BE
C                             DETERMINED (E.G. SINGULARITY,
C                             DISCONTINUITY WITHIN THE INTERVAL) ONE
C                             WILL PROBABLY GAIN FROM SPLITTING UP THE
C                             INTERVAL AT THIS POINT AND CALLING THE
C                             INTEGRATOR ON THE SUBRANGES. IF POSSIBLE,
C                             AN APPROPRIATE SPECIAL-PURPOSE INTEGRATOR
C                             SHOULD BE USED, WHICH IS DESIGNED FOR
C                             HANDLING THE TYPE OF DIFFICULTY INVOLVED.
C                         = 2 THE OCCURRENCE OF ROUNDOFF ERROR IS DETEC-
C                             TED, WHICH PREVENTS THE REQUESTED
C                             TOLERANCE FROM BEING ACHIEVED.
C                             THE ERROR MAY BE UNDER-ESTIMATED.
C                         = 3 EXTREMELY BAD INTEGRAND BEHAVIOUR
C                             OCCURS AT SOME POINTS OF THE INTEGRATION
C                             INTERVAL.
C                         = 4 THE ALGORITHM DOES NOT CONVERGE.
C                             ROUNDOFF ERROR IS DETECTED IN THE
C                             EXTRAPOLATION TABLE. IT IS PRESUMED THAT
C                             THE REQUESTED TOLERANCE CANNOT BE
C                             ACHIEVED, AND THAT THE RETURNED RESULT IS
C                             THE BEST WHICH CAN BE OBTAINED.
C                         = 5 THE INTEGRAL IS PROBABLY DIVERGENT, OR
C                             SLOWLY CONVERGENT. IT MUST BE NOTED THAT
C                             DIVERGENCE CAN OCCUR WITH ANY OTHER VALUE
C                             OF IER.
C                         = 6 THE INPUT IS INVALID, BECAUSE
C                             (EPSABS.LE.0 AND
C                              EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28)
C                             OR LIMIT.LT.1 OR LENW.LT.LIMIT*46 OR
C                                LENIW .LT. LIMIT*3
C                             RESULT, ABSERR, LAST ARE SET TO
C                             ZERO.EXCEPT WHEN LIMIT OR LENW OR LENIW
C                             IS INVALID, IWORK(1), WORK(LIMIT*2+1) AND
C                             WORK(LIMIT*3+1) ARE SET TO ZERO, WORK(1)
C                             IS SET TO A AND WORK(LIMIT+1) TO B.
C
C         DIMENSIONING PARAMETERS
C            LIMIT - INTEGER
C                    LIMIT DETERMINES THE MAXIMUM NUMBER OF SUBINTERVALS
C                    IN THE PARTITION OF THE GIVEN INTEGRATION INTERVAL
C                    (A,B), LIMIT.GE.1.
C                    IF LIMIT.LT.1, THE ROUTINE WILL END WITH IER = 6.
C
C            LENW  - INTEGER
C                    DIMENSIONING PARAMETER FOR WORK
C                    LENW MUST BE AT LEAST LIMIT*46
C                    IF LENW.LT.LIMIT*46, THE ROUTINE WILL END
C                    WITH IER = 6.
C
C            LENIW - INTEGER
C                    DIMENSIONING PARAMETER FOR IWORK
C                    LENIW MUST BE AT LEAST LIMIT*3
C                    IF LENW.LT.LIMIT*3, THE ROUTINE WILL END
C                    WITH IER = 6.
C
C            LAST  - INTEGER
C                    ON RETURN, LAST EQUALS THE NUMBER OF SUBINTERVALS
C                    PRODUCED IN THE SUBDIVISION PROCESS, DETERMINES THE
C                    NUMBER OF SIGNIFICANT ELEMENTS ACTUALLY IN THE WORK
C                    ARRAYS.
C
C         WORK ARRAYS
C            IWORK - INTEGER
C                    VECTOR OF DIMENSION AT LEAST 3*LIMIT, THE FIRST K
C                    ELEMENTS OF WHICH CONTAIN POINTERS
C                    TO THE ERROR ESTIMATES OVER THE SUBINTERVALS
C                    SUCH THAT WORK(LIMIT*3+IWORK(1)),... ,
C                    WORK(LIMIT*3+IWORK(K)) FORM A DECREASING
C                    SEQUENCE, WITH K = LAST IF LAST.LE.(LIMIT/2+2),
C                    AND K = LIMIT+1-LAST OTHERWISE
C
C            WORK  - DOUBLE PRECISION
C                    VECTOR OF DIMENSION AT LEAST LENW
C                    ON RETURN
C                    WORK(1), ..., WORK(LAST) CONTAIN THE LEFT
C                     END-POINTS OF THE SUBINTERVALS IN THE
C                     PARTITION OF (A,B),
C                    WORK(LIMIT+1), ..., WORK(LIMIT+LAST) CONTAIN
C                     THE RIGHT END-POINTS,
C                    WORK(LIMIT*2+1), ..., WORK(LIMIT*2+LAST) CONTAIN
C                     THE INTEGRAL APPROXIMATIONS OVER THE SUBINTERVALS,
C                    WORK(LIMIT*3+1), ..., WORK(LIMIT*3+LAST)
C                     CONTAIN THE ERROR ESTIMATES.
C                    WORK(LIMIT*4+1), ... IS THE AREA RESERVED TO STORE
C                     FUNCTIONAL VALUES.
C***ROUTINES CALLED  DQXGE,XERROR
C
C THE CALL TO XERROR HAS BEEN COMMENTED TO GET A STAND ALONE ROUTINE
C IT CAN BE ADAPTED TO FIT THE LOCAL ERROR HANDLING PROCEDURES
C
C
      DOUBLE PRECISION A,ABSERR,B,EPSABS,EPSREL,F,RESULT,WORK
      INTEGER KEY,IER,IWORK,LAST,LENIW,LENW,LIMIT,LVL,L1,L2,L3,L4,L5
C
      DIMENSION IWORK(LENIW),WORK(LENW)
C
      EXTERNAL F
C
C         CHECK VALIDITY OF LIMIT,LENIW AND LENW.
C
C***FIRST EXECUTABLE STATEMENT  DQXG
      IER = 6
      LAST = 0
      RESULT = 0.0D+00
      ABSERR = 0.0D+00
      IF(LIMIT.LT.1.OR.LENIW.LT.LIMIT*3.OR.LENW.LT.LIMIT*46) GO TO 10
C
C         PREPARE CALL FOR DQXGE.
C
      L1 = LIMIT+1
      L2 = LIMIT+L1
      L3 = LIMIT+L2
      L4 = LIMIT+L3
      L5 = 21*LIMIT+L3
C
      CALL DQXGE(F,A,B,EPSABS,EPSREL,KEY,LIMIT,RESULT,ABSERR,
     *  IER,WORK(1),WORK(L1),WORK(L2),WORK(L3),IWORK,LAST,
     *  WORK(L4),WORK(L5),IWORK(L1),IWORK(L2))
C
C         CALL ERROR HANDLER IF NECESSARY.
C
      LVL = 0
10    IF(IER.EQ.6) LVL = 1
C
C THE CALL TO XERROR HAS BEEN COMMENTED TO GET A STAND ALONE ROUTINE
C IT CAN BE ADAPTED TO FIT THE LOCAL ERROR HANDLING PROCEDURES
C
C     IF(IER.NE.0)CALL XERROR('ABNORMAL RETURN FROM DQXG   ',28,IER,LVL)
      RETURN
      END
*DECK DQXGE
      SUBROUTINE DQXGE(F,A,B,EPSABS,EPSREL,KEY,LIMIT,RESULT,ABSERR,
     *   IER,ALIST,BLIST,RLIST,ELIST,IORD,LAST,VALP,VALN,LP,LN)
C
C
C            THE ROUTINE CALCULATES AN APPROXIMATION RESULT TO A GIVEN
C            DEFINITE INTEGRAL   I = INTEGRAL OF F OVER (A,B),
C            HOPEFULLY SATISFYING FOLLOWING CLAIM FOR ACCURACY
C            ABS(I-RESLT).LE.MAX(EPSABS,EPSREL*ABS(I)).
C
C        PARAMETERS
C         ON ENTRY
C            F      - DOUBLE PRECISION
C                     FUNCTION SUBPROGRAM DEFINING THE INTEGRAND
C                     FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE
C                     DECLARED E X T E R N A L IN THE DRIVER PROGRAM.
C
C            A      - DOUBLE PRECISION
C                     LOWER LIMIT OF INTEGRATION
C
C            B      - DOUBLE PRECISION
C                     UPPER LIMIT OF INTEGRATION
C
C            EPSABS - DOUBLE PRECISION
C                     ABSOLUTE ACCURACY REQUESTED
C
C            EPSREL - DOUBLE PRECISION
C                     RELATIVE ACCURACY REQUESTED
C                     IF  EPSABS.LE.0
C                     AND EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28),
C                     THE ROUTINE WILL END WITH IER = 6.
C
C            KEY    - INTEGER
C                     KEY FOR CHOICE OF LOCAL INTEGRATION RULE
C                     RMS FORMULAS ARE USED WITH
C                      13 - 19               POINTS IF KEY.LT.1,
C                      13 - 19 - (27)        POINTS IF KEY = 1,
C                      13 - 19 - (27) - (41) POINTS IF KEY = 2,
C                           19 -  27  - (41) POINTS IF KEY = 3,
C                                 27  -  41  POINTS IF KEY.GT.3.
C
C                         (RULES) USED IF THE FUNCTION APPEARS
C                         ENOUGH REGULAR IN THE SUBINTERVAL
C
C            LIMIT  - INTEGER
C                     GIVES AN UPPERBOUND ON THE NUMBER OF SUBINTERVALS
C                     IN THE PARTITION OF (A,B), LIMIT.GE.1.
C
C         ON RETURN
C            RESULT - DOUBLE PRECISION
C                     APPROXIMATION TO THE INTEGRAL
C
C            ABSERR - DOUBLE PRECISION
C                     ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR,
C                     WHICH SHOULD EQUAL OR EXCEED ABS(I-RESULT)
C
C            IER    - INTEGER
C                     IER = 0 NORMAL AND RELIABLE TERMINATION OF THE
C                             ROUTINE. IT IS ASSUMED THAT THE REQUESTED
C                             ACCURACY HAS BEEN ACHIEVED.
C                     IER.GT.0 ABNORMAL TERMINATION OF THE ROUTINE
C                             THE ESTIMATES FOR RESULT AND ERROR ARE
C                             LESS RELIABLE. IT IS ASSUMED THAT THE
C                             REQUESTED ACCURACY HAS NOT BEEN ACHIEVED.
C            ERROR MESSAGES
C                     IER = 1 MAXIMUM NUMBER OF SUBDIVISIONS ALLOWED
C                             HAS BEEN ACHIEVED. ONE CAN ALLOW MORE
C                             SUBDIVISIONS BY INCREASING THE VALUE
C                             OF LIMIT.
C                             HOWEVER, IF THIS YIELDS NO IMPROVEMENT IT
C                             IS RATHER ADVISED TO ANALYZE THE INTEGRAND
C                             IN ORDER TO DETERMINE THE INTEGRATION
C                             DIFFICULTIES. IF THE POSITION OF A LOCAL
C                             DIFFICULTY CAN BE DETERMINED(E.G.
C                             SINGULARITY, DISCONTINUITY WITHIN THE
C                             INTERVAL) ONE WILL PROBABLY GAIN FROM
C                             SPLITTING UP THE INTERVAL AT THIS POINT
C                             AND CALLING THE INTEGRATOR ON THE
C                             SUBRANGES. IF POSSIBLE, AN APPROPRIATE
C                             SPECIAL-PURPOSE INTEGRATOR SHOULD BE USED
C                             WHICH IS DESIGNED FOR HANDLING THE TYPE OF
C                             DIFFICULTY INVOLVED.
C                         = 2 THE OCCURRENCE OF ROUNDOFF ERROR IS
C                             DETECTED, WHICH PREVENTS THE REQUESTED
C                             TOLERANCE FROM BEING ACHIEVED.
C                         = 3 EXTREMELY BAD INTEGRAND BEHAVIOUR OCCURS
C                             AT SOME POINTS OF THE INTEGRATION
C                             INTERVAL.
C                         = 6 THE INPUT IS INVALID, BECAUSE
C                             (EPSABS.LE.0 AND
C                              EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28),
C                             RESULT, ABSERR, LAST, RLIST(1) ,
C                             ELIST(1) AND IORD(1) ARE SET TO ZERO.
C                             ALIST(1) AND BLIST(1) ARE SET TO A AND B
C                             RESPECTIVELY.
C
C            ALIST   - DOUBLE PRECISION
C                      VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
C                       LAST  ELEMENTS OF WHICH ARE THE LEFT
C                      END POINTS OF THE SUBINTERVALS IN THE PARTITION
C                      OF THE GIVEN INTEGRATION RANGE (A,B)
C
C            BLIST   - DOUBLE PRECISION
C                      VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
C                       LAST  ELEMENTS OF WHICH ARE THE RIGHT
C                      END POINTS OF THE SUBINTERVALS IN THE PARTITION
C                      OF THE GIVEN INTEGRATION RANGE (A,B)
C
C            RLIST   - DOUBLE PRECISION
C                      VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
C                       LAST  ELEMENTS OF WHICH ARE THE
C                      INTEGRAL APPROXIMATIONS ON THE SUBINTERVALS
C
C            ELIST   - DOUBLE PRECISION
C                      VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
C                       LAST  ELEMENTS OF WHICH ARE THE MODULI OF THE
C                      ABSOLUTE ERROR ESTIMATES ON THE SUBINTERVALS
C
C            IORD    - INTEGER
C                      VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST K
C                      ELEMENTS OF WHICH ARE POINTERS TO THE
C                      ERROR ESTIMATES OVER THE SUBINTERVALS,
C                      SUCH THAT ELIST(IORD(1)), ...,
C                      ELIST(IORD(K)) FORM A DECREASING SEQUENCE,
C                      WITH K = LAST IF LAST.LE.(LIMIT/2+2), AND
C                      K = LIMIT+1-LAST OTHERWISE
C
C            LAST    - INTEGER
C                      NUMBER OF SUBINTERVALS ACTUALLY PRODUCED IN THE
C                      SUBDIVISION PROCESS
C
C            VALP   - DOUBLE PRECISION
C            VALN     ARRAYS OF DIMENSION AT LEAST (21,LIMIT) USED TO
C                     SAVE THE FUNCTIONAL VALUES
C
C            LP     - INTEGER
C            LN       VECTORS OF DIMENSION AT LEAST LIMIT, USED TO
C                     STORE THE ACTUAL NUMBER OF FUNCTIONAL VALUES
C                     SAVED IN THE CORRESPONDING COLUMN
C                     OF VALP,VALN
C
C***ROUTINES CALLED  F,D1MACH,DQXLQM,DQXRRD,DQPSRT,DQXCPY
C
      DOUBLE PRECISION A,ABSERR,ALIST,AREA,AREA1,AREA12,AREA2,A1,A2,B,
     *  BLIST,B1,B2,DABS,DEFABS,DEFAB1,DEFAB2,DMAX1,D1MACH,ELIST,EPMACH,
     *  EPSABS,EPSREL,ERRBND,ERRMAX,ERROR1,ERROR2,ERRO12,ERRSUM,F,
     *  RESABS,RESULT,RLIST,UFLOW,VALP,VALN,VP1,VP2,VN1,VN2
      INTEGER KEY,IER,IORD,IROFF1,IROFF2,K,LAST,LIMIT,MAXERR,
     *  NRMAX,LP,LN,LP1,LP2,LN1,LN2
C
      DIMENSION ALIST(LIMIT),BLIST(LIMIT),ELIST(LIMIT),IORD(LIMIT),
     *  RLIST(LIMIT),VALP(21,LIMIT),VALN(21,LIMIT),LP(LIMIT),LN(LIMIT),
     * VP1(21),VP2(21),VN1(21),VN2(21)
C
      EXTERNAL F
C
C            MACHINE DEPENDENT CONSTANTS
C            ---------------------------
C
C          EPMACH IS THE LARGEST RELATIVE SPACING.
C          UFLOW IS THE SMALLEST POSITIVE MAGNITUDE.
C
C***FIRST EXECUTABLE STATEMENT  DQXGE
      EPMACH = D1MACH(4)
      UFLOW = D1MACH(1)
C
C           TEST ON VALIDITY OF PARAMETERS
C           ------------------------------
C
      IER = 0
      LAST = 0
      RESULT = 0.0D+00
      ABSERR = 0.0D+00
      ALIST(1) = A
      BLIST(1) = B
      RLIST(1) = 0.0D+00
      ELIST(1) = 0.0D+00
      IORD(1) = 0
      IF(EPSABS.LE.0.0D+00.AND.
     *  EPSREL.LT.DMAX1(0.5D+02*EPMACH,0.5D-28)) IER = 6
      IF(IER.EQ.6) GO TO 999
C
C           FIRST APPROXIMATION TO THE INTEGRAL
C           -----------------------------------
C
      LP(1)=1
      LN(1)=1
      VALP(1,1)=F((A+B)*0.5D0)
      VALN(1,1)=VALP(1,1)
      CALL DQXLQM(F,A,B,RESULT,ABSERR,DEFABS,RESABS,
     *           VALP(1,1),VALN(1,1),LP(1),LN(1),KEY)
      LAST = 1
      RLIST(1) = RESULT
      ELIST(1) = ABSERR
      IORD(1) = 1
C
C           TEST ON ACCURACY.
C
      ERRBND = DMAX1(EPSABS,EPSREL*DABS(RESULT))
      IF(ABSERR.LE.0.5D+02*EPMACH*DEFABS.AND.ABSERR.GT.ERRBND) IER = 2
      IF(LIMIT.EQ.1) IER = 1
      IF(IER.NE.0.OR.(ABSERR.LE.ERRBND.AND.ABSERR.NE.RESABS)
     *  .OR.ABSERR.EQ.0.0D+00) GO TO 999
C
C           INITIALIZATION
C           --------------
C
C
      ERRMAX = ABSERR
      MAXERR = 1
      AREA = RESULT
      ERRSUM = ABSERR
      NRMAX = 1
      IROFF1 = 0
      IROFF2 = 0
C
C           MAIN DO-LOOP
C           ------------
C
      DO 30 LAST = 2,LIMIT
C
C           BISECT THE SUBINTERVAL WITH THE LARGEST ERROR ESTIMATE.
C
        A1 = ALIST(MAXERR)
        B1 = 0.5D+00*(ALIST(MAXERR)+BLIST(MAXERR))
        A2 = B1
        B2 = BLIST(MAXERR)
        CALL DQXRRD(F,VALN(1,MAXERR),LN(MAXERR),B1,A1,VN1,VP1,LN1,LP1)
        CALL DQXLQM(F,A1,B1,AREA1,ERROR1,RESABS,DEFAB1,VP1,VN1,LP1,LN1,
     *              KEY)
        CALL DQXRRD(F,VALP(1,MAXERR),LP(MAXERR),A2,B2,VP2,VN2,LP2,LN2)
        CALL DQXLQM(F,A2,B2,AREA2,ERROR2,RESABS,DEFAB2,VP2,VN2,LP2,LN2,
     *              KEY)
C
C           IMPROVE PREVIOUS APPROXIMATIONS TO INTEGRAL
C           AND ERROR AND TEST FOR ACCURACY.
C
        AREA12 = AREA1+AREA2
        ERRO12 = ERROR1+ERROR2
        ERRSUM = ERRSUM+ERRO12-ERRMAX
        AREA = AREA+AREA12-RLIST(MAXERR)
        IF(DEFAB1.EQ.ERROR1.OR.DEFAB2.EQ.ERROR2) GO TO 5
        IF(DABS(RLIST(MAXERR)-AREA12).LE.0.1D-04*DABS(AREA12)
     *  .AND.ERRO12.GE.0.99D+00*ERRMAX) IROFF1 = IROFF1+1
        IF(LAST.GT.10.AND.ERRO12.GT.ERRMAX) IROFF2 = IROFF2+1
    5   RLIST(MAXERR) = AREA1
        RLIST(LAST) = AREA2
        ERRBND = DMAX1(EPSABS,EPSREL*DABS(AREA))
        IF(ERRSUM.LE.ERRBND) GO TO 8
C
C           TEST FOR ROUNDOFF ERROR AND EVENTUALLY SET ERROR FLAG.
C
        IF(IROFF1.GE.6.OR.IROFF2.GE.20) IER = 2
C
C           SET ERROR FLAG IN THE CASE THAT THE NUMBER OF SUBINTERVALS
C           EQUALS LIMIT.
C
        IF(LAST.EQ.LIMIT) IER = 1
C
C           SET ERROR FLAG IN THE CASE OF BAD INTEGRAND BEHAVIOUR
C           AT A POINT OF THE INTEGRATION RANGE.
C
        IF(DMAX1(DABS(A1),DABS(B2)).LE.(0.1D+01+0.1D+03*
     *  EPMACH)*(DABS(A2)+0.1D+04*UFLOW)) IER = 3
C
C           APPEND THE NEWLY-CREATED INTERVALS TO THE LIST.
C
    8   IF(ERROR2.GT.ERROR1) GO TO 10
        ALIST(LAST) = A2
        BLIST(MAXERR) = B1
        BLIST(LAST) = B2
        ELIST(MAXERR) = ERROR1
        ELIST(LAST) = ERROR2
        CALL DQXCPY(VALP(1,MAXERR),VP1,LP1)
        LP(MAXERR)=LP1
        CALL DQXCPY(VALN(1,MAXERR),VN1,LN1)
        LN(MAXERR)=LN1
        CALL DQXCPY(VALP(1,LAST),VP2,LP2)
        LP(LAST)=LP2
        CALL DQXCPY(VALN(1,LAST),VN2,LN2)
        LN(LAST)=LN2
        GO TO 20
   10   ALIST(MAXERR) = A2
        ALIST(LAST) = A1
        BLIST(LAST) = B1
        RLIST(MAXERR) = AREA2
        RLIST(LAST) = AREA1
        ELIST(MAXERR) = ERROR2
        ELIST(LAST) = ERROR1
        CALL DQXCPY(VALP(1,MAXERR),VP2,LP2)
        LP(MAXERR)=LP2
        CALL DQXCPY(VALN(1,MAXERR),VN2,LN2)
        LN(MAXERR)=LN2
        CALL DQXCPY(VALP(1,LAST),VP1,LP1)
        LP(LAST)=LP1
        CALL DQXCPY(VALN(1,LAST),VN1,LN1)
        LN(LAST)=LN1
C
C           CALL SUBROUTINE DQPSRT TO MAINTAIN THE DESCENDING ORDERING
C           IN THE LIST OF ERROR ESTIMATES AND SELECT THE SUBINTERVAL
C           WITH THE LARGEST ERROR ESTIMATE (TO BE BISECTED NEXT).
C
   20   CALL DQPSRT(LIMIT,LAST,MAXERR,ERRMAX,ELIST,IORD,NRMAX)
C ***JUMP OUT OF DO-LOOP
        IF(IER.NE.0.OR.ERRSUM.LE.ERRBND) GO TO 40
   30 CONTINUE
C
C           COMPUTE FINAL RESULT.
C           ---------------------
C
   40 RESULT = 0.0D+00
      DO 50 K=1,LAST
        RESULT = RESULT+RLIST(K)
   50 CONTINUE
      ABSERR = ERRSUM
  999 RETURN
      END
*DECK DQXLQM
      SUBROUTINE DQXLQM(F,A,B,RESULT,ABSERR,RESABS,RESASC,VR,VS,LR,LS,
     *                  KEY)
C
C            TO COMPUTE I = INTEGRAL OF F OVER (A,B), WITH ERROR
C                           ESTIMATE
C                       J = INTEGRAL OF ABS(F) OVER (A,B)
C
C           PARAMETERS
C            ON ENTRY
C              F      - DOUBLE PRECISION
C                       FUNCTION SUBPROGRAM DEFINING THE INTEGRAND
C                       FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE
C                       DECLARED E X T E R N A L IN THE DRIVER PROGRAM.
C
C              A      - DOUBLE PRECISION
C                       LOWER LIMIT OF INTEGRATION
C
C              B      - DOUBLE PRECISION
C                       UPPER LIMIT OF INTEGRATION
C
C              VR     - DOUBLE PRECISION
C                       VECTOR OF LENGTH LR CONTAINING THE
C                       SAVED  FUNCTIONAL VALUES OF POSITIVE ABSCISSAS
C
C              VS     - DOUBLE PRECISION
C                       VECTOR OF LENGTH LS CONTAINING THE
C                       SAVED  FUNCTIONAL VALUES OF NEGATIVE ABSCISSAS
C
C              LR     - INTEGER
C              LS       NUMBER OF ELEMENTS IN
C                       VR,VS RESPECTIVELY
C
C            KEY    - INTEGER
C                     KEY FOR CHOICE OF LOCAL INTEGRATION RULE
C                     RMS FORMULAS ARE USED WITH
C                      13 - 19               POINTS IF KEY.LT.1,
C                      13 - 19 - (27)        POINTS IF KEY = 1,
C                      13 - 19 - (27) - (41) POINTS IF KEY = 2,
C                           19 -  27  - (41) POINTS IF KEY = 3,
C                                 27  -  41  POINTS IF KEY.GT.3.
C
C                         (RULES) USED IF THE FUNCTION APPEARS
C                         ENOUGH REGULAR
C
C            ON RETURN
C              RESULT - DOUBLE PRECISION
C                       APPROXIMATION TO THE INTEGRAL I
C
C              ABSERR - DOUBLE PRECISION
C                       ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR,
C                       WHICH SHOULD NOT EXCEED ABS(I-RESULT)
C
C              RESABS - DOUBLE PRECISION
C                       APPROXIMATION TO THE INTEGRAL J
C
C              RESASC - DOUBLE PRECISION
C                       APPROXIMATION TO THE INTEGRAL OF ABS(F-I/(B-A))
C                       OVER (A,B)
C
C              VR     - DOUBLE PRECISION
C                       VECTOR OF LENGTH LR CONTAINING THE
C                       SAVED  FUNCTIONAL VALUES OF POSITIVE ABSCISSAS
C
C              VS     - DOUBLE PRECISION
C                       VECTOR OF LENGTH LS CONTAINING THE
C                       SAVED  FUNCTIONAL VALUES OF NEGATIVE ABSCISSAS
C
C              LR     - INTEGER
C              LS       NUMBER OF ELEMENTS IN
C                       VR,VS RESPECTIVELY
C
C***ROUTINES CALLED  D1MACH,DQXRUL
C
      DOUBLE PRECISION F,A,B,RESULT,ABSERR,RESABS,RESASC,
     *          D1MACH,EPMACH,RESG,RESK,UFLOW,ERROLD,VR(21),VS(21)
      INTEGER K,K0,K1,K2,KEY,KEY1,LR,LS
      EXTERNAL F
C
C            MACHINE DEPENDENT CONSTANTS
C            ---------------------------
C
C          EPMACH IS THE LARGEST RELATIVE SPACING.
C          UFLOW IS THE SMALLEST POSITIVE MAGNITUDE.
C          ERROLD IS THE LARGEST MAGNITUDE.
C
C***FIRST EXECUTABLE STATEMENT DQXLQM
      EPMACH = D1MACH(4)
      UFLOW  = D1MACH(1)
      ERROLD = D1MACH(2)
C
      KEY1 = MAX(KEY ,  0)
      KEY1 = MIN(KEY1,  4)
      K0   = MAX(KEY1-2,0)
      K1   = K0+1
      K2   = MIN(KEY1+1,3)
C
      CALL DQXRUL(F,A,B,RESG,RESABS,RESASC,K0,K1,VR,VS,LR,LS)
      DO 99 K=K1,K2
        CALL DQXRUL(F,A,B,RESK,RESABS,RESASC,K,K1,VR,VS,LR,LS)
        RESULT=RESK
        ABSERR = DABS((RESK-RESG))
        IF(RESASC.NE.0.0D+00.AND.ABSERR.NE.0.0D+00)
     *  ABSERR = RESASC*DMIN1(1.D0,(200D0*ABSERR/RESASC)**1.5D+00)
        IF(RESABS.GT.UFLOW/(10D0*EPMACH)) ABSERR = DMAX1
     *  ((EPMACH*10D0)*RESABS,ABSERR)
        RESG=RESK
        IF(ABSERR.GT.ERROLD*0.16)GOTO 3000
        IF(ABSERR.LT.1000*EPMACH*RESABS)GOTO 3000
        ERROLD=ABSERR
99    CONTINUE
3000  CONTINUE
      RETURN
      END
*DECK DQXRUL
      SUBROUTINE DQXRUL(F,XL,XU,Y,YA,YM,KE,K1,FV1,FV2,L1,L2)
C
C            TO COMPUTE I = INTEGRAL OF F OVER (A,B), WITH ERROR
C                           ESTIMATE
C            AND CONDITIONALLY COMPUTE
C                       J = INTEGRAL OF ABS(F) OVER (A,B)
C                       BY USING AN  RMS RULE
C           PARAMETERS
C            ON ENTRY
C              F      - DOUBLE PRECISION
C                       FUNCTION SUBPROGRAM DEFINING THE INTEGRAND
C                       FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE
C                       DECLARED E X T E R N A L IN THE DRIVER PROGRAM.
C
C              XL     - DOUBLE PRECISION
C                       LOWER LIMIT OF INTEGRATION
C
C              XU     - DOUBLE PRECISION
C                       UPPER LIMIT OF INTEGRATION
C
C
C              KE     - INTEGER
C                     KEY FOR CHOICE OF LOCAL INTEGRATION RULE
C                     AN RMS RULE IS USED WITH
C                         13      POINTS IF KE  = 2,
C                         19      POINTS IF KE  = 3,
C                         27      POINTS IF KE  = 4,
C                         42      POINTS IF KE  = 5
C
C              K1     INTEGER
C                     VALUE OF KEY FOR WHICH THE ADDITIONAL ESTIMATES
C                     YA, YM ARE TO BE COMPUTED
C
C              FV1    - DOUBLE PRECISION
C                       VECTOR CONTAINING L1
C                       SAVED  FUNCTIONAL VALUES OF POSITIVE ABSCISSAS
C
C              FV2    - DOUBLE PRECISION
C                       VECTOR CONTAINING L2
C                       SAVED  FUNCTIONAL VALUES OF NEGATIVE ABSCISSAS
C
C              L1     - INTEGER
C              L2       NUMBER OF ELEMENTS IN FV1,FV2  RESPECTIVELY
C
C            ON RETURN
C              Y      - DOUBLE PRECISION
C                       APPROXIMATION TO THE INTEGRAL I
C                       RESULT IS COMPUTED BY APPLYING THE
C                       REQUESTED RMS RULE
C
C              YA     - DOUBLE PRECISION
C                       IF KEY = K1  APPROXIMATION TO THE INTEGRAL J
C                       ELSE UNCHANGED
C
C              YM     - DOUBLE PRECISION
C                       IF KEY = K1  APPROXIMATION TO THE INTEGRAL OF
C                                      ABS(F-I/(XU-XL)   OVER (XL,XU)
C                       ELSE UNCHANGED
C
C              FV1    - DOUBLE PRECISION
C                       VECTOR L1 CONTAINING L1
C                       SAVED  FUNCTIONAL VALUES OF POSITIVE ABSCISSAS
C
C              FV2    - DOUBLE PRECISION
C                       VECTOR CONTAINING L2
C                       SAVED  FUNCTIONAL VALUES OF NEGATIVE ABSCISSAS
C
C              L1     - INTEGER
C              L2       NUMBER OF ELEMENTS IN FV1,FV2  RESPECTIVELY
C
C***ROUTINES CALLED  F
C
      DOUBLE PRECISION F,XL,XU,LDL,Y,YA,YM,Y2,XX(41),WW(52),
     *                 FV1(21),FV2(21),AA,BB,C
      EXTERNAL F
      INTEGER ISTART(4),LEN(4),KE,K1,L1,L2
      SAVE ISTART,LEN,XX,WW
      DATA ISTART/0, 7, 17, 31/,LEN/7, 10, 14, 21/
      DATA XX(  1)/.0                       /
      DATA XX(  2)/.25000000000000000000D+00/
      DATA XX(  3)/.50000000000000000000D+00/
      DATA XX(  4)/.75000000000000000000D+00/
      DATA XX(  5)/.87500000000000000000D+00/
      DATA XX(  6)/.93750000000000000000D+00/
      DATA XX(  7)/.10000000000000000000D+01/
      DATA XX(  8)/.37500000000000000000D+00/
      DATA XX(  9)/.62500000000000000000D+00/
      DATA XX( 10)/.96875000000000000000D+00/
      DATA XX( 11)/.12500000000000000000D+00/
      DATA XX( 12)/.68750000000000000000D+00/
      DATA XX( 13)/.81250000000000000000D+00/
      DATA XX( 14)/.98437500000000000000D+00/
      DATA XX( 15)/.18750000000000000000D+00/
      DATA XX( 16)/.31250000000000000000D+00/
      DATA XX( 17)/.43750000000000000000D+00/
      DATA XX( 18)/.56250000000000000000D+00/
      DATA XX( 19)/.84375000000000000000D+00/
      DATA XX( 20)/.90625000000000000000D+00/
      DATA XX( 21)/.99218750000000000000D+00/
C   NUMBER OF NODES 13
      DATA WW(1)/1.303262173284849021810473057638590518409112513421D-1/
      DATA WW(2)/2.390632866847646220320329836544615917290026806242D-1/
      DATA WW(3)/2.630626354774670227333506083741355715758124943143D-1/
      DATA WW(4)/2.186819313830574175167853094864355208948886875898D-1/
      DATA WW(5)/2.757897646642836865859601197607471574336674206700D-2/
      DATA WW(6)/1.055750100538458443365034879086669791305550493830D-1/
      DATA WW(7)/1.571194260595182254168429283636656908546309467968D-2/
C   NUMBER OF NODES 19
      DATA WW(8)/1.298751627936015783241173611320651866834051160074D-1/
      DATA WW(9)/2.249996826462523640447834514709508786970828213187D-1/
      DATA WW(15)/5.542699233295875168406783695143646338274805359780D-2/
      DATA WW(10)/1.680415725925575286319046726692683040162290325505D-1/
      DATA WW(16)/9.986735247403367525720377847755415293097913496236D-2/
      DATA WW(11)/1.415567675701225879892811622832845252125600939627D-1/
      DATA WW(12)/1.006482260551160175038684459742336605269707889822D-1/
      DATA WW(13)/2.510604860724282479058338820428989444699235030871D-2/
      DATA WW(17)/4.507523056810492466415880450799432587809828791196D-2/
      DATA WW(14)/9.402964360009747110031098328922608224934320397592D-3/
C   NUMBER OF NODES 27
      DATA WW(18)/6.300942249647773931746170540321811473310938661469D-2/
      DATA WW(28)/1.239572396231834242194189674243818619042280816640D-1/
      DATA WW(19)/1.261383225537664703012999637242003647020326905948D-1/
      DATA WW(25)/1.235837891364555000245004813294817451524633100256D-1/
      DATA WW(20)/1.273864433581028272878709981850307363453523117880D-1/
      DATA WW(26)/1.148933497158144016800199601785309838604146040215D-1/
      DATA WW(29)/2.501306413750310579525950767549691151739047969345D-2/
      DATA WW(21)/8.576500414311820514214087864326799153427368592787D-2/
      DATA WW(30)/4.915957918146130094258849161350510503556792927578D-2/
      DATA WW(22)/7.102884842310253397447305465997026228407227220665D-2/
      DATA WW(23)/5.026383572857942403759829860675892897279675661654D-2/
      DATA WW(27)/1.252575774226122633391477702593585307254527198070D-2/
      DATA WW(31)/2.259167374956474713302030584548274729936249753832D-2/
      DATA WW(24)/4.683670010609093810432609684738393586390722052124D-3/
C   NUMBER OF NODES 41
      DATA WW(32)/6.362762978782724559269342300509058175967124446839D-2/
      DATA WW(42)/1.187141856692283347609436153545356484256869129472D-1/
      DATA WW(46)/1.533126874056586959338368742803997744815413565014D-2/
      DATA WW(33)/9.950065827346794643193261975720606296171462239514D-2/
      DATA WW(47)/3.527159369750123100455704702965541866345781113903D-2/
      DATA WW(39)/8.140326425945938045967829319725797511040878579808D-2/
      DATA WW(48)/5.000556431653955124212795201196389006184693561679D-2/
      DATA WW(34)/7.048220002718565366098742295389607994441704889441D-2/
      DATA WW(49)/5.744164831179720106340717579281831675999717767532D-2/
      DATA WW(40)/6.583213447600552906273539578430361199084485578379D-2/
      DATA WW(43)/5.999947605385971985589674757013565610751028128731D-2/
      DATA WW(35)/6.512297339398335645872697307762912795346716454337D-2/
      DATA WW(44)/5.500937980198041736910257988346101839062581489820D-2/
      DATA WW(50)/1.598823797283813438301248206397233634639162043386D-2/
      DATA WW(36)/3.998229150313659724790527138690215186863915308702D-2/
      DATA WW(51)/2.635660410220884993472478832884065450876913559421D-2/
      DATA WW(37)/3.456512257080287509832054272964315588028252136044D-2/
      DATA WW(41)/2.592913726450792546064232192976262988065252032902D-2/
      DATA WW(45)/5.264422421764655969760271538981443718440340270116D-3/
      DATA WW(52)/1.196003937945541091670106760660561117114584656319D-2/
      DATA WW(38)/2.212167975884114432760321569298651047876071264944D-3/
C
C***FIRST EXECUTABLE STATEMENT DQXRUL
      K=KE+1
      IS=ISTART(K)
      KS=LEN(K)
      LDL= XU-XL
      BB= LDL*0.5D0
      AA= XL+BB
      Y =0.D0
      DO 100 I=1,KS
         IF(I.GT.L1.OR.I.GT.L2)   C=BB*XX(I)
         IF(I.GT.L1)              FV1(I)=F(AA+C)
         IF(I.GT.L2)              FV2(I)=F(AA-C)
100      Y=Y+(FV1(I)+FV2(I))*WW(IS+I)
      Y2=Y
      Y=Y*BB
      IF(L1.LT.KS)L1=KS
      IF(L2.LT.KS)L2=KS
      IF(KE.NE.K1)RETURN
      YA=0.D0
      DO 25 I=1,KS
  25    YA=YA+(DABS(FV1(I))+DABS(FV2(I)))*WW(IS+I)
      Y2=Y2*0.5D0
      YM=0.D0
      YA=YA*DABS(BB)
      DO 27 I=1,KS
   27 YM=YM+(DABS(FV1(I)-Y2)+DABS(FV2(I)-Y2))*WW(IS+I)
      YM=YM*DABS(BB)
      RETURN
      END
*DECK DQXRRD
      SUBROUTINE DQXRRD(F,Z,LZ,XL,XU,R,S,LR,LS)
C
C            TO REORDER THE COMPUTED FUNCTIONAL VALUES BEFORE
C            THE BISECTION OF AN INTERVAL
C           PARAMETERS
C            ON ENTRY
C              F      - DOUBLE PRECISION
C                       FUNCTION SUBPROGRAM DEFINING THE INTEGRAND
C                       FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE
C                       DECLARED E X T E R N A L IN THE DRIVER PROGRAM.
C
C              XL     - DOUBLE PRECISION
C                       LOWER LIMIT OF INTERVAL
C
C              XU     - DOUBLE PRECISION
C                       UPPER LIMIT OF INTERVAL
C
C              Z      - DOUBLE PRECISION
C                       VECTOR CONTAINING LZ
C                       SAVED  FUNCTIONAL VALUES
C
C              LZ     - INTEGER
C                       NUMBER OF ELEMENTS IN LZ
C
C
C            ON RETURN
C              R      - DOUBLE PRECISION
C              S        VECTORS CONTAINING LR, LS
C                       SAVED  FUNCTIONAL VALUES FOR THE NEW INTERVALS
C
C              LR     - INTEGER
C              LS       NUMBER OF ELEMENTES IN R,S RESPECTIVELY
C
C***ROUTINES CALLED  F
C
      DOUBLE PRECISION F,R,S,Z,XU,XL,DLEN,CENTR
      INTEGER LR,LS,LZ
      DIMENSION R(21),S(21),Z(21)
C***FIRST EXECUTABLE STATEMENT DQXRRD
      DLEN= (XU-XL)*0.5D0
      CENTR= XL+DLEN
      R(1)=  Z(3)
      R(2)=  Z(9)
      R(3)=  Z(4)
      R(4)=  Z(5)
      R(5)=  Z(6)
      R(6)=  Z(10)
      R(7)=  Z(7)
      S(1)=  Z(3)
      S(2)=  Z(8)
      S(3)=  Z(2)
      S(7)=  Z(1)
      IF(LZ.GT.11)GOTO 2
      R(8)=  F(CENTR+DLEN*0.37500000D0)
      R(9)=  F(CENTR+DLEN*0.62500000D0)
      R(10)=  F(CENTR+DLEN*0.96875000D0)
      LR=  10
      IF(LZ.NE.11)S(4)=  F(CENTR-DLEN*0.75000000D0)
      IF(LZ.EQ.11)S(4)=  Z(11)
      S(5)=  F(CENTR-DLEN*0.87500000D0)
      S(6)=  F(CENTR-DLEN*0.93750000D0)
      S(8)=  F(CENTR-DLEN*0.37500000D0)
      S(9)=  F(CENTR-DLEN*0.62500000D0)
      S(10)=  F(CENTR-DLEN*0.96875000D0)
      LS=  10
      GOTO 3000
2     R(8)= Z(12)
      R(9)= Z(13)
      R(10)= Z(14)
      LR=  10
      S(4)= Z(11)
      S(5)= F(CENTR-DLEN*0.87500000D0)
      S(6)= F(CENTR-DLEN*0.93750000D0)
      IF(LZ.GT.14)GOTO3
      S(8)= F(CENTR-DLEN*0.37500000D0)
      S(9)= F(CENTR-DLEN*0.62500000D0)
      S(10)= F(CENTR-DLEN*0.96875000D0)
      LS=  10
      GOTO 3000
3     R(11)= Z(18)
      R(12)= Z(19)
      R(13)= Z(20)
      R(14)= Z(21)
      LR=  14
      S(8)= Z(16)
      S(9)= Z(15)
      S(10)= F(CENTR-DLEN*0.96875000D0)
      S(11)= Z(17)
      LS=  11
3000  RETURN
      END
*DECK DQXCPY
      SUBROUTINE DQXCPY(A,B,L)
C
C  TO COPY THE DOUBLE PRECISION VECTOR B OF LENGTH L   I N T O
C          THE DOUBLE PRECISION VECTOR A OF LENGTH L
C
C***REMARK  THIS ROUTINE CAN BE IMPROVED, BY CODING IT IN THE
C           ASSEMBLER LANGUAGE OF THE USED MACHINE
C***ROUTINES CALLED  (NONE)
C
      INTEGER L
      DOUBLE PRECISION A(L),B(L)
C***FIRST EXECUTABLE STATEMENT DQXCPY
      DO 1 I=1,L
1     A(I)=B(I)
      RETURN
      END
*DECK DTEST  *** TEST PROGRAM ***
C  THIS DRIVER CALLS VARIOUS INTEGRATORS WITH  VARIOUS RELATIVE
C  TOLERACES TO INTEGRATE THE FUNCTION
C     F(X) = DSQRT(X/(1-X)) * DLOG(X)   OVER (0, 1)
C
C THE EXACT INTEGRAL IS                           -0.606789763508705...
C
      DOUBLE PRECISION F,A,B,ATOL,RTOL,TEE,Y,WORK(4600)
      INTEGER NCALLS,IER,NEVAL,LAST,IWORK(300),KKK,KEY,LIMIT,LENW,LENIW
      EXTERNAL F
      COMMON /CTEST/NCALLS
      DATA LIMIT/100/,LENIW/300/,LENW/4600/
      WRITE(6,600)
600   FORMAT(32H1 COMPUTATION OF THE INTEGRAL OF,
     *      /38H       F(X) = DSQRT(X/(1-X)) * DLOG(X),
     *      /13H  OVER (0, 1),
     *      /23H  THE EXACT INTEGRAL IS,30X,
     *       22H -0.606789763508705...)
      A=0.0D0
      B=1.0D0
      WRITE(6,61)
 61   FORMAT(//22H TEST OF DQXGS               //
     *  53H RELATIVE TOLERANCE FUNCTION CALLS ERROR ESTIMATE IER,
     *  13H       RESULT )
      DO 1 I=3,13,2
         ATOL=1D-20
         RTOL=10D0**(-I)
         NCALLS=0
         CALL DQXGS (F,A,B,ATOL,RTOL,Y,TEE,IER,LIMIT,
     *               LENIW,LENW,LAST,IWORK,WORK)
         WRITE(6,60)I,NCALLS,TEE,IER,Y
60       FORMAT(4X,6H10**(-,I2,1H),I17,9X,D7.2,4X,I2,F20.15)
 1    CONTINUE
      DO 2 KKK = 1,5
       KEY = KKK - 1
       WRITE(6,62) KEY
 62    FORMAT(//23H TEST OF DQXG , KEY =    ,I4   //
     *  53H RELATIVE TOLERANCE FUNCTION CALLS ERROR ESTIMATE IER,
     *  13H       RESULT )
       DO 2 I=3,13,2
         ATOL=1D-20
         RTOL=10D0**(-I)
         NCALLS=0
         CALL DQXG  (F,A,B,ATOL,RTOL,KEY,Y,TEE,IER,LIMIT,
     *               LENIW,LENW,LAST,IWORK,WORK)
         WRITE(6,60)I,NCALLS,TEE,IER,Y
 2    CONTINUE
      STOP
      END
*DECK F
      DOUBLE PRECISION FUNCTION F(X)
      DOUBLE PRECISION X
C        INTEGRAND FUNCTION, IN THE INTERVAL (0,1) THE EXACT INTEGRAL
C   TO 15 DIGITS IS
C          -0.60678 97635 08705D0
C
      COMMON /CTEST/NCALLS
      NCALLS=NCALLS+1
      IF (X.EQ.0D0) GO TO  10
      IF (X.EQ.1D0) GO TO  10
      F=DSQRT(X/(1-X))*DLOG(X)
      RETURN
10    F=0.0
      RETURN
      END
CC+------------------------------------------------------------------+CC
CC+                                                                  +CC
CC+                LAWRENCE LIVERMORE NATIONAL LABORATORY            +CC
CC+                                                                  +CC
CC+        MATHEMATICS AND STATISTICS DIVISION SOFTWARE LIBRARY      +CC
CC+                                                                  +CC
CC+------------------------------------------------------------------+CC
CC+                                                                  +CC
CC+  THE SLATEC MATHEMATICAL PROGRAM LIBRARY IS ISSUED BY            +CC
CC+  THE FOLLOWING:                                                  +CC
CC+                                                                  +CC
CC+       AIR FORCE WEAPONS LABORATORY, ALBUQUERQUE                  +CC
CC+       LAWRENCE LIVERMORE NATIONAL LABORATORY, LIVERMORE          +CC
CC+       LOS ALAMOS NATIONAL SCIENTIFIC LABORATORY, LOS ALAMOS      +CC
CC+       NATIONAL BUREAU OF STANDARDS, WASHINGTON                   +CC
CC+       OAK RIDGE NATIONAL LABORATORY, OAK RIDGE                   +CC
CC+       SANDIA NATIONAL LABORATORIES, ALBUQUERQUE                  +CC
CC+       SANDIA NATIONAL LABORATORIES, LIVERMORE                    +CC
CC+                                                                  +CC
CC+  PLEASE REFER QUESTIONS REGARDING THIS SOFTWARE TO THE           +CC
CC+  MSD CONSULTING OFFICE, EXTENSION 3-2976.                        +CC
CC+                                                                  +CC
CC+                                                                  +CC
CC+              * * * * * N O T I C E * * * * *                     +CC
CC+                                                                  +CC
CC+  THIS MATERIAL WAS PREPARED AS AN ACCOUNT OF WORK SPONSORED      +CC
CC+  BY THE UNITED STATES GOVERNMENT.  NEITHER THE UNITED            +CC
CC+  STATES GOVERNMENT, NOR THE DEPARTMENT OF DEFENSE, NOR           +CC
CC+  ANY OF THEIR EMPLOYEES, MAKES ANY WARRANTY, EXPRESSED           +CC
CC+  OR IMPLIED, OR ASSUMES ANY LEGAL LIABILITY OR                   +CC
CC+  RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR               +CC
CC+  USEFULNESS OF ANY INFORMATION, APPARATUS, PRODUCT,              +CC
CC+  OR PROCESS DISCLOSED, OR REPRESENTS ITS USE WOULD NOT           +CC
CC+  INFRINGE UPON PRIVATELY OWNED RIGHTS.                           +CC
CC+                                                                  +CC
CC+------------------------------------------------------------------+CC
*DECK QPSRT
      SUBROUTINE QPSRT(LIMIT,LAST,MAXERR,ERMAX,ELIST,IORD,NRMAX)
C***BEGIN PROLOGUE  QPSRT
C***REFER TO  QAGE,QAGIE,QAGPE,QAGSE,QAWCE,QAWOE,QAWSE
C***ROUTINES CALLED  (NONE)
C***KEYWORDS  SEQUENTIAL SORTING
C***DESCRIPTION
C
C 1.        QPSRT
C           ORDERING ROUTINE
C              STANDARD FORTRAN SUBROUTINE
C              REAL VERSION
C
C 2.        PURPOSE
C              THIS ROUTINE MAINTAINS THE DESCENDING ORDERING
C              IN THE LIST OF THE LOCAL ERROR ESTIMATES RESULTING FROM
C              THE INTERVAL SUBDIVISION PROCESS. AT EACH CALL TWO ERROR
C              ESTIMATES ARE INSERTED USING THE SEQUENTIAL SEARCH
C              METHOD, TOP-DOWN FOR THE LARGEST ERROR ESTIMATE
C              AND BOTTOM-UP FOR THE SMALLEST ERROR ESTIMATE.
C
C 3.        CALLING SEQUENCE
C              CALL QPSRT(LIMIT,LAST,MAXERR,ERMAX,ELIST,IORD,NRMAX)
C
C           PARAMETERS (MEANING AT OUTPUT)
C              LIMIT  - INTEGER
C                       MAXIMUM NUMBER OF ERROR ESTIMATES THE LIST
C                       CAN CONTAIN
C
C              LAST   - INTEGER
C                       NUMBER OF ERROR ESTIMATES CURRENTLY
C                       IN THE LIST
C
C              MAXERR - INTEGER
C                       MAXERR POINTS TO THE NRMAX-TH LARGEST ERROR
C                       ESTIMATE CURRENTLY IN THE LIST
C
C              ERMAX  - REAL
C                       NRMAX-TH LARGEST ERROR ESTIMATE
C                       ERMAX = ELIST(MAXERR)
C
C              ELIST  - REAL
C                       VECTOR OF DIMENSION LAST CONTAINING
C                       THE ERROR ESTIMATES
C
C              IORD   - INTEGER
C                       VECTOR OF DIMENSION LAST, THE FIRST K
C                       ELEMENTS OF WHICH CONTAIN POINTERS
C                       TO THE ERROR ESTIMATES, SUCH THAT
C                       ELIST(IORD(1)),... , ELIST(IORD(K))
C                       FORM A DECREASING SEQUENCE, WITH
C                       K = LAST IF LAST.LE.(LIMIT/2+2), AND
C                       K = LIMIT+1-LAST OTHERWISE
C
C              NRMAX  - INTEGER
C                       MAXERR = IORD(NRMAX)
C
C 4.        NO SUBROUTINES OR FUNCTIONS NEEDED
C***END PROLOGUE  QPSRT
C
      REAL ELIST,ERMAX,ERRMAX,ERRMIN
      INTEGER I,IBEG,IDO,IORD,ISUCC,J,JBND,JUPBN,K,LAST,LIMIT,MAXERR,
     1  NRMAX
      DIMENSION ELIST(LAST),IORD(LAST)
C
C           CHECK WHETHER THE LIST CONTAINS MORE THAN
C           TWO ERROR ESTIMATES.
C
C***FIRST EXECUTABLE STATEMENT  QPSRT
      IF(LAST.GT.2) GO TO 10
      IORD(1) = 1
      IORD(2) = 2
      GO TO 90
C
C           THIS PART OF THE ROUTINE IS ONLY EXECUTED
C           IF, DUE TO A DIFFICULT INTEGRAND, SUBDIVISION
C           INCREASED THE ERROR ESTIMATE. IN THE NORMAL CASE
C           THE INSERT PROCEDURE SHOULD START AFTER THE
C           NRMAX-TH LARGEST ERROR ESTIMATE.
C
   10 ERRMAX = ELIST(MAXERR)
      IF(NRMAX.EQ.1) GO TO 30
      IDO = NRMAX-1
      DO 20 I = 1,IDO
        ISUCC = IORD(NRMAX-1)
C ***JUMP OUT OF DO-LOOP
        IF(ERRMAX.LE.ELIST(ISUCC)) GO TO 30
        IORD(NRMAX) = ISUCC
        NRMAX = NRMAX-1
   20    CONTINUE
C
C           COMPUTE THE NUMBER OF ELEMENTS IN THE LIST TO
C           BE MAINTAINED IN DESCENDING ORDER. THIS NUMBER
C           DEPENDS ON THE NUMBER OF SUBDIVISIONS STILL
C           ALLOWED.
C
   30 JUPBN = LAST
      IF(LAST.GT.(LIMIT/2+2)) JUPBN = LIMIT+3-LAST
      ERRMIN = ELIST(LAST)
C
C           INSERT ERRMAX BY TRAVERSING THE LIST TOP-DOWN,
C           STARTING COMPARISON FROM THE ELEMENT ELIST(IORD(NRMAX+1)).
C
      JBND = JUPBN-1
      IBEG = NRMAX+1
      IF(IBEG.GT.JBND) GO TO 50
      DO 40 I=IBEG,JBND
        ISUCC = IORD(I)
C ***JUMP OUT OF DO-LOOP
        IF(ERRMAX.GE.ELIST(ISUCC)) GO TO 60
        IORD(I-1) = ISUCC
   40 CONTINUE
   50 IORD(JBND) = MAXERR
      IORD(JUPBN) = LAST
      GO TO 90
C
C           INSERT ERRMIN BY TRAVERSING THE LIST BOTTOM-UP.
C
   60 IORD(I-1) = MAXERR
      K = JBND
      DO 70 J=I,JBND
        ISUCC = IORD(K)
C ***JUMP OUT OF DO-LOOP
        IF(ERRMIN.LT.ELIST(ISUCC)) GO TO 80
        IORD(K+1) = ISUCC
        K = K-1
   70 CONTINUE
      IORD(I) = LAST
      GO TO 90
   80 IORD(K+1) = LAST
C
C           SET MAXERR AND ERMAX.
C
   90 MAXERR = IORD(NRMAX)
      ERMAX = ELIST(MAXERR)
      RETURN
      END
*DECK QELG
      SUBROUTINE QELG(N,EPSTAB,RESULT,ABSERR,RES3LA,NRES)
C***BEGIN PROLOGUE  QELG
C***REFER TO  QAGIE,QAGOE,QAGPE,QAGSE
C***ROUTINES CALLED  R1MACH
C***REVISION DATE  830518   (YYMMDD)
C***KEYWORDS  CONVERGENCE ACCELERATION,EPSILON ALGORITHM,EXTRAPOLATION
C***AUTHOR  PIESSENS, ROBERT, APPLIED MATH. AND PROGR. DIV. -
C             K. U. LEUVEN
C           DE DONCKER, ELISE, APPLIED MATH. AND PROGR. DIV. -
C             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           REAL 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 - REAL
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 - REAL
C                       RESULTING APPROXIMATION TO THE INTEGRAL
C
C              ABSERR - REAL
C                       ESTIMATE OF THE ABSOLUTE ERROR COMPUTED FROM
C                       RESULT AND THE 3 PREVIOUS RESULTS
C
C              RES3LA - REAL
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***END PROLOGUE  QELG
C
      REAL ABSERR,DELTA1,DELTA2,DELTA3,R1MACH,
     1  EPMACH,EPSINF,EPSTAB,ERROR,ERR1,ERR2,ERR3,E0,E1,E1ABS,E2,E3,
     2  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
C           E1       COMPUTATION OF A NEW ELEMENT IN
C           E2       THE EPSILON TABLE IS BASED
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  QELG
      EPMACH = R1MACH(4)
      OFLOW = R1MACH(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 = ABS(E1)
        DELTA2 = E2-E1
        ERR2 = ABS(DELTA2)
        TOL2 = AMAX1(ABS(E2),E1ABS)*EPMACH
        DELTA3 = E1-E0
        ERR3 = ABS(DELTA3)
        TOL3 = AMAX1(E1ABS,ABS(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 = ABS(DELTA1)
        TOL1 = AMAX1(E1ABS,ABS(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.1E+01/DELTA1+0.1E+01/DELTA2-0.1E+01/DELTA3
        EPSINF = ABS(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.1E-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.1E+01/SS
        EPSTAB(K1) = RES
        K1 = K1-2
        ERROR = ERR2+ABS(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 = ABS(RESULT-RES3LA(3))+ABS(RESULT-RES3LA(2))
     1  +ABS(RESULT-RES3LA(1))
      RES3LA(1) = RES3LA(2)
      RES3LA(2) = RES3LA(3)
      RES3LA(3) = RESULT
  100 ABSERR = AMAX1(ABSERR,0.5E+01*EPMACH*ABS(RESULT))
      RETURN
      END
*DECK R1MACH
      REAL FUNCTION R1MACH(I)
C***BEGIN PROLOGUE  R1MACH
C***DATE WRITTEN   790101   (YYMMDD)
C***REVISION DATE  890213   (YYMMDD)
C***CATEGORY NO.  R1
C***KEYWORDS  LIBRARY=SLATEC,TYPE=SINGLE PRECISION(R1MACH-S D1MACH-D),
C             MACHINE CONSTANTS
C***AUTHOR  FOX, P. A., (BELL LABS)
C           HALL, A. D., (BELL LABS)
C           SCHRYER, N. L., (BELL LABS)
C***PURPOSE  RETURNS SINGLE PRECISION MACHINE DEPENDENT CONSTANTS
C***DESCRIPTION
C
C   R1MACH CAN BE USED TO OBTAIN MACHINE-DEPENDENT PARAMETERS
C   FOR THE LOCAL MACHINE ENVIRONMENT.  IT IS A FUNCTION
C   SUBROUTINE WITH ONE (INPUT) ARGUMENT, AND CAN BE CALLED
C   AS FOLLOWS, FOR EXAMPLE
C
C        A = R1MACH(I)
C
C   WHERE I=1,...,5.  THE (OUTPUT) VALUE OF A ABOVE IS
C   DETERMINED BY THE (INPUT) VALUE OF I.  THE RESULTS FOR
C   VARIOUS VALUES OF I ARE DISCUSSED BELOW.
C
C   SINGLE-PRECISION MACHINE CONSTANTS
C   R1MACH(1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE.
C   R1MACH(2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE.
C   R1MACH(3) = B**(-T), THE SMALLEST RELATIVE SPACING.
C   R1MACH(4) = B**(1-T), THE LARGEST RELATIVE SPACING.
C   R1MACH(5) = LOG10(B)
C
C   ASSUME SINGLE PRECISION NUMBERS ARE REPRESENTED IN THE T-DIGIT,
C   BASE-B FORM
C
C              SIGN (B**E)*( (X(1)/B) + ... + (X(T)/B**T) )
C
C   WHERE 0 .LE. X(I) .LT. B FOR I=1,...,T, 0 .LT. X(1), AND
C   EMIN .LE. E .LE. EMAX.
C
C   THE VALUES OF B, T, EMIN AND EMAX ARE PROVIDED IN I1MACH AS
C   FOLLOWS:
C   I1MACH(10) = B, THE BASE.
C   I1MACH(11) = T, THE NUMBER OF BASE-B DIGITS.
C   I1MACH(12) = EMIN, THE SMALLEST EXPONENT E.
C   I1MACH(13) = EMAX, THE LARGEST EXPONENT E.
C
C   TO ALTER THIS FUNCTION FOR A PARTICULAR ENVIRONMENT,
C   THE DESIRED SET OF DATA STATEMENTS SHOULD BE ACTIVATED BY
C   REMOVING THE C FROM COLUMN 1.  ALSO, THE VALUES OF
C   R1MACH(1) - R1MACH(4) SHOULD BE CHECKED FOR CONSISTENCY
C   WITH THE LOCAL OPERATING SYSTEM.
C
C***REFERENCES  FOX, P.A., HALL, A.D., SCHRYER, N.L, *FRAMEWORK FOR
C                 A PORTABLE LIBRARY*, ACM TRANSACTIONS ON MATHE-
C                 MATICAL SOFTWARE, VOL. 4, NO. 2, JUNE 1978,
C                 PP. 177-188.
C***ROUTINES CALLED  XERROR
C
C THE CALL TO XERROR HAS BEEN COMMENTED TO GET A STAND ALONE ROUTINE
C IT CAN BE ADAPTED TO FIT THE LOCAL ERROR HANDLING PROCEDURES
C        NOTE ADDED BY F. ROMANI  7/11/89
C
C***END PROLOGUE  R1MACH
C
      INTEGER SMALL(2)
      INTEGER LARGE(2)
      INTEGER RIGHT(2)
      INTEGER DIVER(2)
      INTEGER LOG10(2)
C
      REAL RMACH(5)
      SAVE RMACH
C
      EQUIVALENCE (RMACH(1),SMALL(1))
      EQUIVALENCE (RMACH(2),LARGE(1))
      EQUIVALENCE (RMACH(3),RIGHT(1))
      EQUIVALENCE (RMACH(4),DIVER(1))
      EQUIVALENCE (RMACH(5),LOG10(1))
C
C     MACHINE CONSTANTS FOR THE AMIGA
C     ABSOFT FORTRAN COMPILER USING THE 68020/68881 COMPILER OPTION
C
C     DATA SMALL(1) / Z'00800000' /
C     DATA LARGE(1) / Z'7F7FFFFF' /
C     DATA RIGHT(1) / Z'33800000' /
C     DATA DIVER(1) / Z'34000000' /
C     DATA LOG10(1) / Z'3E9A209B' /
C
C     MACHINE CONSTANTS FOR THE AMIGA
C     ABSOFT FORTRAN COMPILER USING SOFTWARE FLOATING POINT
C
C     DATA SMALL(1) / Z'00800000' /
C     DATA LARGE(1) / Z'7EFFFFFF' /
C     DATA RIGHT(1) / Z'33800000' /
C     DATA DIVER(1) / Z'34000000' /
C     DATA LOG10(1) / Z'3E9A209B' /
C
C     MACHINE CONSTANTS FOR THE APOLLO
C
C     DATA SMALL(1) / 16#00800000 /
C     DATA LARGE(1) / 16#7FFFFFFF /
C     DATA RIGHT(1) / 16#33800000 /
C     DATA DIVER(1) / 16#34000000 /
C     DATA LOG10(1) / 16#3E9A209B /
C
C     MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM
C
C     DATA RMACH(1) / Z400800000 /
C     DATA RMACH(2) / Z5FFFFFFFF /
C     DATA RMACH(3) / Z4E9800000 /
C     DATA RMACH(4) / Z4EA800000 /
C     DATA RMACH(5) / Z500E730E8 /
C
C     MACHINE CONSTANTS FOR THE BURROUGHS 5700/6700/7700 SYSTEMS
C
C     DATA RMACH(1) / O1771000000000000 /
C     DATA RMACH(2) / O0777777777777777 /
C     DATA RMACH(3) / O1311000000000000 /
C     DATA RMACH(4) / O1301000000000000 /
C     DATA RMACH(5) / O1157163034761675 /
C
C     MACHINE CONSTANTS FOR THE CDC 170/180 SERIES USING NOS/VE
C
C     DATA RMACH(1) / Z"3001800000000000" /
C     DATA RMACH(2) / Z"4FFEFFFFFFFFFFFE" /
C     DATA RMACH(3) / Z"3FD2800000000000" /
C     DATA RMACH(4) / Z"3FD3800000000000" /
C     DATA RMACH(5) / Z"3FFF9A209A84FBCF" /
C
C     MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES
C
C     DATA RMACH(1) / 00564000000000000000B /
C     DATA RMACH(2) / 37767777777777777776B /
C     DATA RMACH(3) / 16414000000000000000B /
C     DATA RMACH(4) / 16424000000000000000B /
C     DATA RMACH(5) / 17164642023241175720B /
C
C     MACHINE CONSTANTS FOR THE CELERITY C1260
C
C     DATA SMALL(1) / Z'00800000' /
C     DATA LARGE(1) / Z'7F7FFFFF' /
C     DATA RIGHT(1) / Z'33800000' /
C     DATA DIVER(1) / Z'34000000' /
C     DATA LOG10(1) / Z'3E9A209B' /
C
C     MACHINE CONSTANTS FOR THE CONVEX C-1
C
C     DATA SMALL(1) / '00800000'X /
C     DATA LARGE(1) / '7FFFFFFF'X /
C     DATA RIGHT(1) / '34800000'X /
C     DATA DIVER(1) / '35000000'X /
C     DATA LOG10(1) / '3F9A209B'X /
C
C     MACHINE CONSTANTS FOR THE CRAY-1
C
C     DATA RMACH(1) / 200034000000000000000B /
C     DATA RMACH(2) / 577767777777777777776B /
C     DATA RMACH(3) / 377224000000000000000B /
C     DATA RMACH(4) / 377234000000000000000B /
C     DATA RMACH(5) / 377774642023241175720B /
C
C     MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200
C
C     NOTE - IT MAY BE APPROPRIATE TO INCLUDE THE FOLLOWING CARD -
C     STATIC RMACH(5)
C
C     DATA SMALL /    20K,       0 /
C     DATA LARGE / 77777K, 177777K /
C     DATA RIGHT / 35420K,       0 /
C     DATA DIVER / 36020K,       0 /
C     DATA LOG10 / 40423K,  42023K /
C
C     MACHINE CONSTANTS FOR THE HARRIS 220
C
C     DATA SMALL(1), SMALL(2) / '20000000, '00000201 /
C     DATA LARGE(1), LARGE(2) / '37777777, '00000177 /
C     DATA RIGHT(1), RIGHT(2) / '20000000, '00000352 /
C     DATA DIVER(1), DIVER(2) / '20000000, '00000353 /
C     DATA LOG10(1), LOG10(2) / '23210115, '00000377 /
C
C     MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 SERIES
C
C     DATA RMACH(1) / O402400000000 /
C     DATA RMACH(2) / O376777777777 /
C     DATA RMACH(3) / O714400000000 /
C     DATA RMACH(4) / O716400000000 /
C     DATA RMACH(5) / O776464202324 /
C
C     MACHINE CONSTANTS FOR THE HP 2100
C     3 WORD DOUBLE PRECISION WITH FTN4
C
C     DATA SMALL(1), SMALL(2) / 40000B,       1 /
C     DATA LARGE(1), LARGE(2) / 77777B, 177776B /
C     DATA RIGHT(1), RIGHT(2) / 40000B,    325B /
C     DATA DIVER(1), DIVER(2) / 40000B,    327B /
C     DATA LOG10(1), LOG10(2) / 46420B,  46777B /
C
C     MACHINE CONSTANTS FOR THE HP 2100
C     4 WORD DOUBLE PRECISION WITH FTN4
C
C     DATA SMALL(1), SMALL(2) / 40000B,       1 /
C     DATA LARGE91), LARGE(2) / 77777B, 177776B /
C     DATA RIGHT(1), RIGHT(2) / 40000B,    325B /
C     DATA DIVER(1), DIVER(2) / 40000B,    327B /
C     DATA LOG10(1), LOG10(2) / 46420B,  46777B /
C
C     MACHINE CONSTANTS FOR THE HP 9000
C
C     DATA SMALL(1) / 00004000000B /
C     DATA LARGE(1) / 17677777777B /
C     DATA RIGHT(1) / 06340000000B /
C     DATA DIVER(1) / 06400000000B /
C     DATA LOG10(1) / 07646420233B /
C
C     MACHINE CONSTANTS FOR THE ELXSI 6400
C       ASSUMING REAL*4 IS THE DEFAULT REAL
C
C     DATA SMALL(1) / '00800000'X /
C     DATA LARGE(1) / '7F7FFFFF'X /
C     DATA RIGHT(1) / '33800000'X /
C     DATA DIVER(1) / '34000000'X /
C     DATA LOG10(1) / '3E9A209B'X /
C
C     MACHINE CONSTANTS FOR THE IBM 360/370 SERIES,
C     THE XEROX SIGMA 5/7/9, THE SEL SYSTEMS 85/86  AND
C     THE PERKIN ELMER (INTERDATA) 7/32.
C
C     DATA RMACH(1) / Z00100000 /
C     DATA RMACH(2) / Z7FFFFFFF /
C     DATA RMACH(3) / Z3B100000 /
C     DATA RMACH(4) / Z3C100000 /
C     DATA RMACH(5) / Z41134413 /
C
C     MACHINE CONSTANTS FOR THE IBM PC
C
C     DATA SMALL(1) / 1.18E-38      /
C     DATA LARGE(1) / 3.40E+38      /
C     DATA RIGHT(1) / 0.595E-07     /
C     DATA DIVER(1) / 1.19E-07      /
C     DATA LOG10(1) / 0.30102999566 /
C
C     MACHINE CONSTANTS FOR THE PDP-10 (KA OR KI PROCESSOR)
C
C     DATA RMACH(1) / "000400000000 /
C     DATA RMACH(2) / "377777777777 /
C     DATA RMACH(3) / "146400000000 /
C     DATA RMACH(4) / "147400000000 /
C     DATA RMACH(5) / "177464202324 /
C
C     MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING
C     32-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL).
C
C     DATA SMALL(1) /    8388608 /
C     DATA LARGE(1) / 2147483647 /
C     DATA RIGHT(1) /  880803840 /
C     DATA DIVER(1) /  889192448 /
C     DATA LOG10(1) / 1067065499 /
C
C     DATA RMACH(1) / O00040000000 /
C     DATA RMACH(2) / O17777777777 /
C     DATA RMACH(3) / O06440000000 /
C     DATA RMACH(4) / O06500000000 /
C     DATA RMACH(5) / O07746420233 /
C
C     MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING
C     16-BIT INTEGERS  (EXPRESSED IN INTEGER AND OCTAL).
C
C     DATA SMALL(1), SMALL(2) /   128,     0 /
C     DATA LARGE(1), LARGE(2) / 32767,    -1 /
C     DATA RIGHT(1), RIGHT(2) / 13440,     0 /
C     DATA DIVER(1), DIVER(2) / 13568,     0 /
C     DATA LOG10(1), LOG10(2) / 16282,  8347 /
C
C     DATA SMALL(1), SMALL(2) / O000200, O000000 /
C     DATA LARGE(1), LARGE(2) / O077777, O177777 /
C     DATA RIGHT(1), RIGHT(2) / O032200, O000000 /
C     DATA DIVER(1), DIVER(2) / O032400, O000000 /
C     DATA LOG10(1), LOG10(2) / O037632, O020233 /
C
C     MACHINE CONSTANTS FOR THE SUN
C
C     DATA SMALL(1) / Z'00800000' /
C     DATA LARGE(1) / Z'7F7FFFFF' /
C     DATA RIGHT(1) / Z'33800000' /
C     DATA DIVER(1) / Z'34000000' /
C     DATA LOG10(1) / Z'3E9A209B' /
C
C     MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES
C
C     DATA RMACH(1) / O000400000000 /
C     DATA RMACH(2) / O377777777777 /
C     DATA RMACH(3) / O146400000000 /
C     DATA RMACH(4) / O147400000000 /
C     DATA RMACH(5) / O177464202324 /
C
C     MACHINE CONSTANTS FOR THE VAX 11/780
C     (EXPRESSED IN INTEGER AND HEXADECIMAL)
C     THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSTEMS
C     THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS
C
C     DATA SMALL(1) /       128 /
C     DATA LARGE(1) /    -32769 /
C     DATA RIGHT(1) /     13440 /
C     DATA DIVER(1) /     13568 /
C     DATA LOG10(1) / 547045274 /
C
C     DATA SMALL(1) / Z00000080 /
C     DATA LARGE(1) / ZFFFF7FFF /
C     DATA RIGHT(1) / Z00003480 /
C     DATA DIVER(1) / Z00003500 /
C     DATA LOG10(1) / Z209B3F9A /
C
C     MACHINE CONSTANTS FOR THE Z80 MICROPROCESSOR
C
C     DATA SMALL(1), SMALL(2) /     0,    256/
C     DATA LARGE(1), LARGE(2) /    -1,   -129/
C     DATA RIGHT(1), RIGHT(2) /     0,  26880/
C     DATA DIVER(1), DIVER(2) /     0,  27136/
C     DATA LOG10(1), LOG10(2) /  8347,  32538/
C
C
C***FIRST EXECUTABLE STATEMENT  R1MACH
C
C THE CALL TO XERROR HAS BEEN COMMENTED TO GET A STAND ALONE ROUTINE
C IT CAN BE ADAPTED TO FIT THE LOCAL ERROR HANDLING PROCEDURES
C        NOTE ADDED BY F. ROMANI  7/11/89
C
C     IF (I .LT. 1  .OR.  I .GT. 5)
C    1   CALL XERROR ('R1MACH -- I OUT OF BOUNDS', 25, 1, 2)
C
      R1MACH = RMACH(I)
      RETURN
C
      END
*DECK QXGS
      SUBROUTINE QXGS(F,A,B,EPSABS,EPSREL,RESULT,ABSERR,IER,
     1   LIMIT,LENIW,LENW,LAST,IWORK,WORK)
C
C            THE ROUTINE CALCULATES AN APPROXIMATION RESULT TO A GIVEN
C            DEFINITE INTEGRAL  I = INTEGRAL OF F OVER (A,B),
C            HOPEFULLY SATISFYING FOLLOWING CLAIM FOR ACCURACY
C            ABS(I-RESULT).LE.MAX(EPSABS,EPSREL*ABS(I)).
C
C        PARAMETERS
C         ON ENTRY
C            F      - REAL
C                     FUNCTION SUBPROGRAM DEFINING THE INTEGRAND
C                     FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE
C                     DECLARED E X T E R N A L IN THE DRIVER PROGRAM.
C
C            A      - REAL
C                     LOWER LIMIT OF INTEGRATION
C
C            B      - REAL
C                     UPPER LIMIT OF INTEGRATION
C
C            EPSABS - REAL
C                     ABSOLUTE ACCURACY REQUESTED
C            EPSREL - REAL
C                     RELATIVE ACCURACY REQUESTED
C                     IF  EPSABS.LE.0
C                     AND EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5E-28),
C                     THE ROUTINE WILL END WITH IER = 6.
C
C         ON RETURN
C            RESULT - REAL
C                     APPROXIMATION TO THE INTEGRAL
C
C            ABSERR - REAL
C                     ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR,
C                     WHICH SHOULD EQUAL OR EXCEED ABS(I-RESULT)
C
C            IER    - INTEGER
C                     IER = 0 NORMAL AND RELIABLE TERMINATION OF THE
C                             ROUTINE. IT IS ASSUMED THAT THE REQUESTED
C                             ACCURACY HAS BEEN ACHIEVED.
C                     IER.GT.0 ABNORMAL TERMINATION OF THE ROUTINE
C                             THE ESTIMATES FOR INTEGRAL AND ERROR ARE
C                             LESS RELIABLE. IT IS ASSUMED THAT THE
C                             REQUESTED ACCURACY HAS NOT BEEN ACHIEVED.
C            ERROR MESSAGES
C                     IER = 1 MAXIMUM NUMBER OF SUBDIVISIONS ALLOWED
C                             HAS BEEN ACHIEVED. ONE CAN ALLOW MORE SUB-
C                             DIVISIONS BY INCREASING THE VALUE OF LIMIT
C                             (AND TAKING THE ACCORDING DIMENSION
C                             ADJUSTMENTS INTO ACCOUNT. HOWEVER, IF
C                             THIS YIELDS NO IMPROVEMENT IT IS ADVISED
C                             TO ANALYZE THE INTEGRAND IN ORDER TO
C                             DETERMINE THE INTEGRATION DIFFICULTIES. IF
C                             THE POSITION OF A LOCAL DIFFICULTY CAN BE
C                             DETERMINED (E.G. SINGULARITY,
C                             DISCONTINUITY WITHIN THE INTERVAL) ONE
C                             WILL PROBABLY GAIN FROM SPLITTING UP THE
C                             INTERVAL AT THIS POINT AND CALLING THE
C                             INTEGRATOR ON THE SUBRANGES. IF POSSIBLE,
C                             AN APPROPRIATE SPECIAL-PURPOSE INTEGRATOR
C                             SHOULD BE USED, WHICH IS DESIGNED FOR
C                             HANDLING THE TYPE OF DIFFICULTY INVOLVED.
C                         = 2 THE OCCURRENCE OF ROUNDOFF ERROR IS DETEC-
C                             TED, WHICH PREVENTS THE REQUESTED
C                             TOLERANCE FROM BEING ACHIEVED.
C                             THE ERROR MAY BE UNDER-ESTIMATED.
C                         = 3 EXTREMELY BAD INTEGRAND BEHAVIOUR
C                             OCCURS AT SOME POINTS OF THE INTEGRATION
C                             INTERVAL.
C                         = 4 THE ALGORITHM DOES NOT CONVERGE.
C                             ROUNDOFF ERROR IS DETECTED IN THE
C                             EXTRAPOLATION TABLE. IT IS PRESUMED THAT
C                             THE REQUESTED TOLERANCE CANNOT BE
C                             ACHIEVED, AND THAT THE RETURNED RESULT IS
C                             THE BEST WHICH CAN BE OBTAINED.
C                         = 5 THE INTEGRAL IS PROBABLY DIVERGENT, OR
C                             SLOWLY CONVERGENT. IT MUST BE NOTED THAT
C                             DIVERGENCE CAN OCCUR WITH ANY OTHER VALUE
C                             OF IER.
C                         = 6 THE INPUT IS INVALID, BECAUSE
C                             (EPSABS.LE.0 AND
C                              EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5E-28)
C                             OR LIMIT.LT.1 OR LENW.LT.LIMIT*46 OR
C                                LENIW .LT. LIMIT*3
C                             RESULT, ABSERR, LAST ARE SET TO
C                             ZERO.EXCEPT WHEN LIMIT OR LENW OR LENIW
C                             IS INVALID, IWORK(1), WORK(LIMIT*2+1) AND
C                             WORK(LIMIT*3+1) ARE SET TO ZERO, WORK(1)
C                             IS SET TO A AND WORK(LIMIT+1) TO B.
C
C         DIMENSIONING PARAMETERS
C            LIMIT - INTEGER
C                    LIMIT DETERMINES THE MAXIMUM NUMBER OF SUBINTERVALS
C                    IN THE PARTITION OF THE GIVEN INTEGRATION INTERVAL
C                    (A,B), LIMIT.GE.1.
C                    IF LIMIT.LT.1, THE ROUTINE WILL END WITH IER = 6.
C
C            LENW  - INTEGER
C                    DIMENSIONING PARAMETER FOR WORK
C                    LENW MUST BE AT LEAST LIMIT*46
C                    IF LENW.LT.LIMIT*46, THE ROUTINE WILL END
C                    WITH IER = 6.
C
C            LENIW - INTEGER
C                    DIMENSIONING PARAMETER FOR IWORK
C                    LENIW MUST BE AT LEAST LIMIT*3
C                    IF LENW.LT.LIMIT*3, THE ROUTINE WILL END
C                    WITH IER = 6.
C
C            LAST  - INTEGER
C                    ON RETURN, LAST EQUALS THE NUMBER OF SUBINTERVALS
C                    PRODUCED IN THE SUBDIVISION PROCESS, DETERMINES THE
C                    NUMBER OF SIGNIFICANT ELEMENTS ACTUALLY IN THE WORK
C                    ARRAYS.
C
C         WORK ARRAYS
C            IWORK - INTEGER
C                    VECTOR OF DIMENSION AT LEAST 3*LIMIT, THE FIRST K
C                    ELEMENTS OF WHICH CONTAIN POINTERS
C                    TO THE ERROR ESTIMATES OVER THE SUBINTERVALS
C                    SUCH THAT WORK(LIMIT*3+IWORK(1)),... ,
C                    WORK(LIMIT*3+IWORK(K)) FORM A DECREASING
C                    SEQUENCE, WITH K = LAST IF LAST.LE.(LIMIT/2+2),
C                    AND K = LIMIT+1-LAST OTHERWISE
C
C            WORK  - REAL
C                    VECTOR OF DIMENSION AT LEAST LENW
C                    ON RETURN
C                    WORK(1), ..., WORK(LAST) CONTAIN THE LEFT
C                     END-POINTS OF THE SUBINTERVALS IN THE
C                     PARTITION OF (A,B),
C                    WORK(LIMIT+1), ..., WORK(LIMIT+LAST) CONTAIN
C                     THE RIGHT END-POINTS,
C                    WORK(LIMIT*2+1), ..., WORK(LIMIT*2+LAST) CONTAIN
C                     THE INTEGRAL APPROXIMATIONS OVER THE SUBINTERVALS,
C                    WORK(LIMIT*3+1), ..., WORK(LIMIT*3+LAST)
C                     CONTAIN THE ERROR ESTIMATES.
C                    WORK(LIMIT*4+1), ... IS THE AREA RESERVED TO STORE
C                     FUNCTIONAL VALUES.
C***ROUTINES CALLED  QXGSE,XERROR
C
C THE CALL TO XERROR HAS BEEN COMMENTED TO GET A STAND ALONE ROUTINE
C IT CAN BE ADAPTED TO FIT THE LOCAL ERROR HANDLING PROCEDURES
C
      REAL A,ABSERR,B,EPSABS,EPSREL,F,RESULT,WORK
      INTEGER IER,IWORK,LAST,LENIW,LENW,LIMIT,LVL,L1,L2,L3,L4,L5
C
      DIMENSION IWORK(LENIW),WORK(LENW)
C
      EXTERNAL F
C
C         CHECK VALIDITY OF LIMIT,LENIW AND LENW.
C
C***FIRST EXECUTABLE STATEMENT  QXGS
      IER = 6
      LAST = 0
      RESULT = 0.0E+00
      ABSERR = 0.0E+00
      IF(LIMIT.LT.1.OR.LENIW.LT.LIMIT*3.OR.LENW.LT.LIMIT*46) GO TO 10
C
C         PREPARE CALL FOR QXGSE.
C
      L1 = LIMIT+1
      L2 = LIMIT+L1
      L3 = LIMIT+L2
      L4 = LIMIT+L3
      L5 = 21*LIMIT+L3
C
      CALL QXGSE(F,A,B,EPSABS,EPSREL,LIMIT,RESULT,ABSERR,
     *  IER,WORK(1),WORK(L1),WORK(L2),WORK(L3),IWORK,LAST,
     *  WORK(L4),WORK(L5),IWORK(L1),IWORK(L2))
C
C         CALL ERROR HANDLER IF NECESSARY.
C
      LVL = 0
10    IF(IER.EQ.6) LVL = 1
C
C THE CALL TO XERROR HAS BEEN COMMENTED TO GET A STAND ALONE ROUTINE
C IT CAN BE ADAPTED TO FIT THE LOCAL ERROR HANDLING PROCEDURES
C
C     IF(IER.NE.0)CALL XERROR('ABNORMAL RETURN FROM QXGS  ',28,IER,LVL)
      RETURN
      END
*DECK QXGSE
      SUBROUTINE QXGSE(F,A,B,EPSABS,EPSREL,LIMIT,RESULT,ABSERR,
     *   IER,ALIST,BLIST,RLIST,ELIST,IORD,LAST,VALP,VALN,LP,LN)
C
C            THE ROUTINE CALCULATES AN APPROXIMATION RESULT TO A GIVEN
C            DEFINITE INTEGRAL I = INTEGRAL OF F OVER (A,B),
C            HOPEFULLY SATISFYING FOLLOWING CLAIM FOR ACCURACY
C            ABS(I-RESULT).LE.MAX(EPSABS,EPSREL*ABS(I)).
C
C        PARAMETERS
C         ON ENTRY
C            F      - REAL
C                     FUNCTION SUBPROGRAM DEFINING THE INTEGRAND
C                     FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE
C                     DECLARED E X T E R N A L IN THE DRIVER PROGRAM.
C
C            A      - REAL
C                     LOWER LIMIT OF INTEGRATION
C
C            B      - REAL
C                     UPPER LIMIT OF INTEGRATION
C
C            EPSABS - REAL
C                     ABSOLUTE ACCURACY REQUESTED
C            EPSREL - REAL
C                     RELATIVE ACCURACY REQUESTED
C                     IF  EPSABS.LE.0
C                     AND EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5E-28),
C                     THE ROUTINE WILL END WITH IER = 6.
C
C            LIMIT  - INTEGER
C                     GIVES AN UPPERBOUND ON THE NUMBER OF SUBINTERVALS
C                     IN THE PARTITION OF (A,B)
C
C         ON RETURN
C            RESULT - REAL
C                     APPROXIMATION TO THE INTEGRAL
C
C            ABSERR - REAL
C                     ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR,
C                     WHICH SHOULD EQUAL OR EXCEED ABS(I-RESULT)
C
C            IER    - INTEGER
C                     IER = 0 NORMAL AND RELIABLE TERMINATION OF THE
C                             ROUTINE. IT IS ASSUMED THAT THE REQUESTED
C                             ACCURACY HAS BEEN ACHIEVED.
C                     IER.GT.0 ABNORMAL TERMINATION OF THE ROUTINE
C                             THE ESTIMATES FOR INTEGRAL AND ERROR ARE
C                             LESS RELIABLE. IT IS ASSUMED THAT THE
C                             REQUESTED ACCURACY HAS NOT BEEN ACHIEVED.
C            ERROR MESSAGES
C                         = 1 MAXIMUM NUMBER OF SUBDIVISIONS ALLOWED
C                             HAS BEEN ACHIEVED. ONE CAN ALLOW MORE SUB-
C                             DIVISIONS BY INCREASING THE VALUE OF LIMIT
C                             (AND TAKING THE ACCORDING DIMENSION
C                             ADJUSTMENTS INTO ACCOUNT). HOWEVER, IF
C                             THIS YIELDS NO IMPROVEMENT IT IS ADVISED
C                             TO ANALYZE THE INTEGRAND IN ORDER TO
C                             DETERMINE THE INTEGRATION DIFFICULTIES. IF
C                             THE POSITION OF A LOCAL DIFFICULTY CAN BE
C                             DETERMINED (E.G. SINGULARITY,
C                             DISCONTINUITY WITHIN THE INTERVAL) ONE
C                             WILL PROBABLY GAIN FROM SPLITTING UP THE
C                             INTERVAL AT THIS POINT AND CALLING THE
C                             INTEGRATOR ON THE SUBRANGES. IF POSSIBLE,
C                             AN APPROPRIATE SPECIAL-PURPOSE INTEGRATOR
C                             SHOULD BE USED, WHICH IS DESIGNED FOR
C                             HANDLING THE TYPE OF DIFFICULTY INVOLVED.
C                         = 2 THE OCCURRENCE OF ROUNDOFF ERROR IS DETEC-
C                             TED, WHICH PREVENTS THE REQUESTED
C                             TOLERANCE FROM BEING ACHIEVED.
C                             THE ERROR MAY BE UNDER-ESTIMATED.
C                         = 3 EXTREMELY BAD INTEGRAND BEHAVIOUR
C                             OCCURS AT SOME POINTS OF THE INTEGRATION
C                             INTERVAL.
C                         = 4 THE ALGORITHM DOES NOT CONVERGE.
C                             ROUNDOFF ERROR IS DETECTED IN THE
C                             EXTRAPOLATION TABLE.
C                             IT IS PRESUMED THAT THE REQUESTED
C                             TOLERANCE CANNOT BE ACHIEVED, AND THAT THE
C                             RETURNED RESULT IS THE BEST WHICH CAN BE
C                             OBTAINED.
C                         = 5 THE INTEGRAL IS PROBABLY DIVERGENT, OR
C                             SLOWLY CONVERGENT. IT MUST BE NOTED THAT
C                             DIVERGENCE CAN OCCUR WITH ANY OTHER VALUE
C                             OF IER.
C                         = 6 THE INPUT IS INVALID, BECAUSE
C                             EPSABS.LE.0 AND
C                             EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5E-28).
C                             RESULT, ABSERR, LAST, RLIST(1),
C                             IORD(1) AND ELIST(1) ARE SET TO ZERO.
C                             ALIST(1) AND BLIST(1) ARE SET TO A AND B
C                             RESPECTIVELY.
C
C            ALIST  - REAL
C                     VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
C                      LAST  ELEMENTS OF WHICH ARE THE LEFT END POINTS
C                     OF THE SUBINTERVALS IN THE PARTITION OF THE
C                     GIVEN INTEGRATION RANGE (A,B)
C
C            BLIST  - REAL
C                     VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
C                      LAST  ELEMENTS OF WHICH ARE THE RIGHT END POINTS
C                     OF THE SUBINTERVALS IN THE PARTITION OF THE GIVEN
C                     INTEGRATION RANGE (A,B)
C
C            RLIST  - REAL
C                     VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
C                      LAST  ELEMENTS OF WHICH ARE THE INTEGRAL
C                     APPROXIMATIONS ON THE SUBINTERVALS
C
C            ELIST  - REAL
C                     VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
C                      LAST  ELEMENTS OF WHICH ARE THE MODULI OF THE
C                     ABSOLUTE ERROR ESTIMATES ON THE SUBINTERVALS
C
C            IORD   - INTEGER
C                     VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST K
C                     ELEMENTS OF WHICH ARE POINTERS TO THE
C                     ERROR ESTIMATES OVER THE SUBINTERVALS,
C                     SUCH THAT ELIST(IORD(1)), ..., ELIST(IORD(K))
C                     FORM A DECREASING SEQUENCE, WITH K = LAST
C                     IF LAST.LE.(LIMIT/2+2), AND K = LIMIT+1-LAST
C                     OTHERWISE
C
C            LAST   - INTEGER
C                     NUMBER OF SUBINTERVALS ACTUALLY PRODUCED IN THE
C                     SUBDIVISION PROCESS
C
C            VALP   - REAL
C            VALN     ARRAYS OF DIMENSION AT LEAST (21,LIMIT) USED TO
C                     SAVE THE FUNCTIONAL VALUES
C
C            LP     - INTEGER
C            LN       VECTORS OF DIMENSION AT LEAST LIMIT, USED TO
C                     STORE THE ACTUAL NUMBER OF FUNCTIONAL VALUES
C                     SAVED IN THE CORRESPONDING COLUMN
C                     OF VALP,VALN
C
C***ROUTINES CALLED  F,R1MACH,QELG,QXLQM,QPSRT,QXRRD,QXCPY
C
      REAL A,ABSEPS,ABSERR,ALIST,AREA,AREA1,AREA12,AREA2,A1,
     *  A2,B,BLIST,B1,B2,CORREC,ABS,DEFABS,DEFAB1,DEFAB2,R1MACH,AMAX1,
     *  DRES,ELIST,EPMACH,EPSABS,EPSREL,ERLARG,ERLAST,ERRBND,ERRMAX,
     *  ERROR1,ERROR2,ERRO12,ERRSUM,ERTEST,F,OFLOW,RESABS,RESEPS,RESULT,
     *  RES3LA,RLIST,RLIST2,SMALL,UFLOW,
     *  VALP,VALN,VP1,VP2,VN1,VN2
      INTEGER ID,IER,IERRO,IORD,IROFF1,IROFF2,IROFF3,JUPBND,K,KSGN,
     *  KTMIN,LAST,LIMIT,MAXERR,NRES,NRMAX,NUMRL2,
     *  LP,LN,LP1,LP2,LN1,LN2
      LOGICAL EXTRAP,NOEXT
      DIMENSION ALIST(LIMIT),BLIST(LIMIT),ELIST(LIMIT),IORD(LIMIT),
     * RES3LA(3),RLIST(LIMIT),RLIST2(52),
     * VALP(21,LIMIT),VALN(21,LIMIT),LP(LIMIT),LN(LIMIT),
     * VP1(21),VP2(21),VN1(21),VN2(21)
      EXTERNAL F
C
C            MACHINE DEPENDENT CONSTANTS
C            ---------------------------
C
C          EPMACH IS THE LARGEST RELATIVE SPACING.
C          UFLOW IS THE SMALLEST POSITIVE MAGNITUDE.
C          OFLOW IS THE LARGEST POSITIVE MAGNITUDE.
C
C***FIRST EXECUTABLE STATEMENT  QXGSE
      EPMACH = R1MACH(4)
C
C            TEST ON VALIDITY OF PARAMETERS
C            ------------------------------
      IER = 0
      LAST = 0
      RESULT = 0.0E+00
      ABSERR = 0.0E+00
      ALIST(1) = A
      BLIST(1) = B
      RLIST(1) = 0.0E+00
      ELIST(1) = 0.0E+00
      IF(EPSABS.LE.0.0E+00.AND.EPSREL.LT.AMAX1(0.5E+02*EPMACH,0.5E-28))
     *   IER = 6
      IF(IER.EQ.6) GO TO 999
C
C           FIRST APPROXIMATION TO THE INTEGRAL
C           -----------------------------------
C
      UFLOW = R1MACH(1)
      OFLOW = R1MACH(2)
      IERRO = 0
      LP(1)=1
      LN(1)=1
      VALP(1,1)=F((A+B)*0.5E0)
      VALN(1,1)=VALP(1,1)
      CALL QXLQM(F,A,B,RESULT,ABSERR,DEFABS,RESABS,
     *           VALP(1,1),VALN(1,1),LP(1),LN(1), 2  )
C
C           TEST ON ACCURACY.
C
      DRES = ABS(RESULT)
      ERRBND = AMAX1(EPSABS,EPSREL*DRES)
      LAST = 1
      RLIST(1) = RESULT
      ELIST(1) = ABSERR
      IORD(1) = 1
      IF(ABSERR.LE.1.0E+02*EPMACH*DEFABS.AND.ABSERR.GT.ERRBND) IER = 2
      IF(LIMIT.EQ.1) IER = 1
      IF(IER.NE.0.OR.(ABSERR.LE.ERRBND.AND.ABSERR.NE.RESABS).OR.
     *  ABSERR.EQ.0.0E+00) GO TO 999
C
C           INITIALIZATION
C           --------------
C
      RLIST2(1) = RESULT
      ERRMAX = ABSERR
      MAXERR = 1
      AREA = RESULT
      ERRSUM = ABSERR
      ABSERR = OFLOW
      NRMAX = 1
      NRES = 0
      NUMRL2 = 2
      KTMIN = 0
      EXTRAP = .FALSE.
      NOEXT = .FALSE.
      IROFF1 = 0
      IROFF2 = 0
      IROFF3 = 0
      KSGN = -1
      IF(DRES.GE.(0.1E+01-0.5E+02*EPMACH)*DEFABS) KSGN = 1
C
C           MAIN DO-LOOP
C           ------------
C
      DO 90 LAST = 2,LIMIT
C
C           BISECT THE SUBINTERVAL WITH THE NRMAX-TH LARGEST ERROR
C           ESTIMATE.
C
        A1 = ALIST(MAXERR)
        B1 = 0.5E+00*(ALIST(MAXERR)+BLIST(MAXERR))
        A2 = B1
        B2 = BLIST(MAXERR)
        ERLAST = ERRMAX
        CALL QXRRD(F,VALN(1,MAXERR),LN(MAXERR),B1,A1,VN1,VP1,LN1,LP1)
        CALL QXLQM(F,A1,B1,AREA1,ERROR1,RESABS,DEFAB1,VP1,VN1,LP1,LN1,
     *              2)
        CALL QXRRD(F,VALP(1,MAXERR),LP(MAXERR),A2,B2,VP2,VN2,LP2,LN2)
        CALL QXLQM(F,A2,B2,AREA2,ERROR2,RESABS,DEFAB2,VP2,VN2,LP2,LN2,
     *              2)
C
C           IMPROVE PREVIOUS APPROXIMATIONS TO INTEGRAL
C           AND ERROR AND TEST FOR ACCURACY.
C
        AREA12 = AREA1+AREA2
        ERRO12 = ERROR1+ERROR2
        ERRSUM = ERRSUM+ERRO12-ERRMAX
        AREA = AREA+AREA12-RLIST(MAXERR)
        IF(DEFAB1.EQ.ERROR1.OR.DEFAB2.EQ.ERROR2) GO TO 15
        IF(ABS(RLIST(MAXERR)-AREA12).GT.0.1E-04*ABS(AREA12)
     *  .OR.ERRO12.LT.0.99E+00*ERRMAX) GO TO 10
        IF(EXTRAP) IROFF2 = IROFF2+1
        IF(.NOT.EXTRAP) IROFF1 = IROFF1+1
   10   IF(LAST.GT.10.AND.ERRO12.GT.ERRMAX) IROFF3 = IROFF3+1
   15   RLIST(MAXERR) = AREA1
        RLIST(LAST) = AREA2
        ERRBND = AMAX1(EPSABS,EPSREL*ABS(AREA))
C
C           TEST FOR ROUNDOFF ERROR AND EVENTUALLY SET ERROR FLAG.
C
        IF(IROFF1+IROFF2.GE.10.OR.IROFF3.GE.20) IER = 2
        IF(IROFF2.GE.5) IERRO = 3
C
C           SET ERROR FLAG IN THE CASE THAT THE NUMBER OF SUBINTERVALS
C           EQUALS LIMIT.
C
        IF(LAST.EQ.LIMIT) IER = 1
C
C           SET ERROR FLAG IN THE CASE OF BAD INTEGRAND BEHAVIOUR
C           AT A POINT OF THE INTEGRATION RANGE.
C
        IF(AMAX1(ABS(A1),ABS(B2)).LE.(0.1E+01+0.1E+03*EPMACH)*
     *  (ABS(A2)+0.1E+04*UFLOW)) IER = 4
C
C           APPEND THE NEWLY-CREATED INTERVALS TO THE LIST.
C
        IF(ERROR2.GT.ERROR1) GO TO 20
        ALIST(LAST) = A2
        BLIST(MAXERR) = B1
        BLIST(LAST) = B2
        ELIST(MAXERR) = ERROR1
        ELIST(LAST) = ERROR2
        CALL QXCPY(VALP(1,MAXERR),VP1,LP1)
        LP(MAXERR)=LP1
        CALL QXCPY(VALN(1,MAXERR),VN1,LN1)
        LN(MAXERR)=LN1
        CALL QXCPY(VALP(1,LAST),VP2,LP2)
        LP(LAST)=LP2
        CALL QXCPY(VALN(1,LAST),VN2,LN2)
        LN(LAST)=LN2
        GO TO 30
   20   ALIST(MAXERR) = A2
        ALIST(LAST) = A1
        BLIST(LAST) = B1
        RLIST(MAXERR) = AREA2
        RLIST(LAST) = AREA1
        ELIST(MAXERR) = ERROR2
        ELIST(LAST) = ERROR1
        CALL QXCPY(VALP(1,MAXERR),VP2,LP2)
        LP(MAXERR)=LP2
        CALL QXCPY(VALN(1,MAXERR),VN2,LN2)
        LN(MAXERR)=LN2
        CALL QXCPY(VALP(1,LAST),VP1,LP1)
        LP(LAST)=LP1
        CALL QXCPY(VALN(1,LAST),VN1,LN1)
        LN(LAST)=LN1
C
C           CALL SUBROUTINE QPSRT TO MAINTAIN THE DESCENDING ORDERING
C           IN THE LIST OF ERROR ESTIMATES AND SELECT THE SUBINTERVAL
C           WITH NRMAX-TH LARGEST ERROR ESTIMATE (TO BE BISECTED NEXT).
C
   30   CALL QPSRT(LIMIT,LAST,MAXERR,ERRMAX,ELIST,IORD,NRMAX)
C ***JUMP OUT OF DO-LOOP
        IF(ERRSUM.LE.ERRBND) GO TO 115
C ***JUMP OUT OF DO-LOOP
        IF(IER.NE.0) GO TO 100
        IF(LAST.EQ.2) GO TO 80
        IF(NOEXT) GO TO 90
        ERLARG = ERLARG-ERLAST
        IF(ABS(B1-A1).GT.SMALL) ERLARG = ERLARG+ERRO12
        IF(EXTRAP) GO TO 40
C
C           TEST WHETHER THE INTERVAL TO BE BISECTED NEXT IS THE
C           SMALLEST INTERVAL.
C
        IF(ABS(BLIST(MAXERR)-ALIST(MAXERR)).GT.SMALL) GO TO 90
        EXTRAP = .TRUE.
        NRMAX = 2
C
C           THE BOUND 0.3*ERTEST HAS BEEN INTRODUCED TO PERFORM A
C           MORE CAUTIOUS EXTRAPOLATION THAN IN THE ORIGINAL DQAGSE
C           ROUTINE
C
   40   IF(IERRO.EQ.3.OR.ERLARG.LE.0.3E0*ERTEST) GO TO 60
C
C           THE SMALLEST INTERVAL HAS THE LARGEST ERROR.
C           BEFORE BISECTING DECREASE THE SUM OF THE ERRORS OVER THE
C           LARGER INTERVALS (ERLARG) AND PERFORM EXTRAPOLATION.
C
        ID = NRMAX
        JUPBND = LAST
        IF(LAST.GT.(2+LIMIT/2)) JUPBND = LIMIT+3-LAST
        DO 50 K = ID,JUPBND
          MAXERR = IORD(NRMAX)
          ERRMAX = ELIST(MAXERR)
C ***JUMP OUT OF DO-LOOP
          IF(ABS(BLIST(MAXERR)-ALIST(MAXERR)).GT.SMALL) GO TO 90
          NRMAX = NRMAX+1
   50   CONTINUE
C
C           PERFORM EXTRAPOLATION.
C
   60   NUMRL2 = NUMRL2+1
        RLIST2(NUMRL2) = AREA
        CALL QELG(NUMRL2,RLIST2,RESEPS,ABSEPS,RES3LA,NRES)
        KTMIN = KTMIN+1
        IF(KTMIN.GT.5.AND.ABSERR.LT.0.1E-02*ERRSUM) IER = 5
        IF(ABSEPS.GE.ABSERR) GO TO 70
        KTMIN = 0
        ABSERR = ABSEPS
        RESULT = RESEPS
        CORREC = ERLARG
        ERTEST = AMAX1(EPSABS,EPSREL*ABS(RESEPS))
C ***JUMP OUT OF DO-LOOP
        IF(ABSERR.LE.ERTEST) GO TO 100
C
C           PREPARE BISECTION OF THE SMALLEST INTERVAL.
C
   70   IF(NUMRL2.EQ.1) NOEXT = .TRUE.
        IF(IER.EQ.5) GO TO 100
        MAXERR = IORD(1)
        ERRMAX = ELIST(MAXERR)
        NRMAX = 1
        EXTRAP = .FALSE.
        SMALL = SMALL*0.5E+00
        ERLARG = ERRSUM
        GO TO 90
   80   SMALL = ABS(B-A)*0.375E+00
        ERLARG = ERRSUM
        ERTEST = ERRBND
        RLIST2(2) = AREA
   90 CONTINUE
C
C           SET FINAL RESULT AND ERROR ESTIMATE.
C           ------------------------------------
C
  100 IF(ABSERR.EQ.OFLOW) GO TO 115
      IF(IER+IERRO.EQ.0) GO TO 110
      IF(IERRO.EQ.3) ABSERR = ABSERR+CORREC
      IF(IER.EQ.0) IER = 3
      IF(RESULT.NE.0.0E+00.AND.AREA.NE.0.0E+00) GO TO 105
      IF(ABSERR.GT.ERRSUM) GO TO 115
      IF(AREA.EQ.0.0E+00) GO TO 130
      GO TO 110
  105 IF(ABSERR/ABS(RESULT).GT.ERRSUM/ABS(AREA)) GO TO 115
C
C           TEST ON DIVERGENCE.
C
  110 IF(KSGN.EQ.(-1).AND.AMAX1(ABS(RESULT),ABS(AREA)).LE.
     * DEFABS*0.1E-01) GO TO 130
      IF(0.1E-01.GT.(RESULT/AREA).OR.(RESULT/AREA).GT.0.1E+03
     * .OR.ERRSUM.GT.ABS(AREA)) IER = 6
      GO TO 130
C
C           COMPUTE GLOBAL INTEGRAL SUM.
C
  115 RESULT = 0.0E+00
      DO 120 K = 1,LAST
         RESULT = RESULT+RLIST(K)
  120 CONTINUE
      ABSERR = ERRSUM
  130 IF(IER.GT.2) IER = IER-1
  999 RETURN
      END
*DECK QXG
      SUBROUTINE QXG(F,A,B,EPSABS,EPSREL,KEY,RESULT,ABSERR,IER,
     1   LIMIT,LENIW,LENW,LAST,IWORK,WORK)
C
C            THE ROUTINE CALCULATES AN APPROXIMATION RESULT TO A GIVEN
C            DEFINITE INTEGRAL  I = INTEGRAL OF F OVER (A,B),
C            HOPEFULLY SATISFYING FOLLOWING CLAIM FOR ACCURACY
C            ABS(I-RESULT).LE.MAX(EPSABS,EPSREL*ABS(I)).
C
C        PARAMETERS
C         ON ENTRY
C            F      - REAL
C                     FUNCTION SUBPROGRAM DEFINING THE INTEGRAND
C                     FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE
C                     DECLARED E X T E R N A L IN THE DRIVER PROGRAM.
C
C            A      - REAL
C                     LOWER LIMIT OF INTEGRATION
C
C            B      - REAL
C                     UPPER LIMIT OF INTEGRATION
C
C            EPSABS - REAL
C                     ABSOLUTE ACCURACY REQUESTED
C            EPSREL - REAL
C                     RELATIVE ACCURACY REQUESTED
C                     IF  EPSABS.LE.0
C                     AND EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5E-28),
C                     THE ROUTINE WILL END WITH IER = 6.
C
C            KEY    - INTEGER
C                     KEY FOR CHOICE OF LOCAL INTEGRATION RULE
C                     RMS FORMULAS ARE USED WITH
C                      13 - 19               POINTS IF KEY.LT.1,
C                      13 - 19 - (27)        POINTS IF KEY = 1,
C                      13 - 19 - (27) - (41) POINTS IF KEY = 2,
C                           19 -  27  - (41) POINTS IF KEY = 3,
C                                 27  -  41  POINTS IF KEY.GT.3.
C
C                         (RULES) USED IF THE FUNCTION APPEARS
C                         ENOUGH REGULAR IN THE SUBINTERVAL
C
C         ON RETURN
C            RESULT - REAL
C                     APPROXIMATION TO THE INTEGRAL
C
C            ABSERR - REAL
C                     ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR,
C                     WHICH SHOULD EQUAL OR EXCEED ABS(I-RESULT)
C
C            IER    - INTEGER
C                     IER = 0 NORMAL AND RELIABLE TERMINATION OF THE
C                             ROUTINE. IT IS ASSUMED THAT THE REQUESTED
C                             ACCURACY HAS BEEN ACHIEVED.
C                     IER.GT.0 ABNORMAL TERMINATION OF THE ROUTINE
C                             THE ESTIMATES FOR RESULT AND ERROR ARE
C                             LESS RELIABLE. IT IS ASSUMED THAT THE
C                             REQUESTED ACCURACY HAS NOT BEEN ACHIEVED.
C            ERROR MESSAGES
C                     IER = 1 MAXIMUM NUMBER OF SUBDIVISIONS ALLOWED
C                             HAS BEEN ACHIEVED. ONE CAN ALLOW MORE SUB-
C                             DIVISIONS BY INCREASING THE VALUE OF LIMIT
C                             (AND TAKING THE ACCORDING DIMENSION
C                             ADJUSTMENTS INTO ACCOUNT. HOWEVER, IF
C                             THIS YIELDS NO IMPROVEMENT IT IS ADVISED
C                             TO ANALYZE THE INTEGRAND IN ORDER TO
C                             DETERMINE THE INTEGRATION DIFFICULTIES. IF
C                             THE POSITION OF A LOCAL DIFFICULTY CAN BE
C                             DETERMINED (E.G. SINGULARITY,
C                             DISCONTINUITY WITHIN THE INTERVAL) ONE
C                             WILL PROBABLY GAIN FROM SPLITTING UP THE
C                             INTERVAL AT THIS POINT AND CALLING THE
C                             INTEGRATOR ON THE SUBRANGES. IF POSSIBLE,
C                             AN APPROPRIATE SPECIAL-PURPOSE INTEGRATOR
C                             SHOULD BE USED, WHICH IS DESIGNED FOR
C                             HANDLING THE TYPE OF DIFFICULTY INVOLVED.
C                         = 2 THE OCCURRENCE OF ROUNDOFF ERROR IS DETEC-
C                             TED, WHICH PREVENTS THE REQUESTED
C                             TOLERANCE FROM BEING ACHIEVED.
C                             THE ERROR MAY BE UNDER-ESTIMATED.
C                         = 3 EXTREMELY BAD INTEGRAND BEHAVIOUR
C                             OCCURS AT SOME POINTS OF THE INTEGRATION
C                             INTERVAL.
C                         = 4 THE ALGORITHM DOES NOT CONVERGE.
C                             ROUNDOFF ERROR IS DETECTED IN THE
C                             EXTRAPOLATION TABLE. IT IS PRESUMED THAT
C                             THE REQUESTED TOLERANCE CANNOT BE
C                             ACHIEVED, AND THAT THE RETURNED RESULT IS
C                             THE BEST WHICH CAN BE OBTAINED.
C                         = 5 THE INTEGRAL IS PROBABLY DIVERGENT, OR
C                             SLOWLY CONVERGENT. IT MUST BE NOTED THAT
C                             DIVERGENCE CAN OCCUR WITH ANY OTHER VALUE
C                             OF IER.
C                         = 6 THE INPUT IS INVALID, BECAUSE
C                             (EPSABS.LE.0 AND
C                              EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5E-28)
C                             OR LIMIT.LT.1 OR LENW.LT.LIMIT*46 OR
C                                LENIW .LT. LIMIT*3
C                             RESULT, ABSERR, LAST ARE SET TO
C                             ZERO.EXCEPT WHEN LIMIT OR LENW OR LENIW
C                             IS INVALID, IWORK(1), WORK(LIMIT*2+1) AND
C                             WORK(LIMIT*3+1) ARE SET TO ZERO, WORK(1)
C                             IS SET TO A AND WORK(LIMIT+1) TO B.
C
C         DIMENSIONING PARAMETERS
C            LIMIT - INTEGER
C                    LIMIT DETERMINES THE MAXIMUM NUMBER OF SUBINTERVALS
C                    IN THE PARTITION OF THE GIVEN INTEGRATION INTERVAL
C                    (A,B), LIMIT.GE.1.
C                    IF LIMIT.LT.1, THE ROUTINE WILL END WITH IER = 6.
C
C            LENW  - INTEGER
C                    DIMENSIONING PARAMETER FOR WORK
C                    LENW MUST BE AT LEAST LIMIT*46
C                    IF LENW.LT.LIMIT*46, THE ROUTINE WILL END
C                    WITH IER = 6.
C
C            LENIW - INTEGER
C                    DIMENSIONING PARAMETER FOR IWORK
C                    LENIW MUST BE AT LEAST LIMIT*3
C                    IF LENW.LT.LIMIT*3, THE ROUTINE WILL END
C                    WITH IER = 6.
C
C            LAST  - INTEGER
C                    ON RETURN, LAST EQUALS THE NUMBER OF SUBINTERVALS
C                    PRODUCED IN THE SUBDIVISION PROCESS, DETERMINES THE
C                    NUMBER OF SIGNIFICANT ELEMENTS ACTUALLY IN THE WORK
C                    ARRAYS.
C
C         WORK ARRAYS
C            IWORK - INTEGER
C                    VECTOR OF DIMENSION AT LEAST 3*LIMIT, THE FIRST K
C                    ELEMENTS OF WHICH CONTAIN POINTERS
C                    TO THE ERROR ESTIMATES OVER THE SUBINTERVALS
C                    SUCH THAT WORK(LIMIT*3+IWORK(1)),... ,
C                    WORK(LIMIT*3+IWORK(K)) FORM A DECREASING
C                    SEQUENCE, WITH K = LAST IF LAST.LE.(LIMIT/2+2),
C                    AND K = LIMIT+1-LAST OTHERWISE
C
C            WORK  - REAL
C                    VECTOR OF DIMENSION AT LEAST LENW
C                    ON RETURN
C                    WORK(1), ..., WORK(LAST) CONTAIN THE LEFT
C                     END-POINTS OF THE SUBINTERVALS IN THE
C                     PARTITION OF (A,B),
C                    WORK(LIMIT+1), ..., WORK(LIMIT+LAST) CONTAIN
C                     THE RIGHT END-POINTS,
C                    WORK(LIMIT*2+1), ..., WORK(LIMIT*2+LAST) CONTAIN
C                     THE INTEGRAL APPROXIMATIONS OVER THE SUBINTERVALS,
C                    WORK(LIMIT*3+1), ..., WORK(LIMIT*3+LAST)
C                     CONTAIN THE ERROR ESTIMATES.
C                    WORK(LIMIT*4+1), ... IS THE AREA RESERVED TO STORE
C                     FUNCTIONAL VALUES.
C***ROUTINES CALLED  QXGE,XERROR
C
C THE CALL TO XERROR HAS BEEN COMMENTED TO GET A STAND ALONE ROUTINE
C IT CAN BE ADAPTED TO FIT THE LOCAL ERROR HANDLING PROCEDURES
C
C
      REAL A,ABSERR,B,EPSABS,EPSREL,F,RESULT,WORK
      INTEGER KEY,IER,IWORK,LAST,LENIW,LENW,LIMIT,LVL,L1,L2,L3,L4,L5
C
      DIMENSION IWORK(LENIW),WORK(LENW)
C
      EXTERNAL F
C
C         CHECK VALIDITY OF LIMIT,LENIW AND LENW.
C
C***FIRST EXECUTABLE STATEMENT  QXG
      IER = 6
      LAST = 0
      RESULT = 0.0E+00
      ABSERR = 0.0E+00
      IF(LIMIT.LT.1.OR.LENIW.LT.LIMIT*3.OR.LENW.LT.LIMIT*46) GO TO 10
C
C         PREPARE CALL FOR QXGE.
C
      L1 = LIMIT+1
      L2 = LIMIT+L1
      L3 = LIMIT+L2
      L4 = LIMIT+L3
      L5 = 21*LIMIT+L3
C
      CALL QXGE(F,A,B,EPSABS,EPSREL,KEY,LIMIT,RESULT,ABSERR,
     *  IER,WORK(1),WORK(L1),WORK(L2),WORK(L3),IWORK,LAST,
     *  WORK(L4),WORK(L5),IWORK(L1),IWORK(L2))
C
C         CALL ERROR HANDLER IF NECESSARY.
C
      LVL = 0
10    IF(IER.EQ.6) LVL = 1
C
C THE CALL TO XERROR HAS BEEN COMMENTED TO GET A STAND ALONE ROUTINE
C IT CAN BE ADAPTED TO FIT THE LOCAL ERROR HANDLING PROCEDURES
C
C     IF(IER.NE.0)CALL XERROR('ABNORMAL RETURN FROM QXG   ',28,IER,LVL)
      RETURN
      END
*DECK QXGE
      SUBROUTINE QXGE(F,A,B,EPSABS,EPSREL,KEY,LIMIT,RESULT,ABSERR,
     *   IER,ALIST,BLIST,RLIST,ELIST,IORD,LAST,VALP,VALN,LP,LN)
C
C
C            THE ROUTINE CALCULATES AN APPROXIMATION RESULT TO A GIVEN
C            DEFINITE INTEGRAL   I = INTEGRAL OF F OVER (A,B),
C            HOPEFULLY SATISFYING FOLLOWING CLAIM FOR ACCURACY
C            ABS(I-RESLT).LE.MAX(EPSABS,EPSREL*ABS(I)).
C
C        PARAMETERS
C         ON ENTRY
C            F      - REAL
C                     FUNCTION SUBPROGRAM DEFINING THE INTEGRAND
C                     FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE
C                     DECLARED E X T E R N A L IN THE DRIVER PROGRAM.
C
C            A      - REAL
C                     LOWER LIMIT OF INTEGRATION
C
C            B      - REAL
C                     UPPER LIMIT OF INTEGRATION
C
C            EPSABS - REAL
C                     ABSOLUTE ACCURACY REQUESTED
C
C            EPSREL - REAL
C                     RELATIVE ACCURACY REQUESTED
C                     IF  EPSABS.LE.0
C                     AND EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5E-28),
C                     THE ROUTINE WILL END WITH IER = 6.
C
C            KEY    - INTEGER
C                     KEY FOR CHOICE OF LOCAL INTEGRATION RULE
C                     RMS FORMULAS ARE USED WITH
C                      13 - 19               POINTS IF KEY.LT.1,
C                      13 - 19 - (27)        POINTS IF KEY = 1,
C                      13 - 19 - (27) - (41) POINTS IF KEY = 2,
C                           19 -  27  - (41) POINTS IF KEY = 3,
C                                 27  -  41  POINTS IF KEY.GT.3.
C
C                         (RULES) USED IF THE FUNCTION APPEARS
C                         ENOUGH REGULAR IN THE SUBINTERVAL
C
C            LIMIT  - INTEGER
C                     GIVES AN UPPERBOUND ON THE NUMBER OF SUBINTERVALS
C                     IN THE PARTITION OF (A,B), LIMIT.GE.1.
C
C         ON RETURN
C            RESULT - REAL
C                     APPROXIMATION TO THE INTEGRAL
C
C            ABSERR - REAL
C                     ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR,
C                     WHICH SHOULD EQUAL OR EXCEED ABS(I-RESULT)
C
C            IER    - INTEGER
C                     IER = 0 NORMAL AND RELIABLE TERMINATION OF THE
C                             ROUTINE. IT IS ASSUMED THAT THE REQUESTED
C                             ACCURACY HAS BEEN ACHIEVED.
C                     IER.GT.0 ABNORMAL TERMINATION OF THE ROUTINE
C                             THE ESTIMATES FOR RESULT AND ERROR ARE
C                             LESS RELIABLE. IT IS ASSUMED THAT THE
C                             REQUESTED ACCURACY HAS NOT BEEN ACHIEVED.
C            ERROR MESSAGES
C                     IER = 1 MAXIMUM NUMBER OF SUBDIVISIONS ALLOWED
C                             HAS BEEN ACHIEVED. ONE CAN ALLOW MORE
C                             SUBDIVISIONS BY INCREASING THE VALUE
C                             OF LIMIT.
C                             HOWEVER, IF THIS YIELDS NO IMPROVEMENT IT
C                             IS RATHER ADVISED TO ANALYZE THE INTEGRAND
C                             IN ORDER TO DETERMINE THE INTEGRATION
C                             DIFFICULTIES. IF THE POSITION OF A LOCAL
C                             DIFFICULTY CAN BE DETERMINED(E.G.
C                             SINGULARITY, DISCONTINUITY WITHIN THE
C                             INTERVAL) ONE WILL PROBABLY GAIN FROM
C                             SPLITTING UP THE INTERVAL AT THIS POINT
C                             AND CALLING THE INTEGRATOR ON THE
C                             SUBRANGES. IF POSSIBLE, AN APPROPRIATE
C                             SPECIAL-PURPOSE INTEGRATOR SHOULD BE USED
C                             WHICH IS DESIGNED FOR HANDLING THE TYPE OF
C                             DIFFICULTY INVOLVED.
C                         = 2 THE OCCURRENCE OF ROUNDOFF ERROR IS
C                             DETECTED, WHICH PREVENTS THE REQUESTED
C                             TOLERANCE FROM BEING ACHIEVED.
C                         = 3 EXTREMELY BAD INTEGRAND BEHAVIOUR OCCURS
C                             AT SOME POINTS OF THE INTEGRATION
C                             INTERVAL.
C                         = 6 THE INPUT IS INVALID, BECAUSE
C                             (EPSABS.LE.0 AND
C                              EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5E-28),
C                             RESULT, ABSERR, LAST, RLIST(1) ,
C                             ELIST(1) AND IORD(1) ARE SET TO ZERO.
C                             ALIST(1) AND BLIST(1) ARE SET TO A AND B
C                             RESPECTIVELY.
C
C            ALIST   - REAL
C                      VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
C                       LAST  ELEMENTS OF WHICH ARE THE LEFT
C                      END POINTS OF THE SUBINTERVALS IN THE PARTITION
C                      OF THE GIVEN INTEGRATION RANGE (A,B)
C
C            BLIST   - REAL
C                      VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
C                       LAST  ELEMENTS OF WHICH ARE THE RIGHT
C                      END POINTS OF THE SUBINTERVALS IN THE PARTITION
C                      OF THE GIVEN INTEGRATION RANGE (A,B)
C
C            RLIST   - REAL
C                      VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
C                       LAST  ELEMENTS OF WHICH ARE THE
C                      INTEGRAL APPROXIMATIONS ON THE SUBINTERVALS
C
C            ELIST   - REAL
C                      VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST
C                       LAST  ELEMENTS OF WHICH ARE THE MODULI OF THE
C                      ABSOLUTE ERROR ESTIMATES ON THE SUBINTERVALS
C
C            IORD    - INTEGER
C                      VECTOR OF DIMENSION AT LEAST LIMIT, THE FIRST K
C                      ELEMENTS OF WHICH ARE POINTERS TO THE
C                      ERROR ESTIMATES OVER THE SUBINTERVALS,
C                      SUCH THAT ELIST(IORD(1)), ...,
C                      ELIST(IORD(K)) FORM A DECREASING SEQUENCE,
C                      WITH K = LAST IF LAST.LE.(LIMIT/2+2), AND
C                      K = LIMIT+1-LAST OTHERWISE
C
C            LAST    - INTEGER
C                      NUMBER OF SUBINTERVALS ACTUALLY PRODUCED IN THE
C                      SUBDIVISION PROCESS
C
C            VALP   - REAL
C            VALN     ARRAYS OF DIMENSION AT LEAST (21,LIMIT) USED TO
C                     SAVE THE FUNCTIONAL VALUES
C
C            LP     - INTEGER
C            LN       VECTORS OF DIMENSION AT LEAST LIMIT, USED TO
C                     STORE THE ACTUAL NUMBER OF FUNCTIONAL VALUES
C                     SAVED IN THE CORRESPONDING COLUMN
C                     OF VALP,VALN
C
C***ROUTINES CALLED  F,R1MACH,QXLQM,QXRRD,QPSRT,QXCPY
C
      REAL A,ABSERR,ALIST,AREA,AREA1,AREA12,AREA2,A1,A2,B,
     *  BLIST,B1,B2,ABS,DEFABS,DEFAB1,DEFAB2,AMAX1,R1MACH,ELIST,EPMACH,
     *  EPSABS,EPSREL,ERRBND,ERRMAX,ERROR1,ERROR2,ERRO12,ERRSUM,F,
     *  RESABS,RESULT,RLIST,UFLOW,VALP,VALN,VP1,VP2,VN1,VN2
      INTEGER KEY,IER,IORD,IROFF1,IROFF2,K,LAST,LIMIT,MAXERR,
     *  NRMAX,LP,LN,LP1,LP2,LN1,LN2
C
      DIMENSION ALIST(LIMIT),BLIST(LIMIT),ELIST(LIMIT),IORD(LIMIT),
     *  RLIST(LIMIT),VALP(21,LIMIT),VALN(21,LIMIT),LP(LIMIT),LN(LIMIT),
     * VP1(21),VP2(21),VN1(21),VN2(21)
C
      EXTERNAL F
C
C            MACHINE DEPENDENT CONSTANTS
C            ---------------------------
C
C          EPMACH IS THE LARGEST RELATIVE SPACING.
C          UFLOW IS THE SMALLEST POSITIVE MAGNITUDE.
C
C***FIRST EXECUTABLE STATEMENT  QXGE
      EPMACH = R1MACH(4)
      UFLOW = R1MACH(1)
C
C           TEST ON VALIDITY OF PARAMETERS
C           ------------------------------
C
      IER = 0
      LAST = 0
      RESULT = 0.0E+00
      ABSERR = 0.0E+00
      ALIST(1) = A
      BLIST(1) = B
      RLIST(1) = 0.0E+00
      ELIST(1) = 0.0E+00
      IORD(1) = 0
      IF(EPSABS.LE.0.0E+00.AND.
     *  EPSREL.LT.AMAX1(0.5E+02*EPMACH,0.5E-28)) IER = 6
      IF(IER.EQ.6) GO TO 999
C
C           FIRST APPROXIMATION TO THE INTEGRAL
C           -----------------------------------
C
      LP(1)=1
      LN(1)=1
      VALP(1,1)=F((A+B)*0.5E0)
      VALN(1,1)=VALP(1,1)
      CALL QXLQM(F,A,B,RESULT,ABSERR,DEFABS,RESABS,
     *           VALP(1,1),VALN(1,1),LP(1),LN(1),KEY)
      LAST = 1
      RLIST(1) = RESULT
      ELIST(1) = ABSERR
      IORD(1) = 1
C
C           TEST ON ACCURACY.
C
      ERRBND = AMAX1(EPSABS,EPSREL*ABS(RESULT))
      IF(ABSERR.LE.0.5E+02*EPMACH*DEFABS.AND.ABSERR.GT.ERRBND) IER = 2
      IF(LIMIT.EQ.1) IER = 1
      IF(IER.NE.0.OR.(ABSERR.LE.ERRBND.AND.ABSERR.NE.RESABS)
     *  .OR.ABSERR.EQ.0.0E+00) GO TO 999
C
C           INITIALIZATION
C           --------------
C
C
      ERRMAX = ABSERR
      MAXERR = 1
      AREA = RESULT
      ERRSUM = ABSERR
      NRMAX = 1
      IROFF1 = 0
      IROFF2 = 0
C
C           MAIN DO-LOOP
C           ------------
C
      DO 30 LAST = 2,LIMIT
C
C           BISECT THE SUBINTERVAL WITH THE LARGEST ERROR ESTIMATE.
C
        A1 = ALIST(MAXERR)
        B1 = 0.5E+00*(ALIST(MAXERR)+BLIST(MAXERR))
        A2 = B1
        B2 = BLIST(MAXERR)
        CALL QXRRD(F,VALN(1,MAXERR),LN(MAXERR),B1,A1,VN1,VP1,LN1,LP1)
        CALL QXLQM(F,A1,B1,AREA1,ERROR1,RESABS,DEFAB1,VP1,VN1,LP1,LN1,
     *              KEY)
        CALL QXRRD(F,VALP(1,MAXERR),LP(MAXERR),A2,B2,VP2,VN2,LP2,LN2)
        CALL QXLQM(F,A2,B2,AREA2,ERROR2,RESABS,DEFAB2,VP2,VN2,LP2,LN2,
     *              KEY)
C
C           IMPROVE PREVIOUS APPROXIMATIONS TO INTEGRAL
C           AND ERROR AND TEST FOR ACCURACY.
C
        AREA12 = AREA1+AREA2
        ERRO12 = ERROR1+ERROR2
        ERRSUM = ERRSUM+ERRO12-ERRMAX
        AREA = AREA+AREA12-RLIST(MAXERR)
        IF(DEFAB1.EQ.ERROR1.OR.DEFAB2.EQ.ERROR2) GO TO 5
        IF(ABS(RLIST(MAXERR)-AREA12).LE.0.1E-04*ABS(AREA12)
     *  .AND.ERRO12.GE.0.99E+00*ERRMAX) IROFF1 = IROFF1+1
        IF(LAST.GT.10.AND.ERRO12.GT.ERRMAX) IROFF2 = IROFF2+1
    5   RLIST(MAXERR) = AREA1
        RLIST(LAST) = AREA2
        ERRBND = AMAX1(EPSABS,EPSREL*ABS(AREA))
        IF(ERRSUM.LE.ERRBND) GO TO 8
C
C           TEST FOR ROUNDOFF ERROR AND EVENTUALLY SET ERROR FLAG.
C
        IF(IROFF1.GE.6.OR.IROFF2.GE.20) IER = 2
C
C           SET ERROR FLAG IN THE CASE THAT THE NUMBER OF SUBINTERVALS
C           EQUALS LIMIT.
C
        IF(LAST.EQ.LIMIT) IER = 1
C
C           SET ERROR FLAG IN THE CASE OF BAD INTEGRAND BEHAVIOUR
C           AT A POINT OF THE INTEGRATION RANGE.
C
        IF(AMAX1(ABS(A1),ABS(B2)).LE.(0.1E+01+0.1E+03*
     *  EPMACH)*(ABS(A2)+0.1E+04*UFLOW)) IER = 3
C
C           APPEND THE NEWLY-CREATED INTERVALS TO THE LIST.
C
    8   IF(ERROR2.GT.ERROR1) GO TO 10
        ALIST(LAST) = A2
        BLIST(MAXERR) = B1
        BLIST(LAST) = B2
        ELIST(MAXERR) = ERROR1
        ELIST(LAST) = ERROR2
        CALL QXCPY(VALP(1,MAXERR),VP1,LP1)
        LP(MAXERR)=LP1
        CALL QXCPY(VALN(1,MAXERR),VN1,LN1)
        LN(MAXERR)=LN1
        CALL QXCPY(VALP(1,LAST),VP2,LP2)
        LP(LAST)=LP2
        CALL QXCPY(VALN(1,LAST),VN2,LN2)
        LN(LAST)=LN2
        GO TO 20
   10   ALIST(MAXERR) = A2
        ALIST(LAST) = A1
        BLIST(LAST) = B1
        RLIST(MAXERR) = AREA2
        RLIST(LAST) = AREA1
        ELIST(MAXERR) = ERROR2
        ELIST(LAST) = ERROR1
        CALL QXCPY(VALP(1,MAXERR),VP2,LP2)
        LP(MAXERR)=LP2
        CALL QXCPY(VALN(1,MAXERR),VN2,LN2)
        LN(MAXERR)=LN2
        CALL QXCPY(VALP(1,LAST),VP1,LP1)
        LP(LAST)=LP1
        CALL QXCPY(VALN(1,LAST),VN1,LN1)
        LN(LAST)=LN1
C
C           CALL SUBROUTINE QPSRT TO MAINTAIN THE DESCENDING ORDERING
C           IN THE LIST OF ERROR ESTIMATES AND SELECT THE SUBINTERVAL
C           WITH THE LARGEST ERROR ESTIMATE (TO BE BISECTED NEXT).
C
   20   CALL QPSRT(LIMIT,LAST,MAXERR,ERRMAX,ELIST,IORD,NRMAX)
C ***JUMP OUT OF DO-LOOP
        IF(IER.NE.0.OR.ERRSUM.LE.ERRBND) GO TO 40
   30 CONTINUE
C
C           COMPUTE FINAL RESULT.
C           ---------------------
C
   40 RESULT = 0.0E+00
      DO 50 K=1,LAST
        RESULT = RESULT+RLIST(K)
   50 CONTINUE
      ABSERR = ERRSUM
  999 RETURN
      END
*DECK QXLQM
      SUBROUTINE QXLQM(F,A,B,RESULT,ABSERR,RESABS,RESASC,VR,VS,LR,LS,
     *                  KEY)
C
C            TO COMPUTE I = INTEGRAL OF F OVER (A,B), WITH ERROR
C                           ESTIMATE
C                       J = INTEGRAL OF ABS(F) OVER (A,B)
C
C           PARAMETERS
C            ON ENTRY
C              F      - REAL
C                       FUNCTION SUBPROGRAM DEFINING THE INTEGRAND
C                       FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE
C                       DECLARED E X T E R N A L IN THE DRIVER PROGRAM.
C
C              A      - REAL
C                       LOWER LIMIT OF INTEGRATION
C
C              B      - REAL
C                       UPPER LIMIT OF INTEGRATION
C
C              VR     - REAL
C                       VECTOR OF LENGTH LR CONTAINING THE
C                       SAVED  FUNCTIONAL VALUES OF POSITIVE ABSCISSAS
C
C              VS     - REAL
C                       VECTOR OF LENGTH LS CONTAINING THE
C                       SAVED  FUNCTIONAL VALUES OF NEGATIVE ABSCISSAS
C
C              LR     - INTEGER
C              LS       NUMBER OF ELEMENTS IN
C                       VR,VS RESPECTIVELY
C
C            KEY    - INTEGER
C                     KEY FOR CHOICE OF LOCAL INTEGRATION RULE
C                     RMS FORMULAS ARE USED WITH
C                      13 - 19               POINTS IF KEY.LT.1,
C                      13 - 19 - (27)        POINTS IF KEY = 1,
C                      13 - 19 - (27) - (41) POINTS IF KEY = 2,
C                           19 -  27  - (41) POINTS IF KEY = 3,
C                                 27  -  41  POINTS IF KEY.GT.3.
C
C                         (RULES) USED IF THE FUNCTION APPEARS
C                         ENOUGH REGULAR
C
C            ON RETURN
C              RESULT - REAL
C                       APPROXIMATION TO THE INTEGRAL I
C
C              ABSERR - REAL
C                       ESTIMATE OF THE MODULUS OF THE ABSOLUTE ERROR,
C                       WHICH SHOULD NOT EXCEED ABS(I-RESULT)
C
C              RESABS - REAL
C                       APPROXIMATION TO THE INTEGRAL J
C
C              RESASC - REAL
C                       APPROXIMATION TO THE INTEGRAL OF ABS(F-I/(B-A))
C                       OVER (A,B)
C
C              VR     - REAL
C                       VECTOR OF LENGTH LR CONTAINING THE
C                       SAVED  FUNCTIONAL VALUES OF POSITIVE ABSCISSAS
C
C              VS     - REAL
C                       VECTOR OF LENGTH LS CONTAINING THE
C                       SAVED  FUNCTIONAL VALUES OF NEGATIVE ABSCISSAS
C
C              LR     - INTEGER
C              LS       NUMBER OF ELEMENTS IN
C                       VR,VS RESPECTIVELY
C
C***ROUTINES CALLED  R1MACH,QXRUL
C
      REAL F,A,B,RESULT,ABSERR,RESABS,RESASC,
     *          R1MACH,EPMACH,RESG,RESK,UFLOW,ERROLD,VR(21),VS(21)
      INTEGER K,K0,K1,K2,KEY,KEY1,LR,LS
      EXTERNAL F
C
C            MACHINE DEPENDENT CONSTANTS
C            ---------------------------
C
C          EPMACH IS THE LARGEST RELATIVE SPACING.
C          UFLOW IS THE SMALLEST POSITIVE MAGNITUDE.
C          ERROLD IS THE LARGEST MAGNITUDE.
C
C***FIRST EXECUTABLE STATEMENT QXLQM
      EPMACH = R1MACH(4)
      UFLOW  = R1MACH(1)
      ERROLD = R1MACH(2)
C
      KEY1 = MAX(KEY ,  0)
      KEY1 = MIN(KEY1,  4)
      K0   = MAX(KEY1-2,0)
      K1   = K0+1
      K2   = MIN(KEY1+1,3)
C
      CALL QXRUL(F,A,B,RESG,RESABS,RESASC,K0,K1,VR,VS,LR,LS)
      DO 99 K=K1,K2
        CALL QXRUL(F,A,B,RESK,RESABS,RESASC,K,K1,VR,VS,LR,LS)
        RESULT=RESK
        ABSERR = ABS((RESK-RESG))
        IF(RESASC.NE.0.0E+00.AND.ABSERR.NE.0.0E+00)
     *  ABSERR = RESASC*AMIN1(1.E0,(200E0*ABSERR/RESASC)**1.5E+00)
        IF(RESABS.GT.UFLOW/(10E0*EPMACH)) ABSERR = AMAX1
     *  ((EPMACH*10E0)*RESABS,ABSERR)
        RESG=RESK
        IF(ABSERR.GT.ERROLD*0.16)GOTO 3000
        IF(ABSERR.LT.1000*EPMACH*RESABS)GOTO 3000
        ERROLD=ABSERR
99    CONTINUE
3000  CONTINUE
      RETURN
      END
*DECK QXRUL
      SUBROUTINE QXRUL(F,XL,XU,Y,YA,YM,KE,K1,FV1,FV2,L1,L2)
C
C            TO COMPUTE I = INTEGRAL OF F OVER (A,B), WITH ERROR
C                           ESTIMATE
C            AND CONDITIONALLY COMPUTE
C                       J = INTEGRAL OF ABS(F) OVER (A,B)
C                       BY USING AN  RMS RULE
C           PARAMETERS
C            ON ENTRY
C              F      - REAL
C                       FUNCTION SUBPROGRAM DEFINING THE INTEGRAND
C                       FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE
C                       DECLARED E X T E R N A L IN THE DRIVER PROGRAM.
C
C              XL     - REAL
C                       LOWER LIMIT OF INTEGRATION
C
C              XU     - REAL
C                       UPPER LIMIT OF INTEGRATION
C
C
C              KE     - INTEGER
C                     KEY FOR CHOICE OF LOCAL INTEGRATION RULE
C                     AN RMS RULE IS USED WITH
C                         13      POINTS IF KE  = 2,
C                         19      POINTS IF KE  = 3,
C                         27      POINTS IF KE  = 4,
C                         42      POINTS IF KE  = 5
C
C              K1     INTEGER
C                     VALUE OF KEY FOR WHICH THE ADDITIONAL ESTIMATES
C                     YA, YM ARE TO BE COMPUTED
C
C              FV1    - REAL
C                       VECTOR CONTAINING L1
C                       SAVED  FUNCTIONAL VALUES OF POSITIVE ABSCISSAS
C
C              FV2    - REAL
C                       VECTOR CONTAINING L2
C                       SAVED  FUNCTIONAL VALUES OF NEGATIVE ABSCISSAS
C
C              L1     - INTEGER
C              L2       NUMBER OF ELEMENTS IN FV1,FV2  RESPECTIVELY
C
C            ON RETURN
C              Y      - REAL
C                       APPROXIMATION TO THE INTEGRAL I
C                       RESULT IS COMPUTED BY APPLYING THE
C                       REQUESTED RMS RULE
C
C              YA     - REAL
C                       IF KEY = K1  APPROXIMATION TO THE INTEGRAL J
C                       ELSE UNCHANGED
C
C              YM     - REAL
C                       IF KEY = K1  APPROXIMATION TO THE INTEGRAL OF
C                                      ABS(F-I/(XU-XL)   OVER (XL,XU)
C                       ELSE UNCHANGED
C
C              FV1    - REAL
C                       VECTOR L1 CONTAINING L1
C                       SAVED  FUNCTIONAL VALUES OF POSITIVE ABSCISSAS
C
C              FV2    - REAL
C                       VECTOR CONTAINING L2
C                       SAVED  FUNCTIONAL VALUES OF NEGATIVE ABSCISSAS
C
C              L1     - INTEGER
C              L2       NUMBER OF ELEMENTS IN FV1,FV2  RESPECTIVELY
C
C***ROUTINES CALLED  F
C
      REAL F,XL,XU,LDL,Y,YA,YM,Y2,XX(41),WW(52),
     *                 FV1(21),FV2(21),AA,BB,C
      EXTERNAL F
      INTEGER ISTART(4),LEN(4),KE,K1,L1,L2
      SAVE ISTART,LEN,XX,WW
      DATA ISTART/0, 7, 17, 31/,LEN/7, 10, 14, 21/
      DATA XX(  1)/.0                       /
      DATA XX(  2)/.25000000000000000000E+00/
      DATA XX(  3)/.50000000000000000000E+00/
      DATA XX(  4)/.75000000000000000000E+00/
      DATA XX(  5)/.87500000000000000000E+00/
      DATA XX(  6)/.93750000000000000000E+00/
      DATA XX(  7)/.10000000000000000000E+01/
      DATA XX(  8)/.37500000000000000000E+00/
      DATA XX(  9)/.62500000000000000000E+00/
      DATA XX( 10)/.96875000000000000000E+00/
      DATA XX( 11)/.12500000000000000000E+00/
      DATA XX( 12)/.68750000000000000000E+00/
      DATA XX( 13)/.81250000000000000000E+00/
      DATA XX( 14)/.98437500000000000000E+00/
      DATA XX( 15)/.18750000000000000000E+00/
      DATA XX( 16)/.31250000000000000000E+00/
      DATA XX( 17)/.43750000000000000000E+00/
      DATA XX( 18)/.56250000000000000000E+00/
      DATA XX( 19)/.84375000000000000000E+00/
      DATA XX( 20)/.90625000000000000000E+00/
      DATA XX( 21)/.99218750000000000000E+00/
C   NUMBER OF NODES 13
      DATA WW(1)/1.303262173284849021810473057638590518409112513421E-1/
      DATA WW(2)/2.390632866847646220320329836544615917290026806242E-1/
      DATA WW(3)/2.630626354774670227333506083741355715758124943143E-1/
      DATA WW(4)/2.186819313830574175167853094864355208948886875898E-1/
      DATA WW(5)/2.757897646642836865859601197607471574336674206700E-2/
      DATA WW(6)/1.055750100538458443365034879086669791305550493830E-1/
      DATA WW(7)/1.571194260595182254168429283636656908546309467968E-2/
C   NUMBER OF NODES 19
      DATA WW(8)/1.298751627936015783241173611320651866834051160074E-1/
      DATA WW(9)/2.249996826462523640447834514709508786970828213187E-1/
      DATA WW(15)/5.542699233295875168406783695143646338274805359780E-2/
      DATA WW(10)/1.680415725925575286319046726692683040162290325505E-1/
      DATA WW(16)/9.986735247403367525720377847755415293097913496236E-2/
      DATA WW(11)/1.415567675701225879892811622832845252125600939627E-1/
      DATA WW(12)/1.006482260551160175038684459742336605269707889822E-1/
      DATA WW(13)/2.510604860724282479058338820428989444699235030871E-2/
      DATA WW(17)/4.507523056810492466415880450799432587809828791196E-2/
      DATA WW(14)/9.402964360009747110031098328922608224934320397592E-3/
C   NUMBER OF NODES 27
      DATA WW(18)/6.300942249647773931746170540321811473310938661469E-2/
      DATA WW(28)/1.239572396231834242194189674243818619042280816640E-1/
      DATA WW(19)/1.261383225537664703012999637242003647020326905948E-1/
      DATA WW(25)/1.235837891364555000245004813294817451524633100256E-1/
      DATA WW(20)/1.273864433581028272878709981850307363453523117880E-1/
      DATA WW(26)/1.148933497158144016800199601785309838604146040215E-1/
      DATA WW(29)/2.501306413750310579525950767549691151739047969345D-2/
      DATA WW(21)/8.576500414311820514214087864326799153427368592787E-2/
      DATA WW(30)/4.915957918146130094258849161350510503556792927578E-2/
      DATA WW(22)/7.102884842310253397447305465997026228407227220665E-2/
      DATA WW(23)/5.026383572857942403759829860675892897279675661654E-2/
      DATA WW(27)/1.252575774226122633391477702593585307254527198070E-2/
      DATA WW(31)/2.259167374956474713302030584548274729936249753832E-2/
      DATA WW(24)/4.683670010609093810432609684738393586390722052124E-3/
C   NUMBER OF NODES 41
      DATA WW(32)/6.362762978782724559269342300509058175967124446839E-2/
      DATA WW(42)/1.187141856692283347609436153545356484256869129472E-1/
      DATA WW(46)/1.533126874056586959338368742803997744815413565014E-2/
      DATA WW(33)/9.950065827346794643193261975720606296171462239514E-2/
      DATA WW(47)/3.527159369750123100455704702965541866345781113903E-2/
      DATA WW(39)/8.140326425945938045967829319725797511040878579808E-2/
      DATA WW(48)/5.000556431653955124212795201196389006184693561679E-2/
      DATA WW(34)/7.048220002718565366098742295389607994441704889441E-2/
      DATA WW(49)/5.744164831179720106340717579281831675999717767532E-2/
      DATA WW(40)/6.583213447600552906273539578430361199084485578379E-2/
      DATA WW(43)/5.999947605385971985589674757013565610751028128731E-2/
      DATA WW(35)/6.512297339398335645872697307762912795346716454337E-2/
      DATA WW(44)/5.500937980198041736910257988346101839062581489820E-2/
      DATA WW(50)/1.598823797283813438301248206397233634639162043386E-2/
      DATA WW(36)/3.998229150313659724790527138690215186863915308702E-2/
      DATA WW(51)/2.635660410220884993472478832884065450876913559421E-2/
      DATA WW(37)/3.456512257080287509832054272964315588028252136044E-2/
      DATA WW(41)/2.592913726450792546064232192976262988065252032902E-2/
      DATA WW(45)/5.264422421764655969760271538981443718440340270116E-3/
      DATA WW(52)/1.196003937945541091670106760660561117114584656319E-2/
      DATA WW(38)/2.212167975884114432760321569298651047876071264944E-3/
C
C***FIRST EXECUTABLE STATEMENT QXRUL
      K=KE+1
      IS=ISTART(K)
      KS=LEN(K)
      LDL= XU-XL
      BB= LDL*0.5E0
      AA= XL+BB
      Y =0.E0
      DO 100 I=1,KS
         IF(I.GT.L1.OR.I.GT.L2)   C=BB*XX(I)
         IF(I.GT.L1)              FV1(I)=F(AA+C)
         IF(I.GT.L2)              FV2(I)=F(AA-C)
100      Y=Y+(FV1(I)+FV2(I))*WW(IS+I)
      Y2=Y
      Y=Y*BB
      IF(L1.LT.KS)L1=KS
      IF(L2.LT.KS)L2=KS
      IF(KE.NE.K1)RETURN
      YA=0.E0
      DO 25 I=1,KS
  25    YA=YA+(ABS(FV1(I))+ABS(FV2(I)))*WW(IS+I)
      Y2=Y2*0.5E0
      YM=0.E0
      YA=YA*ABS(BB)
      DO 27 I=1,KS
   27 YM=YM+(ABS(FV1(I)-Y2)+ABS(FV2(I)-Y2))*WW(IS+I)
      YM=YM*ABS(BB)
      RETURN
      END
*DECK QXRRD
      SUBROUTINE QXRRD(F,Z,LZ,XL,XU,R,S,LR,LS)
C
C            TO REORDER THE COMPUTED FUNCTIONAL VALUES BEFORE
C            THE BISECTION OF AN INTERVAL
C           PARAMETERS
C            ON ENTRY
C              F      - REAL
C                       FUNCTION SUBPROGRAM DEFINING THE INTEGRAND
C                       FUNCTION F(X). THE ACTUAL NAME FOR F NEEDS TO BE
C                       DECLARED E X T E R N A L IN THE DRIVER PROGRAM.
C
C              XL     - REAL
C                       LOWER LIMIT OF INTERVAL
C
C              XU     - REAL
C                       UPPER LIMIT OF INTERVAL
C
C              Z      - REAL
C                       VECTOR CONTAINING LZ
C                       SAVED  FUNCTIONAL VALUES
C
C              LZ     - INTEGER
C                       NUMBER OF ELEMENTS IN LZ
C
C
C            ON RETURN
C              R      - REAL
C              S        VECTORS CONTAINING LR, LS
C                       SAVED  FUNCTIONAL VALUES FOR THE NEW INTERVALS
C
C              LR     - INTEGER
C              LS       NUMBER OF ELEMENTES IN R,S RESPECTIVELY
C
C***ROUTINES CALLED  F
C
      REAL F,R,S,Z,XU,XL,DLEN,CENTR
      INTEGER LR,LS,LZ
      DIMENSION R(21),S(21),Z(21)
C***FIRST EXECUTABLE STATEMENT QXRRD
      DLEN= (XU-XL)*0.5E0
      CENTR= XL+DLEN
      R(1)=  Z(3)
      R(2)=  Z(9)
      R(3)=  Z(4)
      R(4)=  Z(5)
      R(5)=  Z(6)
      R(6)=  Z(10)
      R(7)=  Z(7)
      S(1)=  Z(3)
      S(2)=  Z(8)
      S(3)=  Z(2)
      S(7)=  Z(1)
      IF(LZ.GT.11)GOTO 2
      R(8)=  F(CENTR+DLEN*0.37500000E0)
      R(9)=  F(CENTR+DLEN*0.62500000E0)
      R(10)=  F(CENTR+DLEN*0.96875000E0)
      LR=  10
      IF(LZ.NE.11)S(4)=  F(CENTR-DLEN*0.75000000E0)
      IF(LZ.EQ.11)S(4)=  Z(11)
      S(5)=  F(CENTR-DLEN*0.87500000E0)
      S(6)=  F(CENTR-DLEN*0.93750000E0)
      S(8)=  F(CENTR-DLEN*0.37500000E0)
      S(9)=  F(CENTR-DLEN*0.62500000E0)
      S(10)=  F(CENTR-DLEN*0.96875000E0)
      LS=  10
      GOTO 3000
2     R(8)= Z(12)
      R(9)= Z(13)
      R(10)= Z(14)
      LR=  10
      S(4)= Z(11)
      S(5)= F(CENTR-DLEN*0.87500000E0)
      S(6)= F(CENTR-DLEN*0.93750000E0)
      IF(LZ.GT.14)GOTO3
      S(8)= F(CENTR-DLEN*0.37500000E0)
      S(9)= F(CENTR-DLEN*0.62500000E0)
      S(10)= F(CENTR-DLEN*0.96875000E0)
      LS=  10
      GOTO 3000
3     R(11)= Z(18)
      R(12)= Z(19)
      R(13)= Z(20)
      R(14)= Z(21)
      LR=  14
      S(8)= Z(16)
      S(9)= Z(15)
      S(10)= F(CENTR-DLEN*0.96875000E0)
      S(11)= Z(17)
      LS=  11
3000  RETURN
      END
*DECK QXCPY
      SUBROUTINE QXCPY(A,B,L)
C
C  TO COPY THE REAL VECTOR B OF LENGTH L   I N T O
C          THE REAL VECTOR A OF LENGTH L
C
C***REMARK  THIS ROUTINE CAN BE IMPROVED, BY CODING IT IN THE
C           ASSEMBLER LANGUAGE OF THE USED MACHINE
C***ROUTINES CALLED  (NONE)
C
      INTEGER L
      REAL A(L),B(L)
C***FIRST EXECUTABLE STATEMENT QXCPY
      DO 1 I=1,L
1     A(I)=B(I)
      RETURN
      END
*DECK MAIN *** TEST PROGRAM ***
C  THIS DRIVER CALLS VARIOUS INTEGRATORS WITH  VARIOUS RELATIVE
C  TOLERACES TO INTEGRATE THE FUNCTION
C     F(X) = SQRT(X/(1-X)) * LOG(X)   OVER (0, 1)
C
C THE EXACT INTEGRAL IS                           -0.606789763508705...
C
      REAL F,A,B,ATOL,RTOL,TEE,Y,WORK(4600)
      INTEGER NCALLS,IER,NEVAL,LAST,IWORK(300),KKK,KEY,LIMIT,LENW,LENIW
      EXTERNAL F
      COMMON /CTEST/NCALLS
      DATA LIMIT/100/,LENIW/300/,LENW/4600/
      WRITE(6,600)
600   FORMAT(32H1 COMPUTATION OF THE INTEGRAL OF,
     *      /38H       F(X) = SQRT(X/(1-X)) * LOG(X),
     *      /13H  OVER (0, 1),
     *      /23H  THE EXACT INTEGRAL IS,30X,
     *       22H -0.606789763508705...)
      A=0.0
      B=1.0
      WRITE(6,61)
 61   FORMAT(//22H TEST OF QXGS                //
     *  53H RELATIVE TOLERANCE FUNCTION CALLS ERROR ESTIMATE IER,
     *  13H       RESULT )
      DO 1 I=1,4
         ATOL=1E-8
         RTOL=10E0**(-I)
         NCALLS=0
         CALL  QXGS (F,A,B,ATOL,RTOL,Y,TEE,IER,LIMIT,
     *               LENIW,LENW,LAST,IWORK,WORK)
         WRITE(6,60)I,NCALLS,TEE,IER,Y
60       FORMAT(4X,6H10**(-,I2,1H),I17,9X,E7.2,4X,I2,F20.8)
 1    CONTINUE
      DO 2 KKK = 1,5
       KEY = KKK - 1
       WRITE(6,62) KEY
 62    FORMAT(//23H TEST OF  QXG , KEY =    ,I4   //
     *  53H RELATIVE TOLERANCE FUNCTION CALLS ERROR ESTIMATE IER,
     *  13H       RESULT )
      DO 2 I=1,4
         ATOL=1E-8
         RTOL=10E0**(-I)
         NCALLS=0
         CALL  QXG  (F,A,B,ATOL,RTOL,KEY,Y,TEE,IER,LIMIT,
     *               LENIW,LENW,LAST,IWORK,WORK)
         WRITE(6,60)I,NCALLS,TEE,IER,Y
 2    CONTINUE
      STOP
      END
*DECK F
      REAL FUNCTION F(X)
      REAL X
C        INTEGRAND FUNCTION, IN THE INTERVAL (0,1) THE EXACT INTEGRAL
C   TO 15 DIGITS IS
C          -0.60678 97635 08705D0
C
      COMMON /CTEST/NCALLS
      NCALLS=NCALLS+1
      IF (X.EQ.0.0) GO TO  10
      IF (X.EQ.1.0) GO TO  10
      F=SQRT(X/(1-X))*LOG(X)
      RETURN
10    F=0.0
      RETURN
      END
