/* nag_ode_ivp_bdf_gen (d02ejc) Example Program.
 *
 * Copyright 2017 Numerical Algorithms Group.
 *
 * Mark 26.1, 2017.
 *
 */

#include <nag.h>
#include <math.h>
#include <stdio.h>
#include <nag_stdlib.h>
#include <nagd02.h>

#ifdef __cplusplus
extern "C"
{
#endif
  static void NAG_CALL fcn(Integer neq, double x, const double y[],
                           double f[], Nag_User *comm);
  static void NAG_CALL pederv(Integer neq, double x, const double y[],
                              double pw[], Nag_User *comm);
  static double NAG_CALL g(Integer neq, double x, const double y[],
                           Nag_User *comm);
  static void NAG_CALL out(Integer neq, double *tsol, const double y[],
                           Nag_User *comm);
#ifdef __cplusplus
}
#endif

struct user
{
  double xend, h;
  Integer k;
  Integer *use_comm;
};

#define NEQ 3
int main(void)
{
  static Integer use_comm[4] = { 1, 1, 1, 1 };
  Integer exit_status = 0, j, neq;
  NagError fail;
  Nag_User comm;
  double tol, x, *y = 0;
  struct user s;

  INIT_FAIL(fail);

  printf("nag_ode_ivp_bdf_gen (d02ejc) Example Program Results\n");

  /* For communication with user-supplied functions
   * assign address of user defined structure
   * to comm.p.
   */
  s.use_comm = use_comm;
  comm.p = (Pointer) &s;

  neq = NEQ;
  if (neq >= 1) {
    if (!(y = NAG_ALLOC(neq, double)))
    {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }
  }
  else {
    exit_status = 1;
    return exit_status;
  }
  s.xend = 10.0;
  printf("\nCase 1: calculating Jacobian internally\n");
  printf(" intermediate output, root-finding\n\n");

  for (j = 3; j <= 4; ++j) {
    tol = pow(10.0, -(double) j);
    printf("\n Calculation with tol = %10.1e\n", tol);
    x = 0.0;
    y[0] = 1.0;
    y[1] = 0.0;
    y[2] = 0.0;
    s.k = 4;
    s.h = (s.xend - x) / (double) (s.k + 1);
    printf("     X         Y(1)         Y(2)         Y(3)\n");
    /* nag_ode_ivp_bdf_gen (d02ejc).
     * Ordinary differential equations solver, stiff, initial
     * value problems using the Backward Differentiation
     * Formulae
     */
    nag_ode_ivp_bdf_gen(neq, fcn, NULLFN, &x, y, s.xend, tol, Nag_Relative,
                        out, g, &comm, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_ode_ivp_bdf_gen (d02ejc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }
    printf("  Root of Y(1)-0.9 at %5.3f\n", x);
    printf("  Solution is ");
    printf("%7.4f %8.5f %7.4f\n", y[0], y[1], y[2]);
  }
  printf("\nCase 2: calculating Jacobian by pederv\n");
  printf(" intermediate output, root-finding\n\n");

  for (j = 3; j <= 4; ++j) {
    tol = pow(10.0, -(double) j);
    printf("\n Calculation with tol = %10.1e\n", tol);
    x = 0.0;
    y[0] = 1.0;
    y[1] = 0.0;
    y[2] = 0.0;
    s.k = 4;
    s.h = (s.xend - x) / (double) (s.k + 1);
    printf("     X         Y(1)         Y(2)         Y(3)\n");
    /* nag_ode_ivp_bdf_gen (d02ejc), see above. */
    nag_ode_ivp_bdf_gen(neq, fcn, pederv, &x, y, s.xend, tol, Nag_Relative,
                        out, g, &comm, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_ode_ivp_bdf_gen (d02ejc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }
    printf("  Root of Y(1)-0.9 at %5.3f\n", x);
    printf("  Solution is ");
    printf("%7.4f %8.5f %7.4f\n", y[0], y[1], y[2]);
  }
  printf("\nCase 3: calculating Jacobian internally\n");
  printf(" no intermediate output, root-finding\n\n");
  for (j = 3; j <= 4; ++j) {
    tol = pow(10.0, -(double) j);
    printf("\n Calculation with tol = %10.1e\n", tol);
    x = 0.0;
    y[0] = 1.0;
    y[1] = 0.0;
    y[2] = 0.0;

    /* nag_ode_ivp_bdf_gen (d02ejc), see above. */
    nag_ode_ivp_bdf_gen(neq, fcn, NULLFN, &x, y, s.xend, tol, Nag_Relative,
                        NULLFN, g, &comm, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_ode_ivp_bdf_gen (d02ejc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }

    printf("  Root of Y(1)-0.9 at %5.3f\n", x);
    printf("  Solution is ");
    printf("%7.4f %8.5f %7.4f\n", y[0], y[1], y[2]);
  }
  printf("\nCase 4: calculating Jacobian internally\n");
  printf(" intermediate output, no root-finding\n\n");

  for (j = 3; j <= 4; ++j) {
    tol = pow(10.0, -(double) j);
    printf("\n Calculation with tol = %10.1e\n", tol);
    x = 0.0;
    y[0] = 1.0;
    y[1] = 0.0;
    y[2] = 0.0;
    s.k = 4;
    s.h = (s.xend - x) / (double) (s.k + 1);
    printf("     X         Y(1)         Y(2)         Y(3)\n");
    /* nag_ode_ivp_bdf_gen (d02ejc), see above. */
    nag_ode_ivp_bdf_gen(neq, fcn, NULLFN, &x, y, s.xend, tol, Nag_Relative,
                        out, NULLDFN, &comm, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_ode_ivp_bdf_gen (d02ejc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }
    printf("%8.2f", x);
    printf("%13.4f%13.5f%13.4f\n", y[0], y[1], y[2]);
  }

  printf("\nCase 5: calculating Jacobian internally\n");
  printf(" no intermediate output, no root-finding (integrate to xend)\n\n");

  for (j = 3; j <= 4; ++j) {
    tol = pow(10.0, -(double) j);
    printf("\n Calculation with tol = %10.1e\n", tol);
    x = 0.0;
    y[0] = 1.0;
    y[1] = 0.0;
    y[2] = 0.0;
    printf("     X         Y(1)         Y(2)         Y(3)\n");
    printf("%8.2f", x);
    printf("%13.4f%13.5f%13.4f\n", y[0], y[1], y[2]);
    /* nag_ode_ivp_bdf_gen (d02ejc), see above. */
    nag_ode_ivp_bdf_gen(neq, fcn, NULLFN, &x, y, s.xend, tol, Nag_Relative,
                        NULLFN, NULLDFN, &comm, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_ode_ivp_bdf_gen (d02ejc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }
    printf("%8.2f", x);
    printf("%13.4f%13.5f%13.4f\n", y[0], y[1], y[2]);
  }
END:
  NAG_FREE(y);
  return exit_status;
}

static void NAG_CALL fcn(Integer neq, double x, const double y[], double f[],
                         Nag_User *comm)
{
  struct user *s = (struct user *) comm->p;

  if (s->use_comm[0]) {
    printf("(User-supplied callback fcn, first invocation.)\n");
    s->use_comm[0] = 0;
  }

  f[0] = y[0] * -0.04 + y[1] * 1e4 * y[2];
  f[1] = y[0] * 0.04 - y[1] * 1e4 * y[2] - y[1] * 3e7 * y[1];
  f[2] = y[1] * 3e7 * y[1];
}

static void NAG_CALL pederv(Integer neq, double x, const double y[],
                            double pw[], Nag_User *comm)
{
#define PW(I, J) pw[((I) -1)*neq + (J) -1]
  struct user *s = (struct user *) comm->p;

  if (s->use_comm[1]) {
    printf("(User-supplied callback pederv, first invocation.)\n");
    s->use_comm[1] = 0;
  }

  PW(1, 1) = -0.04;
  PW(1, 2) = y[2] * 1e4;
  PW(1, 3) = y[1] * 1e4;
  PW(2, 1) = 0.04;
  PW(2, 2) = y[2] * -1e4 - y[1] * 6e7;
  PW(2, 3) = y[1] * -1e4;
  PW(3, 1) = 0.0;
  PW(3, 2) = y[1] * 6e7;
  PW(3, 3) = 0.0;
}

static double NAG_CALL g(Integer neq, double x, const double y[],
                         Nag_User *comm)
{
  struct user *s = (struct user *) comm->p;

  if (s->use_comm[2]) {
    printf("(User-supplied callback g, first invocation.)\n");
    s->use_comm[2] = 0;
  }

  return y[0] - 0.9;
}

static void NAG_CALL out(Integer neq, double *xsol, const double y[],
                         Nag_User *comm)
{
  struct user *s = (struct user *) comm->p;

  printf("%8.2f", *xsol);
  printf("%13.4f%13.5f%13.4f\n", y[0], y[1], y[2]);

  *xsol = s->xend - (double) s->k * s->h;
  s->k--;
}