NAG Library Manual, Mark 27.3
/* nag_rand_field_2d_user_setup (g05zqc) Example Program.
*
* Copyright 2021 Numerical Algorithms Group.
*
* Mark 27.3, 2021.
*/
#include <math.h>
#include <nag.h>

#ifdef __cplusplus
extern "C" {
#endif
static void NAG_CALL cov2(double t1, double t2, double *gamma, Nag_Comm *comm);
#ifdef __cplusplus
}
#endif

static void display_results(Integer approx, Integer *m, double rho, double *eig,
Integer icount, double *lam);
static void read_input_data(Nag_NormType *norm, double *l1, double *l2,
double *nu, double *var, double *xmin, double *xmax,
double *ymin, double *ymax, Integer *ns,
Integer *maxm, Nag_EmbedScale *corr,
int main(void) {
/*  Scalars */
Integer exit_status = 0;
double l1, l2, nu, rho, var, xmax, xmin, ymax, ymin;
Integer approx, icount;
/*  Arrays */
double eig[3];
double *lam = 0, *xx = 0, *yy = 0;
Integer m[2], maxm[2], ns[2];
/* Nag types */
Nag_NormType norm;
Nag_EmbedScale corr;
Nag_Parity even = Nag_Even;
Nag_Comm comm;
NagError fail;

INIT_FAIL(fail);

printf("nag_rand_field_2d_user_setup (g05zqc) Example Program Results\n\n");
/* Get problem specifications from data file */
read_input_data(&norm, &l1, &l2, &nu, &var, &xmin, &xmax, &ymin, &ymax, ns,
if (!(lam = NAG_ALLOC(maxm[0] * maxm[1], double)) ||
!(xx = NAG_ALLOC(ns[0], double)) || !(yy = NAG_ALLOC(ns[1], double)) ||
!(comm.iuser = NAG_ALLOC(1, Integer)) ||
!(comm.user = NAG_ALLOC(3, double))) {
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/* Put covariance parameters in communication arrays */
comm.iuser[0] = (Integer)norm;
comm.user[0] = l1;
comm.user[1] = l2;
comm.user[2] = nu;
/* Get square roots of the eigenvalues of the embedding matrix. These are
* obtained from the setup for simulating two-dimensional random fields,
* with a user-defined variogram, by the circulant embedding method using
* nag_rand_field_2d_user_setup (g05zqc).
*/
nag_rand_field_2d_user_setup(ns, xmin, xmax, ymin, ymax, maxm, var, cov2,
even, pad, corr, lam, xx, yy, m, &approx, &rho,
&icount, eig, &comm, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_rand_field_2d_user_setup (g05zqc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
/* Output results */
display_results(approx, m, rho, eig, icount, lam);
END:
NAG_FREE(lam);
NAG_FREE(xx);
NAG_FREE(yy);
NAG_FREE(comm.iuser);
NAG_FREE(comm.user);
return exit_status;
}

void read_input_data(Nag_NormType *norm, double *l1, double *l2, double *nu,
double *var, double *xmin, double *xmax, double *ymin,
double *ymax, Integer *ns, Integer *maxm,
char nag_enum_arg[40];

/* Read in norm type by name and convert to value using
* nag_enum_name_to_value (x04nac).
*/
scanf("%*[^\n] %39s%*[^\n]", nag_enum_arg);
*norm = (Nag_NormType)nag_enum_name_to_value(nag_enum_arg);
/* read in l1, l2 and nu for cov function */
scanf("%lf %lf %lf%*[^\n]", l1, l2, nu);
/* Read in variance of random field */
scanf("%lf%*[^\n]", var);
/* Read in domain endpoints */
scanf("%lf %lf%*[^\n]", xmin, xmax);
scanf("%lf %lf%*[^\n]", ymin, ymax);
/* Read in number of sample points in each direction */
scanf("%" NAG_IFMT " %" NAG_IFMT "%*[^\n]", &ns[0], &ns[1]);
/* Read in maximum size of embedding matrix */
scanf("%" NAG_IFMT " %" NAG_IFMT "%*[^\n]", &maxm[0], &maxm[1]);
/* Read name of scaling in case of approximation and convert to value. */
scanf(" %39s%*[^\n]", nag_enum_arg);
*corr = (Nag_EmbedScale)nag_enum_name_to_value(nag_enum_arg);
/* Read in choice of padding and convert name to value. */
scanf(" %39s%*[^\n]", nag_enum_arg);
}

void display_results(Integer approx, Integer *m, double rho, double *eig,
Integer icount, double *lam) {
/*  Scalars */
Integer i, j;

/* Display size of embedding matrix */
printf("\nSize of embedding matrix = %" NAG_IFMT "\n\n", m[0] * m[1]);
/* Display approximation information if approximation used. */
if (approx == 1) {
printf("Approximation required\n\n");
printf("rho = %10.5f\n", rho);
printf("eig = ");
for (j = 0; j < 3; j++)
printf("%10.5f", eig[j]);
printf("\nicount = %" NAG_IFMT "\n", icount);
} else {
printf("Approximation not required\n");
}
/* Display square roots of the eigenvalues of the embedding matrix. */
printf("\nSquare roots of eigenvalues of embedding matrix:\n\n");
for (i = 0; i < m[0]; i++) {
for (j = 0; j < m[1]; j++) {
printf("%8.4f", lam[i + j * m[0]]);
}
printf("\n");
}
}

static void NAG_CALL cov2(double t1, double t2, double *gamma, Nag_Comm *comm) {
/*  Scalars */
double l1, l2, nu, rnorm, tc1, tc2;
Integer norm;

/* Covariance parameters stored in user array. */
norm = comm->iuser[0];
l1 = comm->user[0];
l2 = comm->user[1];
nu = comm->user[2];

tc1 = fabs(t1) / l1;
tc2 = fabs(t2) / l2;
if (norm == (Integer)Nag_OneNorm) {
rnorm = tc1 + tc2;
} else if (norm == (Integer)Nag_TwoNorm) {
rnorm = sqrt(tc1 * tc1 + tc2 * tc2);
} else
rnorm = 0.0;
*gamma = exp(-(pow(rnorm, nu)));
}