/* nag_zero_sparse_nonlin_eqns_easy (c05qsc) Example Program.
 *
 * NAGPRODCODE Version.
 *
 * Copyright 2016 Numerical Algorithms Group.
 *
 * Mark 26, 2016.
 */

#include <nag.h>
#include <stdio.h>
#include <nag_stdlib.h>
#include <math.h>
#include <nagc05.h>
#include <nagx02.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_zero_sparse_nonlin_eqns_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_zero_sparse_nonlin_eqns_easy (c05qsc).
     * Solution of a sparse system of nonlinear equations using function
     * values only (easy-to-use).
     */
    nag_zero_sparse_nonlin_eqns_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_zero_sparse_nonlin_eqns_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;
}