/* nag_opt_nlin_lsq (e04unc) Example Program.
 *
 * Copyright 2014 Numerical Algorithms Group.
 *
 * Mark 5, 1998.
 *
 * Mark 6 revised, 2000.
 * Mark 7 revised, 2001.
 * Mark 8 revised, 2004.
 *
 */

#include <nag.h>
#include <stdio.h>
#include <string.h>
#include <nag_stdlib.h>
#include <math.h>
#include <nage04.h>

#ifdef __cplusplus
extern "C" {
#endif
static void NAG_CALL objfun(Integer m, Integer n, const double x[], double f[],
                            double fjac[], Integer tdfjac, Nag_Comm *comm);
static void NAG_CALL confun(Integer n, Integer ncnlin, const Integer needc[],
                            const double x[], double conf[], double cjac[],
                            Nag_Comm *comm);
#ifdef __cplusplus
}
#endif

static void NAG_CALL objfun(Integer m, Integer n, const double x[], double f[],
                            double fjac[], Integer tdfjac, Nag_Comm *comm)
{
#define FJAC(I, J) fjac[(I) *tdfjac + (J)]

  /* Initialized data */
  static double a[44] = { 8.0, 8.0, 10.0, 10.0, 10.0, 10.0, 12.0, 12.0, 12.0,
                          12.0, 14.0, 14.0, 14.0, 16.0, 16.0, 16.0, 18.0, 18.0,
                          20.0, 20.0, 20.0, 22.0, 22.0, 22.0, 24.0, 24.0, 24.0,
                          26.0, 26.0, 26.0, 28.0, 28.0, 30.0, 30.0, 30.0, 32.0,
                          32.0, 34.0, 36.0, 36.0, 38.0, 38.0, 40.0, 42.0 };

  /* Local variables */
  double        temp;
  Integer       i;
  double        x0, x1, ai;

  /* Function to evaluate the objective subfunctions
   * and their 1st derivatives.
   */
  x0 = x[0];
  x1 = x[1];
  for (i = 0; i < m; ++i)
    {
      ai = a[i];
      temp = exp(-x1 * (ai - 8.0));
      /* Evaluate objective subfunction f(i+1) only if required */
      if (comm->needf == i+1 || comm->needf == 0)
        f[i] = x0 + (.49 - x0) * temp;
      if (comm->flag == 2)
        {
          FJAC(i, 0) = 1.0 - temp;
          FJAC(i, 1) = -(.49 - x0) * (ai - 8.0) * temp;
        }
    }
} /* objfun */

static void NAG_CALL confun(Integer n, Integer ncnlin, const Integer needc[],
                            const double x[], double conf[], double cjac[],
                            Nag_Comm *comm)
{
#define CJAC(I, J) cjac[(I) *n + (J)]

  /* Function to evaluate the nonlinear constraints and its 1st derivatives. */

  if (comm->first == Nag_TRUE)
    {
      /* First call to confun.  Set all Jacobian elements to zero.
       *  Note that this will only work when  options.obj_deriv = TRUE
       *  (the default).
       */
      CJAC(0, 0) = CJAC(0, 1) = 0.0;
    }

  if (needc[0] > 0)
    {
      conf[0] = -0.09 - x[0]*x[1] + 0.49*x[1];

      if (comm->flag == 2)
        {
          CJAC(0, 0) = -x[1];
          CJAC(0, 1) = -x[0] + 0.49;
        }
    }
} /* confun */

#define A(I, J) a[(I) *tda + J]

int main(void)
{
  Integer     exit_status = 0, i, j, m, n, nbnd, nclin, ncnlin, tda, tdfjac;
  Nag_E04_Opt options;
  Nag_Comm    comm;
  double      *a = 0, *bl = 0, *bu = 0, *f = 0, *fjac = 0, objf, *x = 0;
  double      *y = 0;
  NagError    fail;

  INIT_FAIL(fail);

  printf("nag_opt_nlin_lsq (e04unc) Example Program Results\n");
  fflush(stdout);
  scanf(" %*[^\n]"); /* Skip heading in data file */

  /* Read problem dimensions */
  scanf(" %*[^\n]");
  scanf("%ld%ld%*[^\n]", &m, &n);
  scanf(" %*[^\n]");
  scanf("%ld%ld%*[^\n]", &nclin, &ncnlin);

  if (m > 0 && n > 0 && nclin >= 0 && ncnlin >= 0)
    {
      nbnd = n + nclin + ncnlin;
      if (!(x = NAG_ALLOC(n, double)) ||
          !(a = NAG_ALLOC(nclin*n, double)) ||
          !(f = NAG_ALLOC(m, double)) ||
          !(y = NAG_ALLOC(m, double)) ||
          !(fjac = NAG_ALLOC(m*n, double)) ||
          !(bl = NAG_ALLOC(nbnd, double)) ||
          !(bu = NAG_ALLOC(nbnd, double)))
        {
          printf("Allocation failure\n");
          exit_status = -1;
          goto END;
        }
      tda = n;
      tdfjac = n;
    }
  else
    {
      printf("Invalid m or n nclin or ncnlin.\n");
      exit_status = 1;
      return exit_status;
    }
  /* Read a, y, bl, bu and x from data file */

  if (nclin > 0)
    {
      scanf(" %*[^\n]");
      for (i = 0; i < nclin; ++i)
        for (j = 0; j < n; ++j)
          scanf("%lf", &A(i, j));
    }

  /* Read the y vector of the objective */
  scanf(" %*[^\n]");
  for (i = 0; i < m; ++i)
    scanf("%lf", &y[i]);

  /* Read lower bounds */
  scanf(" %*[^\n]");
  for (i = 0; i < n + nclin + ncnlin; ++i)
    scanf("%lf", &bl[i]);

  /* Read upper bounds */
  scanf(" %*[^\n]");
  for (i = 0; i < n + nclin + ncnlin; ++i)
    scanf("%lf", &bu[i]);

  /* Read the initial point x */
  scanf(" %*[^\n]");
  for (i = 0; i < n; ++i)
    scanf("%lf", &x[i]);

  /* Set an option */
  /* nag_opt_init (e04xxc).
   * Initialization function for option setting
   */
  nag_opt_init(&options);

  /* Solve the problem */
  /* nag_opt_nlin_lsq (e04unc), see above. */
  nag_opt_nlin_lsq(m, n, nclin, ncnlin, a, tda, bl, bu, y, objfun,
                   confun, x, &objf, f, fjac, tdfjac, &options,
                   &comm, &fail);
  if (fail.code != NE_NOERROR)
    {
      printf("Error from nag_opt_nlin_lsq (e04unc).\n%s\n",
              fail.message);
      exit_status = 1;
    }
  /* nag_opt_free (e04xzc).
   * Memory freeing function for use with option setting
   */
  nag_opt_free(&options, "all", &fail);
  if (fail.code != NE_NOERROR)
    {
      printf("Error from nag_opt_free (e04xzc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }

 END:
  NAG_FREE(x);
  NAG_FREE(a);
  NAG_FREE(f);
  NAG_FREE(y);
  NAG_FREE(fjac);
  NAG_FREE(bl);
  NAG_FREE(bu);

  return exit_status;
}