/* nag_fft_multiple_sine(c06hac) 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)*row_len + (J)]

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

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

        Vprintf("\nDiscrete Fourier sine transforms\n\n");
        for (i = 0; i < m; ++i)
          {
            for (j = 0; j < row_len; ++j)
              Vprintf(" %10.4f%s", X(i,j),
                      (j%7==6 && j!=row_len-1 ? "\n" : ""));
            Vprintf("\n");
          }

        c06hac(m, n, x, trig, NAGERR_DEFAULT); /* Compute inverse transform */

        Vprintf("\nOriginal data as restored by inverse transform\n\n");
        for (i = 0; i < m; ++i)
          {
            for (j = 0; j < row_len; ++j)
              Vprintf(" %10.4f%s", X(i,j),
                      (j%7==6 && j!=row_len-1 ? "\n" : ""));
            Vprintf("\n");
          }
      }
    else Vfprintf(stderr,"\nInvalid value of m or n.\n");
  return EXIT_SUCCESS;
}