/* nag_fft_multiple_qtr_cosine(c06hdc) Example Program
 *
 * Copyright 1991 Numerical Algorithms Group.
 *
 * Mark 2, 1991.
 */

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

#define MMAX 5
#define NMAX 20
#define X(I,J) x[(I)*n + (J)]

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

  Vprintf("c06hdc Example Program Results\n");
  Vscanf(" %*[^\n]");    /* Skip heading in data file */
  while (scanf("%ld %ld", &m, &n) != EOF)
    if (m <= MMAX && n <= NMAX)
      {
        Vscanf(" %*[^\n]"); /* Skip text in data file */
        Vscanf(" %*[^\n]");
        for (i = 0; i < m; ++i)
          for (j = 0; j < n; ++j)
            Vscanf("%lf", &X(i,j));
        Vprintf("\nOriginal data values\n\n");
        for (i = 0; i < m; ++i)
          {
            for (j = 0; j < n; ++j)
              Vprintf(" %10.4f%s", X(i,j),
                      (j%7==6 && j!=n-1 ? "\n" : ""));
            Vprintf("\n");
          }
        /* Initialise trig array */
        c06gzc(n, trig, NAGERR_DEFAULT);
        /* Compute transform */
        c06hdc(Nag_ForwardTransform, m, n, x, trig, NAGERR_DEFAULT);

        Vprintf("\nDiscrete quarter-wave Fourier cosine transforms\n\n");
        for (i = 0; i < m; ++i)
          {
            for (j = 0; j < n; ++j)
              Vprintf(" %10.4f%s", X(i,j),
                      (j%7==6 && j!=n-1 ? "\n" : ""));
            Vprintf("\n");
          }
        /* Compute inverse transform */
        c06hdc(Nag_BackwardTransform, m, n, x, trig, NAGERR_DEFAULT);
        Vprintf("\nOriginal data as restored by inverse transform\n\n");
        for (i = 0; i < m; ++i)
          {
            for (j = 0; j < n; ++j)
              Vprintf(" %10.4f%s", X(i,j),
                      (j%7==6 && j!=n-1 ? "\n" : ""));
            Vprintf("\n");
          }
      }
    else Vfprintf(stderr,"\nInvalid value of m or n.\n");
  return EXIT_SUCCESS;
}