/* nag_mldwt_2d (c09ecc) Example Program.
 *
 * Copyright 2011, Numerical Algorithms Group.
 *
 * Mark 23, 2011.
 */

#include <string.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagc09.h>
#include <nagx04.h>

int main(void)
{
  /* Scalars */
  Integer         exit_status = 0;
  Integer         dim1, dim2, i, ilevel, itype_coeffs, j, jstart, lenc, m, n,
                  nf, nwcn, nwct, nwl, pda, pdb;
  /* Arrays */
  char            mode[24], wavnam[20], title[27];
  double          *a = 0, *b = 0, *c = 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_mldwt_2d (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("%s%s%*[^\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);

  /*
   * 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_gen_real_mat_print_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_gen_real_mat_print_comp (x04cbc).\n%s\n",
	     fail.message);
      exit_status = 1;
      goto END;
    }

  /* nag_wfilt_2d (c09abc).
   *  Two-dimensional wavelet filter initialization.
   */
  nag_wfilt_2d(wavnamenum, Nag_MultiLevel, modenum, m, n, &nwl, &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((nwl), Integer)) ||
    !(dwtlvn = NAG_ALLOC((nwl), Integer))
    )
    {
      printf("Allocation failure\n");
      exit_status = -2;
      goto END;
    }

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

  /* Print decomposition */
  printf("\n\n Number of Levels : %ld\n", nwl);
  printf(" Number of coefficients in 1st dimension for each level :\n");
  for (j = 0; j < nwl; j++) 
    printf("%8ld%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("%8ld%s", dwtlvn[j], (j+1)%8?" ":"\n");
  jstart = 0;
  printf("\n\nWavelet coefficients C : \n");
  for (ilevel = nwl; ilevel > 0; ilevel -= 1)
    {
      dim1 = dwtlvm[nwl - ilevel];
      dim2 = dwtlvn[nwl - ilevel];
      for (j = 0; j < 55; j++) printf("-");
      printf("\n Level : %ld; output is %ld by %ld\n",
	     ilevel, dim1, dim2);
      for (j = 0; j < 55; j++) printf("-");
      printf("\n");
      for (itype_coeffs = 1; itype_coeffs <= 4; itype_coeffs++)
        {
          switch (itype_coeffs)
            {
            case 1:
              if (ilevel == nwl) strcpy(title, "Approximation coefficients  ");
              break;
            case 2:
              strcpy(title, "Vertical coefficients       ");
              break;
            case 3:
	      strcpy(title, "Horizontal coefficients     ");
              break;
            case 4:
	      strcpy(title, "Diagonal coefficients       ");
            }
          if (itype_coeffs > 1 || ilevel == nwl)
            {
	      nag_gen_real_mat_print_comp(order, matrix, diag, dim1, dim2,
					  &c[jstart], dim1, "%8.4f", title,
					  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 = 4;
		  goto END;
		}
              jstart = jstart + dim1 * dim2;
            }
        }
    }

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

  /* Print reconstruction */
  printf("\n");
  strcpy(title, "Reconstruction           B :");
  nag_gen_real_mat_print_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_gen_real_mat_print_comp (x04cbc).\n%s\n",
	     fail.message);
      exit_status = 6;
      goto END;
    }

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