/* nag_anova_factorial (g04cac) Example Program.
 *
 * Copyright 1998 Numerical Algorithms Group.
 *
 * Mark 5, 1998.
 *
 * Mark 6 revised, 2000.
 * Mark 7b revised, 2004.
 * Mark 8 revised, 2004.
 */

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

#define MTERM 6

#define LFAC(I) lfac[(I)-1]
#define IWK(I) iwk[(I)-1]
#define IMEAN(I) imean[(I)-1]
#define Y(I) y[(I)-1]
#define TMEAN(I) tmean[(I)-1]
#define SEMEAN(I) semean[(I)-1]
#define R(I) r[(I)-1]
#define E(I) e[(I)-1]
#define BMEAN(I) bmean[(I)-1]
#define TABLE(I,J) table[((I)-1) * ( 5) + ((J)-1)]
int main(void)
{

  Integer c__27=27, exit_status=0, i, *imean=0, inter, irdf, j, k, l, *lfac=0;
  Integer mterm=MTERM, n, nblock, nfac, ntreat, num;
  NagError fail;
  double *bmean=0, *e=0, *r=0, *semean=0, *table=0, *tmean=0, *y=0;
  INIT_FAIL(fail);

  Vprintf("nag_anova_factorial (g04cac) Example Program Results\n\n");

  /* Skip heading in data file */
  Vscanf("%*[^\n]");

  Vscanf("%ld%ld%ld%ld%*[^\n]",&n, &nblock,
         &nfac, &inter);
  if (n>=4 && nfac>=1 && nblock >=1 && !(n%nblock))
    {
      if ( !( r = NAG_ALLOC(n, double)) ||
           !( y = NAG_ALLOC(n, double))  ||
           !( lfac = NAG_ALLOC(nfac, Integer)) )
        {
          Vprintf("Allocation failure\n");
          exit_status = -1;
          goto END;
        }
    }
  else
    {
      Vprintf("Invalid n or nfac or nblock.\n");
      exit_status = 1;
      return exit_status;
    }
  for (j = 0; j < nfac; ++j)
    Vscanf("%ld",&lfac[j]);
  Vscanf("%*[^\n]");

  for (i = 0; i < n; ++i)
    Vscanf("%lf",&y[i]);
  Vscanf("%*[^\n]");

  irdf = 0;
  /* nag_anova_factorial (g04cac).
   * Complete factorial design
   */
  nag_anova_factorial(n, y, nfac, lfac, nblock, inter, irdf, &mterm, &table,
                      &tmean, &c__27, &e, &imean, &semean, &bmean, r,
                      &fail);
  if (fail.code != NE_NOERROR)
    {
      Vprintf("Error from nag_anova_factorial (g04cac).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }

  Vprintf("\n  ANOVA table\n\n");
  Vprintf("   Source      df         SS          MS           F"
          "          Prob\n\n");
  k = 0;
  if (nblock > 1)
    {
      ++k;
      Vprintf("%s   "," Blocks  ");
      Vprintf("%4ld  ",(Integer)TABLE(1,1));
      for (j = 2; j <= 5; ++j)
        Vprintf("%12.2f",TABLE(1,j));
      Vprintf("\n");
    }
  ntreat = mterm - 2 - k;
  for (i = 1; i <= ntreat; ++i)
    {
      Vprintf("%s%2ld "," Effect  ",i);
      Vprintf("%4ld  ",(Integer)TABLE(k+i,1));
      for (j = 2; j <= 5; ++j)
        Vprintf("%12.2f",TABLE(k+i,j));
      Vprintf("\n");
    }
  Vprintf("%s   "," Residual");
  Vprintf("%4ld  ",(Integer)TABLE(mterm-1,1));
  for (j = 2; j <= 3; ++j)
    Vprintf("%12.2f",TABLE(mterm-1,j));
  Vprintf("\n");

  Vprintf("%s   "," Total   ");
  Vprintf("%4ld  ",(Integer)TABLE(mterm,1));
  Vprintf("%12.2f\n\n",TABLE(mterm,2));

  Vprintf("  Treatment Means and Standard Errors \n\n");
  k = 0;
  for (i = 0; i < ntreat; ++i)
    {
      l = imean[i];
      Vprintf("%s%2ld\n\n"," Effect ",i+1);

      num=1;
      for (j = k; j < l; ++j)
        {
          Vprintf("%10.2f%s",tmean[j],num%8?"":"\n");
          num++;
        }

      Vprintf("\n\n%s%10.2f\n\n"," SE of difference in means  = ",semean[i]);
      k = l;
    }
  /* nag_anova_factorial_free (g04czc).
   * Memory freeing function for nag_anova_factorial (g04cac)
   */
  nag_anova_factorial_free(&table, &tmean, &e, &imean, &semean, &bmean);
 END:
  if (r) NAG_FREE(r);
  if (y) NAG_FREE(y);
  if (lfac) NAG_FREE(lfac);
  return exit_status;
}