/* nag_quasi_rand_uniform (g05ymc) Example Program. * * Copyright 2008, Numerical Algorithms Group. * * Mark 9, 2009. */ /* Pre-processor includes */ #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, d, i, j, q_size; Integer *iref = 0; Integer pdquas; /* NAG structures */ NagError fail; /* Double scalar and array declarations */ double sum, tmp, vsbl; double *quas = 0; /* Number of dimensions */ Integer idim = 8; /* Set the sample size */ Integer n = 200; /* Skip the first 1000 variates */ Integer iskip = 1000; /* Use row major order */ Nag_OrderType order = Nag_RowMajor; /* Choose the quasi generator */ Nag_QuasiRandom_Sequence genid = Nag_QuasiRandom_Sobol; /* Initialise the error structure */ INIT_FAIL(fail); printf("nag_quasi_rand_uniform (g05ymc) Example Program Results\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 n quasi-random variates */ nag_quasi_rand_uniform(order, n, quas, pdquas, iref, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_quasi_rand_uniform (g05ymc).\n%s\n", fail.message); exit_status = 1; goto END; } /* Estimate integral by evaluating function at each variate and summing */ sum = 0.0e0; for (i = 0; i < n; i++) { tmp = 1.0e0; for (d = 0; d < idim; d++) tmp *= fabs(4.0e0 * QUAS(i, d) - 2.0e0); sum += tmp; } /* Convert sum to mean value */ vsbl = sum / (double) n; /* Print the estimated value of the integral */ printf("Value of integral = %8.4f\n\n", vsbl); /* Display the first 10 variates used */ printf("First 10 variates\n"); for (i = 0; i < 10; i++) { printf(" %3ld ", i + 1); for (j = 0; j < idim; j++) printf("%8.4f%s", QUAS(i, j), ((j+1)%20)?" ":"\n"); if (idim%20) printf("\n"); } END: if (quas) NAG_FREE(quas); if (iref) NAG_FREE(iref); return exit_status; }