/* nag_ode_bvp_ps_lin_solve (d02uec) Example Program.
 *
 * Copyright 2023 Numerical Algorithms Group.
 *
 * Mark 29.2, 2023.
 */
#include <math.h>
#include <nag.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_math_pi / 2.0, b = nag_math_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:
  NAG_FREE(c);
  NAG_FREE(f0);
  NAG_FREE(u);
  NAG_FREE(uc);
  NAG_FREE(x);
  NAG_FREE(bmat);
  NAG_FREE(bvec);
  NAG_FREE(f);
  NAG_FREE(uerr);
  NAG_FREE(y);
  return exit_status;
}

static double NAG_CALL 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 NAG_CALL 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 NAG_CALL pdedef(Integer m, double f[]) {
  f[0] = 1.0;
  f[1] = 2.0;
  f[2] = 3.0;
  f[3] = 4.0;
}