/* nag_quasi_rand_normal (g05yjc) Example Program. * * Copyright 2009, Numerical Algorithms Group. * * Mark 9, 2009. */ /* Pre-processor includes */ #include #include #include #include #include #include #define QUAS(I, J) quas[(order == Nag_ColMajor)?(J*pdquas + I):(I*pdquas + J)] int main(void) { /* Integer scalar and array declarations */ Integer exit_status = 0; Integer liref, i, j, q_size; Integer *iref = 0; Integer pdquas; /* NAG structures */ NagError fail; /* Double scalar and array declarations */ double *quas = 0; /* Number of dimensions */ Integer idim = 4; /* Mean and standard deviation of the normal distribution */ double xmean[] = { 1.0e0, 2.0e0, 3.0e0, 4.0e0 }; double std[] = { 1.0e0, 1.0e0, 1.0e0, 1.0e0 }; /* Set the sample size */ Integer n = 5; /* Skip the first 1000 variates */ Integer iskip = 1000; /* Use column major order */ Nag_OrderType order = Nag_ColMajor; /* Choose the quasi generator */ Nag_QuasiRandom_Sequence genid = Nag_QuasiRandom_Sobol; /* Initialise the error structure */ INIT_FAIL(fail); printf("nag_quasi_rand_normal (g05yjc) Example Program Results\n\n"); pdquas = (order == Nag_RowMajor)?idim:n; q_size = (order == Nag_RowMajor)?pdquas * n:pdquas * idim; /* Calculate the size of the reference vector */ liref = (genid == Nag_QuasiRandom_Faure)?407:32 * idim + 7; /* Allocate arrays */ if (!(quas = NAG_ALLOC(q_size, double)) || !(iref = NAG_ALLOC(liref, Integer))) { printf("Allocation failure\n"); exit_status = -1; goto END; } /* Initialise the Sobol generator */ nag_quasi_init(genid, idim, iref, liref, iskip, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_quasi_init (g05ylc).\n%s\n", fail.message); exit_status = 1; goto END; } /* Generate a normal quasi-random number sequence */ nag_quasi_rand_normal(order, xmean, std, n, quas, pdquas, iref, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_quasi_rand_normal (g05yjc).\n%s\n", fail.message); exit_status = 1; goto END; } /* Print the estimated value of the integral */ for (i = 0; i < n; i++) { printf(" "); for (j = 0; j < idim; j++) printf("%9.4f%s", QUAS(i, j), ((j+1)%4)?" ":"\n"); if (idim%4) printf("\n"); } END: NAG_FREE(quas); NAG_FREE(iref); return exit_status; }