/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:33:10 */
/*FOR_C Options SET: ftn=u io=c no=p op=aimnv pf=,p_dintmr s=dbov str=l x=f - prototypes */
#include
#include "fcrt.h"
#include
#include
#include "p_dintmr.h"
/* program DRDINTMR
*>> 2003-03-30 DRDINTMR Krogh Changed way of dealing with 0 denominator.
*>> 2001-05-22 DRDINTMR Krogh Minor change for making .f90 version.
*>> 1994-10-19 DRDINTMR Krogh Changes to use M77CON
*>> 1994-08-31 DRDINTMR Snyder Moved formats for C conversion
*>> 1994-08-08 DRDINTMR Snyder Took '0' out of formats for C conversion
*>> 1991-11-20 CLL Edited for Fortran 90.
*>> 1987-12-09 DRDINTMR Snyder Initial Code.
*--D replaces "?": DR?INTMR, ?INTM, ?INTMA, ?INTA
*
* DEMO DRIVER for multi-dimensional quadrature subprogram DINTM.
*
* Compute the integral for Y = 0.0 to Y = PI of
* the integral for X = 0.0 to X = Y of
* X * COS(Y) /(X*X+Y*Y). The ANSWER should be zero.
*
* The integrand and all limits are provided by reverse
* communication.
*
* This sample demonstrates the use of functional transformation of
* the inner integral to reduce the cost of the overall computation.
* The factor "cos y" does not depend on the variable of integration
* for the inner integral, and therefore is factored out. Similarly,
* the term y*y in the denominator of the inner integrand is pre-
* computed when the limits of the inner integral are requested.
*
* This example also demonstrates the necessity to cope with overflow
* of the integrand, even though the integral is well behaved.
* */
/* PARAMETER translations */
#define NDIMI 2
#define NWORK (3*NDIMI + NWPD*(NDIMI - 1))
#define NWPD 217
/* end of PARAMETER translations */
int main( )
{
long int iflag, iopt[10];
double answer, cosy, denom, work[NWORK], xdy, ysq;
/* OFFSET Vectors w/subscript range: 1 to dimension */
long *const Iopt = &iopt[0] - 1;
double *const Work = &work[0] - 1;
/* end of OFFSET VECTORS */
printf("\n DRDINTMR:\n Compute the integral for Y = 0.0 to Y = PI of\n the integral for X = 0.0 to X = Y of\n X * COS(Y) /(X*X+Y*Y). The ANSWER should be zero.\n");
Iopt[2] = 10;
Iopt[3] = 0;
Iopt[4] = 6;
Iopt[5] = 0;
dintm( NDIMI, &answer, work, NWORK, iopt );
L_30:
;
dintma( &answer, work, iopt );
iflag = Iopt[1];
if (iflag == 0)
{
/* IFLAG = 0, compute innermost integrand.
*
* Iterate on the innermost integrand by calling DINTA here.
* The code would be slightly simpler, but slightly slower,
* if control of the iteration were relegated to DINTMA.
* */
L_40:
;
denom = Work[1]*Work[1] + ysq;
if (denom != 0.0e0)
{
/* This test will not detect overflow of the integrand if
* the arithmetic underflows gradually. */
answer = Work[1]/denom;
}
else
{
/* Special care to avoid a zero denominator. */
xdy = Work[1]/Work[2];
answer = xdy/(Work[2]*(1.e0 + SQ(xdy)));
}
dinta( &answer, work, iopt );
if (Iopt[1] == 0)
goto L_40;
if (Iopt[1] > 0)
{
/* Done with the integration */
Iopt[1] = -(Iopt[1] + NDIMI);
goto L_60;
}
}
else if (iflag == 1)
{
/* Compute limits of inner dimension.
*
* Set WORK(1) = COS(WORK(2)) = Partial derivative, with respect to
* the integral over the inner dimension, of the transformation
* applied when integration over the inner dimension is complete.
* This is so dintm knows how much accuracy it needs for this integral.
* */
Work[NDIMI + iflag] = 0.0e0;
Work[2*NDIMI + iflag] = Work[2];
ysq = Work[2]*Work[2];
cosy = cos( Work[2] );
Work[1] = cosy;
}
else if (iflag == 2)
{
/* Compute limits of outer dimension.
* */
Work[NDIMI + iflag] = 0.0e0;
Work[2*NDIMI + iflag] = 4.0e0*atan( 1.0e0 );
}
else
{
if (iflag + NDIMI <= 0)
goto L_60;
/* IFLAG < 0, transform inner integrand.
* */
answer *= cosy;
Work[1] *= cosy;
}
goto L_30;
L_60:
;
printf("\n ANSWER = %15.8g\n ERROR ESTIMATE =%15.8g\n STATUS FLAG = %3ld\n FUNCTION VALUES (INNERMOST INTEGRALS) =%6ld\n",
answer, Work[1], Iopt[1], Iopt[3]);
exit(0);
} /* end of function */