/* nag_quad_1d_gen_vec_multi_rcomm (d01rac) Example Program.
*
* Copyright 2014 Numerical Algorithms Group.
*
* Mark 24, 2013.
*/
#include <stdio.h>
#include <math.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagd01.h>
#include <nagx01.h>
/* Print information on splitting and evaluations over subregions? */
Nag_Boolean disp_integration_info = Nag_TRUE;
static void display_integration_details(const Integer ni,
const Integer iopts[],
const double opts[],
const Integer icom[],
const double com[]);
static void display_option(const char *optstr, const Nag_VariableType optype,
const Integer ivalue, const double rvalue,
const char *cvalue);
int main(void)
{
#define FM(J,I) fm[(I-1)*ldfm + J-1]
/* Scalars */
int exit_status = 0;
Integer len_cvalue;
double a, b, rvalue;
Integer irevcm, ivalue, i, j, lcmax, lcmin, lcom, ldfm, ldfmrq,
lenx, lenxrq, licmax, licmin, licom, liopts, lopts, ni, nx,
sdfm, sdfmrq, sid;
/* Arrays */
char cvalue[17];
double *com = 0, *dinest = 0, *errest = 0, *fm = 0, *opts = 0, *x = 0;
Integer *icom = 0, *iopts = 0, *needi = 0;
/* NAG types */
Nag_VariableType optype;
NagError fail;
printf("nag_quad_1d_gen_vec_multi_rcomm (d01rac) Example Program Results"
"\n\n");
/* Setup phase.*/
/* Set problem parameters. */
ni = 2;
a = 0.0;
b = nag_pi;
liopts = 100;
lopts = 100;
if (
!(opts = NAG_ALLOC((lopts), double))||
!(iopts = NAG_ALLOC((liopts), Integer))
)
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
INIT_FAIL(fail);
/* Initialize option arrays using nag_quad_opt_set (d01zkc). */
nag_quad_opt_set("Initialize = nag_quad_1d_gen_vec_multi_rcomm",
iopts, liopts, opts, lopts, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_quad_opt_set (d01zkc).\n%s\n", fail.message);
exit_status = 1;
goto END;
}
nag_quad_opt_set("Quadrature Rule = gk41", iopts, liopts, opts, lopts, &fail);
nag_quad_opt_set("Absolute Tolerance = 1.0e-7", iopts, liopts, opts, lopts,
&fail);
nag_quad_opt_set("Relative Tolerance = 1.0e-7", iopts, liopts, opts, lopts,
&fail);
/* Determine required array dimensions for
* nag_quad_1d_gen_vec_multi_rcomm (d01rac) using
* nag_quad_1d_gen_vec_multi_dimreq (d01rcc).
*/
nag_quad_1d_gen_vec_multi_dimreq(ni, &lenxrq, &ldfmrq, &sdfmrq,
&licmin, &licmax, &lcmin, &lcmax,
iopts, opts, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_quad_1d_gen_vec_multi_dimreq (d01rcc).\n%s\n",
fail.message);
exit_status = 2;
goto END;
}
ldfm = ldfmrq;
sdfm = sdfmrq;
lenx = lenxrq;
licom = licmax;
lcom = lcmax;
/* Allocate remaining arrays.*/
if (
!(x = NAG_ALLOC((lenx), double))||
!(needi = NAG_ALLOC((ni), Integer))||
!(fm = NAG_ALLOC((ldfm)*(sdfm), double))||
!(dinest = NAG_ALLOC((ni), double))||
!(errest = NAG_ALLOC((ni), double))||
!(com = NAG_ALLOC((lcom), double))||
!(icom = NAG_ALLOC((licom), Integer))
)
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/* Solve phase.*/
/* Use nag_quad_1d_gen_vec_multi_rcomm (d01rac) to evaluate the
* definite integrals of:
* f_1 = (x*sin(2*x))*cos(15*x)
* f_2 = (x*sin(2*x))*(x*cos(50*x))
*/
INIT_FAIL(fail);
/* Set initial irevcm. */
irevcm = 1;
while (irevcm)
{
/* nag_quad_1d_gen_vec_multi_rcomm (d01rac).
* One-dimensional quadrature, adaptive, vectorized, multi-integral,
* reverse communication.
*/
nag_quad_1d_gen_vec_multi_rcomm(&irevcm, ni, a, b,
&sid, needi, x, lenx, &nx, fm, ldfm,
dinest, errest,
iopts, opts, icom, licom, com, lcom,
&fail);
switch (irevcm) {
case 11:
/* Initial returns.
* These will occur during the non-adaptive phase.
* All values must be supplied.
* dinest and errest do not contain approximations over the complete
* interval at this stage.
*/
for (i=1; i<=nx; i++) {
FM(2, i) = x[i-1]*sin(2.0*x[i-1]);
FM(1, i) = FM(2, i)*cos(15.0*x[i-1]);
/* Complete f_2 calculation.*/
FM(2, i) = FM(2, i)*x[i-1]*cos(50.0*x[i-1]);
}
break;
case 12:
/* Intermediate returns.
* These will occur during the adaptive phase.
* All requested values must be supplied.
* dinest and errest contain approximations over the complete
* interval at this stage.
*/
if ( (needi[0]==1) && (needi[1]==1)) {
for (i=1; i<=nx; i++) {
FM(2, i) = x[i-1]*sin(2.0*x[i-1]);
FM(1, i) = FM(2, i)*cos(15.0*x[i-1]);
/* Complete f_2 calculation.*/
FM(2, i) = FM(2, i)*x[i-1]*cos(50.0*x[i-1]);
}
} else if (needi[0]==1) {
/* Only calculation of f_1 is requried.*/
for (i=1; i<=nx; i++)
FM(1, i) = (x[i-1]*sin(2.0*x[i-1]))*(cos(15.0*x[i-1]));
} else if (needi[1]==1) {
/* Only calculation of f_2 is requried.*/
for (i=1; i<=nx; i++)
FM(2, i) = (x[i-1]*sin(2.0*x[i-1]))*(x[i-1]*cos(50.0*x[i-1]));
}
break;
case 0:
/* Final return. Test fail.code.*/
switch (fail.code) {
case NE_NOERROR:
break;
case NE_ACCURACY:
printf("Warning: nag_quad_1d_gen_vec_multi_rcomm (d01rac) has "
"returned at \n least one error estimate exceeding the"
" requested tolerances\n");
break;
case NE_QUAD_BAD_SUBDIV_INT:
/* Useful information has been returned.*/
printf("Warning: nag_quad_1d_gen_vec_multi_rcomm (d01rac) has "
"detected \n extremely bad behaviour for at least"
" one integral\n");
break;
default: ;
/* An unrecoverable error has been detected.*/
printf("Error from nag_quad_1d_gen_vec_multi_rcomm (d01rac).\n%s\n",
fail.message);
exit_status = 3;
goto END;
}
}
}
/* Query some currently set options using nag_quad_opt_get (d01zlc). */
len_cvalue = 17;
nag_quad_opt_get("Quadrature rule", &ivalue, &rvalue, cvalue, len_cvalue,
&optype, iopts, opts, &fail);
display_option("Quadrature rule", optype, ivalue, rvalue, cvalue);
nag_quad_opt_get("Maximum Subdivisions", &ivalue, &rvalue, cvalue, len_cvalue,
&optype, iopts, opts, &fail);
display_option("Maximum Subdivisions", optype, ivalue, rvalue, cvalue);
nag_quad_opt_get("Extrapolation", &ivalue, &rvalue, cvalue, len_cvalue,
&optype, iopts, opts, &fail);
display_option("Extrapolation", optype, ivalue, rvalue, cvalue);
nag_quad_opt_get("Extrapolation Safeguard", &ivalue, &rvalue,
cvalue, len_cvalue, &optype, iopts, opts, &fail);
display_option("Extrapolation safeguard", optype, ivalue, rvalue, cvalue);
/* Print solution.*/
printf("\n Integral | needi | dinest | errest \n");
for ( j=1; j<=ni; j++)
printf(" %9ld %9ld %12.4e %12.4e\n",
j,needi[j-1],dinest[j-1],errest[j-1]);
/* Investigate integration strategy. */
if(disp_integration_info)
display_integration_details(ni,iopts,opts,icom,com);
END:
NAG_FREE(com);
NAG_FREE(dinest);
NAG_FREE(errest);
NAG_FREE(fm);
NAG_FREE(opts);
NAG_FREE(x);
NAG_FREE(icom);
NAG_FREE(iopts);
NAG_FREE(needi);
return exit_status;
}
static void display_integration_details(const Integer ni,
const Integer iopts[],
const double opts[],
const Integer icom[],
const double com[])
{
#define FCOUNT(J) icom[fcp + J-2]
#define EVALS(J,K) icom[evalsp +(K-1)*ldi + J-2]
#define SINFOI(L,K) icom[sinfoip +(K-1)*ldi + L-2]
#define DS(J,K) com[dsp + (K-1)*ldr + J-2]
#define ES(J,K) com[esp + (K-1)*ldr + J-2]
#define SINFOR(L,K) com[sinforp + (K-1)*ldr + L-2]
double lbnd, ubnd,rvalue;
Integer ldi,ldr,sinfoip,sinforp,evalsp,fcp,dsp,esp,nseg,nsdiv;
Integer child1, child2, j, k, level, parent, sid, index,len_cvalue;
char cvalue[17];
NagError fail;
Nag_VariableType optype;
/* Request communication array indices */
INIT_FAIL(fail);
len_cvalue = 17;
nag_quad_opt_get("Index nseg", &index, &rvalue, cvalue, len_cvalue,
&optype, iopts, opts, &fail);
nseg = icom[index-1];
nag_quad_opt_get("Index nsdiv", &index, &rvalue, cvalue, len_cvalue,
&optype, iopts, opts, &fail);
nsdiv = icom[index-1];
nag_quad_opt_get("Index ldi", &index, &rvalue, cvalue, len_cvalue,
&optype, iopts, opts, &fail);
ldi = icom[index-1];
nag_quad_opt_get("Index ldr", &index, &rvalue, cvalue, len_cvalue,
&optype, iopts, opts, &fail);
ldr = icom[index-1];
nag_quad_opt_get("Index fcp", &index, &rvalue, cvalue, len_cvalue,
&optype, iopts, opts, &fail);
fcp = icom[index-1];
nag_quad_opt_get("Index evalsp", &index, &rvalue, cvalue, len_cvalue,
&optype, iopts, opts, &fail);
evalsp = icom[index-1];
nag_quad_opt_get("Index sinfoip", &index, &rvalue, cvalue, len_cvalue,
&optype, iopts, opts, &fail);
sinfoip = icom[index-1];
nag_quad_opt_get("Index dsp", &index, &rvalue, cvalue, len_cvalue,
&optype, iopts, opts, &fail);
dsp = icom[index-1];
nag_quad_opt_get("Index esp", &index, &rvalue, cvalue, len_cvalue,
&optype, iopts, opts, &fail);
esp = icom[index-1];
nag_quad_opt_get("Index sinforp", &index, &rvalue, cvalue, len_cvalue,
&optype, iopts, opts, &fail);
sinforp = icom[index-1];
printf("\n Information on integration:\n ni = %3"NAG_IFMT
", nseg = %3ld, nsdiv = %3ld.",ni,nseg,nsdiv);
for(j=1;j<=ni;j++)
printf("\n Integral %2ld total approximations: %3ld.",
j,FCOUNT(j));
printf("\n\n Information on subdivision and evaluations over segments.\n");
for ( k=1; k<=nseg; k++) {
printf("\n");
sid = SINFOI(1, k);
parent = SINFOI(2, k);
child1 = SINFOI(3, k);
child2 = SINFOI(4, k);
level = SINFOI(5, k);
lbnd = SINFOR(1, k);
ubnd = SINFOR(2, k);
printf(" Segment %3ld.\n Sid = %3ld, Parent = "
"%3ld, Level = %3ld.\n", k, sid, parent, level);
if ( child1>0)
printf(" Children = (%3ld,%3ld).\n",child1,child2);
printf(" Bounds (%11.4e,%11.4e).\n", lbnd,ubnd);
for ( j=1; j<=ni; j++) {
if (EVALS(j, k)>0 && EVALS(j,k)<5) {
printf(" Integral %2ld approximation : %11.4e.\n", j,
DS(j,k));
printf(" Integral %2ld error estimate: %11.4e.\n", j,
ES(j,k));
if ( EVALS(j, k) == 3)
printf(" Integral %2ld evaluation has been superseded by "
"descendants.\n",j);
}
}
}
fflush(stdout);
}
static void display_option(const char *optstr, const Nag_VariableType optype,
const Integer ivalue, const double rvalue,
const char *cvalue)
{
/* Query optype and print the appropriate option values. */
switch (optype) {
case Nag_Integer:
printf(" %30s : %13ld\n", optstr, ivalue);
break;
case Nag_Real:
printf(" %30s : %13.4e\n", optstr, rvalue);
break;
case Nag_Character:
printf(" %30s : %16s\n", optstr, cvalue);
break;
case Nag_Integer_Additional:
printf(" %30s : %3ld %16s\n", optstr, ivalue, cvalue);
break;
case Nag_Real_Additional:
printf(" %30s : %13.4e %16s\n", optstr, rvalue, cvalue);
break;
default: ;
}
fflush(stdout);
}