NAG Library Manual, Mark 27.3
/* nag_rand_field_1d_user_setup (g05zmc) 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 cov1(double x, double *gamma, Nag_Comm *comm);
#ifdef __cplusplus
}
#endif
static void read_input_data(double *l, double *nu, double *var, double *xmin,
double *xmax, Integer *ns, Integer *maxm,
static void display_results(Integer approx, Integer m, double rho, double *eig,
Integer icount, double *lam);
int main(void) {
Integer exit_status = 0;
/* Scalars */
double c, nu, rho, var, xmax, xmin;
Integer approx, icount, m, maxm, ns;
/* Arrays */
double eig[3], *lam = 0, *xx = 0;
/* Nag types */
Nag_EmbedScale corr;
Nag_Comm comm;
NagError fail;

INIT_FAIL(fail);

printf("nag_rand_field_1d_user_setup (g05zmc) Example Program Results\n\n");

/* Get problem specifications from data file */
if (!(lam = NAG_ALLOC(maxm, double)) || !(xx = NAG_ALLOC(ns, double)) ||
!(comm.user = NAG_ALLOC(2, double))) {
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/* Put covariance parameters in communication array */
comm.user[0] = c;
comm.user[1] = nu;
/* Get square roots of the eigenvalues of the embedding matrix. These are
* obtained from the setup for simulating one-dimensional random fields,
* with a user-defined variogram, by the circulant embedding method using
* nag_rand_field_1d_user_setup (g05zmc).
*/
nag_rand_field_1d_user_setup(ns, xmin, xmax, maxm, var, cov1, pad, corr, lam,
xx, &m, &approx, &rho, &icount, eig, &comm,
&fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_rand_field_1d_user_setup (g05zmc).\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(comm.user);
return exit_status;
}

void read_input_data(double *l, double *nu, double *var, double *xmin,
double *xmax, Integer *ns, Integer *maxm,
/* Arrays */
char nag_enum_arg[40];
/* Skip heading and get l and nu for cov1 function. */
scanf("%*[^\n] %lf %lf%*[^\n]", l, nu);
/* Read in variance of random field. */
scanf("%lf%*[^\n]", var);
/* Read in domain endpoints. */
scanf("%lf %lf%*[^\n]", xmin, xmax);
/* Read in number of sample points. */
scanf("%" NAG_IFMT "%*[^\n]", ns);
/* Read in maximum size of embedding matrix. */
scanf("%" NAG_IFMT "%*[^\n]", maxm);
/* Read in choice of scaling in case of approximation. */
scanf(" %39s%*[^\n]", nag_enum_arg);
/* nag_enum_name_to_value (x04nac).
* Converts NAG enum member name to value
*/
*corr = (Nag_EmbedScale)nag_enum_name_to_value(nag_enum_arg);
/* Read in choice of padding and convert to enum name to value. */
scanf(" %39s%*[^\n]", nag_enum_arg);
}

void display_results(Integer approx, Integer m, double rho, double *eig,
Integer icount, double *lam) {
Integer j;
/* Display size of embedding matrix */
printf("\nSize of embedding matrix = %" NAG_IFMT "\n\n", m);
/* 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 (j = 0; j < m; j++)
printf("%10.5f%s", lam[j], j % 4 == 3 ? "\n" : "");
printf("\n");
}

static void NAG_CALL cov1(double x, double *gamma, Nag_Comm *comm) {
/* Scalars */
double dummy, l, nu;

/* Correlation length and exponent in comm->ruser. */
l = comm->user[0];
nu = comm->user[1];
if (x == 0.0) {
*gamma = 1.0;
} else {
dummy = pow(x / l, nu);
*gamma = exp(-dummy);
}
}