/* nag_fft_multiple_hermitian(c06fqc) Example Program
 *
 * Copyright 1990 Numerical Algorithms Group.
 *
 * Mark 1, 1990.
 */

#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 u[MMAX*NMAX], v[MMAX*NMAX];
  double x[MMAX*NMAX];

  Vprintf("c06fqc Example Program Results\n");
  /* Skip heading in data file */
  Vscanf("%*[^\n]");
  while (scanf("%ld%ld", &m, &n)!=EOF)
    if (m<=MMAX && n<=NMAX)
      {
        Vprintf("\n\nm = %2ld  n = %2ld\n", m, n);
        /* Read in data and print out. */
        for (j = 0; j<m; ++j)
          for (i = 0; i<n; ++i)
            Vscanf("%lf", &x[j*n + i]);
        Vprintf("\nOriginal data values\n\n");
        for (j = 0; j<m; ++j)
          {
            Vprintf("    ");
            for (i = 0; i<n; ++i)
              Vprintf("%10.4f%s", x[j*n + i],
                      (i%6==5 && i!=n-1 ? "\n     " : ""));
            Vprintf("\n");
          }
        /* Calculate full complex form of Hermitian data sequences */
        c06gsc(m, n, x, u, v, NAGERR_DEFAULT);
        Vprintf("\nOriginal data written in full complex form\n\n");
        for (j = 0; j<m; ++j)
          {
            Vprintf("Real");
            for (i = 0; i<n; ++i)
              Vprintf("%10.4f%s", u[j*n + i],
                      (i%6==5 && i!=n-1 ? "\n     " : ""));
            Vprintf("\nImag");
            for (i = 0; i<n; ++i)
              Vprintf("%10.4f%s", v[j*n + i],
                      (i%6==5 && i!=n-1 ? "\n     " : ""));
            Vprintf("\n\n");
          }
        /* Initialise trig array */
        c06gzc(n, trig, NAGERR_DEFAULT);
        /* Calculate transforms */
        c06fqc(m, n, x, trig, NAGERR_DEFAULT);
        Vprintf("\nDiscrete Fourier transforms (real values)\n\n");
        for (j = 0; j<m; ++j)
          {
            Vprintf("    ");
            for (i = 0; i<n; ++i)
              Vprintf("%10.4f%s", x[j*n + i],
                      (i%6==5 && i!=n-1 ? "\n     " : ""));
            Vprintf("\n");
          }
        /* Calculate inverse transforms */
        c06fpc(m, n, x, trig, NAGERR_DEFAULT);
        c06gqc(m, n, x, NAGERR_DEFAULT);
        Vprintf("\nOriginal data as restored by inverse transform\n\n");
        for (j = 0; j<m; ++j)
          {
            Vprintf("    ");
            for (i = 0; i<n; ++i)
              Vprintf("%10.4f%s", x[j*n + i],
                      (i%6==5 && i!=n-1 ? "\n     " : ""));
            Vprintf("\n");
          }
      }
    else
      Vfprintf(stderr,"\nInvalid value of m or n.\n");
  return EXIT_SUCCESS;
}