/* nag_ode_bvp_ps_lin_solve (d02uec) Example Program.
 *
 * Copyright 2011, Numerical Algorithms Group.
 *
 * Mark 23, 2011.
 */
#include <math.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagd02.h>
#include <nagx01.h>
#include <nagx02.h>

#ifdef __cplusplus
extern "C" {
#endif
  static double NAG_CALL exact(double x, Integer q);
  static void   NAG_CALL bndary(Integer m, double a, double b, double y[],
				double bmat[], double bvec[]);
  static void   NAG_CALL pdedef(Integer m, double f[]);
#ifdef __cplusplus
}
#endif

int main(void)
{
  /* Scalars */
  Integer  exit_status = 0;
  double   a = -nag_pi/2.0, b = nag_pi/2.0, resid;
  Integer  i, j, n, m = 3;
  double   teneps = 10.0 * nag_machine_precision;
  /* Arrays */
  double   *bmat = 0, *bvec = 0, *f = 0, *uerr = 0, *y = 0, *c = 0, *f0 = 0,
           *u = 0, *uc = 0, *x = 0;
  /* NAG types */
  Nag_Boolean reqerr = Nag_FALSE;
  NagError    fail;

  INIT_FAIL(fail);

  printf("nag_ode_bvp_ps_lin_solve (d02uec) Example Program Results\n\n");

  /* Skip heading in data file */
  scanf("%*[^\n] ");

  scanf("%"NAG_IFMT "%*[^\n] ", &n);
  if (
    !(u = NAG_ALLOC((n + 1)*(m + 1), double)) ||
    !(f0 = NAG_ALLOC((n + 1), double)) ||
    !(c = NAG_ALLOC((n + 1), double)) ||
    !(uc = NAG_ALLOC((n + 1)*(m + 1), double)) ||
    !(x = NAG_ALLOC((n + 1), double)) ||
    !(bmat = NAG_ALLOC(m*(m + 1), double)) ||
    !(bvec = NAG_ALLOC(m, double)) ||
    !(f = NAG_ALLOC((m + 1), double)) ||
    !(uerr = NAG_ALLOC((m + 1), double)) ||
    !(y = NAG_ALLOC(m, double))
    )
    {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }

  /* Set up domain, boundary conditions and definition. */
  bndary(m, a, b, y, bmat, bvec);
  pdedef(m, f);

  /* nag_ode_bvp_ps_lin_cgl_grid (d02ucc).
   * Generate Chebyshev Gauss-Lobatto solution grid.
   */
  nag_ode_bvp_ps_lin_cgl_grid(n, a, b, x, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_ode_bvp_ps_lin_cgl_grid (d02ucc).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }

  /* Set up problem right-hand sides for grid and transform. */
  for (i = 0; i < n + 1; i++) f0[i] = 2.0 * sin(x[i]) - 2.0 * cos(x[i]);

  /* nag_ode_bvp_ps_lin_coeffs (d02uac).
   * Coefficients of Chebyshev interpolating polynomial from function values f0
   * on Chebyshev grid.
   */
  nag_ode_bvp_ps_lin_coeffs(n, f0, c, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_ode_bvp_ps_lin_coeffs (d02uac).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }

  /* nag_ode_bvp_ps_lin_solve (d02uec).
   * Solve given boundary value problem on Chebyshev grid, in coefficient space
   * using an integral formulation of the pseudospectral method.
   */
  nag_ode_bvp_ps_lin_solve(n, a, b, m, c, bmat, y, bvec, f, uc, &resid, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_ode_bvp_ps_lin_solve (d02uec).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }

  /* nag_ode_bvp_ps_lin_cgl_vals (d02ubc).
   * Obtain function values from coefficients of Chebyshev polynomial.
   * Also obtain first- to third-derivative values.
   */
  for (i = 0; i < m + 1; i++) {
    nag_ode_bvp_ps_lin_cgl_vals(n, a, b, i, &uc[(n+1)*i], &u[(n+1)*i], &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_ode_bvp_ps_lin_cgl_vals (d02ubc).\n%s\n",
             fail.message);
      exit_status = 1;
      goto END;
    }
  }
  /* Print solution. */
  printf("Numerical Solution U and its first three derivatives\n\n");
  printf("%7s%12s%12s%11s%11s\n","x","U","Ux","Uxx","Uxxx");
  for (i = 0; i < n + 1; i++)
    printf("%10.4f %10.4f %10.4f %10.4f %10.4f\n", x[i], u[i], u[(n+1)+i],
	   u[(n+1)*2+i], u[(n+1)*3+i]);

  if (reqerr) {
    for (i = 0; i < m + 1; i++) uerr[i] = 0.0;

    for (i = 0; i < n + 1; i++)
      for (j = 0; j <= m; j++)
        uerr[j] = MAX(uerr[j], fabs(u[(n+1)*j+i] - exact(x[i], j)));

    for (i = 0; i <= m; i++) {
      printf("Error in the order %1"NAG_IFMT " derivative of U is < %8"
	     NAG_IFMT " * machine precision.\n", i,
	     10 * ((Integer) (uerr[i]/teneps) + 1));
    }
  }
 END:
  if (c) NAG_FREE(c);
  if (f0) NAG_FREE(f0);
  if (u) NAG_FREE(u);
  if (uc) NAG_FREE(uc);
  if (x) NAG_FREE(x);
  if (bmat) NAG_FREE(bmat);
  if (bvec) NAG_FREE(bvec);
  if (f) NAG_FREE(f);
  if (uerr) NAG_FREE(uerr);
  if (y) NAG_FREE(y);
  return exit_status;
}

static double exact(double x, Integer q)
{
  switch (q) {
    case 0:
      return cos(x);
      break;
    case 1:
      return -sin(x);
      break;
    case 2:
      return -cos(x);
      break;
    case 3:
      return sin(x);
      break;
  }
  return 0.0;
}
static void bndary(Integer m, double a, double b, double y[], double bmat[],
                   double bvec[])
{
  Integer i;
  /* Boundary condition on left side of domain. */
  for (i = 0; i < 2; i++) y[i] = a;

  y[2] = b;
  /* Set up Dirichlet condition using exact solution at x = a. */
  for (i = 0; i < m*(m+1); i++) bmat[i] = 0.0;
  for (i = 0; i < 3; i++) bmat[i] = 1.0;
  for (i = 1; i < 3; i++) bmat[m+i] = 2.0;
  for (i = 1; i < 3; i++) bmat[m*2+i] = 3.0;

  bvec[0] = 0.0;
  bvec[1] = 2.0;
  bvec[2] = -2.0;
}
static void pdedef(Integer m, double f[])
{
  f[0] = 1.0;
  f[1] = 2.0;
  f[2] = 3.0;
  f[3] = 4.0;
}