/* nag_dwt_3d (c09fac) Example Program.
 *
 * Copyright 2017 Numerical Algorithms Group.
 *
 * Mark 26.1, 2017.
 */
#include <stdio.h>
#include <math.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagc09.h>

#define A(I,J,K) a[I-1 + (J-1)* lda + (K-1)* lda * sda]
#define B(I,J,K) b[I-1 + (J-1)* ldb + (K-1)* ldb * sdb]
#define D(I,J,K) d[I-1 + (J-1)* nwcm + (K-1)* nwcm * nwcn]

int main(void)
{
  /* Scalars */
  Integer exit_status = 0, zero = 0;
  Integer cindex, i, j, k, lda, ldb, lenc;
  Integer m, n, fr, nf, nwcfr, nwcm, nwcn, nwct, nwl, sda, sdb;
  /* Arrays */
  char mode[25], wavnam[25];
  double *a = 0, *b = 0, *c = 0, *d = 0;
  Integer icomm[260];
  /* Nag Types */
  Nag_Wavelet wavnamenum;
  Nag_WaveletMode modenum;
  Nag_MatrixType matrix = Nag_GeneralMatrix;
  Nag_OrderType order = Nag_ColMajor;
  Nag_DiagType diag = Nag_NonUnitDiag;
  NagError fail;

  INIT_FAIL(fail);

  printf("nag_dwt_3d (c09fac) Example Program Results\n\n");
  fflush(stdout);

  /* Skip heading in data file and read problem parameters */
  scanf("%*[^\n] %" NAG_IFMT "%" NAG_IFMT "%" NAG_IFMT "%*[^\n]", &m, &n,
        &fr);
  lda = m;
  ldb = m;
  sda = n;
  sdb = n;
  scanf("%24s%24s%*[^\n]\n", wavnam, mode);

  if (!(a = NAG_ALLOC((lda) * (sda) * (fr), double)) ||
      !(b = NAG_ALLOC((ldb) * (sdb) * (fr), double)))
  {
    printf("Allocation failure\n");
    exit_status = 1;
    goto END;
  }

  printf("Parameters read from file :: \n");
  printf("DWT :: Wavelet  : %s\n", wavnam);
  printf("       End mode : %s\n", mode);
  printf("       m  : %4" NAG_IFMT "\n", m);
  printf("       n  : %4" NAG_IFMT "\n", n);
  printf("       fr : %4" NAG_IFMT "\n\n", fr);

  /* nag_enum_name_to_value (x04nac).
   * Converts NAG enum member name to value
   */
  wavnamenum = (Nag_Wavelet) nag_enum_name_to_value(wavnam);
  modenum = (Nag_WaveletMode) nag_enum_name_to_value(mode);

  /* Read data array */
  for (k = 1; k <= fr; k++) {
    for (i = 1; i <= m; i++) {
      for (j = 1; j <= n; j++)
        scanf("%lf", &A(i, j, k));
    }
    scanf("%*[^\n] ");
  }

  /* Print out the input data */
  printf("Input Data :\n");
  fflush(stdout);
  for (k = 1; k <= fr; k++) {
    /* nag_gen_real_mat_print_comp (x04cbc).
     * Prints out a matrix.
     */
    nag_gen_real_mat_print_comp(order, matrix, diag, m, n, &A(1, 1, k), lda,
                                "%8.4f", " ", Nag_NoLabels, 0, Nag_NoLabels,
                                0, 80, 0, 0, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_gen_real_mat_print_comp (x04cbc).\n%s\n",
             fail.message);
      exit_status = 2;
      goto END;
    }
    printf("\n");
    fflush(stdout);
  }

  /* nag_wfilt_3d (c09acc).
   * Three-dimensional wavelet filter initialization
   */
  nag_wfilt_3d(wavnamenum, Nag_SingleLevel, modenum, m, n, fr, &nwl, &nf,
               &nwct, &nwcn, &nwcfr, icomm, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_wfilt_3d (c09acc).\n%s\n", fail.message);
    exit_status = 3;
    goto END;
  }

  /* Calculate the number of wavelet coefficients in 
   * the first dimension, nwcm.
   */
  nwcm = nwct / (8 * nwcn * nwcfr);
  lenc = nwct;

  /* Allocate space for the coefficients array, C */
  if (!(c = NAG_ALLOC((lenc), double)))
  {
    printf("Allocation failure\n");
    exit_status = 4;
    goto END;
  }

  /* nag_dwt_3d (c09fac).
   * Three-dimensional discrete wavelet transform
   */
  nag_dwt_3d(m, n, fr, a, lda, sda, lenc, c, icomm, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_dwt_3d (c09fac).\n%s\n", fail.message);
    exit_status = 5;
    goto END;
  }

  /* Allocate space for extraction of coefficients of a single type */
  if (!(d = NAG_ALLOC((nwcm) * (nwcn) * (nwcfr), double)))
  {
    printf("Allocation failure\n");
    exit_status = 6;
    goto END;
  }

  for (cindex = 0; cindex <= 7; cindex++) {
    /* Use the extraction routine c09fyc to retrieve the required
     * coefficients.
     */

    /* nag_wav_3d_coeff_ext (c09fyc).
     * Extract the nominated coefficients.
     */
    nag_wav_3d_coeff_ext(zero, cindex, lenc, c, d, nwcm, nwcn, icomm, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_wav_3d_coeff_ext (c09fyc).\n%s\n", fail.message);
      exit_status = 7;
      goto END;
    }

    /* Print out the extracted coefficients */
    switch (cindex) {
    case 0:
      printf("Approximation coefficients (LLL)\n");
      break;
    case 1:
      printf("Detail coefficients (LLH)\n");
      break;
    case 2:
      printf("Detail coefficients (LHL)\n");
      break;
    case 3:
      printf("Detail coefficients (LHH)\n");
      break;
    case 4:
      printf("Detail coefficients (HLL)\n");
      break;
    case 5:
      printf("Detail coefficients (HLH)\n");
      break;
    case 6:
      printf("Detail coefficients (HHL)\n");
      break;
    case 7:
      printf("Detail coefficients (HHH)\n");
      break;
    }

    for (i = 1; i <= nwcm; i++) {
      if (i == 1) {
        printf("Coefficients   ");
        for (k = 1; k <= nwcfr; k++) {
          printf("Frame %4" NAG_IFMT, k);
          for (j = 1; j <= 9 * nwcn - 8; j++)
            printf(" ");
        }
        printf("\n");
      }

      for (k = 1; k <= nwcfr; k++) {
        if (k == 1 && i == 1)
          printf("%5" NAG_IFMT "%8s", cindex, " ");
        else if (k == 1)
          printf("%13s", " ");
        else
          printf("%2s", " ");
        for (j = 1; j <= nwcn; j++) {
          printf("%8.4f ", D(i, j, k));
        }
      }
      printf("\n");
    }
    printf("\n");
  }
  fflush(stdout);

  /* nag_idwt_3d (c09fbc).
   * Three-dimensional inverse discrete wavelet transform
   */
  nag_idwt_3d(m, n, fr, lenc, c, b, ldb, sdb, icomm, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_idwt_3d (c09fbc).\n%s\n", fail.message);
    exit_status = 8;
    goto END;
  }

  printf("Output Data :\n");
  fflush(stdout);
  for (k = 1; k <= fr; k++) {
    /* nag_gen_real_mat_print_comp (x04cbc).
     * Prints out a matrix.
     */
    nag_gen_real_mat_print_comp(order, matrix, diag, m, n, &B(1, 1, k), ldb,
                                "%8.4f", " ", Nag_NoLabels, 0, Nag_NoLabels,
                                0, 80, 0, 0, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_gen_real_mat_print_comp (x04cbc).\n%s\n",
             fail.message);
      exit_status = 9;
      goto END;
    }
    printf("\n");
    fflush(stdout);
  }

END:
  NAG_FREE(a);
  NAG_FREE(b);
  NAG_FREE(c);
  NAG_FREE(d);

  return exit_status;
}