/* nag_mv_procustes (g03bcc) Example Program.
 *
 * Copyright 1998 Numerical Algorithms Group.
 *
 * Mark 5, 1998.
 *
 */

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

#define NMAX 3
#define MMAX 2

int main(void)
{
  double r[MMAX][MMAX], res[NMAX],
  x[NMAX][MMAX], y[NMAX][MMAX], yhat[NMAX][MMAX];
  double alpha;
  double  rss;

  Integer i, j, m, n;
  Integer tdx = MMAX, tdr = MMAX, tdy = MMAX;

  char char_scale[2];
  char char_stand[2];

  Nag_TransNorm stand;
  Nag_RotationScale scale;

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

  /* Skip heading in data file */
  Vscanf("%*[^\n]");
  Vscanf("%ld",&n);
  Vscanf("%ld",&m);
  Vscanf("%s",char_stand);
  Vscanf("%s",char_scale);

  if (n <= NMAX && m <= MMAX)
    {
      for (i = 0; i < n; ++i)
        {
          for (j = 0; j < m; ++j)
            Vscanf("%lf",&x[i][j]);
        }
      for (i = 0; i < n; ++i)
        {
          for (j = 0; j < m; ++j)
            Vscanf("%lf",&y[i][j]);
        }

      if (*char_stand == 'N')
        {
          stand = Nag_NoTransNorm;
        }
      else if (*char_stand == 'Z')
        {
          stand = Nag_Orig;
        }
      else if (*char_stand == 'C')
        {
          stand = Nag_OrigCentroid;
        }
      else if (*char_stand == 'U')
        {
          stand = Nag_Norm;
        }
      else if (*char_stand == 'S')
        {
          stand = Nag_OrigNorm;
        }
      else if (*char_stand == 'M')
        {
          stand = Nag_OrigNormCentroid;
        }

      if (*char_scale == 'S')
        {
          scale = Nag_LsqScale;
        }
      else if (*char_scale == 'U')
        {
          scale = Nag_NotLsqScale;
        }

      g03bcc(stand, scale, n, m, &x[0][0], tdx, &y[0][0], tdy,
             &yhat[0][0], &r[0][0], tdr,
             &alpha, &rss, res, NAGERR_DEFAULT);

      Vprintf("\n          Rotation Matrix\n\n");
      for (i = 0; i < m; ++i)
        {
          for (j = 0; j < m; ++j)
            Vprintf("  %7.3f  ",r[i][j]);
          Vprintf("\n");
        }
      if (*char_scale == 'S' || *char_scale == 's')
        {
          Vprintf("\n%s%10.3f\n"," Scale factor = ",alpha);
        }
      Vprintf("\n         Target Matrix \n\n");
      for (i = 0; i < n; ++i)
        {
          for (j = 0; j < m; ++j)
            Vprintf("  %7.3f  ",y[i][j]);
          Vprintf("\n");
        }
      Vprintf("\n         Fitted Matrix\n\n");
      for (i = 0; i < n; ++i)
        {
          for (j = 0; j < m; ++j)
            Vprintf("  %7.3f  ",yhat[i][j]);
          Vprintf("\n");
        }
      Vprintf("\n%s%10.3f\n","RSS = ",rss);
      return EXIT_SUCCESS;
    }
  else
    {
      Vprintf("Incorrect input value of n or m.\n");
      return EXIT_FAILURE;
    }
}