/* nag_anova_random}(g04bbc) Example Program.
 *
 * Copyright 2023 Numerical Algorithms Group.
 *
 * Mark 29.1, 2023.
 *
 */

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

#define C(I, J) c[(I)*tdc + J]
#define TABLE(I, J) table[(I)*tdtable + J]
int main(void) {

  Integer exit_status = 0, i, irdf, j, n, nblock, nt, tdc, tdtable;
  Integer *irep = 0, *it = 0;
  NagError fail;
  Nag_Blocks blocks;
  double gmean, tol;
  double *bmean = 0, *c = 0, *ef = 0, *r = 0, *table = 0, *tmean = 0;
  double *y = 0;

  INIT_FAIL(fail);

  printf("nag_anova_random (g04bbc) Example Program Results\n");

  /* Skip heading in data file */
  scanf("%*[^\n] ");
  scanf("%" NAG_IFMT "%" NAG_IFMT "%" NAG_IFMT "%*[^\n] ", &n, &nt, &nblock);

  if ((nblock >= 2 && !(n % nblock)) || (nblock == 0)) {
    if (nblock != 0) {
      if (!(bmean = NAG_ALLOC(nblock, double))) {
        printf("Allocation failure\n");
        exit_status = -1;
        goto END;
      }
      blocks = Nag_SerialBlocks;
    } else {
      blocks = Nag_NoBlocks;
    }

    if (!(c = NAG_ALLOC((nt) * (nt), double)) ||
        !(ef = NAG_ALLOC(nt, double)) || !(r = NAG_ALLOC(n, double)) ||
        !(table = NAG_ALLOC((4) * (5), double)) ||
        !(tmean = NAG_ALLOC(nt, double)) || !(y = NAG_ALLOC(n, double)) ||
        !(irep = NAG_ALLOC(nt, Integer)) || !(it = NAG_ALLOC(n, Integer))) {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }
    tdc = nt;
    tdtable = 5;
  } else {
    printf("Invalid nblock or n.\n");
    exit_status = 1;
    return exit_status;
  }
  for (i = 0; i < n; ++i)
    scanf("%lf", &y[i]);
  scanf("%*[^\n]");
  for (i = 0; i < n; ++i)
    scanf("%" NAG_IFMT "", &it[i]);
  scanf("%*[^\n]");

  tol = 5e-6;
  irdf = 0;

  /* nag_anova_random (g04bbc).
   * General block design or completely randomized design
   */
  nag_anova_random(n, y, blocks, nblock, nt, it, &gmean, bmean, tmean, table, c,
                   tdc, irep, r, ef, tol, irdf, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_anova_random (g04bbc).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }

  printf("\nANOVA table\n\n");
  printf("   Source        df         SS          MS          F"
         "          Prob\n\n");
  printf(" Blocks      ");
  for (j = 0; j < 5; ++j)
    printf("%10.4f  ", TABLE(0, j));
  printf("\n");

  printf(" Treatments  ");
  for (j = 0; j < 5; ++j)
    printf("%10.4f  ", TABLE(1, j));
  printf("\n");

  printf(" Residual    ");
  for (j = 0; j < 3; ++j)
    printf("%10.4f  ", TABLE(2, j));
  printf("\n");

  printf(" Total       ");
  for (j = 1; j <= 2; ++j)
    printf("%10.4f  ", TABLE(3, j - 1));
  printf("\n");

  printf("\nEfficiency Factors\n\n");
  for (i = 0; i < nt; ++i)
    printf("%10.5f", ef[i]);
  printf("\n");

  printf("\n%s%10.5f\n", "  Grand Mean", gmean);

  printf("\nTreatment Means\n\n");
  for (i = 1; i <= nt; ++i)
    printf("%10.5f", tmean[i - 1]);
  printf("\n");

  printf("\nStandard errors of differences between means\n\n");
  for (i = 1; i < nt; ++i) {
    for (j = 0; j < i; ++j)
      printf("%10.5f", C(i, j));
    printf("\n");
  }
END:
  NAG_FREE(bmean);
  NAG_FREE(c);
  NAG_FREE(ef);
  NAG_FREE(r);
  NAG_FREE(table);
  NAG_FREE(tmean);
  NAG_FREE(y);
  NAG_FREE(irep);
  NAG_FREE(it);
  return exit_status;
}