/* nag_quad_md_sgq_multi_vec (d01esc) Example Program.
*
* Copyright 2014 Numerical Algorithms Group.
*
* Mark 25, 2014.
*/
#include <stdio.h>
#include <math.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagd01.h>
#include <nagx06.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 *private;
} 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.private = 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.private[thdnum].s = NAG_ALLOC(maxnx, double)) ||
!(parcom.private[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;
goto END;
}
/* 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|%8ld\n",
i, dinest[i], errest[i], ivalid[i]);
END:
NAG_FREE(maxdlv);
NAG_FREE(dinest);
NAG_FREE(errest);
NAG_FREE(ivalid);
NAG_FREE(iopts);
NAG_FREE(opts);
for (thdnum = 0; thdnum < smpthd; ++thdnum) {
NAG_FREE(parcom.private[thdnum].s);
NAG_FREE(parcom.private[thdnum].logs);
}
NAG_FREE(parcom.private);
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) = Σ k*x_i(k).
* k=1
*
* Split the S expression into two components, one involving only the
* 'trivial' value xtr:
*
* ndim ndim
* S(i) = Σ (k*xtr) + Σ (k*(x_i(k)-xtr))
* k=1 k=1
*
* ndim*(ndim+1) ndim
* = xtr * ------------- + Σ (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 = comm->p;
s = (*parcom).private[tid].s;
logs = (*parcom).private[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];
}