/* nag_quad_md_sgq_multi_vec (d01esc) Example Program.
*
* Copyright 2022 Numerical Algorithms Group.
*
* Mark 28.7, 2022.
*/
#include <math.h>
#include <nag.h>
#include <stdio.h>
#ifdef __cplusplus
extern "C" {
#endif
static void NAG_CALL f(Integer ni, Integer ndim, Integer nx, double xtr,
Integer nntr, const Integer *icolzp,
const Integer *irowix, const double *xs,
const Integer *qs, double *fm, Integer *iflag,
Nag_Comm *comm);
#ifdef __cplusplus
}
#endif
/* We define some structures to serve as a demonstration of safely operating
* with the NAG communication structure comm when running in parallel.
*/
/* par_sh: any information to be shared between threads in the function f. */
typedef struct {
double s_tr;
} par_sh;
/* par_pr: any private workspace that a single thread will require in the
* execution of the function f.
*/
typedef struct {
double *s;
double *logs;
} par_pr;
/* parallel_comm: a container for par_sh and par_pr. */
typedef struct {
par_sh shared;
par_pr *tprivate;
} parallel_comm;
int main(void) {
Integer exit_status = 0;
Integer ndim, ni;
Integer maxnx, smpthd, lcvalue;
double rvalue;
char cvalue[16];
Integer *ivalid, *iopts, *maxdlv;
double *dinest, *errest, *opts;
parallel_comm parcom;
int i, thdnum;
/* Nag Types */
Nag_VariableType optype;
Nag_Comm comm;
NagError fail;
INIT_FAIL(fail);
printf("nag_quad_md_sgq_multi_vec (d01esc) Example Program Results\n");
ni = 10;
ndim = 4;
lcvalue = 16;
if (!(iopts = NAG_ALLOC(100, Integer)) || !(opts = NAG_ALLOC(100, double)) ||
!(maxdlv = NAG_ALLOC(ndim, Integer)) ||
!(dinest = NAG_ALLOC(ni, double)) || !(errest = NAG_ALLOC(ni, double)) ||
!(ivalid = NAG_ALLOC(ni, Integer))) {
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/* Initialize option arrays. */
nag_quad_opt_set("Initialize = d01esc", iopts, 100, opts, 100, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_quad_opt_set (d01zkc).\n%s\n", fail.message);
exit_status = 1;
goto END;
}
/* Set any required options. */
nag_quad_opt_set("Absolute Tolerance = 0.0", iopts, 100, opts, 100, &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("Relative Tolerance = 1.0e-3", iopts, 100, opts, 100, &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("Maximum Level = 6", iopts, 100, opts, 100, &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("Index Level = 5", iopts, 100, opts, 100, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_quad_opt_set (d01zkc).\n%s\n", fail.message);
exit_status = 1;
goto END;
}
/* Set any required maximum dimension levels. */
for (i = 0; i < ndim; ++i)
maxdlv[i] = 0;
/* As a demonstration of safely operating with the communication structure
* comm when running in parallel, we will create an instance of our
* parallel_comm structure with fields indexed by the current thread number.
* The size of the array subfields in this structure are a function of
* Maximum Nx.
*/
nag_quad_opt_get("Maximum Nx", &maxnx, &rvalue, cvalue, lcvalue, &optype,
iopts, opts, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_quad_opt_get (d01zlc).\n%s\n", fail.message);
exit_status = 1;
goto END;
}
smpthd = nag_omp_get_max_threads();
/* Allocate an array of smpthd pointers to private structures. */
if (!(parcom.tprivate = NAG_ALLOC(smpthd, par_pr))) {
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/* For each thread, allocate a double array of size maxnx in the
* private component of the parallel structure.
*/
for (thdnum = 0; thdnum < smpthd; thdnum++) {
if (!(parcom.tprivate[thdnum].s = NAG_ALLOC(maxnx, double)) ||
!(parcom.tprivate[thdnum].logs = NAG_ALLOC(maxnx, double))) {
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
}
/* Finally, store the parallel structure in comm for communication to f. */
comm.p = &parcom;
/* Approximate the integrals. */
nag_quad_md_sgq_multi_vec(ni, ndim, f, maxdlv, dinest, errest, ivalid, iopts,
opts, &comm, &fail);
if (fail.code != NE_NOERROR && fail.code != NE_ACCURACY &&
fail.code != NE_TOTAL_PRECISION_LOSS && fail.code != NE_USER_STOP) {
/* If internal memory allocation failed consider reducing the options
* 'Maximum Nx' and 'Index Level', or run with fewer threads.
*/
printf("Error from nag_quad_md_sgq_multi_vec (d01esc).\n%s\n",
fail.message);
exit_status = 2;
} else {
/* NE_NOERROR:
* The result returned satisfies the requested accuracy requirements.
* NE_ACCURACY, NE_TOTAL_PRECISION_LOSS:
* The result returned is inaccurate for at least one integral.
* NE_USER_STOP:
* Exit was requested by setting iflag negative in f.
* A result will be returned if at least one call to f was successful.
*/
printf("Integral # | Estimated value | Error estimate |"
" Final state of integral\n");
for (i = 0; i < ni; ++i)
printf("%11d|%17.5e|%16.5e|%8" NAG_IFMT "\n", i, dinest[i], errest[i],
ivalid[i]);
}
for (thdnum = 0; thdnum < smpthd; ++thdnum) {
NAG_FREE(parcom.tprivate[thdnum].s);
NAG_FREE(parcom.tprivate[thdnum].logs);
}
END:
NAG_FREE(maxdlv);
NAG_FREE(dinest);
NAG_FREE(errest);
NAG_FREE(ivalid);
NAG_FREE(iopts);
NAG_FREE(opts);
NAG_FREE(parcom.tprivate);
return exit_status;
}
static void NAG_CALL f(Integer ni, Integer ndim, Integer nx, double xtr,
Integer nntr, const Integer *icolzp,
const Integer *irowix, const double *xs,
const Integer *qs, double *fm, Integer *iflag,
Nag_Comm *comm) {
Integer i, j, k, tid;
double s_tr, s_ntr;
double *s, *logs;
parallel_comm *parcom;
/* For each evaluation point x_i, i = 1, ..., nx, return in fm the computed
* values of the ni integrals f_j, j = 1, ..., ni defined by
*
* fm(j,i) = f_j(x_i)
* ndim
* = sin(j + S(i))*log(S(i)), where S(i) = Sum k*x_i(k).
* k=1
*
* Split the S expression into two components, one involving only the
* 'trivial' value xtr:
*
* ndim ndim
* S(i) = Sum (k*xtr) + Sum (k*(x_i(k)-xtr))
* k=1 k=1
*
* ndim*(ndim+1) ndim
* = xtr * ------------- + Sum (k*(x_i(k)-xtr))
* 2 k=1
*
* := s_tr + s_ntr(i)
*
* By definition the summands in the s_ntr(i) term on the right-hand side
* are zero for those k outside the range of indices defined in irowix.
*/
/* As a demonstration of safely operating with the communication structure
* comm when running in parallel, the p field of comm is itself (a pointer
* to) an instance of our parallel_comm structure 'partitioned' by the current
* thread number. Store some of the s_tr and s_ntr computations in these.
*/
/* The thread number. */
tid = nag_omp_get_thread_num();
/* The S and log(S) terms from above, extracted from comm. */
parcom = (parallel_comm *)comm->p;
s = (*parcom).tprivate[tid].s;
logs = (*parcom).tprivate[tid].logs;
if (*iflag == 0) {
/* First call: nx=1, no non-trivial dimensions.
* The constant s_tr can be reused by all subsequent calculations.
*/
s_tr = 0.5 * xtr * ((double)(ndim * (ndim + 1)));
parcom->shared.s_tr = s_tr;
s[0] = s_tr;
logs[0] = log(s_tr);
} else {
/* Calculate S(i) = s_tr + s_ntr(i). */
s_tr = parcom->shared.s_tr;
for (i = 0; i < nx; ++i) {
s_ntr = 0.0;
for (k = icolzp[i]; k < icolzp[i + 1]; ++k)
s_ntr = s_ntr + ((double)irowix[k - 1]) * (xs[k - 1] - xtr);
s[i] = s_tr + s_ntr;
logs[i] = log(s[i]);
}
}
/* Finally we obtain fm(j,:) = sin(j+S(:))*log(S(:)). */
for (j = 0; j < ni; ++j)
for (i = 0; i < nx; ++i)
fm[i * ni + j] = sin(((double)j + 1) + s[i]) * logs[i];
}