/* nag_wav_dim2_multi_fwd (c09ecc) Example Program.
 *
 * Copyright 2022 Numerical Algorithms Group.
 *
 * Mark 28.3, 2022.
 */

#include <nag.h>
#include <string.h>

int main(void) {
  /* Scalars */
  Integer exit_status = 0;
  Integer nwcm, i, ilevel, itype_coeffs, j, lenc, m, n, nf, nwcn, nwct, nwlmax,
      nwl, nwlinv, pda, pdb;
  /* Arrays */
  char mode[24], wavnam[20], title[50];
  double *a = 0, *b = 0, *c = 0, *d = 0;
  Integer *dwtlvm = 0, *dwtlvn = 0;
  Integer icomm[180];
  /* 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);

  /* Output preamble */
  printf("nag_wav_dim2_multi_fwd (c09ecc) Example Program Results\n\n");

  /* Skip heading in data file and read problem parameters */
  scanf("%*[^\n] %" NAG_IFMT "%" NAG_IFMT "%*[^\n] ", &m, &n);
  scanf("%19s%23s%*[^\n] ", wavnam, mode);
  pda = m;
  pdb = m;
  if (!(a = NAG_ALLOC(pda * n, double)) || !(b = NAG_ALLOC(pdb * n, double))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }
  printf(" Parameters read from file :: \n");
  printf(" MLDWT :: Wavelet  : %s\n", wavnam);
  printf("          End mode : %s\n", mode);
  printf("          m        : %" NAG_IFMT "\n", m);
  printf("          n        : %" NAG_IFMT "\n\n", n);
  fflush(stdout);

  /* 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 and write it out */
#define A(I, J) a[(J - 1) * pda + I - 1]
  for (i = 1; i <= m; i++)
    for (j = 1; j <= n; j++)
      scanf("%lf", &A(i, j));

  nag_file_print_matrix_real_gen_comp(
      order, matrix, diag, m, n, a, pda, "%8.4f",
      "Input Data     A :", Nag_NoLabels, 0, Nag_NoLabels, 0, 80, 0, 0, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_file_print_matrix_real_gen_comp (x04cbc).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }

  /* nag_wav_dim2_init (c09abc).
   *  Two-dimensional wavelet filter initialization.
   */
  nag_wav_dim2_init(wavnamenum, Nag_MultiLevel, modenum, m, n, &nwlmax, &nf,
                    &nwct, &nwcn, icomm, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_nfilt_2d (c09abc).\n%s\n", fail.message);
    exit_status = 2;
    goto END;
  }
  lenc = nwct;
  if (!(c = NAG_ALLOC(lenc, double)) ||
      !(dwtlvm = NAG_ALLOC(nwlmax, Integer)) ||
      !(dwtlvn = NAG_ALLOC(nwlmax, Integer))) {
    printf("Allocation failure\n");
    exit_status = -2;
    goto END;
  }

  nwl = nwlmax;

  /* nag_wav_dim2_multi_fwd (c09ecc).
   * Two-dimensional multi-level discrete wavelet transform
   */
  nag_wav_dim2_multi_fwd(m, n, a, pda, lenc, c, nwl, dwtlvm, dwtlvn, icomm,
                         &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_wav_dim2_multi_fwd (c09ecc).\n%s\n", fail.message);
    exit_status = 3;
    goto END;
  }

  /* Print decomposition */
  printf("\n Number of Levels : %" NAG_IFMT "\n", nwl);
  printf(" Number of coefficients in 1st dimension for each level :\n");
  for (j = 0; j < nwl; j++)
    printf("%8" NAG_IFMT "%s", dwtlvm[j], (j + 1) % 8 ? " " : "\n");

  printf("\n Number of coefficients in 2nd dimension for each level :\n");
  for (j = 0; j < nwl; j++)
    printf("%8" NAG_IFMT "%s", dwtlvn[j], (j + 1) % 8 ? " " : "\n");
  printf("\n\nWavelet coefficients C : \n");
  for (ilevel = nwl; ilevel > 0; ilevel -= 1) {
    nwcm = dwtlvm[nwl - ilevel];
    nwcn = dwtlvn[nwl - ilevel];
    if (!(d = NAG_ALLOC(nwcm * nwcn, double))) {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }
    for (j = 0; j < 55; j++)
      printf("-");
    printf("\n Level : %" NAG_IFMT "; output is %" NAG_IFMT " by %" NAG_IFMT
           "\n",
           ilevel, nwcm, nwcn);
    for (j = 0; j < 55; j++)
      printf("-");
    printf("\n");
    fflush(stdout);
    for (itype_coeffs = 0; itype_coeffs <= 3; itype_coeffs++) {
      switch (itype_coeffs) {
      case 0:
        if (ilevel == nwl)
          strcpy(title, "Approximation coefficients  ");
        break;
      case 1:
        strcpy(title, "Vertical coefficients       ");
        break;
      case 2:
        strcpy(title, "Horizontal coefficients     ");
        break;
      case 3:
        strcpy(title, "Diagonal coefficients       ");
      }
      if (itype_coeffs > 0 || ilevel == nwl) {
        /* nag_wav_dim2_coeff_ext (c09aec).
         * Call the 2D extraction routine c09eac
         */
        nag_wav_dim2_coeff_ext(ilevel, itype_coeffs, lenc, c, d, nwcm, icomm,
                               &fail);
        nag_file_print_matrix_real_gen_comp(
            order, matrix, diag, nwcm, nwcn, d, nwcm, "%8.4f", title,
            Nag_NoLabels, 0, Nag_NoLabels, 0, 80, 0, 0, &fail);
        if (fail.code != NE_NOERROR) {
          printf("Error from nag_file_print_matrix_real_gen_comp (x04cbc)."
                 "\n%s\n",
                 fail.message);
          exit_status = 4;
          goto END;
        }
      }
    }
    NAG_FREE(d);
  }

  nwlinv = nwl;

  /* nag_wav_dim2_multi_inv (c09edc).
   * Two-dimensional inverse multi-level discrete wavelet transform
   */
  nag_wav_dim2_multi_inv(nwlinv, lenc, c, m, n, b, pdb, icomm, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_wav_dim2_multi_inv (c09edc).\n%s\n", fail.message);
    exit_status = 5;
    goto END;
  }

  /* Print reconstruction */
  printf("\n");
  fflush(stdout);
  strcpy(title, "Reconstruction           B :");
  nag_file_print_matrix_real_gen_comp(order, matrix, diag, m, n, b, pdb,
                                      "%8.4f", title, Nag_NoLabels, 0,
                                      Nag_NoLabels, 0, 80, 0, 0, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_file_print_matrix_real_gen_comp (x04cbc).\n%s\n",
           fail.message);
    exit_status = 6;
    goto END;
  }

END:
  NAG_FREE(a);
  NAG_FREE(b);
  NAG_FREE(c);
  NAG_FREE(d);
  NAG_FREE(dwtlvm);
  NAG_FREE(dwtlvn);
  return exit_status;
}