/* nag_dae_ivp_dassl_setup (d02mwc) Example Program.
 *
 * Copyright 2009, Numerical Algorithms Group.
 *
 * Mark 9, 2009.
 *
 */

/* Pre-processor includes */
#include <stdio.h>
#include <math.h>
#include <nag.h>
#include <nagx04.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);

#ifdef __cplusplus
}
#endif
int main(int argc, char *argv[])
{
  FILE        *fpout;
  /*Integer scalar and array declarations */
  Integer     exit_status = 0;
  Integer     neq, maxord, licom, lcom;
  Integer     i, itask;
  Integer     *icom = 0;
  Nag_Boolean vector_tol;
  Nag_Comm    comm;
  /*Double scalar and array declarations */
  double      h0, hmax, g1, g2, t, tout;
  double      *atol = 0, *com = 0, *rtol = 0, *y = 0, *ydot = 0;
  NagError    fail;

  INIT_FAIL(fail);

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

  fprintf(fpout,
          "nag_dae_ivp_dassl_setup (d02mwc) Example Program Results\n\n");
  t = 0.00e0;
  tout = 1.00e0;
  vector_tol = Nag_TRUE;
  maxord = 5;
  neq = 5;
  licom = 50+neq;
  lcom = 40+(maxord+4)*neq+(int) pow(neq, 2);
  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)) ||
      !(icom = NAG_ALLOC(licom, Integer)))
    {
      fprintf(fpout, "Allocation failure\n");
      exit_status = -1;
      goto END;
    }
  for (i = 0; i < neq; i++)
    {
      rtol[i] = 1.00e-8;
      atol[i] = 1.00e-8;
    }
  /* Set initial values*/
  y[0] = 1.00e0;
  y[1] = 0.00e0;
  y[2] = 0.00e0;
  y[3] = 1.00e0;
  y[4] = 1.00e0;
  h0 = 0.0;
  hmax = 0.0;
  /*
   * nag_dae_ivp_dassl_setup (d02mwc)
   * Implicit DAE/ODEs, stiff IVP, setup for d02nec
   */
  nag_dae_ivp_dassl_setup(neq, maxord, Nag_AnalyticalJacobian, hmax, h0,
                          vector_tol,
                          icom, licom, com, lcom,
                          &fail);
  if (fail.code != NE_NOERROR)
    {
      fprintf(fpout, "Error from nag_dae_ivp_dassl_setup (d02mwc).\n%s\n",
              fail.message);
      exit_status = 1;
      goto END;
    }
  fprintf(fpout, "%7s%12s%12s%12s%12s%12s\n", "t", "y(1)", "y(2)", "y(3)",
          "y(4)", "y(5)");
  fprintf(fpout, "   %6.4f", t);
  for (i = 0; i < neq; i++)
    fprintf(fpout, "%11.6f%s", y[i], (i+1)%5?" ":"\n");

  for (i = 0; i < neq; i++)
    ydot[i] = 0.00e0;
  itask = 0;
 L20:;
  /*
   * 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)
    {
      fprintf(fpout, "Error from nag_dae_ivp_dassl_gen (d02nec).\n%s\n",
              fail.message);
      exit_status = 1;
      goto END;
    }
  fprintf(fpout, "   %6.4f", t);
  for (i = 0; i < neq; i++)
    fprintf(fpout, "%11.6f%s", y[i], (i+1)%5?" ":"\n");
  fprintf(fpout, "\n");
  fprintf(fpout, "  d02nec  returned with ITASK = %4ld\n\n", itask);
  if ((itask >= 0) && (itask <= 3))
    {
      if (t < tout)
        goto L20;
      g1 = pow(y[0], 2)+pow(y[1], 2.0)-1;
      g2 = y[0]*y[2]+y[1]*y[3];
      fprintf(fpout, "  The position-level constraint G1 =  %13.4e\n", g1);
      fprintf(fpout, "  The velocity-level constraint G2 =  %13.4e\n", g2);
    }

 END:
  if (fpout != stdout) fclose(fpout);
  if (atol) NAG_FREE(atol);
  if (com) NAG_FREE(com);
  if (rtol) NAG_FREE(rtol);
  if (y) NAG_FREE(y);
  if (ydot) NAG_FREE(ydot);
  if (icom) 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)
{
  r[0] = y[2]-ydot[0];
  r[1] = y[3]-ydot[1];
  r[2] = (-(y[4]*y[0]))-ydot[2];
  r[3] = (-(y[4]*y[1]))-1.0e0-ydot[3];
  r[4] = pow(y[2], 2)+pow(y[3], 2)-y[4]-y[1];
  return;
}
static void NAG_CALL jac(Integer neq, double t, const double y[],
                         const double ydot[], double *pd, double cj,
                         Nag_Comm *comm)
{
  Integer pdpd;
  pdpd = neq;
#define PD(I, J) pd[(J-1)*pdpd + I-1]
  PD(1, 1) = (-(cj));
  PD(1, 3) = 1.00e0;
  PD(2, 2) = (-(cj));
  PD(2, 4) = 1.00e0;
  PD(3, 1) = (-(y[4]));
  PD(3, 3) = (-(cj));
  PD(3, 5) = (-(y[0]));
  PD(4, 2) = (-(y[4]));
  PD(4, 4) = (-(cj));
  PD(4, 5) = (-(y[1]));
  PD(5, 2) = (-(1.00e0));
  PD(5, 3) = 2.00e0*y[2];
  PD(5, 4) = 2.00e0*y[3];
  PD(5, 5) = (-(1.00e0));
  return;
}