/* nag_sum_fast_gauss (c06sac) Example Program.
 *
 * Copyright 2023 Numerical Algorithms Group.
 *
 * Mark 29.3, 2023.
 */

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

int main(void) {

  /* Scalars */
  Integer d, n, m, p1, p2, k, i, j;
  Integer exit_status = 0;
  double tol;

  /* Arrays */
  double *srcs = 0, *trgs = 0, *q = 0, *hin = 0, *v = 0, *term = 0;

  /* Nag Types */
  NagError fail;

  INIT_FAIL(fail);

  printf("nag_sum_fast_gauss (c06sac) Example Program Results\n");

  /* Read dimensions of arrays from data file. */
  scanf("%*[^\n] ");
  scanf("%" NAG_IFMT "%" NAG_IFMT "%" NAG_IFMT "%*[^\n]", &d, &n, &m);
  scanf("%" NAG_IFMT "%" NAG_IFMT "%" NAG_IFMT "%*[^\n]", &p1, &p2, &k);
  scanf("%lf %*[^\n]", &tol);

  /* Allocate arrays accordingly */
  if (!(srcs = NAG_ALLOC((d * n), double)) ||
      !(trgs = NAG_ALLOC((d * m), double)) || !(q = NAG_ALLOC((n), double)) ||
      !(hin = NAG_ALLOC((n), double)) || !(v = NAG_ALLOC((m), double)) ||
      !(term = NAG_ALLOC((m), double))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  /* Read array values from data file */
  for (i = 0; i < n; i++) {
    scanf("%lf", &q[i]);
  }
  scanf("%*[^\n]");

  for (i = 0; i < n; i++) {
    scanf("%lf", &hin[i]);
  }
  scanf("%*[^\n]");

  /* Data in file is stored by column but it is read into rows here
   * First for sources
   */
#define SRCS(I, J) srcs[(J - 1) * d + I - 1]
  for (j = 1; j <= n; j++) {
    for (i = 1; i <= d; i++) {
      scanf("%lf", &SRCS(i, j));
    }
  }
  scanf("%*[^\n]");

#define TRGS(I, J) trgs[(J - 1) * d + I - 1]
  /* And then for targets */
  for (j = 1; j <= m; j++) {
    for (i = 1; i <= d; i++) {
      scanf("%lf", &TRGS(i, j));
    }
  }

  nag_sum_fast_gauss(d, srcs, n, trgs, m, q, &p1, &p2, &k, hin, n, tol, v, term,
                     &fail);

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

  printf("\n      y        FGT(y)     term\n\n");

  for (i = 0; i < m; i++) {
    for (j = 0; j < d; j++)
      printf(" %4.1f", trgs[d * i + j]);
    printf("   %8.3f    %8.6f\n", v[i], term[i]);
  }

END:
  NAG_FREE(srcs);
  NAG_FREE(trgs);
  NAG_FREE(q);
  NAG_FREE(hin);
  NAG_FREE(v);
  NAG_FREE(term);
  return exit_status;
}