/* nag_ode_ivp_bdf_gen (d02ejc) Example Program.
 *
 * Copyright 1992 Numerical Algorithms Group.
 *
 * Mark 3, 1992.
 * Mark 7 revised, 2001.
 * Mark 8 revised, 2004.
 *
 */

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


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

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

#define NEQ 3
int main(void)
{
  Integer exit_status=0, i, j, neq;
  NagError fail;
  Nag_User comm;
  double tol, x, *y=0;
  struct user s;

  INIT_FAIL(fail);
  Vprintf("nag_ode_ivp_bdf_gen (d02ejc) Example Program Results\n");

  /* For communication with function out()
   * assign address of user defined structure
   * to comm.p.
   */
  comm.p = (Pointer)&s;

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

  for (j=3; j<=4; ++j)
    {
      tol = pow(10.0, -(double)j);
      Vprintf("\n Calculation with tol = %3.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);
      Vprintf( "     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)
        {
          Vprintf("Error from nag_ode_ivp_bdf_gen (d02ejc).\n%s\n",
                  fail.message);
          exit_status = 1;
          goto END;
        }
      Vprintf("  Root of Y(1)-0.9 at %5.3f\n", x);
      Vprintf("  Solution is ");
      for (i=0; i<3; ++i)
        Vprintf("%7.5f ", y[i]);
      Vprintf("\n");
    }
  Vprintf("\nCase 2: calculating Jacobian by pederv\n");
  Vprintf(" intermediate output, root-finding\n\n");

  for (j=3; j<=4; ++j)
    {
      tol = pow(10.0, -(double)j);
      Vprintf("\n Calculation with tol = %3.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);
      Vprintf( "     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)
        {
          Vprintf("Error from nag_ode_ivp_bdf_gen (d02ejc).\n%s\n",
                  fail.message);
          exit_status = 1;
          goto END;
        }
      Vprintf("  Root of Y(1)-0.9 at %5.3f\n", x);
      Vprintf("  Solution is ");
      for (i=0; i<3; ++i)
        Vprintf("%7.5f ", y[i]);
      Vprintf("\n");
    }
  Vprintf("\nCase 3: calculating Jacobian internally\n");
  Vprintf(" no intermediate output, root-finding\n\n");
  for (j=3; j<=4; ++j)
    {
      tol = pow(10.0, -(double)j);
      Vprintf("\n Calculation with tol = %3.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)
        {
          Vprintf("Error from nag_ode_ivp_bdf_gen (d02ejc).\n%s\n",
                  fail.message);
          exit_status = 1;
          goto END;
        }

      Vprintf("  Root of Y(1)-0.9 at %5.3f\n", x);
      Vprintf("  Solution is ");
      for (i=0; i<3; ++i)
        Vprintf("%7.5f ", y[i]);
      Vprintf("\n");
    }
  Vprintf("\nCase 4: calculating Jacobian internally\n");
  Vprintf(" intermediate output, no root-finding\n\n");

  for (j=3; j<=4; ++j)
    {
      tol = pow(10.0, -(double)j);
      Vprintf("\n Calculation with tol = %3.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);
      Vprintf( "     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)
        {
          Vprintf("Error from nag_ode_ivp_bdf_gen (d02ejc).\n%s\n",
                  fail.message);
          exit_status = 1;
          goto END;
        }
      Vprintf("%8.2f", x);
      for (i=0; i<3; ++i)
        Vprintf("%13.5f", y[i]);
      Vprintf("\n");
    }

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

  for (j=3; j<=4; ++j)
    {
      tol = pow(10.0, -(double)j);
      Vprintf("\n Calculation with tol = %3.1e\n", tol);
      x = 0.0;
      y[0] = 1.0;
      y[1] = 0.0;
      y[2] = 0.0;
      Vprintf( "     X         Y(1)         Y(2)         Y(3)\n");
      Vprintf("%8.2f", x);
      for (i=0; i<3; ++i)
        Vprintf("%13.5f", y[i]);
      Vprintf("\n");
      /* 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)
        {
          Vprintf("Error from nag_ode_ivp_bdf_gen (d02ejc).\n%s\n",
                  fail.message);
          exit_status = 1;
          goto END;
        }
      Vprintf("%8.2f", x);
      for (i=0; i<3; ++i)
        Vprintf("%13.5f", y[i]);
      Vprintf("\n");
    }
 END:
  if (y) NAG_FREE(y);
  return exit_status;
}
static void fcn(Integer neq, double x, double y[], double f[], Nag_User *comm)
{
  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 pederv(Integer neq, double x, double y[], double pw[],
                   Nag_User *comm)
{

#define PW(I,J) pw[((I)-1)*neq + (J)-1]

  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 g(Integer neq, double x, double y[], Nag_User *comm)
{
  return y[0]-0.9;
}
static void out(Integer neq, double *xsol, double y[], Nag_User *comm)
{
  Integer j;
  struct user *s = (struct user *)comm->p;

  Vprintf("%8.2f", *xsol);
  for (j=0; j<3; ++j)
    Vprintf("%13.5f", y[j]);
  Vprintf("\n");

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