/* nag_sparse_nherm_precon_bdilu (f11dtc) Example Program.
 *
 * Copyright 2017 Numerical Algorithms Group.
 *
 * Mark 26.1, 2017.
 */

#include <nag.h>
#include <nagf11.h>
#include <nag_stdlib.h>

static void overlap(Integer *n, Integer *nnz, Integer *irow, Integer *icol,
                    Integer *nb, Integer *istb, Integer *indb, Integer *lindb,
                    Integer *nover, Integer *iwork);

int main(void)
{

  /* Scalars */
  double dtolg;
  Integer i, j, k, la, lfillg, lindb, liwork, minval, mb, n, nb, nnz, nnzc,
         nover;
  Integer exit_status = 0, maxval_ret = 9999;
  Nag_SparseNsym_Piv pstrag;
  Nag_SparseNsym_Fact milug;

  /* Arrays */
  char nag_enum_arg[40];
  double *dtol = 0;
  Complex *a = 0;
  Integer *icol = 0, *idiag = 0, *indb = 0, *ipivp = 0, *ipivq = 0, *irow = 0;
  Integer *istb = 0, *istr = 0, *iwork = 0, *lfill = 0, *npivm = 0;
  Nag_SparseNsym_Piv *pstrat;
  Nag_SparseNsym_Fact *milu;

  /* Nag Types */
  NagError fail;

  /* Print example header */
  printf("nag_sparse_nherm_precon_bdilu (f11dtc) Example Program Results\n\n");

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

  /* Get the matrix order and number of nonzero entries. */
  scanf("%" NAG_IFMT " %*[^\n]", &n);
  scanf("%" NAG_IFMT " %*[^\n]", &nnz);

  la = 20 * nnz;
  lindb = 3 * n;
  liwork = 9 * n + 3;

  /* Allocate arrays */
  a = NAG_ALLOC(la, Complex);
  irow = NAG_ALLOC(la, Integer);
  icol = NAG_ALLOC(la, Integer);

  idiag = NAG_ALLOC(lindb, Integer);
  indb = NAG_ALLOC(lindb, Integer);
  ipivp = NAG_ALLOC(lindb, Integer);
  ipivq = NAG_ALLOC(lindb, Integer);
  istr = NAG_ALLOC(lindb + 1, Integer);

  iwork = NAG_ALLOC(liwork, Integer);

  if ((!a) || (!irow) || (!icol) || (!idiag) || (!indb) || (!ipivp) ||
      (!ipivq) || (!istr) || (!iwork)) {
    printf("Allocation failure!\n");
    exit_status = -1;
    goto END;
  }

  /* Initialize arrays */
  for (i = 0; i < la; i++) {
    a[i].re = 0.0;
    a[i].im = 0.0;
    irow[i] = 0;
    icol[i] = 0;
  }

  for (i = 0; i < lindb; i++) {
    indb[i] = 0;
    ipivp[i] = 0;
    ipivq[i] = 0;
    istr[i] = 0;
    idiag[i] = 0;
  }
  istr[lindb] = 0;

  for (i = 0; i < liwork; i++) {
    iwork[i] = 0;
  }

  /* Read the matrix A */
  for (i = 0; i < nnz; i++) {
    scanf(" ( %lf , %lf ) %" NAG_IFMT " %" NAG_IFMT "",
          &a[i].re, &a[i].im, &irow[i], &icol[i]);
  }
  scanf("%*[^\n] ");

  /* Read algorithmic parameters */
  scanf("%" NAG_IFMT " %lf %*[^\n]", &lfillg, &dtolg);

  /* nag_enum_name_to_value (x04nac): Converts NAG enum member name to value */
  scanf("%39s %*[^\n]", nag_enum_arg);
  pstrag = (Nag_SparseNsym_Piv) nag_enum_name_to_value(nag_enum_arg);

  /* nag_enum_name_to_value (x04nac): Converts NAG enum member name to value */
  scanf("%39s %*[^\n]", nag_enum_arg);
  milug = (Nag_SparseNsym_Fact) nag_enum_name_to_value(nag_enum_arg);

  /* Read algorithmic parameters */
  scanf("%" NAG_IFMT " %" NAG_IFMT " %*[^\n]", &nb, &nover);

  if (nb < 1)
    {
      printf("Value read for nb is out of range\n");
      exit_status = -4;
      goto END;
    }

  /* Allocate arrays */
  dtol = NAG_ALLOC(nb, double);
  istb = NAG_ALLOC(nb + 1, Integer);
  lfill = NAG_ALLOC(nb, Integer);
  npivm = NAG_ALLOC(nb, Integer);

  pstrat = (Nag_SparseNsym_Piv *) NAG_ALLOC(nb, Nag_SparseNsym_Piv);
  milu = (Nag_SparseNsym_Fact *) NAG_ALLOC(nb, Nag_SparseNsym_Fact);

  if ((!dtol) || (!istb) || (!lfill) || (!npivm) || (!pstrat) || (!milu)) {
    printf("Allocation failure!\n");
    exit_status = -1;
    goto END;
  }

  /* Initialize arrays */
  for (i = 0; i < nb; i++) {
    dtol[i] = 0.0;
    istb[i] = 0;
    lfill[i] = 0;
    npivm[i] = 0;
  }
  istb[nb] = 0;

  /* Define diagonal block indices.
   * In this example use blocks of MB consecutive rows and initialize
   * assuming no overlap.
   */
  mb = (n + nb - 1) / nb;
  for (k = 0; k < nb; k++) {
    istb[k] = k * mb + 1;
  }
  istb[nb] = n + 1;

  for (i = 0; i < n; i++) {
    indb[i] = i + 1;
  }

  /* Modify INDB and ISTB to account for overlap. */
  overlap(&n, &nnz, irow, icol, &nb, istb, indb, &lindb, &nover, iwork);

  /* Output matrix and blocking details */
  printf(" Original Matrix\n");
  printf(" n     = %4" NAG_IFMT "\n", n);
  printf(" nnz   = %4" NAG_IFMT "\n", nnz);
  printf(" nb    = %4" NAG_IFMT "\n", nb);

  for (k = 0; k < nb; k++) {
    printf(" Block = %4" NAG_IFMT ",%12s = %4" NAG_IFMT ",", k + 1, "order",
           istb[k + 1] - istb[k]);
    minval = indb[istb[k] - 1];
    for (j = istb[k]; j < istb[k + 1] - 1; j++) {
      minval = MIN(minval, indb[j]);
    }
    printf("%13s = %4" NAG_IFMT "\n", "start row", minval);
  }
  printf("\n");

  /* Set algorithmic parameters for each block from global values */
  for (k = 0; k < nb; k++) {
    lfill[k] = lfillg;
    dtol[k] = dtolg;
    pstrat[k] = pstrag;
    milu[k] = milug;
  }

  /* Initialize fail */
  INIT_FAIL(fail);

  /* Calculate factorization
   * 
   * nag_sparse_nherm_precon_bdilu (f11dtc). Calculates incomplete LU
   * factorization of local or overlapping diagonal blocks, mostly used
   * as incomplete LU preconditioner for complex sparse matrix.
   */

  nag_sparse_nherm_precon_bdilu(n, nnz, a, la, irow, icol, nb, istb, indb,
                                lindb, lfill, dtol, pstrat, milu, ipivp,
                                ipivq, istr, idiag, &nnzc, npivm, &fail);

  if (fail.code != NE_NOERROR) {
    printf("Error from nag_sparse_nherm_precon_bdilu (f11dtc).\n%s\n\n",
           fail.message);
    exit_status = -2;
    goto END;
  }

  /* Output details of the factorization */
  printf(" Factorization\n");
  printf(" nnzc  = %4" NAG_IFMT "\n\n", nnzc);
  printf(" Elements of factorization\n");
  printf("        i    j           c(i,  j)            Index\n");

  for (k = 0; k < nb; k++) {
    printf("  C_%1" NAG_IFMT
           "  -------------------------------------------\n", k + 1);

    /* Elements of the k-th block */
    for (i = istr[istb[k] - 1] - 1; i < istr[istb[k + 1] - 1] - 1; i++) {
      printf("%9" NAG_IFMT " %4" NAG_IFMT " %13.5e %13.5e %7" NAG_IFMT "\n",
             irow[i], icol[i], a[i].re, a[i].im, i + 1);
    }
  }

  k = 0;
  maxval_ret = npivm[k];
  for (k = 1; k < nb; k++) {
    maxval_ret = MAX(maxval_ret, npivm[k]);
  }

  printf("\n Details of factorized blocks\n");
  if (maxval_ret > 0) {

    /* Including pivoting details. */
    printf("  k   i    istr(i)   idiag(i)    indb(i)  ipivp(i)  ipivq(i)\n");
    for (k = 0; k < nb; k++) {
      i = istb[k] - 1;
      printf("%3" NAG_IFMT " %3" NAG_IFMT " %10" NAG_IFMT " ", k + 1, i + 1,
             istr[i]);
      printf("%10" NAG_IFMT " %10" NAG_IFMT " %10" NAG_IFMT " %10" NAG_IFMT
             "\n", idiag[i], indb[i], ipivp[i], ipivq[i]);

      for (i = istb[k]; i < istb[k + 1] - 1; i++) {
        printf("%3" NAG_IFMT " %10" NAG_IFMT " %10" NAG_IFMT " ",
               i + 1, istr[i], idiag[i]);
        printf("%10" NAG_IFMT " %10" NAG_IFMT " %10" NAG_IFMT "\n",
               indb[i], ipivp[i], ipivq[i]);
      }
      printf(" -----------------------------------------------------\n");
    }
  }
  else {

    /* No pivoting on any block. */
    printf("  k   i     istr(i)   idiag(i)    indb(i)\n");
    for (k = 0; k < nb; k++) {
      i = istb[k] - 1;
      printf("%3" NAG_IFMT " %3" NAG_IFMT " %10" NAG_IFMT " ", k + 1, i + 1,
             istr[i]);
      printf("%10" NAG_IFMT " %10" NAG_IFMT "\n", idiag[i], indb[i]);

      for (i = istb[k]; i < istb[k + 1] - 1; i++) {
        printf("%7" NAG_IFMT " %10" NAG_IFMT " %10" NAG_IFMT " %10" NAG_IFMT
               "\n", i + 1, istr[i], idiag[i], indb[i]);
      }
      printf("  --------------------------------------\n");
    }
  }

END:
  NAG_FREE(a);
  NAG_FREE(irow);
  NAG_FREE(icol);
  NAG_FREE(idiag);
  NAG_FREE(indb);
  NAG_FREE(ipivp);
  NAG_FREE(ipivq);
  NAG_FREE(istr);
  NAG_FREE(dtol);
  NAG_FREE(istb);
  NAG_FREE(lfill);
  NAG_FREE(npivm);
  NAG_FREE(pstrat);
  NAG_FREE(milu);
  NAG_FREE(iwork);

  return exit_status;
}

/* ************************************************************************** */

static void overlap(Integer *n, Integer *nnz, Integer *irow, Integer *icol,
                    Integer *nb, Integer *istb, Integer *indb, Integer *lindb,
                    Integer *nover, Integer *iwork)
{
  /* Purpose
   * =======
   *
   * This routine takes a set of row indices INDB defining the diagonal blocks
   * to be used in nag_sparse_nherm_precon_bdilu (f11dtc) to define a block
   * Jacobi or additive Schwarz preconditioner, and expands them to allow for
   * NOVER levels of overlap.
   *
   * The pointer array ISTB is also updated accordingly, so that the returned
   * values of ISTB and INDB can be passed on to 
   * nag_sparse_nherm_precon_bdilu (f11dtc) to define overlapping diagonal
   * blocks.
   *
   * ----------------------------------------------------------------------- */

  /* Scalars */
  Integer i, ik, ind, iover, j, k, l, n21, nadd, row;

  /* Find the number of nonzero elements in each row of the matrix A, and start
   * address of each row. Store the start addresses in iwork(n,...,2*n-1).
   */

  for (i = 0; i < (*n); i++) {
    iwork[i] = 0;
  }

  for (i = 0; i < (*nnz); i++) {
    iwork[irow[i] - 1] = iwork[irow[i] - 1] + 1;
  }
  iwork[(*n)] = 1;

  for (i = 0; i < (*n); i++) {
    iwork[(*n) + i + 1] = iwork[(*n) + i] + iwork[i];
  }

  /* Loop over blocks. */
  for (k = 0; k < (*nb); k++) {

    /* Initialize marker array. */
    for (j = 0; j < (*n); j++) {
      iwork[j] = 0;
    }

    /* Mark the rows already in block K in the workspace array. */
    for (l = istb[k]; l < istb[k + 1]; l++) {
      iwork[indb[l - 1] - 1] = 1;
    }

    /* Loop over levels of overlap. */
    for (iover = 1; iover <= (*nover); iover++) {

      /* Initialize counter of new row indices to be added. */
      ind = 0;

      /* Loop over the rows currently in the diagonal block. */
      for (l = istb[k]; l < istb[k + 1]; l++) {
        row = indb[l - 1];

        /* Loop over nonzero elements in row ROW. */
        for (i = iwork[(*n) + row - 1]; i < iwork[(*n) + row]; i++) {

          /* If the column index of the nonzero element is not in the
           * existing set for this block, store it to be added later, and
           * mark it in the marker array.
           */
          if (iwork[icol[i - 1] - 1] == 0) {
            iwork[icol[i - 1] - 1] = 1;
            iwork[2 * (*n) + 1 + ind] = icol[i - 1];
            ind = ind + 1;
          }
        }
      }

      /* Shift the indices in INDB and add the new entries for block K.
       * Change ISTB accordingly.
       */
      nadd = ind;
      if (istb[(*nb)] + nadd - 1 > (*lindb)) {
        printf("**** lindb too small, lindb = %" NAG_IFMT " ****\n", *lindb);
        exit(-1);
      }

      for (i = istb[(*nb)] - 1; i >= istb[k + 1]; i--) {
        indb[i + nadd - 1] = indb[i - 1];
      }

      n21 = 2 * (*n) + 1;
      ik = istb[k + 1] - 1;

      for (j = 0; j < nadd; j++) {
        indb[ik + j] = iwork[n21 + j];
      }

      for (j = k + 1; j < (*nb) + 1; j++) {
        istb[j] = istb[j] + nadd;
      }
    }
  }
  return;
}