/* nag_mv_fac_score (g03ccc) 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>
#include <math.h>

#define NMAX 20
#define MMAX 10

int main(void)
{

  double com[MMAX], e[MMAX], fl[MMAX][MMAX], fs[MMAX][MMAX],
  psi[MMAX], r[MMAX][MMAX], stat[4],
  wt[NMAX], x[NMAX][MMAX];
  double *wtptr=0;
  double eps;

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

  char char_method[2], char_weight[2], char_matrix[2];

  Nag_FacMat matrix;
  Nag_E04_Opt options;
  Nag_FacScoreMethod  method;

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

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

  Vscanf("%s",char_matrix);
  Vscanf("%s",char_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 (*char_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-3;
      eps = 1e-5;
      g03cac(matrix, n, m, &x[0][0], tdx, nvar, isx, nfac, wtptr, e,
             stat, com, psi, &r[0][0], &fl[0][0], tdfl, &options, eps, NAGERR_DEFAULT);
      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]);
        }
      Vscanf("%s",char_method);
      if (*char_method == 'R')
        {
          method = Nag_FacScoreRegsn;
        }
      else if (*char_method == 'B')
        {
          method = Nag_FacScoreBart;
        }

      g03ccc(method, Nag_FacNoRotate, nvar, nfac, &fl[0][0], tdfl, psi, e,
             &r[0][0], tdr, &fs[0][0], tdfs, NAGERR_DEFAULT);
      Vprintf("\nFactor score coefficients\n\n");
      for (i = 0; i < nvar; ++i)
        {
          for (j = 0; j < nfac; ++j)
            Vprintf("  %8.3f",fs[i][j]);
          Vprintf("\n");
        }
  return EXIT_SUCCESS;
    }
  else
    {
      Vprintf("Incorrect input value of n or m.\n");
      return EXIT_FAILURE;
    }
}