/* nag_fft_multiple_qtr_cosine (c06hdc) Example Program.

 *

 * Copyright 1991 Numerical Algorithms Group.

 *

 * Mark 2, 1991.

 * Mark 8 revised, 2004.

 */



#include <nag.h>

#include <stdio.h>

#include <nag_stdlib.h>

#include <nagc06.h>



#define X(I,J) x[(I)*n + (J)]



int main(void)

{

  Integer exit_status=0, i, j, m, n;

  NagError fail;

  double *trig=0, *x=0;



  INIT_FAIL(fail);

  Vprintf("nag_fft_multiple_qtr_cosine (c06hdc) Example Program Results\n");

  Vscanf(" %*[^\n]");    /* Skip heading in data file */

  while (scanf("%ld %ld", &m, &n) != EOF)

    {

      if (m>=1 && n>=1)

        {

          if ( !( trig = NAG_ALLOC(2*n, double)) ||

               !( x = NAG_ALLOC(m*n, double)) )

            {

              Vprintf("Allocation failure\n");

              exit_status = -1;

              goto END;

            }

        }

      else

        {

          Vprintf("Invalid m or n.\n");

          exit_status = 1;

          return exit_status;

        }

      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 */

      /* nag_fft_init_trig (c06gzc).

       * Initialization function for other c06 functions

       */

      nag_fft_init_trig(n, trig, &fail);

      if (fail.code != NE_NOERROR)

        {

          Vprintf("Error from nag_fft_init_trig (c06gzc).\n%s\n", fail.message);

          exit_status = 1;

          goto END;

        }



      /* Compute transform */

      /* nag_fft_multiple_qtr_cosine (c06hdc).

       * Discrete quarter-wave cosine transform

       */

      nag_fft_multiple_qtr_cosine(Nag_ForwardTransform, m, n, x, trig, &fail);

      if (fail.code != NE_NOERROR)

        {

          Vprintf("Error from nag_fft_multiple_qtr_cosine (c06hdc).\n%s\n",

                  fail.message);

          exit_status = 1;

          goto END;

        }





      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 */

      /* nag_fft_multiple_qtr_cosine (c06hdc), see above. */

      nag_fft_multiple_qtr_cosine(Nag_BackwardTransform, m, n, x, trig, &fail);

      if (fail.code != NE_NOERROR)

        {

          Vprintf("Error from nag_fft_multiple_qtr_cosine (c06hdc).\n%s\n",

                  fail.message);

          exit_status = 1;

          goto END;

        }



      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");

        }

    }

 END:

  if (trig) NAG_FREE(trig);

  if (x) NAG_FREE(x);

  return exit_status;

}