c program DRDINTMR
c>> 2003-03-30 DRDINTMR Krogh Changed way of dealing with 0 denominator.
c>> 2001-05-22 DRDINTMR Krogh Minor change for making .f90 version.
c>> 1994-10-19 DRDINTMR Krogh Changes to use M77CON
c>> 1994-08-31 DRDINTMR Snyder Moved formats for C conversion
c>> 1994-08-08 DRDINTMR Snyder Took '0' out of formats for C conversion
c>> 1991-11-20 CLL Edited for Fortran 90.
c>> 1987-12-09 DRDINTMR Snyder Initial Code.
c--D replaces "?": DR?INTMR, ?INTM, ?INTMA, ?INTA
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 reverse
c communication.
c
c This sample demonstrates the use of functional transformation of
c the inner integral to reduce the cost of the overall computation.
c The factor "cos y" does not depend on the variable of integration
c for the inner integral, and therefore is factored out. Similarly,
c the term y*y in the denominator of the inner integrand is pre-
c computed when the limits of the inner 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
integer IFLAG, IOPT(10), NDIMI, NWPD, NWORK
parameter (NDIMI=2, NWPD=217, NWORK=3*NDIMI+NWPD*(NDIMI-1))
double precision ANSWER, WORK(NWORK)
double precision COSY, DENOM, XDY, YSQ
c
10 format (/' DRDINTMR:'/
1' Compute the integral for Y = 0.0 to Y = PI of'/
2' the integral for X = 0.0 to X = Y of'/
3' 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) = 6
IOPT(5) = 0
call DINTM (NDIMI, ANSWER, WORK, NWORK, IOPT)
30 continue
call DINTMA (ANSWER, WORK, IOPT)
IFLAG = IOPT(1)
if (IFLAG .eq. 0) then
c
c IFLAG = 0, compute innermost integrand.
c
c Iterate on the innermost integrand by calling DINTA here.
c The code would be slightly simpler, but slightly slower,
c if control of the iteration were relegated to DINTMA.
c
40 continue
DENOM = WORK(1) * WORK(1) + YSQ
if (DENOM .ne. 0.0D0) then
c This test will not detect overflow of the integrand if
c the 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
call DINTA (ANSWER, WORK, IOPT)
if(IOPT(1) .eq. 0) goto 40
c
if (IOPT(1) .gt. 0) then
c Done with the integration
IOPT(1) = -(IOPT(1)+NDIMI)
go to 60
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
if (IFLAG+NDIMI .le. 0) goto 60
c
c IFLAG < 0, transform inner integrand.
c
ANSWER = ANSWER * COSY
WORK(1) = WORK(1) * COSY
c
end if
go to 30
c
60 continue
print 20, ANSWER, WORK(1), IOPT(1), IOPT(3)
stop
end