/* nag_fft_init_trig(c06gzc) 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("c06gzc Example Program Results\n");
  Vscanf("%*[^\n]");   /* Skip heading in data file */
  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");
          }
        c06gzc(n, trig, NAGERR_DEFAULT);  /* Initialize trig array */
        /* Calculate transform */
        c06fpc(m, n, x, trig, NAGERR_DEFAULT);
        Vprintf("\nDiscrete Fourier transforms in Hermitian format\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");
          }
        /* Convert Hermitian form to full complex */
        c06gsc(m, n, x, u, v, NAGERR_DEFAULT);
        Vprintf("\nFourier transforms 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");
          }
      }
    else
      {
        Vfprintf(stderr,"\nInvalid value of m or n.\n");
        return EXIT_FAILURE;
      }
  return EXIT_SUCCESS;
}