/* nag_roots_sys_func_aa_rcomm (c05mdc) Example Program.
 *
 * Copyright 2023 Numerical Algorithms Group.
 *
 * Mark 29.3, 2023.
 */

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

int main(void) {
  Integer exit_status = 0, i, n = 4, irevcm = 0, m = 2, astart = 0, icount = 0,
          imax = 50, exit_loop = 0;
  double *fvec = 0, *x = 0, *rwsav = 0;
  Integer *iwsav = 0;
  double cndtol = 0.0, rtol, atol;
  /* Nag Types */
  NagError fail;

  INIT_FAIL(fail);

  printf("nag_roots_sys_func_aa_rcomm (c05mdc) Example Program Results\n");

  if (!(fvec = NAG_ALLOC(n, double)) || !(x = NAG_ALLOC(n, double)) ||
      !(iwsav = NAG_ALLOC(14 + m, Integer)) ||
      !(rwsav = NAG_ALLOC(2 * m * n + m * m + m + 2 * n + 1 +
                              MIN(m, 1) * MAX(n, 3 * m),
                          double))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  /* The following starting values provide a rough solution. */
  x[0] = 2.0;
  x[1] = 0.5;
  x[2] = 2.0;
  x[3] = 0.5;

  /* nag_machine_precision (x02ajc).
   * The machine precision
   */
  atol = sqrt(nag_machine_precision);
  rtol = sqrt(nag_machine_precision);

  /* nag_roots_sys_func_aa_rcomm (c05mdc).
   * Solution of a system of nonlinear equations using Anderson acceleration
   * (reverse communication)
   */

  do {
    nag_roots_sys_func_aa_rcomm(&irevcm, n, x, fvec, atol, rtol, m, cndtol,
                                astart, iwsav, rwsav, &fail);

    switch (irevcm) {
    case 1:
      if (icount == imax) {
        printf("Exiting after the maximum number of iterations\n\n");
        exit_loop = 1;
      } else {
        icount++;
      }
      /* x and fvec are available for printing */
      break;
    case 2:
      fvec[0] = cos(x[2]) - x[0];
      fvec[1] = sqrt(1.0 - x[3] * x[3]) - x[1];
      fvec[2] = sin(x[0]) - x[2];
      fvec[3] = x[1] * x[1] - x[3];
      break;
    }
  } while (irevcm != 0 && exit_loop == 0);

  if (fail.code != NE_NOERROR) {
    printf("Error from nag_roots_sys_func_aa_rcomm (c05mdc).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }

  printf("Final approximate solution after %" NAG_IFMT " iterations\n\n",
         icount);
  for (i = 0; i < n; i++)
    printf("%12.4f  ", x[i]);

  printf("\n");

  if (fail.code != NE_NOERROR)
    exit_status = 2;

END:
  NAG_FREE(fvec);
  NAG_FREE(x);
  NAG_FREE(iwsav);
  NAG_FREE(rwsav);
  return exit_status;
}