      SUBROUTINE DQAWSE(F,A,B,ALFA,BETA,INTEGR,EPSABS,EPSREL,LIMIT,
     1   RESULT,ABSERR,NEVAL,IER,ALIST,BLIST,RLIST,ELIST,IORD,LAST)
C***BEGIN PROLOGUE  DQAWSE
C***DATE WRITTEN   800101   (YYMMDD)
C***REVISION DATE  830518   (YYMMDD)
C***REVISION HISTORY (YYMMDD)
C   000601   Changed DMAX1/DABS to generic MAX/ABS
C***CATEGORY NO.  H2A2A1
C***KEYWORDS  ALGEBRAICO-LOGARITHMIC END POINT SINGULARITIES,
C             AUTOMATIC INTEGRATOR,CLENSHAW-CURTIS METHOD,
C             SPECIAL-PURPOSE
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 calculates an approximation result to a given
C            definite integral I = Integral of F*W over (A,B),
C            (where W shows a singular behaviour at the end points,
C            see parameter INTEGR).
C            Hopefully satisfying following claim for accuracy
C            ABS(I-RESULT).LE.MAX(EPSABS,EPSREL*ABS(I)).
C***DESCRIPTION
C
C        Integration of functions having algebraico-logarithmic
C        end point singularities
C        Standard fortran subroutine
C        Double precision version
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, B.GT.A
C                     If B.LE.A, the routine will end with IER = 6.
C
C            ALFA   - Double precision
C                     Parameter in the WEIGHT function, ALFA.GT.(-1)
C                     If ALFA.LE.(-1), the routine will end with
C                     IER = 6.
C
C            BETA   - Double precision
C                     Parameter in the WEIGHT function, BETA.GT.(-1)
C                     If BETA.LE.(-1), the routine will end with
C                     IER = 6.
C
C            INTEGR - Integer
C                     Indicates which WEIGHT function is to be used
C                     = 1  (X-A)**ALFA*(B-X)**BETA
C                     = 2  (X-A)**ALFA*(B-X)**BETA*LOG(X-A)
C                     = 3  (X-A)**ALFA*(B-X)**BETA*LOG(B-X)
C                     = 4  (X-A)**ALFA*(B-X)**BETA*LOG(X-A)*LOG(B-X)
C                     If INTEGR.LT.1 or INTEGR.GT.4, the routine
C                     will end with IER = 6.
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 upper bound on the number of subintervals
C                     in the partition of (A,B), LIMIT.GE.2
C                     If LIMIT.LT.2, 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            NEVAL  - Integer
C                     Number of integrand evaluations
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 the integral and error
C                             are 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
C                             subdivisions by increasing the value of
C                             LIMIT. However, if this yields no
C                             improvement, it is advised to analyze the
C                             integrand in order to determine the
C                             integration difficulties which prevent the
C                             requested tolerance from being achieved.
C                             In case of a jump DISCONTINUITY or a local
C                             SINGULARITY of algebraico-logarithmic type
C                             at one or more interior points of the
C                             integration range, one should proceed by
C                             splitting up the interval at these
C                             points and calling the integrator on the
C                             subranges.
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                             B.LE.A or ALFA.LE.(-1) or BETA.LE.(-1), or
C                             INTEGR.LT.1 or INTEGR.GT.4, or
C                             (EPSABS.LE.0 and
C                              EPSREL.LT.MAX(50*REL.MACH.ACC.,0.5D-28),
C                             or LIMIT.LT.2.
C                             RESULT, ABSERR, NEVAL, RLIST(1), ELIST(1),
C                             IORD(1) and LAST are set to zero. ALIST(1)
C                             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 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                     of which are pointers to the error
C                     estimates over the subintervals, so that
C                     ELIST(IORD(1)), ..., ELIST(IORD(K)) with K = LAST
C                     If LAST.LE.(LIMIT/2+2), and K = LIMIT+1-LAST
C                     otherwise form a decreasing sequence
C
C            LAST   - Integer
C                     Number of subintervals actually produced in
C                     the subdivision process
C***REFERENCES  (NONE)
C***ROUTINES CALLED  D1MACH,DQC25S,DQMOMO,DQPSRT
C***END PROLOGUE  DQAWSE
C
      DOUBLE PRECISION A,ABSERR,ALFA,ALIST,AREA,AREA1,AREA12,AREA2,A1,
     1  A2,B,BETA,BLIST,B1,B2,CENTRE,D1MACH,ELIST,EPMACH,
     2  EPSABS,EPSREL,ERRBND,ERRMAX,ERROR1,ERRO12,ERROR2,ERRSUM,F,
     3  RESAS1,RESAS2,RESULT,RG,RH,RI,RJ,RLIST,UFLOW
      INTEGER IER,INTEGR,IORD,IROFF1,IROFF2,K,LAST,LIMIT,MAXERR,NEV,
     1  NEVAL,NRMAX
C
      EXTERNAL F
C
      DIMENSION ALIST(LIMIT),BLIST(LIMIT),RLIST(LIMIT),ELIST(LIMIT),
     1  IORD(LIMIT),RI(25),RJ(25),RH(25),RG(25)
C
C            LIST OF MAJOR VARIABLES
C            -----------------------
C
C           ALIST     - LIST OF LEFT END POINTS OF ALL SUBINTERVALS
C                       CONSIDERED UP TO NOW
C           BLIST     - LIST OF RIGHT END POINTS OF ALL SUBINTERVALS
C                       CONSIDERED UP TO NOW
C           RLIST(I)  - APPROXIMATION TO THE INTEGRAL OVER
C                       (ALIST(I),BLIST(I))
C           ELIST(I)  - ERROR ESTIMATE APPLYING TO RLIST(I)
C           MAXERR    - POINTER TO THE INTERVAL WITH LARGEST
C                       ERROR ESTIMATE
C           ERRMAX    - ELIST(MAXERR)
C           AREA      - SUM OF THE INTEGRALS OVER THE SUBINTERVALS
C           ERRSUM    - SUM OF THE ERRORS OVER THE SUBINTERVALS
C           ERRBND    - REQUESTED ACCURACY MAX(EPSABS,EPSREL*
C                       ABS(RESULT))
C           *****1    - VARIABLE FOR THE LEFT SUBINTERVAL
C           *****2    - VARIABLE FOR THE RIGHT SUBINTERVAL
C           LAST      - INDEX FOR SUBDIVISION
C
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  DQAWSE
      EPMACH = D1MACH(4)
      UFLOW = D1MACH(1)
C
C           TEST ON VALIDITY OF PARAMETERS
C           ------------------------------
C
      IER = 6
      NEVAL = 0
      LAST = 0
      RLIST(1) = 0.0D+00
      ELIST(1) = 0.0D+00
      IORD(1) = 0
      RESULT = 0.0D+00
      ABSERR = 0.0D+00
      IF(B.LE.A.OR.(EPSABS.EQ.0.0D+00.AND.
     1  EPSREL.LT.MAX(0.5D+02*EPMACH,0.5D-28)).OR.ALFA.LE.(-0.1D+01).
     2  OR.BETA.LE.(-0.1D+01).OR.INTEGR.LT.1.OR.INTEGR.GT.4.OR.
     3  LIMIT.LT.2) GO TO 999
      IER = 0
C
C           COMPUTE THE MODIFIED CHEBYSHEV MOMENTS.
C
      CALL DQMOMO(ALFA,BETA,RI,RJ,RG,RH,INTEGR)
C
C           INTEGRATE OVER THE INTERVALS (A,(A+B)/2) AND ((A+B)/2,B).
C
      CENTRE = 0.5D+00*(B+A)
      CALL DQC25S(F,A,B,A,CENTRE,ALFA,BETA,RI,RJ,RG,RH,AREA1,
     1  ERROR1,RESAS1,INTEGR,NEV)
      NEVAL = NEV
      CALL DQC25S(F,A,B,CENTRE,B,ALFA,BETA,RI,RJ,RG,RH,AREA2,
     1  ERROR2,RESAS2,INTEGR,NEV)
      LAST = 2
      NEVAL = NEVAL+NEV
      RESULT = AREA1+AREA2
      ABSERR = ERROR1+ERROR2
C
C           TEST ON ACCURACY.
C
      ERRBND = MAX(EPSABS,EPSREL*ABS(RESULT))
C
C           INITIALIZATION
C           --------------
C
      IF(ERROR2.GT.ERROR1) GO TO 10
      ALIST(1) = A
      ALIST(2) = CENTRE
      BLIST(1) = CENTRE
      BLIST(2) = B
      RLIST(1) = AREA1
      RLIST(2) = AREA2
      ELIST(1) = ERROR1
      ELIST(2) = ERROR2
      GO TO 20
   10 ALIST(1) = CENTRE
      ALIST(2) = A
      BLIST(1) = B
      BLIST(2) = CENTRE
      RLIST(1) = AREA2
      RLIST(2) = AREA1
      ELIST(1) = ERROR2
      ELIST(2) = ERROR1
   20 IORD(1) = 1
      IORD(2) = 2
      IF(LIMIT.EQ.2) IER = 1
      IF(ABSERR.LE.ERRBND.OR.IER.EQ.1) GO TO 999
      ERRMAX = ELIST(1)
      MAXERR = 1
      NRMAX = 1
      AREA = RESULT
      ERRSUM = ABSERR
      IROFF1 = 0
      IROFF2 = 0
C
C            MAIN DO-LOOP
C            ------------
C
      DO 60 LAST = 3,LIMIT
C
C           BISECT THE SUBINTERVAL WITH LARGEST ERROR ESTIMATE.
C
        A1 = ALIST(MAXERR)
        B1 = 0.5D+00*(ALIST(MAXERR)+BLIST(MAXERR))
        A2 = B1
        B2 = BLIST(MAXERR)
C
        CALL DQC25S(F,A,B,A1,B1,ALFA,BETA,RI,RJ,RG,RH,AREA1,
     1  ERROR1,RESAS1,INTEGR,NEV)
        NEVAL = NEVAL+NEV
        CALL DQC25S(F,A,B,A2,B2,ALFA,BETA,RI,RJ,RG,RH,AREA2,
     1  ERROR2,RESAS2,INTEGR,NEV)
        NEVAL = NEVAL+NEV
C
C           IMPROVE PREVIOUS APPROXIMATIONS INTEGRAL AND ERROR
C           AND TEST FOR ACCURACY.
C
        AREA12 = AREA1+AREA2
        ERRO12 = ERROR1+ERROR2
        ERRSUM = ERRSUM+ERRO12-ERRMAX
        AREA = AREA+AREA12-RLIST(MAXERR)
        IF(A.EQ.A1.OR.B.EQ.B2) GO TO 30
        IF(RESAS1.EQ.ERROR1.OR.RESAS2.EQ.ERROR2) GO TO 30
C
C           TEST FOR ROUNDOFF ERROR.
C
        IF(ABS(RLIST(MAXERR)-AREA12).LT.0.1D-04*ABS(AREA12)
     1  .AND.ERRO12.GE.0.99D+00*ERRMAX) IROFF1 = IROFF1+1
        IF(LAST.GT.10.AND.ERRO12.GT.ERRMAX) IROFF2 = IROFF2+1
   30   RLIST(MAXERR) = AREA1
        RLIST(LAST) = AREA2
C
C           TEST ON ACCURACY.
C
        ERRBND = MAX(EPSABS,EPSREL*ABS(AREA))
        IF(ERRSUM.LE.ERRBND) GO TO 35
C
C           SET ERROR FLAG IN THE CASE THAT THE NUMBER OF INTERVAL
C           BISECTIONS EXCEEDS LIMIT.
C
        IF(LAST.EQ.LIMIT) IER = 1
C
C
C           SET ERROR FLAG IN THE CASE OF ROUNDOFF ERROR.
C
        IF(IROFF1.GE.6.OR.IROFF2.GE.20) IER = 2
C
C           SET ERROR FLAG IN THE CASE OF BAD INTEGRAND BEHAVIOUR
C           AT INTERIOR POINTS OF INTEGRATION RANGE.
C
        IF(MAX(ABS(A1),ABS(B2)).LE.(0.1D+01+0.1D+03*EPMACH)*
     1  (ABS(A2)+0.1D+04*UFLOW)) IER = 3
C
C           APPEND THE NEWLY-CREATED INTERVALS TO THE LIST.
C
   35   IF(ERROR2.GT.ERROR1) GO TO 40
        ALIST(LAST) = A2
        BLIST(MAXERR) = B1
        BLIST(LAST) = B2
        ELIST(MAXERR) = ERROR1
        ELIST(LAST) = ERROR2
        GO TO 50
   40   ALIST(MAXERR) = A2
        ALIST(LAST) = A1
        BLIST(LAST) = B1
        RLIST(MAXERR) = AREA2
        RLIST(LAST) = AREA1
        ELIST(MAXERR) = ERROR2
        ELIST(LAST) = ERROR1
C
C           CALL SUBROUTINE DQPSRT TO MAINTAIN THE DESCENDING ORDERING
C           IN THE LIST OF ERROR ESTIMATES AND SELECT THE SUBINTERVAL
C           WITH LARGEST ERROR ESTIMATE (TO BE BISECTED NEXT).
C
   50   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 70
   60 CONTINUE
C
C           COMPUTE FINAL RESULT.
C           ---------------------
C
   70 RESULT = 0.0D+00
      DO 80 K=1,LAST
        RESULT = RESULT+RLIST(K)
   80 CONTINUE
      ABSERR = ERRSUM
  999 RETURN
      END
      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
      DOUBLE PRECISION FUNCTION D1MACH(I)
C***BEGIN PROLOGUE  D1MACH
C***DATE WRITTEN   750101   (YYMMDD)
C***REVISION DATE  880330   (YYMMDD)
C***CATEGORY NO.  R1
C***KEYWORDS  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     This is the CMLIB version of D1MACH, the double precision machine
C     constants subroutine originally developed for the PORT library.
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  Double-precision machine constants
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***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***END PROLOGUE  D1MACH
C
      INTEGER SMALL(4)
      INTEGER LARGE(4)
      INTEGER RIGHT(4)
      INTEGER DIVER(4)
      INTEGER LOG10(4)
C
      DOUBLE PRECISION DMACH(5)
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 IEEE ARITHMETIC MACHINES, SUCH AS THE AT&T
C     3B SERIES AND MOTOROLA 68000 BASED MACHINES (E.G. SUN 3 AND AT&T
C     PC 7300), SUN SPARCSTATIONS, SILICON GRAPHCS WORKSTATIONS, HP
C     9000 WORKSTATIONS, IBM RS/6000 WORKSTATIONS IN WHICH THE MOST 
C     SIGNIFICANT BYTE IS STORED FIRST.
C
C === MACHINE = ATT.3B
C === MACHINE = ATT.7300
C === MACHINE = HP.9000
C === MACHINE = IBM.RS6000
C === MACHINE = IEEE.MOST-SIG-BYTE-FIRST
C === MACHINE = SGI
C === MACHINE = SUN
C === MACHINE = 68000
C      DATA SMALL(1),SMALL(2) /    1048576,          0 /
C      DATA LARGE(1),LARGE(2) / 2146435071,         -1 /
C      DATA RIGHT(1),RIGHT(2) / 1017118720,          0 /
C      DATA DIVER(1),DIVER(2) / 1018167296,          0 /
C      DATA LOG10(1),LOG10(2) / 1070810131, 1352628735 /
C
C     MACHINE CONSTANTS FOR IEEE ARITHMETIC MACHINES AND 8087-BASED
C     MICROS, SUCH AS THE IBM PC AND AT&T 6300, IN WHICH THE LEAST
C     SIGNIFICANT BYTE IS STORED FIRST.
C
C === MACHINE = IEEE.LEAST-SIG-BYTE-FIRST
C === MACHINE = 8087
C === MACHINE = IBM.PC
C === MACHINE = ATT.6300
C      DATA SMALL(1),SMALL(2) /          0,    1048576 /
C      DATA LARGE(1),LARGE(2) /         -1, 2146435071 /
C      DATA RIGHT(1),RIGHT(2) /          0, 1017118720 /
C      DATA DIVER(1),DIVER(2) /          0, 1018167296 /
C      DATA LOG10(1),LOG10(2) / 1352628735, 1070810131 /
C
C     MACHINE CONSTANTS FOR SUN WORKSTATIONS.  f77 WITH -r8 OPTION.
C
C === MACHINE = SUN.F77-WITH-R8-OPTION
C      DATA DMACH(1) / 3.3621031431120935062626778173217526D-4932 /
C      DATA DMACH(2) / 1.1897314953572317650857593266280070D+4932 /
C      DATA DMACH(3) / 9.6296497219361792652798897129246366D-035  /
C      DATA DMACH(4) / 1.9259299443872358530559779425849273D-034  /
C      DATA DMACH(5) / 0.30102999566398119521373889472449302      /
C
C     MACHINE CONSTANTS FOR IBM RS/6000 WORKSTATIONS WITH -qautodbl=dblpad.
C
C === MACHINE = IBM.RS6000.XLF-WITH-AUTODBL-OPTION
C      DATA DMACH(1) / 2.2250738585072E-308 /
C      DATA DMACH(2) / 1.7976931348623E308  /
C      DATA DMACH(3) / 0.48148248609680896326399448564623183E-34  /
C      DATA DMACH(4) / 0.48148248609680896326399448564623183E-34  /
C      DATA DMACH(5) / 0.30102999566398119521373889472449302      /
C
C     MACHINE CONSTANTS FOR SGI Origin 2000 with -r8 -d16 options.
C
C === MACHINE = SGI.ORIGIN.F77-WITH-R8-D16-OPTIONS
C      DATA DMACH(1) / 1.8051943758648295760692620811737E-276 /
C      DATA DMACH(2) / 8.9884656743115795386465259539451E+307  /
C      DATA DMACH(3) / 6.1629758220391547297791294162718E-33  /
C      DATA DMACH(4) / 1.2325951644078309459558258832544E-32  /
C      DATA DMACH(5) / 0.30102999566398119521373889472449302  /
C
C     MACHINE CONSTANTS FOR AMDAHL MACHINES.
C
C === MACHINE = AMDAHL
C      DATA SMALL(1),SMALL(2) /    1048576,          0 /
C      DATA LARGE(1),LARGE(2) / 2147483647,         -1 /
C      DATA RIGHT(1),RIGHT(2) /  856686592,          0 /
C      DATA DIVER(1),DIVER(2) /  873463808,          0 /
C      DATA LOG10(1),LOG10(2) / 1091781651, 1352628735 /
C
C     MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM.
C
C === MACHINE = BURROUGHS.1700
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 === MACHINE = BURROUGHS.5700
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 === MACHINE = BURROUGHS.6700
C === MACHINE = BURROUGHS.7700
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 CONVEX C1, C2, C3 SERIES (NATIVE MODE)
C
C === MACHINE = CONVEX
C      DATA DMACH(1) / 5.562684646268007D-309 /
C      DATA DMACH(2) / 8.988465674311577D+307 /
C      DATA DMACH(3) / 1.110223024625157D-016 /
C      DATA DMACH(4) / 2.220446049250313D-016 /
C      DATA DMACH(5) / 3.010299956639812D-001 /
C
C     MACHINE CONSTANTS FOR THE CONVEX C1, C2, C3 SERIES (NATIVE MODE)
C     WITH -P8 OPTION
C
C === MACHINE = CONVEX.P8
C      DATA DMACH(1) / 8.500000000000000000000000000000000Q-4933 /
C      DATA DMACH(2) / 5.900000000000000000000000000000000Q+4931 /
C      DATA DMACH(3) / 1.925929944387235853055977942584927Q-0034 /
C      DATA DMACH(4) / 3.851859888774471706111955885169854Q-0034 /
C      DATA DMACH(5) / 3.010299956639811952137388947244930Q-0001 /
C
C     MACHINE CONSTANTS FOR THE CONVEX C1, C2, C3 SERIES (IEEE MODE)
C     WITH OR WITHOUT -P8 OPTION
C
C === MACHINE = CONVEX.IEEE
C === MACHINE = CONVEX.IEEE.P8
C      DATA DMACH(1) / 2.225073858507202D-308 /
C      DATA DMACH(2) / 1.797693134862315D+308 /
C      DATA DMACH(3) / 1.110223024625157D-016 /
C      DATA DMACH(4) / 2.220446049250313D-016 /
C      DATA DMACH(5) / 3.010299956639812D-001 /
C
C     MACHINE CONSTANTS FOR THE CYBER 170/180 SERIES USING NOS (FTN5).
C
C === MACHINE = CYBER.170.NOS
C === MACHINE = CYBER.180.NOS
C      DATA SMALL(1) / O"00604000000000000000" /
C      DATA SMALL(2) / O"00000000000000000000" /
C      DATA LARGE(1) / O"37767777777777777777" /
C      DATA LARGE(2) / O"37167777777777777777" /
C      DATA RIGHT(1) / O"15604000000000000000" /
C      DATA RIGHT(2) / O"15000000000000000000" /
C      DATA DIVER(1) / O"15614000000000000000" /
C      DATA DIVER(2) / O"15010000000000000000" /
C      DATA LOG10(1) / O"17164642023241175717" /
C      DATA LOG10(2) / O"16367571421742254654" /
C
C     MACHINE CONSTANTS FOR THE CDC 180 SERIES USING NOS/VE
C
C === MACHINE = CYBER.180.NOS/VE
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 CYBER 205
C
C === MACHINE = CYBER.205
C      DATA SMALL(1) / X'9000400000000000' /
C      DATA SMALL(2) / X'8FD1000000000000' /
C      DATA LARGE(1) / X'6FFF7FFFFFFFFFFF' /
C      DATA LARGE(2) / X'6FD07FFFFFFFFFFF' /
C      DATA RIGHT(1) / X'FF74400000000000' /
C      DATA RIGHT(2) / X'FF45000000000000' /
C      DATA DIVER(1) / X'FF75400000000000' /
C      DATA DIVER(2) / X'FF46000000000000' /
C      DATA LOG10(1) / X'FFD04D104D427DE7' /
C      DATA LOG10(2) / X'FFA17DE623E2566A' /
C
C     MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES.
C
C === MACHINE = CDC.6000
C === MACHINE = CDC.7000
C      DATA SMALL(1) / 00604000000000000000B /
C      DATA SMALL(2) / 00000000000000000000B /
C      DATA LARGE(1) / 37767777777777777777B /
C      DATA LARGE(2) / 37167777777777777777B /
C      DATA RIGHT(1) / 15604000000000000000B /
C      DATA RIGHT(2) / 15000000000000000000B /
C      DATA DIVER(1) / 15614000000000000000B /
C      DATA DIVER(2) / 15010000000000000000B /
C      DATA LOG10(1) / 17164642023241175717B /
C      DATA LOG10(2) / 16367571421742254654B /
C
C     MACHINE CONSTANTS FOR THE CRAY 1, XMP, 2, AND 3.
C
C === MACHINE = CRAY.46-BIT-INTEGER
C === MACHINE = CRAY.64-BIT-INTEGER
C      DATA SMALL(1) / 201354000000000000000B /
C      DATA SMALL(2) / 000000000000000000000B /
C      DATA LARGE(1) / 577767777777777777777B /
C      DATA LARGE(2) / 000007777777777777776B /
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 LINE -
C     STATIC DMACH(5)
C
C === MACHINE = DATA_GENERAL.ECLIPSE.S/200
C      DATA SMALL/20K,3*0/,LARGE/77777K,3*177777K/
C      DATA RIGHT/31420K,3*0/,DIVER/32020K,3*0/
C      DATA LOG10/40423K,42023K,50237K,74776K/
C
C     ELXSI 6400
C
C === MACHINE = ELSXI.6400
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), DIVER(2) / '3FD34413'X,'509F79FF'X /
C
C     MACHINE CONSTANTS FOR THE HARRIS 220
C     MACHINE CONSTANTS FOR THE HARRIS SLASH 6 AND SLASH 7
C
C === MACHINE = HARRIS.220
C === MACHINE = HARRIS.SLASH6
C === MACHINE = HARRIS.SLASH7
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     MACHINE CONSTANTS FOR THE HONEYWELL DPS 8/70 SERIES.
C
C === MACHINE = HONEYWELL.600/6000
C === MACHINE = HONEYWELL.DPS.8/70
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      3 WORD DOUBLE PRECISION OPTION WITH FTN4
C
C === MACHINE = HP.2100.3_WORD_DP
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      4 WORD DOUBLE PRECISION OPTION WITH FTN4
C
C === MACHINE = HP.2100.4_WORD_DP
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 IBM 360/370 SERIES,
C     THE XEROX SIGMA 5/7/9, THE SEL SYSTEMS 85/86, AND
C     THE INTERDATA 3230 AND INTERDATA 7/32.
C
C === MACHINE = IBM.360
C === MACHINE = IBM.370
C === MACHINE = XEROX.SIGMA.5
C === MACHINE = XEROX.SIGMA.7
C === MACHINE = XEROX.SIGMA.9
C === MACHINE = SEL.85
C === MACHINE = SEL.86
C === MACHINE = INTERDATA.3230
C === MACHINE = INTERDATA.7/32
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 INTERDATA 8/32
C     WITH THE UNIX SYSTEM FORTRAN 77 COMPILER.
C
C     FOR THE INTERDATA FORTRAN VII COMPILER REPLACE
C     THE Z'S SPECIFYING HEX CONSTANTS WITH Y'S.
C
C === MACHINE = INTERDATA.8/32.UNIX
C      DATA SMALL(1),SMALL(2) / Z'00100000', Z'00000000' /
C      DATA LARGE(1),LARGE(2) / Z'7EFFFFFF', Z'FFFFFFFF' /
C      DATA RIGHT(1),RIGHT(2) / Z'33100000', Z'00000000' /
C      DATA DIVER(1),DIVER(2) / Z'34100000', Z'00000000' /
C      DATA LOG10(1),LOG10(2) / Z'41134413', Z'509F79FF' /
C
C     MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR).
C
C === MACHINE = PDP-10.KA
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 === MACHINE = PDP-10.KI
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, "047674776746 /
C
C     MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING
C     32-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL).
C
C === MACHINE = PDP-11.32-BIT
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 === MACHINE = PDP-11.16-BIT
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 SEQUENT BALANCE 8000
C
C === MACHINE = SEQUENT.BALANCE.8000
C      DATA SMALL(1),SMALL(2) / $00000000,  $00100000 /
C      DATA LARGE(1),LARGE(2) / $FFFFFFFF,  $7FEFFFFF /
C      DATA RIGHT(1),RIGHT(2) / $00000000,  $3CA00000 /
C      DATA DIVER(1),DIVER(2) / $00000000,  $3CB00000 /
C      DATA LOG10(1),LOG10(2) / $509F79FF,  $3FD34413 /
C
C     MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. FTN COMPILER
C
C === MACHINE = UNIVAC.1100
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 INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS***
C
C === MACHINE = VAX.11/780
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    ***THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSYEMS***
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 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    ***THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSYEMS***
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
      D1MACH = DMACH(I)
      RETURN
C
      END
      SUBROUTINE DQC25S(F,A,B,BL,BR,ALFA,BETA,RI,RJ,RG,RH,RESULT,
     1   ABSERR,RESASC,INTEGR,NEV)
C***BEGIN PROLOGUE  DQC25S
C***DATE WRITTEN   810101   (YYMMDD)
C***REVISION DATE  830518   (YYMMDD)
C***REVISION HISTORY (YYMMDD)
C   000601   Changed DLOG/DABS to generic LOG/ABS
C***CATEGORY NO.  H2A2A2
C***KEYWORDS  25-POINT CLENSHAW-CURTIS INTEGRATION
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  To compute I = Integral of F*W over (BL,BR), with error
C            estimate, where the weight function W has a singular
C            behaviour of ALGEBRAICO-LOGARITHMIC type at the points
C            A and/or B. (BL,BR) is a part of (A,B).
C***DESCRIPTION
C
C        Integration rules for integrands having ALGEBRAICO-LOGARITHMIC
C        end point singularities
C        Standard fortran subroutine
C        Double precision version
C
C        PARAMETERS
C           F      - Double precision
C                    Function subprogram defining the integrand
C                    F(X). The actual name for F needs to be declared
C                    E X T E R N A L  in the driver program.
C
C           A      - Double precision
C                    Left end point of the original interval
C
C           B      - Double precision
C                    Right end point of the original interval, B.GT.A
C
C           BL     - Double precision
C                    Lower limit of integration, BL.GE.A
C
C           BR     - Double precision
C                    Upper limit of integration, BR.LE.B
C
C           ALFA   - Double precision
C                    PARAMETER IN THE WEIGHT FUNCTION
C
C           BETA   - Double precision
C                    Parameter in the weight function
C
C           RI,RJ,RG,RH - Double precision
C                    Modified CHEBYSHEV moments for the application
C                    of the generalized CLENSHAW-CURTIS
C                    method (computed in subroutine DQMOMO)
C
C           RESULT - Double precision
C                    Approximation to the integral
C                    RESULT is computed by using a generalized
C                    CLENSHAW-CURTIS method if B1 = A or BR = B.
C                    in all other cases the 15-POINT KRONROD
C                    RULE is applied, obtained by optimal addition of
C                    Abscissae to the 7-POINT GAUSS RULE.
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           RESASC - Double precision
C                    Approximation to the integral of ABS(F*W-I/(B-A))
C
C           INTEGR - Integer
C                    Which determines the weight function
C                    = 1   W(X) = (X-A)**ALFA*(B-X)**BETA
C                    = 2   W(X) = (X-A)**ALFA*(B-X)**BETA*LOG(X-A)
C                    = 3   W(X) = (X-A)**ALFA*(B-X)**BETA*LOG(B-X)
C                    = 4   W(X) = (X-A)**ALFA*(B-X)**BETA*LOG(X-A)*
C                                 LOG(B-X)
C
C           NEV    - Integer
C                    Number of integrand evaluations
C***REFERENCES  (NONE)
C***ROUTINES CALLED  DQCHEB,DQK15W,DQWGTS
C***END PROLOGUE  DQC25S
C
      DOUBLE PRECISION A,ABSERR,ALFA,B,BETA,BL,BR,CENTR,CHEB12,CHEB24,
     1  DC,F,FACTOR,FIX,FVAL,HLGTH,RESABS,RESASC,RESULT,RES12,
     2  RES24,RG,RH,RI,RJ,U,DQWGTS,X
      INTEGER I,INTEGR,ISYM,NEV
C
      DIMENSION CHEB12(13),CHEB24(25),FVAL(25),RG(25),RH(25),RI(25),
     1  RJ(25),X(11)
C
      EXTERNAL F,DQWGTS
C
C           THE VECTOR X CONTAINS THE VALUES COS(K*PI/24)
C           K = 1, ..., 11, TO BE USED FOR THE COMPUTATION OF THE
C           CHEBYSHEV SERIES EXPANSION OF F.
C
       DATA X(1),X(2),X(3),X(4),X(5),X(6),X(7),X(8),X(9),X(10),X(11)/
     1     0.9914448613738104D+00,     0.9659258262890683D+00,
     2     0.9238795325112868D+00,     0.8660254037844386D+00,
     3     0.7933533402912352D+00,     0.7071067811865475D+00,
     4     0.6087614290087206D+00,     0.5000000000000000D+00,
     5     0.3826834323650898D+00,     0.2588190451025208D+00,
     6     0.1305261922200516D+00/
C
C           LIST OF MAJOR VARIABLES
C           -----------------------
C
C           FVAL   - VALUE OF THE FUNCTION F AT THE POINTS
C                    (BR-BL)*0.5*COS(K*PI/24)+(BR+BL)*0.5
C                    K = 0, ..., 24
C           CHEB12 - COEFFICIENTS OF THE CHEBYSHEV SERIES EXPANSION
C                    OF DEGREE 12, FOR THE FUNCTION F, IN THE
C                    INTERVAL (BL,BR)
C           CHEB24 - COEFFICIENTS OF THE CHEBYSHEV SERIES EXPANSION
C                    OF DEGREE 24, FOR THE FUNCTION F, IN THE
C                    INTERVAL (BL,BR)
C           RES12  - APPROXIMATION TO THE INTEGRAL OBTAINED FROM CHEB12
C           RES24  - APPROXIMATION TO THE INTEGRAL OBTAINED FROM CHEB24
C           DQWGTS - EXTERNAL FUNCTION SUBPROGRAM DEFINING
C                    THE FOUR POSSIBLE WEIGHT FUNCTIONS
C           HLGTH  - HALF-LENGTH OF THE INTERVAL (BL,BR)
C           CENTR  - MID POINT OF THE INTERVAL (BL,BR)
C
C***FIRST EXECUTABLE STATEMENT  DQC25S
      NEV = 25
      IF(BL.EQ.A.AND.(ALFA.NE.0.0D+00.OR.INTEGR.EQ.2.OR.INTEGR.EQ.4))
     1 GO TO 10
      IF(BR.EQ.B.AND.(BETA.NE.0.0D+00.OR.INTEGR.EQ.3.OR.INTEGR.EQ.4))
     1 GO TO 140
C
C           IF A.GT.BL AND B.LT.BR, APPLY THE 15-POINT GAUSS-KRONROD
C           SCHEME.
C
C
      CALL DQK15W(F,DQWGTS,A,B,ALFA,BETA,INTEGR,BL,BR,
     1    RESULT,ABSERR,RESABS,RESASC)
      NEV = 15
      GO TO 270
C
C           THIS PART OF THE PROGRAM IS EXECUTED ONLY IF A = BL.
C           ----------------------------------------------------
C
C           COMPUTE THE CHEBYSHEV SERIES EXPANSION OF THE
C           FOLLOWING FUNCTION
C           F1 = (0.5*(B+B-BR-A)-0.5*(BR-A)*X)**BETA
C                  *F(0.5*(BR-A)*X+0.5*(BR+A))
C
   10 HLGTH = 0.5D+00*(BR-BL)
      CENTR = 0.5D+00*(BR+BL)
      FIX = B-CENTR
      FVAL(1) = 0.5D+00*F(HLGTH+CENTR)*(FIX-HLGTH)**BETA
      FVAL(13) = F(CENTR)*(FIX**BETA)
      FVAL(25) = 0.5D+00*F(CENTR-HLGTH)*(FIX+HLGTH)**BETA
      DO 20 I=2,12
        U = HLGTH*X(I-1)
        ISYM = 26-I
        FVAL(I) = F(U+CENTR)*(FIX-U)**BETA
        FVAL(ISYM) = F(CENTR-U)*(FIX+U)**BETA
   20 CONTINUE
      FACTOR = HLGTH**(ALFA+0.1D+01)
      RESULT = 0.0D+00
      ABSERR = 0.0D+00
      RES12 = 0.0D+00
      RES24 = 0.0D+00
      IF(INTEGR.GT.2) GO TO 70
      CALL DQCHEB(X,FVAL,CHEB12,CHEB24)
C
C           INTEGR = 1  (OR 2)
C
      DO 30 I=1,13
        RES12 = RES12+CHEB12(I)*RI(I)
        RES24 = RES24+CHEB24(I)*RI(I)
   30 CONTINUE
      DO 40 I=14,25
        RES24 = RES24+CHEB24(I)*RI(I)
   40 CONTINUE
      IF(INTEGR.EQ.1) GO TO 130
C
C           INTEGR = 2
C
      DC = LOG(BR-BL)
      RESULT = RES24*DC
      ABSERR = ABS((RES24-RES12)*DC)
      RES12 = 0.0D+00
      RES24 = 0.0D+00
      DO 50 I=1,13
        RES12 = RES12+CHEB12(I)*RG(I)
        RES24 = RES12+CHEB24(I)*RG(I)
   50 CONTINUE
      DO 60 I=14,25
        RES24 = RES24+CHEB24(I)*RG(I)
   60 CONTINUE
      GO TO 130
C
C           COMPUTE THE CHEBYSHEV SERIES EXPANSION OF THE
C           FOLLOWING FUNCTION
C           F4 = F1*LOG(0.5*(B+B-BR-A)-0.5*(BR-A)*X)
C
   70 FVAL(1) = FVAL(1)*LOG(FIX-HLGTH)
      FVAL(13) = FVAL(13)*LOG(FIX)
      FVAL(25) = FVAL(25)*LOG(FIX+HLGTH)
      DO 80 I=2,12
        U = HLGTH*X(I-1)
        ISYM = 26-I
        FVAL(I) = FVAL(I)*LOG(FIX-U)
        FVAL(ISYM) = FVAL(ISYM)*LOG(FIX+U)
   80 CONTINUE
      CALL DQCHEB(X,FVAL,CHEB12,CHEB24)
C
C           INTEGR = 3  (OR 4)
C
      DO 90 I=1,13
        RES12 = RES12+CHEB12(I)*RI(I)
        RES24 = RES24+CHEB24(I)*RI(I)
   90 CONTINUE
      DO 100 I=14,25
        RES24 = RES24+CHEB24(I)*RI(I)
  100 CONTINUE
      IF(INTEGR.EQ.3) GO TO 130
C
C           INTEGR = 4
C
      DC = LOG(BR-BL)
      RESULT = RES24*DC
      ABSERR = ABS((RES24-RES12)*DC)
      RES12 = 0.0D+00
      RES24 = 0.0D+00
      DO 110 I=1,13
        RES12 = RES12+CHEB12(I)*RG(I)
        RES24 = RES24+CHEB24(I)*RG(I)
  110 CONTINUE
      DO 120 I=14,25
        RES24 = RES24+CHEB24(I)*RG(I)
  120 CONTINUE
  130 RESULT = (RESULT+RES24)*FACTOR
      ABSERR = (ABSERR+ABS(RES24-RES12))*FACTOR
      GO TO 270
C
C           THIS PART OF THE PROGRAM IS EXECUTED ONLY IF B = BR.
C           ----------------------------------------------------
C
C           COMPUTE THE CHEBYSHEV SERIES EXPANSION OF THE
C           FOLLOWING FUNCTION
C           F2 = (0.5*(B+BL-A-A)+0.5*(B-BL)*X)**ALFA
C                *F(0.5*(B-BL)*X+0.5*(B+BL))
C
  140 HLGTH = 0.5D+00*(BR-BL)
      CENTR = 0.5D+00*(BR+BL)
      FIX = CENTR-A
      FVAL(1) = 0.5D+00*F(HLGTH+CENTR)*(FIX+HLGTH)**ALFA
      FVAL(13) = F(CENTR)*(FIX**ALFA)
      FVAL(25) = 0.5D+00*F(CENTR-HLGTH)*(FIX-HLGTH)**ALFA
      DO 150 I=2,12
        U = HLGTH*X(I-1)
        ISYM = 26-I
        FVAL(I) = F(U+CENTR)*(FIX+U)**ALFA
        FVAL(ISYM) = F(CENTR-U)*(FIX-U)**ALFA
  150 CONTINUE
      FACTOR = HLGTH**(BETA+0.1D+01)
      RESULT = 0.0D+00
      ABSERR = 0.0D+00
      RES12 = 0.0D+00
      RES24 = 0.0D+00
      IF(INTEGR.EQ.2.OR.INTEGR.EQ.4) GO TO 200
C
C           INTEGR = 1  (OR 3)
C
      CALL DQCHEB(X,FVAL,CHEB12,CHEB24)
      DO 160 I=1,13
        RES12 = RES12+CHEB12(I)*RJ(I)
        RES24 = RES24+CHEB24(I)*RJ(I)
  160 CONTINUE
      DO 170 I=14,25
        RES24 = RES24+CHEB24(I)*RJ(I)
  170 CONTINUE
      IF(INTEGR.EQ.1) GO TO 260
C
C           INTEGR = 3
C
      DC = LOG(BR-BL)
      RESULT = RES24*DC
      ABSERR = ABS((RES24-RES12)*DC)
      RES12 = 0.0D+00
      RES24 = 0.0D+00
      DO 180 I=1,13
        RES12 = RES12+CHEB12(I)*RH(I)
        RES24 = RES24+CHEB24(I)*RH(I)
  180 CONTINUE
      DO 190 I=14,25
        RES24 = RES24+CHEB24(I)*RH(I)
  190 CONTINUE
      GO TO 260
C
C           COMPUTE THE CHEBYSHEV SERIES EXPANSION OF THE
C           FOLLOWING FUNCTION
C           F3 = F2*LOG(0.5*(B-BL)*X+0.5*(B+BL-A-A))
C
  200 FVAL(1) = FVAL(1)*LOG(HLGTH+FIX)
      FVAL(13) = FVAL(13)*LOG(FIX)
      FVAL(25) = FVAL(25)*LOG(FIX-HLGTH)
      DO 210 I=2,12
        U = HLGTH*X(I-1)
        ISYM = 26-I
        FVAL(I) = FVAL(I)*LOG(U+FIX)
        FVAL(ISYM) = FVAL(ISYM)*LOG(FIX-U)
  210 CONTINUE
      CALL DQCHEB(X,FVAL,CHEB12,CHEB24)
C
C           INTEGR = 2  (OR 4)
C
      DO 220 I=1,13
        RES12 = RES12+CHEB12(I)*RJ(I)
        RES24 = RES24+CHEB24(I)*RJ(I)
  220 CONTINUE
      DO 230 I=14,25
        RES24 = RES24+CHEB24(I)*RJ(I)
  230 CONTINUE
      IF(INTEGR.EQ.2) GO TO 260
      DC = LOG(BR-BL)
      RESULT = RES24*DC
      ABSERR = ABS((RES24-RES12)*DC)
      RES12 = 0.0D+00
      RES24 = 0.0D+00
C
C           INTEGR = 4
C
      DO 240 I=1,13
        RES12 = RES12+CHEB12(I)*RH(I)
        RES24 = RES24+CHEB24(I)*RH(I)
  240 CONTINUE
      DO 250 I=14,25
        RES24 = RES24+CHEB24(I)*RH(I)
  250 CONTINUE
  260 RESULT = (RESULT+RES24)*FACTOR
      ABSERR = (ABSERR+ABS(RES24-RES12))*FACTOR
  270 RETURN
      END
      SUBROUTINE DQCHEB(X,FVAL,CHEB12,CHEB24)
C***BEGIN PROLOGUE  DQCHEB
C***REFER TO  DQC25C,DQC25F,DQC25S
C***ROUTINES CALLED  (NONE)
C***REVISION DATE  830518   (YYMMDD)
C***KEYWORDS  CHEBYSHEV SERIES EXPANSION,FAST FOURIER TRANSFORM
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 computes the CHEBYSHEV series expansion
C            of degrees 12 and 24 of a function using A
C            FAST FOURIER TRANSFORM METHOD
C            F(X) = SUM(K=1,..,13) (CHEB12(K)*T(K-1,X)),
C            F(X) = SUM(K=1,..,25) (CHEB24(K)*T(K-1,X)),
C            Where T(K,X) is the CHEBYSHEV POLYNOMIAL OF DEGREE K.
C***DESCRIPTION
C
C        Chebyshev Series Expansion
C        Standard Fortran Subroutine
C        Double precision version
C
C        PARAMETERS
C          ON ENTRY
C           X      - Double precision
C                    Vector of dimension 11 containing the
C                    Values COS(K*PI/24), K = 1, ..., 11
C
C           FVAL   - Double precision
C                    Vector of dimension 25 containing the
C                    function values at the points
C                    (B+A+(B-A)*COS(K*PI/24))/2, K = 0, ...,24,
C                    where (A,B) is the approximation interval.
C                    FVAL(1) and FVAL(25) are divided by two
C                    (these values are destroyed at output).
C
C          ON RETURN
C           CHEB12 - Double precision
C                    Vector of dimension 13 containing the
C                    CHEBYSHEV coefficients for degree 12
C
C           CHEB24 - Double precision
C                    Vector of dimension 25 containing the
C                    CHEBYSHEV Coefficients for degree 24
C***END PROLOGUE  DQCHEB
C
      DOUBLE PRECISION ALAM,ALAM1,ALAM2,CHEB12,CHEB24,FVAL,PART1,PART2,
     1  PART3,V,X
      INTEGER I,J
C
      DIMENSION CHEB12(13),CHEB24(25),FVAL(25),V(12),X(11)
C
C***FIRST EXECUTABLE STATEMENT  DQCHEB
      DO 10 I=1,12
        J = 26-I
        V(I) = FVAL(I)-FVAL(J)
        FVAL(I) = FVAL(I)+FVAL(J)
   10 CONTINUE
      ALAM1 = V(1)-V(9)
      ALAM2 = X(6)*(V(3)-V(7)-V(11))
      CHEB12(4) = ALAM1+ALAM2
      CHEB12(10) = ALAM1-ALAM2
      ALAM1 = V(2)-V(8)-V(10)
      ALAM2 = V(4)-V(6)-V(12)
      ALAM = X(3)*ALAM1+X(9)*ALAM2
      CHEB24(4) = CHEB12(4)+ALAM
      CHEB24(22) = CHEB12(4)-ALAM
      ALAM = X(9)*ALAM1-X(3)*ALAM2
      CHEB24(10) = CHEB12(10)+ALAM
      CHEB24(16) = CHEB12(10)-ALAM
      PART1 = X(4)*V(5)
      PART2 = X(8)*V(9)
      PART3 = X(6)*V(7)
      ALAM1 = V(1)+PART1+PART2
      ALAM2 = X(2)*V(3)+PART3+X(10)*V(11)
      CHEB12(2) = ALAM1+ALAM2
      CHEB12(12) = ALAM1-ALAM2
      ALAM = X(1)*V(2)+X(3)*V(4)+X(5)*V(6)+X(7)*V(8)
     1  +X(9)*V(10)+X(11)*V(12)
      CHEB24(2) = CHEB12(2)+ALAM
      CHEB24(24) = CHEB12(2)-ALAM
      ALAM = X(11)*V(2)-X(9)*V(4)+X(7)*V(6)-X(5)*V(8)
     1  +X(3)*V(10)-X(1)*V(12)
      CHEB24(12) = CHEB12(12)+ALAM
      CHEB24(14) = CHEB12(12)-ALAM
      ALAM1 = V(1)-PART1+PART2
      ALAM2 = X(10)*V(3)-PART3+X(2)*V(11)
      CHEB12(6) = ALAM1+ALAM2
      CHEB12(8) = ALAM1-ALAM2
      ALAM = X(5)*V(2)-X(9)*V(4)-X(1)*V(6)
     1  -X(11)*V(8)+X(3)*V(10)+X(7)*V(12)
      CHEB24(6) = CHEB12(6)+ALAM
      CHEB24(20) = CHEB12(6)-ALAM
      ALAM = X(7)*V(2)-X(3)*V(4)-X(11)*V(6)+X(1)*V(8)
     1  -X(9)*V(10)-X(5)*V(12)
      CHEB24(8) = CHEB12(8)+ALAM
      CHEB24(18) = CHEB12(8)-ALAM
      DO 20 I=1,6
        J = 14-I
        V(I) = FVAL(I)-FVAL(J)
        FVAL(I) = FVAL(I)+FVAL(J)
   20 CONTINUE
      ALAM1 = V(1)+X(8)*V(5)
      ALAM2 = X(4)*V(3)
      CHEB12(3) = ALAM1+ALAM2
      CHEB12(11) = ALAM1-ALAM2
      CHEB12(7) = V(1)-V(5)
      ALAM = X(2)*V(2)+X(6)*V(4)+X(10)*V(6)
      CHEB24(3) = CHEB12(3)+ALAM
      CHEB24(23) = CHEB12(3)-ALAM
      ALAM = X(6)*(V(2)-V(4)-V(6))
      CHEB24(7) = CHEB12(7)+ALAM
      CHEB24(19) = CHEB12(7)-ALAM
      ALAM = X(10)*V(2)-X(6)*V(4)+X(2)*V(6)
      CHEB24(11) = CHEB12(11)+ALAM
      CHEB24(15) = CHEB12(11)-ALAM
      DO 30 I=1,3
        J = 8-I
        V(I) = FVAL(I)-FVAL(J)
        FVAL(I) = FVAL(I)+FVAL(J)
   30 CONTINUE
      CHEB12(5) = V(1)+X(8)*V(3)
      CHEB12(9) = FVAL(1)-X(8)*FVAL(3)
      ALAM = X(4)*V(2)
      CHEB24(5) = CHEB12(5)+ALAM
      CHEB24(21) = CHEB12(5)-ALAM
      ALAM = X(8)*FVAL(2)-FVAL(4)
      CHEB24(9) = CHEB12(9)+ALAM
      CHEB24(17) = CHEB12(9)-ALAM
      CHEB12(1) = FVAL(1)+FVAL(3)
      ALAM = FVAL(2)+FVAL(4)
      CHEB24(1) = CHEB12(1)+ALAM
      CHEB24(25) = CHEB12(1)-ALAM
      CHEB12(13) = V(1)-V(3)
      CHEB24(13) = CHEB12(13)
      ALAM = 0.1D+01/0.6D+01
      DO 40 I=2,12
        CHEB12(I) = CHEB12(I)*ALAM
   40 CONTINUE
      ALAM = 0.5D+00*ALAM
      CHEB12(1) = CHEB12(1)*ALAM
      CHEB12(13) = CHEB12(13)*ALAM
      DO 50 I=2,24
        CHEB24(I) = CHEB24(I)*ALAM
   50 CONTINUE
      CHEB24(1) = 0.5D+00*ALAM*CHEB24(1)
      CHEB24(25) = 0.5D+00*ALAM*CHEB24(25)
      RETURN
      END
      SUBROUTINE DQMOMO(ALFA,BETA,RI,RJ,RG,RH,INTEGR)
C***BEGIN PROLOGUE  DQMOMO
C***DATE WRITTEN   820101   (YYMMDD)
C***REVISION DATE  830518   (YYMMDD)
C***CATEGORY NO.  H2A2A1,C3A2
C***KEYWORDS  MODIFIED CHEBYSHEV MOMENTS
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 computes modified CHEBSYSHEV moments. The K-th
C            modified CHEBYSHEV moment is defined as the integral over
C            (-1,1) of W(X)*T(K,X), where T(K,X) is the CHEBYSHEV
C            POLYNOMIAL of degree K.
C***DESCRIPTION
C
C        MODIFIED CHEBYSHEV MOMENTS
C        STANDARD FORTRAN SUBROUTINE
C        DOUBLE PRECISION VERSION
C
C        PARAMETERS
C           ALFA   - Double precision
C                    Parameter in the weight function W(X), ALFA.GT.(-1)
C
C           BETA   - Double precision
C                    Parameter in the weight functioN W(X), BETA.GT.(-1)
C
C           RI     - Double precision
C                    Vector of dimension 25
C                    RI(K) is the integral over (-1,1) of
C                    (1+X)**ALFA*T(K-1,X), K = 1, ..., 25.
C
C           RJ     - Double precision
C                    Vector of dimension 25
C                    RJ(K) is the integral over (-1,1) of
C                    (1-X)**BETA*T(K-1,X), K = 1, ..., 25.
C
C           RG     - Double precision
C                    Vector of dimension 25
C                    RG(K) is the integral over (-1,1) of
C                    (1+X)**ALFA*LOG((1+X)/2)*T(K-1,X), K = 1, ..., 25.
C
C           RH     - Double precision
C                    Vector of dimension 25
C                    RH(K) is the integral over (-1,1) of
C                    (1-X)**BETA*LOG((1-X)/2)*T(K-1,X), K = 1, ..., 25.
C
C           INTEGR - Integer
C                    Input parameter indicating the modified
C                    Moments to be computed
C                    Integr = 1 compute RI, RJ
C                           = 2 compute RI, RJ, RG
C                           = 3 compute RI, RJ, RH
C                           = 4 compute RI, RJ, RG, RH
C***REFERENCES  (NONE)
C***ROUTINES CALLED  (NONE)
C***END PROLOGUE  DQMOMO
C
      DOUBLE PRECISION ALFA,ALFP1,ALFP2,AN,ANM1,BETA,BETP1,BETP2,RALF,
     1  RBET,RG,RH,RI,RJ
      INTEGER I,IM1,INTEGR
C
      DIMENSION RG(25),RH(25),RI(25),RJ(25)
C
C
C***FIRST EXECUTABLE STATEMENT  DQMOMO
      ALFP1 = ALFA+0.1D+01
      BETP1 = BETA+0.1D+01
      ALFP2 = ALFA+0.2D+01
      BETP2 = BETA+0.2D+01
      RALF = 0.2D+01**ALFP1
      RBET = 0.2D+01**BETP1
C
C           COMPUTE RI, RJ USING A FORWARD RECURRENCE RELATION.
C
      RI(1) = RALF/ALFP1
      RJ(1) = RBET/BETP1
      RI(2) = RI(1)*ALFA/ALFP2
      RJ(2) = RJ(1)*BETA/BETP2
      AN = 0.2D+01
      ANM1 = 0.1D+01
      DO 20 I=3,25
        RI(I) = -(RALF+AN*(AN-ALFP2)*RI(I-1))/(ANM1*(AN+ALFP1))
        RJ(I) = -(RBET+AN*(AN-BETP2)*RJ(I-1))/(ANM1*(AN+BETP1))
        ANM1 = AN
        AN = AN+0.1D+01
   20 CONTINUE
      IF(INTEGR.EQ.1) GO TO 70
      IF(INTEGR.EQ.3) GO TO 40
C
C           COMPUTE RG USING A FORWARD RECURRENCE RELATION.
C
      RG(1) = -RI(1)/ALFP1
      RG(2) = -(RALF+RALF)/(ALFP2*ALFP2)-RG(1)
      AN = 0.2D+01
      ANM1 = 0.1D+01
      IM1 = 2
      DO 30 I=3,25
        RG(I) = -(AN*(AN-ALFP2)*RG(IM1)-AN*RI(IM1)+ANM1*RI(I))/
     1  (ANM1*(AN+ALFP1))
        ANM1 = AN
        AN = AN+0.1D+01
        IM1 = I
   30 CONTINUE
      IF(INTEGR.EQ.2) GO TO 70
C
C           COMPUTE RH USING A FORWARD RECURRENCE RELATION.
C
   40 RH(1) = -RJ(1)/BETP1
      RH(2) = -(RBET+RBET)/(BETP2*BETP2)-RH(1)
      AN = 0.2D+01
      ANM1 = 0.1D+01
      IM1 = 2
      DO 50 I=3,25
        RH(I) = -(AN*(AN-BETP2)*RH(IM1)-AN*RJ(IM1)+
     1  ANM1*RJ(I))/(ANM1*(AN+BETP1))
        ANM1 = AN
        AN = AN+0.1D+01
        IM1 = I
   50 CONTINUE
      DO 60 I=2,25,2
        RH(I) = -RH(I)
   60 CONTINUE
   70 DO 80 I=2,25,2
        RJ(I) = -RJ(I)
   80 CONTINUE
   90 RETURN
      END
      DOUBLE PRECISION FUNCTION DQWGTS(X,A,B,ALFA,BETA,INTEGR)
C***BEGIN PROLOGUE  DQWGTS
C***REFER TO  DQK15W
C***ROUTINES CALLED  (NONE)
C***REVISION DATE  830518   (YYMMDD)
C***REVISION HISTORY (YYMMDD)
C   000601   Changed DLOG to generic LOG
C***KEYWORDS  ALGEBRAICO-LOGARITHMIC,END-POINT SINGULARITIES,
C             WEIGHT FUNCTION
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 function subprogram is used together with the
C            routine DQAWS and defines the WEIGHT function.
C***END PROLOGUE  DQWGTS
C
      DOUBLE PRECISION A,ALFA,B,BETA,BMX,X,XMA
      INTEGER INTEGR
C***FIRST EXECUTABLE STATEMENT  DQWGTS
      XMA = X-A
      BMX = B-X
      DQWGTS = XMA**ALFA*BMX**BETA
      GO TO (40,10,20,30),INTEGR
   10 DQWGTS = DQWGTS*LOG(XMA)
      GO TO 40
   20 DQWGTS = DQWGTS*LOG(BMX)
      GO TO 40
   30 DQWGTS = DQWGTS*LOG(XMA)*LOG(BMX)
   40 RETURN
      END
      SUBROUTINE DQK15W(F,W,P1,P2,P3,P4,KP,A,B,RESULT,ABSERR,RESABS,
     1   RESASC)
C***BEGIN PROLOGUE  DQK15W
C***DATE WRITTEN   810101   (YYMMDD)
C***REVISION DATE  830518   (YYMMDD)
C***REVISION HISTORY (YYMMDD)
C   000601   Changed DMAX1/DMIN1/DABS to generic MAX/MIN/ABS
C***CATEGORY NO.  H2A2A2
C***KEYWORDS  15-POINT GAUSS-KRONROD RULES
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  To compute I = Integral of F*W over (A,B), with error
C                           estimate
C                       J = Integral of ABS(F*W) over (A,B)
C***DESCRIPTION
C
C           Integration rules
C           Standard fortran subroutine
C           Double precision version
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              W      - Double precision
C                       Function subprogram defining the integrand
C                       WEIGHT function W(X). The actual name for W
C                       needs to be declared E X T E R N A L in the
C                       calling program.
C
C              P1, P2, P3, P4 - Double precision
C                       Parameters in the WEIGHT function
C
C              KP     - Integer
C                       Key for indicating the type of WEIGHT function
C
C              A      - Double precision
C                       Lower limit of integration
C
C              B      - Double precision
C                       Upper limit of integration
C
C            ON RETURN
C              RESULT - Double precision
C                       Approximation to the integral I
C                       RESULT is computed by applying the 15-point
C                       Kronrod rule (RESK) obtained by optimal addition
C                       of abscissae to the 7-point Gauss rule (RESG).
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              RESABS - Double precision
C                       Approximation to the integral of ABS(F)
C
C              RESASC - Double precision
C                       Approximation to the integral of ABS(F-I/(B-A))
C***REFERENCES  (NONE)
C***ROUTINES CALLED  D1MACH
C***END PROLOGUE  DQK15W
C
      DOUBLE PRECISION A,ABSC,ABSC1,ABSC2,ABSERR,B,CENTR,DHLGTH,
     1  D1MACH,EPMACH,F,FC,FSUM,FVAL1,FVAL2,FV1,FV2,HLGTH,
     2  P1,P2,P3,P4,RESABS,RESASC,RESG,RESK,RESKH,RESULT,UFLOW,W,WG,WGK,
     3  XGK
      INTEGER J,JTW,JTWM1,KP
      EXTERNAL F,W
C
      DIMENSION FV1(7),FV2(7),XGK(8),WGK(8),WG(4)
C
C           THE ABSCISSAE AND WEIGHTS ARE GIVEN FOR THE INTERVAL (-1,1).
C           BECAUSE OF SYMMETRY ONLY THE POSITIVE ABSCISSAE AND THEIR
C           CORRESPONDING WEIGHTS ARE GIVEN.
C
C           XGK    - ABSCISSAE OF THE 15-POINT GAUSS-KRONROD RULE
C                    XGK(2), XGK(4), ... ABSCISSAE OF THE 7-POINT
C                    GAUSS RULE
C                    XGK(1), XGK(3), ... ABSCISSAE WHICH ARE OPTIMALLY
C                    ADDED TO THE 7-POINT GAUSS RULE
C
C           WGK    - WEIGHTS OF THE 15-POINT GAUSS-KRONROD RULE
C
C           WG     - WEIGHTS OF THE 7-POINT GAUSS RULE
C
      DATA XGK(1),XGK(2),XGK(3),XGK(4),XGK(5),XGK(6),XGK(7),XGK(8)/
     1     0.9914553711208126D+00,     0.9491079123427585D+00,
     2     0.8648644233597691D+00,     0.7415311855993944D+00,
     3     0.5860872354676911D+00,     0.4058451513773972D+00,
     4     0.2077849550078985D+00,     0.0000000000000000D+00/
C
      DATA WGK(1),WGK(2),WGK(3),WGK(4),WGK(5),WGK(6),WGK(7),WGK(8)/
     1     0.2293532201052922D-01,     0.6309209262997855D-01,
     2     0.1047900103222502D+00,     0.1406532597155259D+00,
     3     0.1690047266392679D+00,     0.1903505780647854D+00,
     4     0.2044329400752989D+00,     0.2094821410847278D+00/
C
      DATA WG(1),WG(2),WG(3),WG(4)/
     1     0.1294849661688697D+00,    0.2797053914892767D+00,
     2     0.3818300505051889D+00,    0.4179591836734694D+00/
C
C
C           LIST OF MAJOR VARIABLES
C           -----------------------
C
C           CENTR  - MID POINT OF THE INTERVAL
C           HLGTH  - HALF-LENGTH OF THE INTERVAL
C           ABSC*  - ABSCISSA
C           FVAL*  - FUNCTION VALUE
C           RESG   - RESULT OF THE 7-POINT GAUSS FORMULA
C           RESK   - RESULT OF THE 15-POINT KRONROD FORMULA
C           RESKH  - APPROXIMATION TO THE MEAN VALUE OF F*W OVER (A,B),
C                    I.E. TO I/(B-A)
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  DQK15W
      EPMACH = D1MACH(4)
      UFLOW = D1MACH(1)
C
      CENTR = 0.5D+00*(A+B)
      HLGTH = 0.5D+00*(B-A)
      DHLGTH = ABS(HLGTH)
C
C           COMPUTE THE 15-POINT KRONROD APPROXIMATION TO THE
C           INTEGRAL, AND ESTIMATE THE ERROR.
C
      FC = F(CENTR)*W(CENTR,P1,P2,P3,P4,KP)
      RESG = WG(4)*FC
      RESK = WGK(8)*FC
      RESABS = ABS(RESK)
      DO 10 J=1,3
        JTW = J*2
        ABSC = HLGTH*XGK(JTW)
        ABSC1 = CENTR-ABSC
        ABSC2 = CENTR+ABSC
        FVAL1 = F(ABSC1)*W(ABSC1,P1,P2,P3,P4,KP)
        FVAL2 = F(ABSC2)*W(ABSC2,P1,P2,P3,P4,KP)
        FV1(JTW) = FVAL1
        FV2(JTW) = FVAL2
        FSUM = FVAL1+FVAL2
        RESG = RESG+WG(J)*FSUM
        RESK = RESK+WGK(JTW)*FSUM
        RESABS = RESABS+WGK(JTW)*(ABS(FVAL1)+ABS(FVAL2))
   10 CONTINUE
      DO 15 J=1,4
        JTWM1 = J*2-1
        ABSC = HLGTH*XGK(JTWM1)
        ABSC1 = CENTR-ABSC
        ABSC2 = CENTR+ABSC
        FVAL1 = F(ABSC1)*W(ABSC1,P1,P2,P3,P4,KP)
        FVAL2 = F(ABSC2)*W(ABSC2,P1,P2,P3,P4,KP)
        FV1(JTWM1) = FVAL1
        FV2(JTWM1) = FVAL2
        FSUM = FVAL1+FVAL2
        RESK = RESK+WGK(JTWM1)*FSUM
        RESABS = RESABS+WGK(JTWM1)*(ABS(FVAL1)+ABS(FVAL2))
   15 CONTINUE
      RESKH = RESK*0.5D+00
      RESASC = WGK(8)*ABS(FC-RESKH)
      DO 20 J=1,7
        RESASC = RESASC+WGK(J)*(ABS(FV1(J)-RESKH)+ABS(FV2(J)-RESKH))
   20 CONTINUE
      RESULT = RESK*HLGTH
      RESABS = RESABS*DHLGTH
      RESASC = RESASC*DHLGTH
      ABSERR = ABS((RESK-RESG)*HLGTH)
      IF(RESASC.NE.0.0D+00.AND.ABSERR.NE.0.0D+00)
     1  ABSERR = RESASC*MIN(0.1D+01,(0.2D+03*ABSERR/RESASC)**1.5D+00)
      IF(RESABS.GT.UFLOW/(0.5D+02*EPMACH)) ABSERR = MAX((EPMACH*
     1  0.5D+02)*RESABS,ABSERR)
      RETURN
      END
