/* nag_zero_sparse_nonlin_eqns_easy (c05qsc) Example Program.
 *
 * Copyright 2011 Numerical Algorithms Group.
 *
 * Mark 23, 2011.
 */

#include <nag.h>
#include <stdio.h>
#include <nag_stdlib.h>
#include <math.h>
#include <nagc05.h>
#include <nagx02.h>

static void 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:
  if (fvec) NAG_FREE(fvec);
  if (x) NAG_FREE(x);
  if (rcomm) NAG_FREE(rcomm);
  if (icomm) NAG_FREE(icomm);
  return exit_status;
}