/* nag_numdiff_1d_real_eval (d04bac) Example Program.
 *
 * Copyright 2011, Numerical Algorithms Group.
 *
 * Mark 23, 2011.
 */

#include <stdio.h>
#include <string.h>
#include <math.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagd04.h>
#include <nags.h>
#include <nagx04.h>

int main(void)
{
/* DER_COMP is used to store results in column major format only. */
#define DER_COMP(I, J, K) der_comp[I + J*n_hbase + K*n_hbase*n_der_comp]

  Integer exit_status = 0;
  double  x_0;
  Integer n_der_comp, n_display, n_hbase;
  double  hbase;
  Integer i, j, k;
  char    title[51];
  double  der[14], erest[14], fval[21], xval[21];
  double  *actder=0, *der_comp=0;
  char    *clabs=0, **clabsc=0;
  NagError fail;

  INIT_FAIL(fail);

  printf("nag_numdiff_1d_real_eval (d04bac) Example Program Results\n");

  printf("\n"
         "Find the derivatives of the polygamma (psi) function\n"
         "using function values generated by nag_real_polygamma (s14aec).\n\n"
         "Demonstrate the effect of successively reducing hbase.\n\n");

  n_der_comp = 3;
  n_display = 3;
  n_hbase = 4;

  if ( !(clabs = NAG_ALLOC(n_der_comp*10, char)) ||
       !(clabsc = NAG_ALLOC(n_der_comp, char *)) ||
       !(actder = NAG_ALLOC(n_display, double)) ||
       !(der_comp = NAG_ALLOC(n_hbase*n_der_comp*14, double)))
    {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }

  x_0 = 0.05;

  /* Select an initial separation distance hbase. */
  hbase = 0.0025;

  /* Compute the actual derivatives at target location x_0
   * using nag_real_polygamma (s14aec), for comparison.
   */
  for (j = 0; j < n_display; j++)
    {
      /* nag_real_polygamma (s14aec).
       * Derivative of the psi function psi(x).
       */
      actder[j] = nag_real_polygamma(x_0, j+1, &fail);
      if (fail.code != NE_NOERROR)
        {
          printf("Error from nag_real_polygamma (s14aec).\n%s\n",
                 fail.message);
          exit_status = 1;
          goto END;
        }
    }

  /* Attempt n_hbase approximations, reducing hbase by factor 0.1 each time. */
  for (j = 0; j < n_hbase; j++)
    {
      /* nag_numdiff_1d_real_absci (d04bbc).
       * Generates sample points for function evaluations
       * by nag_numdiff_1d_real_eval (d04bac).
       */
      nag_numdiff_1d_real_absci(x_0, hbase, xval);

      /* Calculate the corresponding objective function values. */
      for (k = 0; k < 21; k++)
        {
          fval[k] = nag_real_polygamma(xval[k], 0, &fail);
          if (fail.code != NE_NOERROR)
            {
              printf("Error from nag_real_polygamma (s14aec).\n%s\n",
                     fail.message);
              exit_status = 1;
              goto END;
            }
        }

      /* Calculate the derivative estimates. */

      /* nag_numdiff_1d_real_eval (d04bac).
       * Numerical differentiation, user-supplied function values,
       * derivatives up to order 14,
       * derivatives with respect to one real variable.
       */
      nag_numdiff_1d_real_eval(xval, fval, der, erest, &fail);
      if (fail.code != NE_NOERROR)
        {
          printf("Error from nag_numdiff_1d_real_eval (d04bac).\n%s\n",
                 fail.message);
          exit_status = 1;
          goto END;
        }

      /* Store results in DER_COMP in column major format. */
      for (i = 0; i < 14; i++)
        {
          DER_COMP(j, 0, i) = hbase;
          DER_COMP(j, 1, i) = der[i];
          DER_COMP(j, 2, i) = erest[i];
        }

      /* Decrease hbase for next loop. */
      hbase = hbase * 0.1;
    }

  /* Display results for first n_display derivatives. */
  for (j = 0; j < n_display; j++)
    {
      printf(" Derivative (%ld) calculated using "
             "nag_real_polygamma (s14aec): %11.4e\n", j+1, actder[j]);

      sprintf(title,
              "Derivative and error estimates for derivative (%ld)",
              j+1);
      sprintf(&clabs[0], "hbase");
      sprintf(&clabs[10], "der[%ld]", j);
      sprintf(&clabs[20], "erest[%ld]", j);
      for(k = 0; k < n_der_comp; k++) clabsc[k] = &clabs[k*10];

      /* nag_gen_real_mat_print_comp (x04cbc).
       * Print real general matrix (comprehensive).
       */
      nag_gen_real_mat_print_comp(Nag_ColMajor, Nag_GeneralMatrix,
                                  Nag_NonUnitDiag, n_hbase, n_der_comp,
                                  &DER_COMP(0, 0, j), n_hbase, NULL, title,
                                  Nag_NoLabels, NULL, Nag_CharacterLabels,
                                  (const char**)clabsc, 0, 0, NULL, &fail);
      if (fail.code != NE_NOERROR)
        {
          printf("Error from nag_gen_real_mat_print_comp (x04cbc).\n%s\n",
                 fail.message);
          exit_status = 1;
          goto END;
        }
      printf("\n");
    }

 END:
  if (actder) NAG_FREE(actder);
  if (der_comp) NAG_FREE(der_comp);
  if (clabs) NAG_FREE(clabs);
  if (clabsc) NAG_FREE(clabsc);

  return exit_status;
}