/* nag_fft_multiple_complex(c06frc) Example Program
 *
 * Copyright 1990 Numerical Algorithms Group.
 *
 * Mark 1, 1990.
 *
 * Mark 3 revised, 1994.
 */

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

#define MMAX  5
#define NMAX 20

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

  /* Skip heading in data file */
  Vscanf("%*[^\n]");
  Vprintf("c06frc Example Program Results\n");
  while (scanf("%ld%ld", &m, &n)!=EOF)
    if (m<=MMAX && n<=NMAX)
      {
        Vprintf("\n\nm = %2ld  n = %2ld\n", m, n);
        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");
          }
        /* Initialise trig array */
        c06gzc(n, trig, NAGERR_DEFAULT);
        /* Compute transforms */
        c06frc(m, n, x, y, trig, NAGERR_DEFAULT);
        Vprintf("\nDiscrete 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 transforms */
        c06gcc(m*n, y, NAGERR_DEFAULT);
        c06frc(m, n, x, y, trig, 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
      {
        Vprintf("\nInvalid value of m or n.\n");
        return EXIT_FAILURE;
      }
  return EXIT_SUCCESS;
}