/* nag_roots_sparsys_func_easy (c05qsc) Example Program.
 *
 * Copyright 2023 Numerical Algorithms Group.
 *
 * Mark 29.2, 2023.
 */

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

static void NAG_CALL fcn(Integer n, Integer lindf, const Integer indf[],
                         const double x[], double fvec[], Nag_Comm *comm,
                         Integer *iflag) {
  double theta;
  Integer i, ind;

  *iflag = 0;
  theta = (double)(comm->iuser[0]) * pow(0.5, 7);
  for (ind = 0; ind < lindf; ind++) {
    i = indf[ind] - 1;
    fvec[i] = (3.0 - (2.0 + theta) * x[i]) * x[i] + 1.0;
    if (i > 0)
      fvec[i] = fvec[i] - x[i - 1];
    if (i < n - 1)
      fvec[i] = fvec[i] - 2.0 * x[i + 1];
  }
}

int main(void) {
  Integer exit_status = 0, n = 9, i, j, licomm, lrcomm;
  double fnorm, xtol;
  Nag_Boolean init;
  Nag_Comm comm;
  Integer iuser[1], *icomm = 0;
  double ruser[1], *rcomm = 0, *fvec = 0, *x = 0;
  NagError fail;

  printf("nag_roots_sparsys_func_easy (c05qsc) Example Program Results\n");
  lrcomm = 12 + 3 * n;
  licomm = 8 * n + 19 + 3 * n;
  if (!(fvec = NAG_ALLOC(n, double)) || !(x = NAG_ALLOC(n, double)) ||
      !(rcomm = NAG_ALLOC(lrcomm, double)) ||
      !(icomm = NAG_ALLOC(licomm, Integer))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }
  comm.iuser = iuser;
  comm.user = ruser;
  xtol = sqrt(nag_machine_precision);
  /* The following starting values provide a rough solution. */
  for (j = 0; j < n; j++)
    x[j] = -1.0E0;
  for (i = 0; i <= 1; i++) {
    INIT_FAIL(fail);
    /* Perturb the system? */
    comm.iuser[0] = i;
    init = (i == 0 ? Nag_TRUE : Nag_FALSE);
    /* nag_roots_sparsys_func_easy (c05qsc).
     * Solution of a sparse system of nonlinear equations using function
     * values only (easy-to-use).
     */
    nag_roots_sparsys_func_easy(fcn, n, x, fvec, xtol, init, rcomm, lrcomm,
                                icomm, licomm, &comm, &fail);
    if (fail.code == NE_NOERROR) {
      /* Compute Euclidean norm. */
      fnorm = 0.0E0;
      for (j = 0; j < n; j++)
        fnorm += pow(fvec[j], 2);
      fnorm = sqrt(fnorm);
      printf("\nFinal 2-norm of the residuals = %12.4e\n", fnorm);
      printf("\nFinal approximate solution\n\n");
      for (j = 0; j < n; j++)
        printf("%12.4f%s", x[j], (j + 1) % 3 ? " " : "\n");
      printf("\n");
    } else {
      printf("Error from nag_roots_sparsys_func_easy (c05qsc).\n%s\n",
             fail.message);
      if (fail.code == NE_TOO_MANY_FEVALS || fail.code == NE_TOO_SMALL ||
          fail.code == NE_NO_IMPROVEMENT) {
        printf("\nApproximate solution\n");
        for (j = 0; j < n; j++)
          printf("%12.4f%s", x[j], (j + 1) % 3 ? " " : "\n");
        printf("\n");
      }
    }
  }
END:
  NAG_FREE(fvec);
  NAG_FREE(x);
  NAG_FREE(rcomm);
  NAG_FREE(icomm);
  return exit_status;
}