/* nag_anova_random}(g04bbc) Example Program.
 *
 * Copyright 1998 Numerical Algorithms Group.
 *
 * Mark 5, 1998.
 *
 * Mark 6 revised, 2000.
 * Mark 8 revised, 2004.
 */

#include <nag.h>
#include <nag_stdlib.h>
#include <stdio.h>
#include <nagg04.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("%ld%ld%ld%*[^\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("%ld", &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:
  if (bmean) NAG_FREE(bmean);
  if (c) NAG_FREE(c);
  if (ef) NAG_FREE(ef);
  if (r) NAG_FREE(r);
  if (table) NAG_FREE(table);
  if (tmean) NAG_FREE(tmean);
  if (y) NAG_FREE(y);
  if (irep) NAG_FREE(irep);
  if (it) NAG_FREE(it);
  return exit_status;
}