c     DRDINTM
c>> 2003-03-30 DRDINTMF Krogh Changed way of dealing with 0 denominator.
c>> 2001-05-22 DRDINTMF Krogh Minor change for making .f90 version.
c>> 1994-11-02 DRDINTMF Krogh  Changes to use M77CON
c>> 1994-08-08 DRDINTMF Snyder  Took '0' out of formats for C conversion
c>> 1993-05-05 DRDINTMF Krogh   Adjusted to simplify conversion to C.
c>> 1987-12-09 DRDINTMF Snyder  Initial Code.
c--D replaces "?": DR?INTM, DR?INTMF, ?INTM, ?INTF
c
c     DEMO DRIVER for multi-dimensional quadrature subprogram DINTM.
c
c     Compute the integral for Y = 0.0 to Y = PI of
c             the integral for X = 0.0 to X = Y of
c             X * COS(Y) /(X*X+Y*Y).  The ANSWER should be zero.
c
c     The integrand and all limits are provided by the subprogram DINTF,
c     although the limits of the outer dimension, and the lower limit of
c     the inner dimension, being constants, could have been provided
c     here.
c
      integer NDIMI, NWPD, NWORK, IOPT(10)
      parameter (NDIMI=2, NWPD=217, NWORK=3*NDIMI+NWPD*(NDIMI-1))
      double precision ANSWER, WORK(NWORK)
10    format (/
     1' DRDINTMF:'/
     2' Compute the integral for Y = 0.0 to Y = PI of'/
     3'         the integral for X = 0.0 to X = Y of'/
     4'         X * COS(Y) /(X*X+Y*Y).  The ANSWER should be zero.')
20    format (/' ANSWER =        ',G15.8/
     1         ' ERROR ESTIMATE =',G15.8/
     2         ' STATUS FLAG =   ',I3/
     3         ' FUNCTION VALUES (INNERMOST INTEGRALS) =',I6)
c
      print 10
      IOPT(2) = 10
      IOPT(3) = 0
      IOPT(4) = 0
      call DINTM (NDIMI, ANSWER, WORK, NWORK, IOPT)
      print 20, ANSWER, WORK(1), IOPT(1), IOPT(3)
      stop
      end
c   End of DRDINTMF

      subroutine DINTF (ANSWER, WORK, IFLAG)
c
c     Subroutine to provide integrand and limits for DINTM.
c
c     This sample subroutine demonstrates the use of functional
c     transformation of the inner integral to reduce the cost of the
c     overall computation.  The factor "cos y" does not depend on the
c     variable of integration for the inner integral, and therefore is
c     factored out.  Similarly, the term y*y in the denominator of the
c     inner integrand is pre-computed when the limits of the inner
c     integral are requested.
c
c     This example also demonstrates the necessity to cope with overflow
c     of the integrand, even though the integral is well behaved.
c
      double precision ANSWER, WORK(*)
      integer IFLAG
      integer NDIMI
      double precision COSY, DENOM, XDY, YSQ
      save COSY, YSQ
      data NDIMI /2/
c
      if (IFLAG .EQ. 0) then
c
c       IFLAG = 0, compute innermost integrand.
c
        DENOM = WORK(1) * WORK(1) + YSQ
        if (DENOM .NE. 0.0D0) then
c        This test will not detect overflow of the integrand if the
c        arithmetic underflows gradually.
           ANSWER = WORK(1) / DENOM
        else
c     Special care to avoid a zero denominator.
           XDY = WORK(1) / WORK(2)
           ANSWER = XDY / (WORK(2) * (1.D0 + XDY**2))
        end if
c
      else if (IFLAG .EQ. 1) then
c
c       Compute limits of inner dimension.
c
c   Set WORK(1) = COS(WORK(2)) = Partial derivative, with respect to
c   the integral over the inner dimension, of the transformation
c   applied when integration over the inner dimension is complete.
c   This is so dintm knows how much accuracy it needs for this integral.
c
        WORK(NDIMI+IFLAG) = 0.0D0
        WORK(2*NDIMI+IFLAG) = WORK(2)
        YSQ = WORK(2) * WORK(2)
        COSY = COS(WORK(2))
        WORK(1) = COSY
c
      else if (IFLAG .EQ. 2) then
c
c       Compute limits of outer dimension.
c
        WORK(NDIMI+IFLAG) = 0.0D0
        WORK(2*NDIMI+IFLAG) = 4.0D0*ATAN(1.0D0)
c
      else
c
c       IFLAG < 0, transform inner integrand.
c
        ANSWER = ANSWER * COSY
        WORK(1) = WORK(1) * COSY
c
      end if
c
      return
c
      end