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

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

#define NMAX 54
#define MAXF 2
#define MAXT 27
#define MTERM 6
#define BMAX 4
#define LDT MTERM

int main(void)
{

  double  r[NMAX], y[NMAX];
  double *bmean=0, *e=0, *semean=0,
  *table=0, *tmean=0;

  Integer c__27 = 27;
  Integer lfac[MAXF];
  Integer *imean=0;
  Integer mterm = MTERM;
  Integer nfac, irdf;
  Integer i, j, k, l, n;
  Integer num;
  Integer inter, nblock;
  Integer itotal, ntreat;

#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)]

  Vprintf("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 <= NMAX && nblock <= BMAX - 1 && nfac <= MAXF)
    {
      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;
      g04cac(n, y, nfac, lfac, nblock, inter, irdf, &mterm, &table,
              &tmean, &c__27, &e, &imean, &semean, &bmean, r,
              NAGERR_DEFAULT);

      itotal = mterm;
      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;
        }
      g04czc(&table, &tmean, &e, &imean, &semean, &bmean);
      return EXIT_SUCCESS;
    }
  else
    {
      Vprintf( "Incorrect input value of n or nblock or nfac.\n");
      return EXIT_FAILURE;
    }
}