/* nag_ode_bvp_coll_nlin_interp (d02tyc) Example Program.
 *
 * Copyright 2014 Numerical Algorithms Group.
 * 
 * Mark 24, 2013.
 */

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

#define Y(I, J)   y[J*neq + I-1]

#ifdef __cplusplus
extern "C" {
#endif
static void NAG_CALL ffun(double x, const double y[], Integer neq, 
                          const Integer m[], double f[], Nag_Comm *comm);
static void NAG_CALL fjac(double x, const double y[], Integer neq, 
                          const Integer m[], double dfdy[], Nag_Comm *comm);
static void NAG_CALL gafun(const double ya[], Integer neq, const Integer m[],
                           Integer nlbc, double ga[], Nag_Comm *comm);
static void NAG_CALL gbfun(const double yb[], Integer neq, const Integer m[],
                           Integer nrbc, double gb[], Nag_Comm *comm);
static void NAG_CALL gajac(const double ya[], Integer neq, const Integer m[],
                           Integer nlbc, double dgady[], Nag_Comm *comm);
static void NAG_CALL gbjac(const double yb[], Integer neq, const Integer m[],
                           Integer nrbc, double dgbdy[], Nag_Comm *comm);
static void NAG_CALL guess(double x, Integer neq, const Integer m[], double y[],
                           double dym[], Nag_Comm *comm);
#ifdef __cplusplus
}
#endif

int main(void)
{
  /* Scalars */
  Integer     exit_status = 0, nmesh_out = 11;
  Integer     neq, mmax, nlbc, nrbc;
  Integer     i, iermx, ijermx, licomm, lrcomm, mxmesh, ncol, nmesh;
  double      a, ainc, ermx, x;
  /* Arrays */
  static Integer iuser[7] = {-1, -1, -1, -1, -1, -1, -1};
  double      *mesh = 0, *rcomm = 0, *tol = 0, *y = 0;
  double      rdum[1], ruser[1];
  Integer     *ipmesh = 0, *icomm = 0, *m = 0;
  Integer     idum[2];
  /* Nag Types */
  Nag_Boolean failed = Nag_FALSE;
  Nag_Comm    comm;
  NagError    fail;

  INIT_FAIL(fail);

  printf("nag_ode_bvp_coll_nlin_interp (d02tyc) Example Program Results\n\n");
  /* Skip heading in data file*/
  scanf("%*[^\n] ");
  scanf("%"NAG_IFMT "%"NAG_IFMT "", &neq, &mmax);
  scanf("%"NAG_IFMT "%"NAG_IFMT "%*[^\n] ", &nlbc, &nrbc);
  scanf("%"NAG_IFMT "%"NAG_IFMT "%"NAG_IFMT "%*[^\n] ", &ncol, &nmesh, &mxmesh);
  if (!(mesh   = NAG_ALLOC(mxmesh, double)) ||
      !(m      = NAG_ALLOC(neq, Integer)) ||
      !(tol    = NAG_ALLOC(neq, double)) ||
      !(y      = NAG_ALLOC(neq*mmax, double)) ||
      !(ipmesh = NAG_ALLOC(mxmesh, Integer)))
    {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }

  comm.user = ruser;
  comm.iuser = iuser;

  for (i = 0; i < neq; i++) {
    scanf("%"NAG_IFMT "", &m[i]);
  } 
  scanf("%*[^\n] ");
  scanf("%lf%*[^\n] ", &a);

  for (i = 0; i < neq; i++) {
    scanf("%lf", &tol[i]);
  }
  scanf("%*[^\n] ");

  ainc = a/(double) (nmesh - 1);
  mesh[0] = 0.0;
  for (i = 1; i < nmesh - 1; i++) {
    mesh[i] = mesh[i - 1] + ainc;
  }
  mesh[nmesh-1] = a;

  ipmesh[0] = 1;
  for (i = 1; i < nmesh - 1; i++) {
    ipmesh[i] = 2;
  }
  ipmesh[nmesh-1] = 1;

  /* For communication with function guess store a in comm.user[0]. */
  ruser[0] = a;

  /* Communication space query to get size of rcomm and icomm 
   * by setting lrcomm=0 in call to
   * nag_ode_bvp_coll_nlin_setup (d02tvc):
   * Ordinary differential equations, general nonlinear boundary value problem,
   * setup for nag_ode_bvp_coll_nlin_solve (d02tlc).
   */
  nag_ode_bvp_coll_nlin_setup(neq, m, nlbc, nrbc, ncol, tol, mxmesh, nmesh,
                              mesh, ipmesh, rdum, 0, idum, 2, &fail);
  if (fail.code == NE_NOERROR) {
    lrcomm = idum[0];
    licomm = idum[1];

    if (!(rcomm = NAG_ALLOC(lrcomm, double)) ||
        !(icomm = NAG_ALLOC(licomm, Integer))) {
      printf("Allocation failure\n");
      exit_status = -2;
      goto END;
    }

    /* Initialize, again using nag_ode_bvp_coll_nlin_setup (d02tvc). */
    nag_ode_bvp_coll_nlin_setup(neq, m, nlbc, nrbc, ncol, tol, mxmesh, nmesh,
                                mesh, ipmesh, rcomm, lrcomm, icomm, licomm,
                                &fail);
  }
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_ode_bvp_coll_nlin_setup (d02tvc).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }

  printf(" Tolerance = %8.1e, a = %8.2f\n", tol[0], a);
  /* Solve*/
  /* nag_ode_bvp_coll_nlin_solve (d02tlc).
   * Ordinary differential equations, general nonlinear boundary value 
   * problem, collocation technique.
   */
  nag_ode_bvp_coll_nlin_solve(ffun, fjac, gafun, gbfun, gajac, gbjac, guess,
                              rcomm, icomm, &comm, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_ode_bvp_coll_nlin_solve (d02tlc).\n%s\n",
           fail.message);
    failed = Nag_TRUE;
  }

  /* Extract mesh.*/
  /* nag_ode_bvp_coll_nlin_diag (d02tzc).
   * Ordinary differential equations, general nonlinear boundary value
   * problem, diagnostics for nag_ode_bvp_coll_nlin_solve (d02tlc).
   */
  nag_ode_bvp_coll_nlin_diag(mxmesh, &nmesh, mesh, ipmesh, &ermx, &iermx,
                             &ijermx, rcomm, icomm, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_ode_bvp_coll_nlin_diag (d02tzc).\n%s\n",
           fail.message);
    exit_status = 2;
    goto END;
  }

  /* Print mesh statistics*/
  printf("\n Used a mesh of %4"NAG_IFMT "  points\n", nmesh);

  printf(" Maximum error = %10.2e  in interval %4ld", ermx, iermx);
  printf("  for component %4"NAG_IFMT " \n\n\n", ijermx);
  printf(" Mesh points:\n");
  for (i = 0; i < nmesh; i++) {
    printf("%4ld(%1ld)%13.4e%s", i+1, ipmesh[i], mesh[i],
           (i+1)%4?"":"\n");
  }
  printf("\n");
  
  if (!failed) {
    /* Print solution on output mesh.*/
    printf("\n\n Computed solution\n     x       solution   derivative\n");
    x = 0.0;
    ainc = a/(double) (nmesh_out - 1);
    for (i = 0; i < nmesh_out; i++) {

      /* nag_ode_bvp_coll_nlin_interp (d02tyc).
       * Ordinary differential equations, general nonlinear boundary value 
       * problem, interpolation for nag_ode_bvp_coll_nlin_solve (d02tlc).
       */
      nag_ode_bvp_coll_nlin_interp(x, y, neq, mmax, rcomm, icomm, &fail);
      if (fail.code != NE_NOERROR) {
        printf("Error from nag_ode_bvp_coll_nlin_interp (d02tyc).\n%s\n",
               fail.message);
        exit_status = 3;
        goto END;
      }

      printf("%8.2f %11.5f %11.5f\n", x, y[0], y[neq]);
      x = x + ainc;
    }
  }

 END:
  NAG_FREE(mesh);
  NAG_FREE(m);
  NAG_FREE(tol);
  NAG_FREE(rcomm);
  NAG_FREE(y);
  NAG_FREE(ipmesh);
  NAG_FREE(icomm);
  return exit_status;
}

static void NAG_CALL ffun(double x, const double y[], Integer neq, 
                          const Integer m[], double f[], Nag_Comm *comm)
{
  if (comm->iuser[0] == -1)
    {
      printf("(User-supplied callback ffun, first invocation.)\n");
      comm->iuser[0] = 0;
    }
  if (y[0] <= 0.0) {
    f[0] = 0.0;
  }
  else {
    f[0] = pow(y[0], 1.5)/sqrt(x);
  }
}

static void NAG_CALL fjac(double x, const double y[], Integer neq, 
                          const Integer m[], double dfdy[], Nag_Comm *comm)
{
  if (comm->iuser[1] == -1)
    {
      printf("(User-supplied callback fjac, first invocation.)\n");
      comm->iuser[1] = 0;
    }
  if (y[0] <= 0.0) {
    dfdy[0] = 0.0;
  }
  else {
    dfdy[0] = 1.5 * sqrt(y[0])/sqrt(x);
  }
}

static void NAG_CALL gafun(const double ya[], Integer neq, const Integer m[],
                           Integer nlbc, double ga[], Nag_Comm *comm)
{
  if (comm->iuser[2] == -1)
    {
      printf("(User-supplied callback gafun, first invocation.)\n");
      comm->iuser[2] = 0;
    }
  ga[0] = ya[0] - 1.0;
}

static void NAG_CALL gbfun(const double yb[], Integer neq, const Integer m[],
                           Integer nrbc, double gb[], Nag_Comm *comm)
{
  if (comm->iuser[3] == -1)
    {
      printf("(User-supplied callback gbfun, first invocation.)\n");
      comm->iuser[3] = 0;
    }
  gb[0] = yb[0];
}

static void NAG_CALL gajac(const double ya[], Integer neq, const Integer m[],
                           Integer nlbc, double dgady[], Nag_Comm *comm)
{
  if (comm->iuser[4] == -1)
    {
      printf("(User-supplied callback gajac, first invocation.)\n");
      comm->iuser[4] = 0;
    }
  dgady[0] = 1.0;
}

static void NAG_CALL gbjac(const double yb[], Integer neq, const Integer m[],
                           Integer nrbc, double dgbdy[], Nag_Comm *comm)
{
  if (comm->iuser[5] == -1)
    {
      printf("(User-supplied callback gbjac, first invocation.)\n");
      comm->iuser[5] = 0;
    }
  dgbdy[0] = 1.0;
}

static void NAG_CALL guess(double x, Integer neq, const Integer m[], double y[],
                           double dym[], Nag_Comm *comm)
{
  double a = comm->user[0];
  if (comm->iuser[6] == -1)
    {
      printf("(User-supplied callback guess, first invocation.)\n");
      comm->iuser[6] = 0;
    }
  y[0] = 1.0 - x/(a);
  y[neq] = -1.0/(a);
  dym[0] = 0.0;
}