/* nag_mv_factor (g03cac) Example Program.
 *
 * Copyright 1998 Numerical Algorithms Group.
 *
 * Mark 5, 1998.
 *
 */

#include <nag.h>
#include <stdio.h>
#include <nag_stdlib.h>
#include <nage04.h>
#include <nagg03.h>

#define NMAX 9
#define MMAX 9

int main(void)
{
  double com[MMAX], e[MMAX], fl[MMAX][MMAX], psi[MMAX],
  res[MMAX*(MMAX-1)/2], stat[4], wt[NMAX],  x[NMAX][MMAX];
  double *wtptr=0;
  double eps;

  Integer nfac, nvar;
  Integer i, j, l, m, n;
  Integer  isx[MMAX];
  Integer tdfl = MMAX, tdx = MMAX;

  char weight[2], char_matrix[2];

  Nag_FacMat matrix;
  Nag_E04_Opt options;

  Vprintf("g03cac Example Program Results\n\n");

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

  Vscanf("%s",char_matrix);
  Vscanf("%s",weight);
  Vscanf("%ld",&n);
  Vscanf("%ld",&m);
  Vscanf("%ld",&nvar);
  Vscanf("%ld",&nfac);

  if (m <= MMAX && (*char_matrix == 'C' || n <= NMAX ))
    {
      if (*char_matrix == 'C')
        {
          for (i = 0; i < m; ++i)
            {
              for (j = 0; j < m; ++j)
                Vscanf("%lf",&x[i][j]);
            }
        }
      else
        {
          if (*weight == 'W' || *weight == 'w')
            {
              for (i = 0; i < n; ++i)
                {
                  for (j = 0; j < m; ++j)
                    Vscanf("%lf",&x[i][j]);
                  Vscanf("%lf",&wt[i]);
                }
              wtptr = wt;
            }
          else
            {
              for (i = 0; i < n; ++i)
                {
                  for (j = 0; j < m; ++j)
                    Vscanf("%lf",&x[i][j]);
                }
            }
        }
      for (j = 0; j < m; ++j)
        Vscanf("%ld",&isx[j]);
      if (*char_matrix == 'D')
        {
          matrix = Nag_DataCorr;
        }
      else if (*char_matrix == 'S')
        {
          matrix = Nag_DataCovar;
        }
      else if (*char_matrix == 'C')
        {
          matrix = Nag_MatCorr_Covar;
        }

      e04xxc(&options);
      options.max_iter = 500;
      options.optim_tol = 1e-2;
      eps = 1e-5;
      g03cac(matrix, n, m, &x[0][0], tdx, nvar, isx, nfac, wtptr, e,
             stat, com, psi, res, &fl[0][0], tdfl, &options, eps, NAGERR_DEFAULT);

      Vprintf("\nEigenvalues\n\n");
      for (j = 0; j < m; ++j)
        {
          Vprintf(" %12.4e%s",e[j], (j+1)%6==0 ? "\n" : "");
        }
      Vprintf("\n\n%s%6.3f\n","    Test Statistic = ",stat[1]);
      Vprintf("%s%6.3f\n","                df = ",stat[2]);
      Vprintf("%s%6.3f\n\n","Significance level = ",stat[3]);
      Vprintf("Residuals\n\n");
      l = 1;
      for (i = 1; i <= nvar-1; ++i)
        {
          for (j = l; j <= l+i-1; ++j)
            Vprintf(" %8.3f",res[j-1]);
          Vprintf("\n");
          l += i;
        }
      Vprintf("\nLoadings, Communalities and PSI\n\n");
      for (i = 0; i < nvar; ++i)
        {
          for (j = 0; j < nfac; ++j)
            Vprintf(" %8.3f",fl[i][j]);
          Vprintf("%8.3f%8.3f\n", com[i], psi[i]);
        }
      return EXIT_SUCCESS;
    }
  else
    {
      Vprintf("Incorrect input value of n or m.\n");
      return EXIT_FAILURE;
    }
}