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