/* 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, *irep=0, *it=0, j, n, nblock, nt, tdc;
  Integer tdtable;
  NagError fail;
  double *bmean=0, *c=0, *ef=0, gmean, *r=0, *table=0, *tmean=0, tol, *y=0;

  INIT_FAIL(fail);
  Vprintf("nag_anova_random (g04bbc) Example Program Results\n");

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

  if (nblock>=2 && !(n%nblock))
    {
      if ( !( bmean = NAG_ALLOC(nblock, double)) ||
           !( 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)) )
        {
          Vprintf("Allocation failure\n");
          exit_status = -1;
          goto END;
        }
      tdc = nt;
      tdtable = 5;
    }
  else
    {
      Vprintf("Invalid nblock or n.\n");
      exit_status = 1;
      return exit_status;
    }
  for (i = 0; i < n; ++i)
    Vscanf("%lf",&y[i]);
  Vscanf("%*[^\n]");
  for (i = 0; i < n; ++i)
    Vscanf("%ld",&it[i]);
  Vscanf("%*[^\n]");

  tol = 5e-6;
  irdf = 0;

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

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

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

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

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

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

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

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

  Vprintf("\nStandard errors of differences between means\n\n");
  for (i = 1; i < nt; ++i)
    {
      for (j = 0; j < i; ++j)
        Vprintf("%10.5f",C(i,j));
      Vprintf("\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;

}