/* nag_dae_ivp_dassl_gen (d02nec) Example Program.
 *
 * NAGPRODCODE Version.
 *
 * Copyright 2016 Numerical Algorithms Group.
 *
 * Mark 26, 2016.
 *
 */

/* Pre-processor includes */
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagd02.h>

#ifdef __cplusplus
extern "C"
{
#endif
  static void NAG_CALL res(Integer neq, double t, const double y[],
                           const double ydot[], double r[], Integer *ires,
                           Nag_Comm *comm);
  static void NAG_CALL jac(Integer neq, double t, const double y[],
                           const double ydot[], double *pd, double cj,
                           Nag_Comm *comm);
  static void NAG_CALL myjac(Integer neq, Integer ml, Integer mu, double t,
                             const double y[], const double ydot[],
                             double *pd, double cj);
#ifdef __cplusplus
}
#endif
int main(void)
{
  /*Integer scalar and array declarations */
  Integer exit_status = 0, maxord = 5;
  Nag_Comm comm;
  Integer neq, licom, mu, ml, lcom;
  Integer i, itask, j;
  Nag_Boolean vector_tol;
  Integer *icom = 0;
  NagError fail;
  /*Double scalar and array declarations */
  double dt, h0, hmax, t, tout;
  double *atol = 0, *com = 0, *rtol = 0, *y = 0, *ydot = 0;
  static double ruser[2] = { -1.0, -1.0 };

  INIT_FAIL(fail);

  printf("nag_dae_ivp_dassl_gen (d02nec) Example Program Results\n\n");

  /* For communication with user-supplied functions: */
  comm.user = ruser;

  /* Set problem parameters required to allocate arrays */
  neq = 3;
  ml = 1;
  mu = 2;
  licom = 50 + neq;
  lcom = 40 + (maxord + 4) * neq + (2 * ml + mu + 1) * neq +
         2 * (neq / (ml + mu + 1) + 1);
  if (!(atol = NAG_ALLOC(neq, double)) || !(com = NAG_ALLOC(lcom, double))
      || !(rtol = NAG_ALLOC(neq, double)) || !(y = NAG_ALLOC(neq, double))
      || !(ydot = NAG_ALLOC(neq, double))
      || !(comm.iuser = NAG_ALLOC(2, Integer))
      || !(icom = NAG_ALLOC(licom, Integer)))
  {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }
  /* Initialize the problem, specifying that the Jacobian is to be */
  /* evaluated analytically using the provided routine jac.        */
  h0 = 0.0;
  hmax = 0.0;
  vector_tol = Nag_TRUE;
  /*
   * nag_dae_ivp_dassl_setup (d02mwc)
   * Implicit DAE/ODEs, stiff ivp, setup for nag_dae_ivp_dassl_gen (d02nec)
   */
  nag_dae_ivp_dassl_setup(neq, maxord, Nag_AnalyticalJacobian, hmax, h0,
                          vector_tol, icom, licom, com, lcom, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_dae_ivp_dassl_setup (d02mwc).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }
  /* Specify that the Jacobian is banded.
   *
   * nag_dae_ivp_dassl_linalg (d02npc)
   * ODE/DAEs, ivp, linear algebra setup routine for
   * nag_dae_ivp_dassl_gen (d02nec)
   */
  nag_dae_ivp_dassl_linalg(neq, ml, mu, icom, licom, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_dae_ivp_dassl_linalg (d02npc).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }

  /* Set initial values */
  t = 0.00e0;
  tout = 0.00e0;
  dt = 0.020e0;
  for (i = 0; i < neq; i++) {
    rtol[i] = 1.00e-3;
    atol[i] = 1.00e-6;
    y[i] = 0.00e0;
    ydot[i] = 0.00e0;
  }
  y[0] = 1.00e0;
  /* Use the comm.iuser array to pass the band dimensions through to jac. */
  /* An alternative would be to hard code values for ml and mu in jac.    */
  comm.iuser[0] = ml;
  comm.iuser[1] = mu;
  printf("    t          y(1)         y(2)         y(3)\n");
  printf("%8.4f", t);
  for (i = 0; i < neq; i++)
    printf("%12.6f%s", y[i], (i + 1) % 3 ? " " : "\n");
  itask = 0;
  /*     Obtain the solution at 5 equally spaced values of t. */
  for (j = 0; j < 5; j++) {
    tout = tout + dt;
    /*
     * nag_dae_ivp_dassl_gen (d02nec)
     * dassl integrator
     */
    nag_dae_ivp_dassl_gen(neq, &t, tout, y, ydot, rtol, atol, &itask, res,
                          jac, icom, com, lcom, &comm, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_dae_ivp_dassl_gen (d02nec).\n%s\n",
             fail.message);
      exit_status = 1;
      goto END;
    }
    printf("%8.4f", t);
    for (i = 0; i < neq; i++)
      printf("%12.6f%s", y[i], (i + 1) % 3 ? " " : "\n");
    /*
     * nag_dae_ivp_dassl_cont (d02mcc)
     * dassl method continuation resetting function
     */
    nag_dae_ivp_dassl_cont(icom);
  }
  printf("\n");
  printf("  The integrator completed task, ITASK = %4" NAG_IFMT "\n", itask);

END:
  NAG_FREE(atol);
  NAG_FREE(com);
  NAG_FREE(rtol);
  NAG_FREE(y);
  NAG_FREE(ydot);
  NAG_FREE(comm.iuser);
  NAG_FREE(icom);

  return exit_status;
}

static void NAG_CALL res(Integer neq, double t, const double y[],
                         const double ydot[], double r[], Integer *ires,
                         Nag_Comm *comm)
{
  if (comm->user[0] == -1.0) {
    printf("(User-supplied callback res, first invocation.)\n");
    comm->user[0] = 0.0;
  }
  r[0] = (-(0.040e0 * y[0])) + 1.00e4 * y[1] * y[2] - ydot[0];
  r[1] = 0.040e0 * y[0] - 1.00e4 * y[1] * y[2] - 3.00e7 * y[1] * y[1] -
         ydot[1];
  r[2] = 3.00e7 * y[1] * y[1] - ydot[2];
  return;
}

static void NAG_CALL jac(Integer neq, double t, const double y[],
                         const double ydot[], double *pd, double cj,
                         Nag_Comm *comm)
{
  Integer ml, mu;
  if (comm->user[1] == -1.0) {
    printf("(User-supplied callback jac, first invocation.)\n");
    comm->user[1] = 0.0;
  }
  ml = comm->iuser[0];
  mu = comm->iuser[1];
  myjac(neq, ml, mu, t, y, ydot, pd, cj);
  return;
}

static void NAG_CALL myjac(Integer neq, Integer ml, Integer mu, double t,
                           const double y[], const double ydot[], double *pd,
                           double cj)
{
  Integer md, ms;
  Integer pdpd;
  pdpd = 2 * ml + mu + 1;
#define PD(I, J) pd[(J-1)*pdpd + I - 1]
  /*     Main diagonal PDFULL(i,i), i=1,neq */
  md = mu + ml + 1;
  PD(md, 1) = -0.040e0 - cj;
  PD(md, 2) = -1.00e4 * y[2] - 6.00e7 * y[1] - cj;
  PD(md, 3) = -cj;
  /*     1 Subdiagonal PDFULL(i+1:i), i=1,neq-1 */
  ms = md + ml;
  PD(ms, 1) = 0.040e0;
  PD(ms, 2) = 6.00e7 * y[1];
  /*     First superdiagonal PDFULL(i-1,i), i=2, neq */
  ms = md - 1;
  PD(ms, 2) = 1.00e4 * y[2];
  PD(ms, 3) = -1.00e4 * y[1];
  /*     Second superdiagonal PDFULL(i-2,i), i=3, neq */
  ms = md - 2;
  PD(ms, 3) = 1.00e4 * y[1];
  return;
}