/* nag_rand_bb_init (g05xac) Example Program.
 *
 * Copyright 2023 Numerical Algorithms Group.
 *
 * Mark 29.1, 2023.
 */
#include <math.h>
#include <nag.h>
#include <stdio.h>
int get_z(Integer nelements, double *z);
void display_results(Nag_OrderType order, Integer npaths, Integer ntimes,
                     Integer d, double *b, Integer pdb);

#define CHECK_FAIL(name, fail)                                                 \
  if (fail.code != NE_NOERROR) {                                               \
    printf("Error calling %s\n%s\n", name, fail.message);                      \
    exit_status = -1;                                                          \
    goto END;                                                                  \
  }

int main(void) {
#define C(I, J) c[(J - 1) * pdc + I - 1]
  Integer exit_status = 0;
  NagError fail;
  /*  Scalars */
  double t0, tend;
  Integer a, d, pdb, pdc, pdz, nmove, npaths, ntimes, i;
  /*  Arrays */
  double *b = 0, *c = 0, *intime = 0, *rcomm = 0, *start = 0, *term = 0,
         *times = 0, *z = 0;
  Integer *move = 0;
  INIT_FAIL(fail);

  /* Parameters which determine the bridge */
  ntimes = 10;
  t0 = 0.0;
  npaths = 2;
  a = 0;
  nmove = 0;
  d = 3;
  pdz = d * (ntimes + 1 - a);
  pdb = d * (ntimes + 1);
  pdc = d;
  /* Allocate memory */
  if (!(intime = NAG_ALLOC((ntimes), double)) ||
      !(times = NAG_ALLOC((ntimes), double)) ||
      !(rcomm = NAG_ALLOC((12 * (ntimes + 1)), double)) ||
      !(start = NAG_ALLOC(d, double)) || !(term = NAG_ALLOC(d, double)) ||
      !(c = NAG_ALLOC(pdc * d, double)) ||
      !(z = NAG_ALLOC(pdz * npaths, double)) ||
      !(b = NAG_ALLOC(pdb * npaths, double)) ||
      !(move = NAG_ALLOC(nmove, Integer))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  /* Fix the time points at which the bridge is required */
  for (i = 0; i < ntimes; i++) {
    intime[i] = t0 + (double)(i + 1);
  }
  tend = t0 + (double)(ntimes + 1);

  /* g05xec.  Creates a Brownian bridge construction order */
  /* out of a set of input times */
  nag_rand_bb_make_bridge_order(Nag_RLRoundDown, t0, tend, ntimes, intime,
                                nmove, move, times, &fail);
  CHECK_FAIL("nag_rand_bb_make_bridge_order", fail);

  /* nag_rand_bb_init (g05xac). Initializes the Brownian bridge generator   */
  nag_rand_bb_init(t0, tend, times, ntimes, rcomm, &fail);
  CHECK_FAIL("nag_rand_bb_init", fail);

  /* We want the following covariance matrix */
  C(1, 1) = 6.0;
  C(2, 1) = 1.0;
  C(3, 1) = -0.2;
  C(1, 2) = 1.0;
  C(2, 2) = 5.0;
  C(3, 2) = 0.3;
  C(1, 3) = -0.2;
  C(2, 3) = 0.3;
  C(3, 3) = 4.0;
  /* nag_rand_bb uses the Cholesky factorization of the covariance matrix C */
  /* f07fdc. Cholesky factorization of real positive definite matrix */
  nag_lapacklin_dpotrf(Nag_ColMajor, Nag_Lower, d, c, pdc, &fail);
  CHECK_FAIL("nag_lapacklin_dpotrf", fail);

  /* Generate the random numbers */
  if (get_z(npaths * d * (ntimes + 1 - a), z) != 0) {
    printf("Error generating random numbers\n");
    exit_status = -1;
    goto END;
  }

  for (i = 0; i < d; i++)
    start[i] = 0.0;
  /* g05xbc. Generate paths for a free or non-free Wiener process using the */
  /* Brownian bridge algorithm   */
  nag_rand_bb(Nag_RowMajor, npaths, d, start, a, term, z, pdz, c, pdc, b, pdb,
              rcomm, &fail);
  CHECK_FAIL("nag_rand_bb", fail);

  /* Display the results */
  display_results(Nag_RowMajor, npaths, ntimes, d, b, pdb);
END:;
  NAG_FREE(b);
  NAG_FREE(c);
  NAG_FREE(intime);
  NAG_FREE(rcomm);
  NAG_FREE(start);
  NAG_FREE(term);
  NAG_FREE(times);
  NAG_FREE(z);
  NAG_FREE(move);
  return exit_status;
}

int get_z(Integer nelements, double *z) {
  NagError fail;
  Integer lseed, lstate, exit_status = 0;
  /*  Arrays */
  Integer seed[1];
  Integer state[80];
  lstate = 80;
  lseed = 1;
  INIT_FAIL(fail);

  /* We now need to generate the input pseudorandom numbers */
  seed[0] = 1023401;
  /* g05kfc. Initializes a pseudorandom number generator */
  /* to give a repeatable sequence */
  nag_rand_init_repeat(Nag_MRG32k3a, 0, seed, lseed, state, &lstate, &fail);
  CHECK_FAIL("nag_rand_init_repeat", fail);

  /* g05skc.  Generates a vector of pseudorandom numbers from */
  /* a Normal distribution */
  nag_rand_dist_normal(nelements, 0.0, 1.0, state, z, &fail);
  CHECK_FAIL("nag_rand_dist_normal", fail);

END:
  return exit_status;
}

void display_results(Nag_OrderType order, Integer npaths, Integer ntimes,
                     Integer d, double *b, Integer pdb) {
#define B(I, J)                                                                \
  (order == Nag_RowMajor ? b[(I - 1) * pdb + J - 1] : b[(J - 1) * pdb + I - 1])

  Integer i, p, k;
  printf("nag_rand_bb_init (g05xac) Example Program Results\n\n");
  for (p = 1; p <= npaths; p++) {
    printf("Wiener Path ");
    printf("%1" NAG_IFMT " ", p);
    printf(",  ");
    printf("%1" NAG_IFMT " ", ntimes + 1);
    printf(" time steps, ");
    printf("%1" NAG_IFMT " ", d);
    printf(" dimensions \n");

    for (k = 1; k <= d; k++) {
      printf("%10" NAG_IFMT " ", k);
    }
    printf("\n");

    for (i = 1; i <= ntimes + 1; i++) {
      printf("%2" NAG_IFMT " ", i);
      for (k = 1; k <= d; k++) {
        printf("%10.4f", B(p, k + (i - 1) * d));
      }
      printf("\n");
    }
    printf("\n");
  }
}