/* nag_opt_qpconvex1_sparse_solve (e04nkc) Example Program.
 *
 * Copyright 2024 Numerical Algorithms Group.
 *
 * Mark 30.3, 2024.
 *
 */

#include <nag.h>
#include <stdio.h>
#include <string.h>

#ifdef __cplusplus
extern "C" {
#endif
static void NAG_CALL qphess(Integer ncolh, const double x[], double hx[],
                            Nag_Comm *comm);
#ifdef __cplusplus
}
#endif

/* Declare a data structure for passing sparse Hessian data to qphess */
typedef struct {
  double *hess;
  Integer *khess;
  Integer *hhess;
} HessianData;

#define NAMES(I, J) names[(I)*9 + J]

int main(void) {
  HessianData hess_data;
  Integer exit_status = 0, *ha = 0, *hhess = 0, i, icol, iobj, j, jcol;
  Integer *ka = 0, *khess = 0, m, n, nbnd, ncolh, ninf, nnz, nnz_hess;
  Nag_Comm comm;
  Nag_E04_Opt options;
  char **crnames = 0, *names = 0;
  double *a = 0, *bl = 0, *bu = 0, *hess = 0, obj, sinf, *x = 0;
  NagError fail;

  INIT_FAIL(fail);

  printf("nag_opt_qpconvex1_sparse_solve (e04nkc) Example Program Results\n");
  fflush(stdout);

  /* Skip heading in data file */
  scanf(" %*[^\n]");
  /* Read the problem dimensions */
  scanf(" %*[^\n]");
  scanf("%" NAG_IFMT "%" NAG_IFMT "", &n, &m);

  /* Read nnz, iobj, ncolh */
  scanf(" %*[^\n]");
  scanf("%" NAG_IFMT "%" NAG_IFMT "%" NAG_IFMT "", &nnz, &iobj, &ncolh);

  if (n >= 1 && m >= 1 && nnz >= 1 && nnz <= n * m) {
    nbnd = n + m;
    if (!(a = NAG_ALLOC(nnz, double)) || !(bl = NAG_ALLOC(nbnd, double)) ||
        !(bu = NAG_ALLOC(nbnd, double)) || !(x = NAG_ALLOC(nbnd, double)) ||
        !(ha = NAG_ALLOC(nnz, Integer)) || !(ka = NAG_ALLOC(n + 1, Integer)) ||
        !(khess = NAG_ALLOC(n + 1, Integer)) ||
        !(crnames = NAG_ALLOC(nbnd, char *)) ||
        !(names = NAG_ALLOC(nbnd * 9, char))) {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }
  } else {
    printf("Invalid n or m or nnz.\n");
    exit_status = 1;
    return exit_status;
  }

  /* Read the matrix and set up ka */
  jcol = 1;
  ka[jcol - 1] = 0;
  scanf(" %*[^\n]");
  for (i = 0; i < nnz; ++i) {
    /* a[i] stores the (ha[i], icol) element of matrix */
    scanf("%lf%" NAG_IFMT "%" NAG_IFMT "", &a[i], &ha[i], &icol);

    /* Check whether we have started a new column */
    if (icol == jcol + 1) {
      ka[icol - 1] = i; /* Start of icol-th column in a */
      jcol = icol;
    } else if (icol > jcol + 1) {
      /* Index in a of the start of the icol-th column
       * equals i, but columns jcol+1, jcol+2, ...,
       * icol-1 are empty.  Set the corresponding elements
       * of ka to i.
       */
      for (j = jcol + 1; j < icol; ++j)
        ka[j - 1] = i;

      ka[icol - 1] = i;
      jcol = icol;
    }
  }
  ka[n] = nnz;

  /* If the last columns are empty, set ka accordingly */
  if (n > icol) {
    for (j = icol; j <= n - 1; ++j)
      ka[j] = nnz;
  }

  /* Read the bounds */
  nbnd = n + m;
  scanf(" %*[^\n]"); /* Skip heading in data file */
  for (i = 0; i < nbnd; ++i)
    scanf("%lf", &bl[i]);
  scanf(" %*[^\n]");
  for (i = 0; i < nbnd; ++i)
    scanf("%lf", &bu[i]);

    /* Read the column and row names */
  scanf(" %*[^\n]"); /* Skip heading in data file */
  scanf(" %*[^']");
  for (i = 0; i < nbnd; ++i) {
    scanf(" '%8c'", &NAMES(i, 0));
    NAMES(i, 8) = '\0';
    crnames[i] = &NAMES(i, 0);
  }

  /* Read the initial estimate of x */
  scanf(" %*[^\n]"); /* Skip heading in data file */
  for (i = 0; i < n; ++i)
    scanf("%lf", &x[i]);

    /* Read nnz_hess */
  scanf(" %*[^\n]");
  scanf("%" NAG_IFMT "", &nnz_hess);

  if (!(hess = NAG_ALLOC(nnz_hess, double)) ||
      !(hhess = NAG_ALLOC(nnz_hess, Integer))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  /* Read the hessian matrix and set up khess */
  jcol = 1;
  khess[jcol - 1] = 0;
  scanf(" %*[^\n]");
  for (i = 0; i < nnz_hess; ++i) {
    /* hess[i] stores the (hhess[i], icol) element of matrix */
    scanf("%lf%" NAG_IFMT "%" NAG_IFMT "", &hess[i], &hhess[i], &icol);

    /* Check whether we have started a new column */
    if (icol == jcol + 1) {
      khess[icol - 1] = i; /* Start of icol-th column in hess */
      jcol = icol;
    } else if (icol > jcol + 1) {
      /* Index in hess of the start of the icol-th column
       * equals i, but columns jcol+1, jcol+2, ...,
       * icol-1 are empty.  Set the corresponding elements
       * of khess to i.
       */
      for (j = jcol + 1; j < icol; ++j)
        khess[j - 1] = i;

      khess[icol - 1] = i;
    }
  }
  khess[ncolh] = nnz_hess;

  /* Initialize options structure */
  /* nag_opt_init (e04xxc).
   * Initialization function for option setting
   */
  nag_opt_init(&options);
  options.crnames = crnames;

  /* Package up Hessian data for communication via comm */
  hess_data.hess = hess;
  hess_data.khess = khess;
  hess_data.hhess = hhess;

  comm.p = (Pointer)&hess_data;

  /* Solve the problem */
  /* nag_opt_qpconvex1_sparse_solve (e04nkc), see above. */
  nag_opt_qpconvex1_sparse_solve(n, m, nnz, iobj, ncolh, qphess, a, ha, ka, bl,
                                 bu, x, &ninf, &sinf, &obj, &options, &comm,
                                 &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_opt_qpconvex1_sparse_solve (e04nkc).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }

  printf("\nPerturb the problem and re-solve with warm start.\n");
  fflush(stdout);
  for (i = 0; i < nnz_hess; ++i)
    hess[i] *= 1.001;

  options.start = Nag_Warm;
  options.print_level = Nag_Soln;
  /* nag_opt_qpconvex1_sparse_solve (e04nkc), see above. */
  nag_opt_qpconvex1_sparse_solve(n, m, nnz, iobj, ncolh, qphess, a, ha, ka, bl,
                                 bu, x, &ninf, &sinf, &obj, &options, &comm,
                                 &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_opt_qpconvex1_sparse_solve (e04nkc).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }

  /* Free memory NAG-allocated to members of options */
  /* nag_opt_free (e04xzc).
   * Memory freeing function for use with option setting
   */
  nag_opt_free(&options, "", &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(a);
  NAG_FREE(bl);
  NAG_FREE(bu);
  NAG_FREE(x);
  NAG_FREE(hess);
  NAG_FREE(ha);
  NAG_FREE(ka);
  NAG_FREE(hhess);
  NAG_FREE(khess);
  NAG_FREE(crnames);
  NAG_FREE(names);

  return exit_status;
}

static void NAG_CALL qphess(Integer ncolh, const double x[], double hx[],
                            Nag_Comm *comm) {
  Integer i, j, jrow;
  HessianData *hd = (HessianData *)(comm->p);
  double *hess = hd->hess;
  Integer *hhess = hd->hhess;
  Integer *khess = hd->khess;

  for (i = 0; i < ncolh; ++i)
    hx[i] = 0.0;

  for (i = 0; i < ncolh; ++i) {
    /* For each column of Hessian */
    for (j = khess[i]; j < khess[i + 1]; ++j) {
      /* For each nonzero in column */

      jrow = hhess[j] - 1;

      /* Using symmetry of hessian, add contribution
       * to hx of hess[j] as a diagonal or upper
       * triangular element of hessian.
       */
      hx[i] += hess[j] * x[jrow];

      /* If hess[j] not a diagonal element add its
       * contribution to hx as a lower triangular
       * element of hessian.
       */
      if (jrow != i)
        hx[jrow] += hess[j] * x[i];
    }
  }
} /* qphess */