/* nag_sum_fast_gauss (c06sac) Example Program.
 *
 * Copyright 2017 Numerical Algorithms Group.
 *
 * Mark 26.1, 2017.
 */

#include <stdio.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagc06.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;
}