/* nag_rand_field_1d_generate (g05zpc) Example Program.
*
* Copyright 2024 Numerical Algorithms Group.
*
* Mark 30.3, 2024.
*/
#include <math.h>
#include <nag.h>
#define NPMAX 4
#define LENST 17
#define LSEED 1
static void read_input_data(Nag_Variogram *cov, Integer *np, double *params,
double *var, double *xmin, double *xmax,
Integer *ns, Integer *maxm, Nag_EmbedScale *corr,
Nag_EmbedPad *pad, Integer *s);
static void display_embedding_results(Integer approx, Integer m, double rho,
double *eig, Integer icount, double *lam);
static void initialize_state(Integer *state);
static void display_realizations(Integer ns, Integer s, double *xx, double *z,
Integer *exit_status);
int main(void) {
/* Scalars */
Integer exit_status = 0;
double rho, var, xmax, xmin;
Integer approx, icount, m, maxm, np, ns, s;
/* Arrays */
double eig[3], params[NPMAX];
double *lam = 0, *xx = 0, *z = 0;
Integer state[LENST];
/* Nag types */
Nag_Variogram cov;
Nag_EmbedPad pad;
Nag_EmbedScale corr;
NagError fail;
INIT_FAIL(fail);
printf("nag_rand_field_1d_generate (g05zpc) Example Program Results\n\n");
fflush(stdout);
/* Get problem specifications from data file */
read_input_data(&cov, &np, params, &var, &xmin, &xmax, &ns, &maxm, &corr,
&pad, &s);
if (!(lam = NAG_ALLOC(maxm, double)) || !(xx = NAG_ALLOC(ns, double))) {
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
/* Get square roots of the eigenvalues of the embedding matrix.
* nag_rand_field_1d_predef_setup (g05znc).
* Setup for simulating one-dimensional random fields, preset variogram,
* circulant embedding method
*/
nag_rand_field_1d_predef_setup(ns, xmin, xmax, maxm, var, cov, np, params,
pad, corr, lam, xx, &m, &approx, &rho, &icount,
eig, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_rand_field_1d_predef_setup (g05znc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}
fflush(stdout);
display_embedding_results(approx, m, rho, eig, icount, lam);
/* Initialize state array */
initialize_state(state);
if (!(z = NAG_ALLOC(ns * s, double))) {
printf("Allocation failure\n");
exit_status = -2;
goto END;
}
/* Compute s random field realizations.
* nag_rand_field_1d_generate (g05zpc).
* Generates s realizations of a one-dimensional random field by the
* circulant embedding method.
*/
nag_rand_field_1d_generate(ns, s, m, lam, rho, state, z, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_rand_field_1d_generate (g05zpc).\n%s\n",
fail.message);
exit_status = 2;
goto END;
}
display_realizations(ns, s, xx, z, &exit_status);
END:
NAG_FREE(lam);
NAG_FREE(xx);
NAG_FREE(z);
return exit_status;
}
void read_input_data(Nag_Variogram *cov, Integer *np, double *params,
double *var, double *xmin, double *xmax, Integer *ns,
Integer *maxm, Nag_EmbedScale *corr, Nag_EmbedPad *pad,
Integer *s) {
Integer j;
char nag_enum_arg[40];
/* Read in covariance function name and convert to value using
* nag_enum_name_to_value (x04nac).
*/
scanf("%*[^\n] %39s%*[^\n]", nag_enum_arg);
*cov = (Nag_Variogram)nag_enum_name_to_value(nag_enum_arg);
/* Read in parameters */
scanf("%" NAG_IFMT "%*[^\n]", np);
for (j = 0; j < *np; j++)
scanf("%lf", ¶ms[j]);
scanf("%*[^\n]");
/* 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 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);
*pad = (Nag_EmbedPad)nag_enum_name_to_value(nag_enum_arg);
/* Read in number of realization samples to be generated */
scanf("%" NAG_IFMT "%*[^\n]", s);
}
void display_embedding_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");
fflush(stdout);
}
void initialize_state(Integer *state) {
/* Scalars */
Integer inseed = 14965, lseed = LSEED, subid = 1;
Integer lstate;
/* Arrays */
Integer seed[LSEED];
/* Nag types */
NagError fail;
INIT_FAIL(fail);
lstate = LENST;
seed[0] = inseed;
/* nag_rand_init_repeat (g05kfc).
* Initializes a pseudorandom number generator to give a repeatable sequence.
*/
nag_rand_init_repeat(Nag_Basic, subid, seed, lseed, state, &lstate, &fail);
}
void display_realizations(Integer ns, Integer s, double *xx, double *z,
Integer *exit_status) {
/* Scalars */
Integer indent = 0, ncols = 80;
Integer i;
/* Arrays */
char **rlabs = 0;
/* Nag types */
NagError fail;
INIT_FAIL(fail);
if (!(rlabs = NAG_ALLOC(ns, char *))) {
printf("Allocation failure\n");
*exit_status = -3;
goto END;
}
/* Set row labels to mesh points (column label is realization number). */
for (i = 0; i < ns; i++) {
if (!(rlabs[i] = NAG_ALLOC(11, char))) {
printf("Allocation failure\n");
*exit_status = -4;
goto END;
}
sprintf(rlabs[i], "%10.5f", xx[i]);
}
printf("\n");
fflush(stdout);
/* Display random field results, z, using the comprehensive real general
* matrix print routine nag_file_print_matrix_real_gen_comp (x04cbc).
*/
nag_file_print_matrix_real_gen_comp(
Nag_ColMajor, Nag_GeneralMatrix, Nag_NonUnitDiag, ns, s, z, ns, "%10.5f",
"Random field realizations:", Nag_CharacterLabels, (const char **)rlabs,
Nag_IntegerLabels, NULL, ncols, indent, 0, &fail);
END:
for (i = 0; i < ns; i++) {
NAG_FREE(rlabs[i]);
}
NAG_FREE(rlabs);
}