/* nag_ftt_2d_complex(c06fuc) Example Program
 *
 * Copyright 1990 Numerical Algorithms Group.
 *
 * Mark 2 revised, 1992.
 */

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

#define MMAX 6
#define NMAX 6
#define MNMAX MMAX*NMAX

int main(void)
{
  Integer i, j, m, n;
  double  trigm[2*MMAX], trign[2*NMAX];
  double  x[MNMAX], y[MNMAX];

  Vprintf("c06fuc Example Program Results\n");
  /* Skip heading in data file */
  Vscanf("%*[^\n]");
  while (scanf("%ld%ld", &m, &n)!=EOF)
    if (m*n>=1 && m*n<=MNMAX)
      {
        Vprintf("\n\nm = %2ld  n = %2ld\n", m, n);
        /* Read in complex data and print out. */
        for (j = 0; j<m; ++j)
          {
            for (i = 0; i<n; ++i)
              Vscanf("%lf", &x[j*n + i]);
            for (i = 0; i<n; ++i)
              Vscanf("%lf", &y[j*n + i]);
          }
        Vprintf("\nOriginal data values\n\n");
        for (j = 0; j<m; ++j)
          {
            Vprintf("Real");
            for (i = 0; i<n; ++i)
              Vprintf("%10.4f%s", x[j*n + i],
                      (i%6==5 && i!=n-1 ? "\n     " : ""));
            Vprintf("\nImag");
            for (i = 0; i<n; ++i)
              Vprintf("%10.4f%s", y[j*n + i],
                      (i%6==5 && i!=n-1 ? "\n     " : ""));
            Vprintf("\n\n");
          }
        /* Initialize trig arrays */
        c06gzc(m, trigm, NAGERR_DEFAULT);
        c06gzc(n, trign, NAGERR_DEFAULT);
        /* Compute transform */
        c06fuc(m, n, x, y, trigm, trign, NAGERR_DEFAULT);
        Vprintf("\nComponents of discrete Fourier transforms\n\n");
        for (j = 0; j<m; ++j)
          {
            Vprintf("Real");
            for (i = 0; i<n; ++i)
              Vprintf("%10.4f%s", x[j*n + i],
                      (i%6==5 && i!=n-1 ? "\n     " : ""));
            Vprintf("\nImag");
            for (i = 0; i<n; ++i)
              Vprintf("%10.4f%s", y[j*n + i],
                      (i%6==5 && i!=n-1 ? "\n     " : ""));
            Vprintf("\n\n");
          }
        /* Compute inverse transform */
        c06gcc(m*n, y, NAGERR_DEFAULT);
        c06fuc(m, n, x, y, trigm, trign, NAGERR_DEFAULT);
        c06gcc(m*n, y, NAGERR_DEFAULT);
        Vprintf("\nOriginal data as restored by inverse transform\n\n");
        for (j = 0; j<m; ++j)
          {
            Vprintf("Real");
            for (i = 0; i<n; ++i)
              Vprintf("%10.4f%s", x[j*n + i],
                      (i%6==5 && i!=n-1 ? "\n     " : ""));
            Vprintf("\nImag");
            for (i = 0; i<n; ++i)
              Vprintf("%10.4f%s", y[j*n + i],
                      (i%6==5 && i!=n-1 ? "\n     " : ""));
            Vprintf("\n\n");
          }
      }
    else
      Vfprintf(stderr,"\nInvalid value of m or n.\n");
  return EXIT_SUCCESS;
}