/* nag_pde_dim1_blackscholes_closed (d03ndc) Example Program.
 *
 * Copyright 2024 Numerical Algorithms Group.
 *
 * Mark 30.1, 2024.
 */

#include <math.h>
#include <nag.h>
#include <stdio.h>
#include <string.h>

#define F(I, J) f[ns * ((J)-1) + (I)-1]
#define THETA(I, J) theta[ns * ((J)-1) + (I)-1]
#define DELTA(I, J) delta[ns * ((J)-1) + (I)-1]
#define GAMMA(I, J) gamma[ns * ((J)-1) + (I)-1]
#define LAMBDA(I, J) lambda[ns * ((J)-1) + (I)-1]
#define RHO(I, J) rho[ns * ((J)-1) + (I)-1]

int main(void) {
  double ds, dt, tmat, x;
  Integer i, igreek, j, ns, nt, exit_status = 0;
  double *delta = 0, *f = 0, *gamma = 0, *lambda = 0, q[3], r[3], *rho = 0,
         *s = 0;
  double sigma[3], *t = 0, *theta = 0, smin, smax, tmin, tmax;
  Nag_Boolean gprnt[5] = {Nag_TRUE, Nag_TRUE, Nag_TRUE, Nag_TRUE, Nag_TRUE};
  Nag_Boolean tdpar[3];
  const char *gname[5] = {"Theta", "Delta", "Gamma", "Lambda", "Rho"};
  NagError fail;

  INIT_FAIL(fail);

  printf(
      "nag_pde_dim1_blackscholes_closed (d03ndc) Example Program Results\n\n");

  /* Skip heading in data file */
  scanf("%*[^\n] ");

  /* Read problem parameters */
  scanf("%lf", &x);
  scanf("%lf", &tmat);
  scanf("%lf", &r[0]);
  scanf("%lf", &q[0]);
  scanf("%lf", &sigma[0]);
  scanf("%" NAG_IFMT "%" NAG_IFMT "", &ns, &nt);
  scanf("%lf%lf", &smin, &smax);
  scanf("%lf%lf", &tmin, &tmax);

  /* Allocate memory */

  if (!(s = NAG_ALLOC(ns, double)) || !(t = NAG_ALLOC(nt, double)) ||
      !(f = NAG_ALLOC(ns * nt, double)) ||
      !(theta = NAG_ALLOC(ns * nt, double)) ||
      !(delta = NAG_ALLOC(ns * nt, double)) ||
      !(gamma = NAG_ALLOC(ns * nt, double)) ||
      !(lambda = NAG_ALLOC(ns * nt, double)) ||
      !(rho = NAG_ALLOC(ns * nt, double))) {
    printf("Allocation failure\n");
    exit_status = 1;
    goto END;
  }

  /* Set up input parameters for nag_pde_dim1_blackscholes_fd (d03ncc) */

  s[0] = smin;
  s[ns - 1] = smax;
  t[0] = tmin;
  t[nt - 1] = tmax;
  tdpar[0] = Nag_FALSE;
  tdpar[1] = Nag_FALSE;
  tdpar[2] = Nag_FALSE;

  ds = (s[ns - 1] - s[0]) / (ns - 1.0);
  dt = (t[nt - 1] - t[0]) / (nt - 1.0);

  /* Loop over times */

  for (j = 1; j <= nt; j++) {
    t[j - 1] = t[0] + (j - 1) * dt;

    /* Loop over stock prices */

    for (i = 1; i <= ns; i++) {
      s[i - 1] = s[0] + (i - 1) * ds;

      /* Evaluate analytic solution of Black-Scholes equation */

      /* nag_pde_dim1_blackscholes_closed (d03ndc).
       * Analytic solution of the Black-Scholes equations
       */
      nag_pde_dim1_blackscholes_closed(Nag_AmericanCall, x, s[i - 1], t[j - 1],
                                       tmat, tdpar, r, q, sigma, &F(i, j),
                                       &THETA(i, j), &DELTA(i, j), &GAMMA(i, j),
                                       &LAMBDA(i, j), &RHO(i, j), &fail);

      if (fail.code != NE_NOERROR) {
        printf("Error from nag_pde_dim1_blackscholes_closed (d03ndc).\n%s\n",
               fail.message);
        exit_status = 1;
        goto END;
      }
    }
  }

  /* Output option values */

  printf("\n");
  printf("Option Values\n");
  printf("-------------\n");
  printf("%14s  |  %s\n", "Stock Price", "Time to Maturity (months)");
  printf("%14s  | ", "");
  for (i = 0; i < nt; i++)
    printf(" %13.4e", 12.0 * (tmat - t[i]));
  printf("\n");
  for (i = 0; i < 74; i++)
    printf("-");
  printf("\n");
  for (i = 1; i <= ns; i++) {
    printf(" %13.4e  | ", s[i - 1]);
    for (j = 1; j <= nt; j++)
      printf(" %13.4e", F(i, j));
    printf("\n");
  }

  for (igreek = 0; igreek < 5; igreek++) {
    if (!gprnt[igreek])
      continue;

    printf("\n");
    printf("%s\n", gname[igreek]);
    for (i = 0; i < (Integer)strlen(gname[igreek]); i++)
      printf("-");
    printf("\n");
    printf("%14s  |  %s\n", "Stock Price", "Time to Maturity (months)");
    printf("%14s  | ", "");
    for (i = 0; i < nt; i++)
      printf(" %13.4e", 12.0 * (tmat - t[i]));
    printf("\n");
    for (i = 0; i < 74; i++)
      printf("-");
    printf("\n");

    for (i = 1; i <= ns; i++) {
      printf(" %13.4e  | ", s[i - 1]);
      switch (igreek) {
      case 0:
        for (j = 1; j <= nt; j++)
          printf(" %13.4e", THETA(i, j));
        break;
      case 1:
        for (j = 1; j <= nt; j++)
          printf(" %13.4e", DELTA(i, j));
        break;
      case 2:
        for (j = 1; j <= nt; j++)
          printf(" %13.4e", GAMMA(i, j));
        break;
      case 3:
        for (j = 1; j <= nt; j++)
          printf(" %13.4e", LAMBDA(i, j));
        break;
      case 4:
        for (j = 1; j <= nt; j++)
          printf(" %13.4e", RHO(i, j));
        break;
      default:
        break;
      }
      printf("\n");
    }
  }
END:
  NAG_FREE(s);
  NAG_FREE(t);
  NAG_FREE(f);
  NAG_FREE(theta);
  NAG_FREE(delta);
  NAG_FREE(gamma);
  NAG_FREE(lambda);
  NAG_FREE(rho);

  return exit_status;
}