/* nag_pde_bs_1d (d03ncc) Example Program.
 *
 * Copyright 2001 Numerical Algorithms Group.
 *
 * Mark 7, 2001.
 */

#include <stdio.h>
#include <string.h>
#include <math.h>
#include <nag.h>
#include <nagx04.h>
#include <nag_stdlib.h>
#include <nagd03.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(int argc, char *argv[])
{
  FILE        *fpin, *fpout;
  double      alpha, x;
  Integer     i, igreek, j, ns, nt, ntkeep, exit_status;
  double      *delta, *f, *gamma, *lambda, q[3], r[3], *rho, *s;
  double      sigma[3], *t, *theta, 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);

  /* Check for command-line IO options */
  fpin = nag_example_file_io(argc, argv, "-data", NULL);
  fpout = nag_example_file_io(argc, argv, "-results", NULL);

  fprintf(fpout, "nag_pde_bs_1d (d03ncc) Example Program Results\n\n");

  /* Skip heading in data file */

  fscanf(fpin, "%*[^\n] ");
  exit_status = 0;

  /* Read problem parameters */

  fscanf(fpin, "%lf", &x);
  fscanf(fpin, "%ld%ld", &ns, &nt);
  fscanf(fpin, "%lf%lf", &smin, &smax);
  fscanf(fpin, "%lf%lf", &tmin, &tmax);
  fscanf(fpin, "%lf", &alpha);
  fscanf(fpin, "%ld", &ntkeep);

  /* Allocate memory */

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

  /* Set up input parameters for nag_pde_bs_1d (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;
  q[0] = 0.0;
  r[0] = 0.10;
  sigma[0] = 0.4;

  /* Call Black-Scholes solver */

  /* nag_pde_bs_1d (d03ncc).
   * Finite difference solution of the Black-Scholes equations
   */
  nag_pde_bs_1d(Nag_AmericanPut, x, Nag_UniformMesh, ns, s,
                nt, t, tdpar, r, q, sigma, alpha, ntkeep, f,
                theta, delta, gamma, lambda, rho, &fail);

  if (fail.code != NE_NOERROR)
    {
      fprintf(fpout, "Error from nag_pde_bs_1d (d03ncc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }

  /* Output option values */

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

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

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

      for (i = 1; i <= ns; i++)
        {
          fprintf(fpout, " %13.4e  | ", s[i-1]);
          switch (igreek)
            {
            case 0:
              for (j = 1; j <= ntkeep; j++)
                fprintf(fpout, " %13.4e", THETA(i, j));
              break;
            case 1:
              for (j = 1; j <= ntkeep; j++)
                fprintf(fpout, " %13.4e", DELTA(i, j));
              break;
            case 2:
              for (j = 1; j <= ntkeep; j++)
                fprintf(fpout, " %13.4e", GAMMA(i, j));
              break;
            case 3:
              for (j = 1; j <= ntkeep; j++)
                fprintf(fpout, " %13.4e", LAMBDA(i, j));
              break;
            case 4:
              for (j = 1; j <= ntkeep; j++)
                fprintf(fpout, " %13.4e", RHO(i, j));
              break;
            default:
              break;
            }
          fprintf(fpout, "\n");
        }
    }

 END:
  if (fpin != stdin) fclose(fpin);
  if (fpout != stdout) fclose(fpout);
  if (s) NAG_FREE(s);
  if (t) NAG_FREE(t);
  if (f) NAG_FREE(f);
  if (theta) NAG_FREE(theta);
  if (delta) NAG_FREE(delta);
  if (gamma) NAG_FREE(gamma);
  if (lambda) NAG_FREE(lambda);
  if (rho) NAG_FREE(rho);

  return exit_status;
}