/* nag_sparse_sym_rcm (f11yec) Example Program.
 *
 * Copyright 2023 Numerical Algorithms Group.
 *
 * Mark 29.3, 2023.
 */
#include <nag.h>
#include <stdio.h>

void plot(const Integer n, const Integer nnz, Integer *perm, Integer *icolzp,
          Integer *irowix);
void uncompress(Integer n, Integer *icolzp, Integer *icol);

int main(void) {
  /* Scalars */
  Integer n, nnz, exit_status = 0, doplot = 0, i;
  /* Arrays */
  Integer *icolzp = 0, *irowix = 0, *mask = 0, *perm = 0;
  Integer info[4];
  Nag_Boolean lopts[5];
  /* Nag Types */
  NagError fail;

  INIT_FAIL(fail);
  printf("nag_sparse_sym_rcm (f11yec) Example Program Results\n");
  /* Skip heading in data file and
   * read Size of the matrix and Number of nonzero elements
   */
  scanf("%*[^\n] ");
  scanf("%" NAG_IFMT "%*[^\n] ", &n);
  scanf("%" NAG_IFMT "%*[^\n] ", &nnz);
  if (!(icolzp = NAG_ALLOC((n + 1), Integer)) ||
      !(irowix = NAG_ALLOC((nnz), Integer)) ||
      !(mask = NAG_ALLOC((n), Integer)) || !(perm = NAG_ALLOC((n), Integer))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }
  /* Read in data */
  for (i = 0; i < nnz; i += 1)
    scanf("%" NAG_IFMT "", &irowix[i]);
  scanf("%*[^\n] ");
  for (i = 0; i < (n + 1); i += 1)
    scanf("%" NAG_IFMT "", &icolzp[i]);
  scanf("%*[^\n] ");
  /* Set options */
  lopts[0 /* Use Mask                 */] = Nag_FALSE;
  lopts[1 /* Don't reverse            */] = Nag_FALSE;
  lopts[2 /* Check symmetry           */] = Nag_TRUE;
  lopts[3 /* Compute bandwidth before */] = Nag_TRUE;
  lopts[4 /* Compute bandwidth after  */] = Nag_TRUE;
  /* nag_sparse_sym_rcm (f11yec).
   * Reverse Cuthill-McKee reordering of a real sparse symmetric
   * matrix in CCS format
   */
  nag_sparse_sym_rcm(n, nnz, icolzp, irowix, lopts, mask, perm, info, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_sparse_sym_rcm (f11yec).\n%s\n", fail.message);
    exit_status = 1;
    goto END;
  }
  /* Print results */
  printf("Permutation (perm):\n");
  for (i = 0; i < n; i += 1) {
    printf("    %3" NAG_IFMT "", perm[i]);
    if (i % 6 == 5)
      printf("\n");
  }
  printf("\n\nStatistics:\n");
  printf(" %s%6" NAG_IFMT "\n", " Before:  Bandwidth = ", info[0]);
  printf(" %s%6" NAG_IFMT "\n", " Before:  Profile   = ", info[1]);
  printf(" %s%6" NAG_IFMT "\n", " After :  Bandwidth = ", info[2]);
  printf(" %s%6" NAG_IFMT "\n", " After :  Profile   = ", info[3]);
  /* Print matrix entries and permuted entries in form suitable
   * for plotting
   */
  if (doplot)
    plot(n, nnz, perm, icolzp, irowix);
END:
  NAG_FREE(icolzp);
  NAG_FREE(irowix);
  NAG_FREE(mask);
  NAG_FREE(perm);

  return exit_status;
}

void uncompress(Integer n, Integer *icolzp, Integer *icol) {
  Integer i, j, col_beg, col_end;
  for (i = 0; i < n; i++) {
    col_end = icolzp[i + 1] - 1;
    col_beg = icolzp[i];
    for (j = col_beg; j <= col_end; j++)
      icol[j - 1] = i + 1;
  }
}

void plot(const Integer n, const Integer nnz, Integer *perm, Integer *icolzp,
          Integer *irowix) {
  /* Put data, suitable for plotting matrix structure, in data file */

  /* Scalars */
  Integer i, nnz2;
  /* Arrays */
  double *a = 0;
  Integer *icolix = 0, *ipcolix = 0, *iperm = 0, *iprowix = 0, *istr = 0;
  /* Nag Types */
  NagError fail;

  INIT_FAIL(fail);
  if (!(icolix = NAG_ALLOC(nnz, Integer)) ||
      !(ipcolix = NAG_ALLOC(nnz, Integer)) ||
      !(iprowix = NAG_ALLOC(nnz, Integer)) ||
      !(iperm = NAG_ALLOC(n, Integer)) || !(a = NAG_ALLOC(nnz, double)) ||
      !(istr = NAG_ALLOC(n + 1, Integer))) {
    printf("Allocation failure\n");
    return;
  }
  /* Decompress icolzp to full set of column indices (icolix)
   * and compute inverse permutation
   */
  uncompress(n, icolzp, icolix);
  for (i = 0; i < n; i++) {
    iperm[perm[i] - 1] = i + 1;
  }
  /* Original matrix structure */
  for (i = 0; i < nnz; i++) {
    a[i] = icolix[i] * .01 + 1.0 * irowix[i];
    printf("%8" NAG_IFMT "    ", irowix[i]);
    printf("%8" NAG_IFMT "    ", icolix[i]);
    printf("%8.2f\n", a[i]);
  }
  printf("\n");
  /* Apply Inverse Permutation */
  for (i = 0; i < nnz; i++) {
    ipcolix[i] = iperm[icolix[i] - 1];
    iprowix[i] = iperm[irowix[i] - 1];
  }
  /* Reorder (in exit: istr contains new CCS icolzp) */
  nnz2 = nnz;
  nag_sparse_real_gen_sort(n, &nnz2, a, ipcolix, iprowix,
                           Nag_SparseNsym_FailDups, Nag_SparseNsym_KeepZeros,
                           istr, &fail);
  /* Permuted matrix structure */
  for (i = 0; i < nnz2; i++) {
    printf("%8" NAG_IFMT "    ", iprowix[i]);
    printf("%8" NAG_IFMT "    ", ipcolix[i]);
    printf("%8.2f\n", a[i]);
  }
  NAG_FREE(icolix);
  NAG_FREE(ipcolix);
  NAG_FREE(iprowix);
  NAG_FREE(iperm);
  NAG_FREE(a);
  NAG_FREE(istr);
}