/* nag_quad_dim1_gen_vec_multi_rcomm (d01rac) Example Program.
*
* Copyright 2024 Numerical Algorithms Group.
*
* Mark 30.1, 2024.
*/
#include <math.h>
#include <nag.h>
#include <stdio.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_dim1_gen_vec_multi_rcomm (d01rac) Example Program Results"
"\n\n");
/* Setup phase. */
/* Set problem parameters. */
ni = 2;
a = 0.0;
b = nag_math_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_dim1_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_dim1_gen_vec_multi_rcomm (d01rac) using
* nag_quad_dim1_gen_vec_multi_dimreq (d01rcc).
*/
nag_quad_dim1_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_dim1_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_dim1_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_dim1_gen_vec_multi_rcomm (d01rac).
* One-dimensional quadrature, adaptive, vectorized, multi-integral,
* reverse communication.
*/
nag_quad_dim1_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_dim1_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_dim1_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_dim1_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(" %9" NAG_IFMT " %9" NAG_IFMT " %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 = %3" NAG_IFMT ", nsdiv = %3" NAG_IFMT ".",
ni, nseg, nsdiv);
for (j = 1; j <= ni; j++)
printf("\n Integral %2" NAG_IFMT " total approximations: %3" NAG_IFMT ".",
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 %3" NAG_IFMT ".\n Sid = %3" NAG_IFMT ", Parent = "
"%3" NAG_IFMT ", Level = %3" NAG_IFMT ".\n",
k, sid, parent, level);
if (child1 > 0)
printf(" Children = (%3" NAG_IFMT ",%3" NAG_IFMT ").\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 %2" NAG_IFMT " approximation : %11.4e.\n", j,
DS(j, k));
printf(" Integral %2" NAG_IFMT " error estimate: %11.4e.\n", j,
ES(j, k));
if (EVALS(j, k) == 3)
printf(" Integral %2" NAG_IFMT " 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 : %13" NAG_IFMT "\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 : %3" NAG_IFMT " %16s\n", optstr, ivalue, cvalue);
break;
case Nag_Real_Additional:
printf(" %30s : %13.4e %16s\n", optstr, rvalue, cvalue);
break;
default:;
}
fflush(stdout);
}