/* nag_opt_nlin_lsq (e04unc) Example Program.
 *
 * Copyright 2022 Numerical Algorithms Group.
 *
 * Mark 28.3, 2022.
 *
 *
 */

#include <math.h>
#include <nag.h>
#include <stdio.h>
#include <string.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] = -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("%" NAG_IFMT "%" NAG_IFMT "%*[^\n]", &m, &n);
  scanf(" %*[^\n]");
  scanf("%" NAG_IFMT "%" NAG_IFMT "%*[^\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;
}