/* nag_eigen_real_gen_sparse_arnoldi (f02ekc) Example Program.
 *
 * Copyright 2017 Numerical Algorithms Group.
 *
 * Mark 26.1, 2017.
 */

#include <math.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <naga02.h>
#include <nagf02.h>
#include <nagf12.h>
#include <nagx02.h>
#include <nagx04.h>

/* User-defined Functions */
#ifdef __cplusplus
extern "C"
{
#endif
  static void NAG_CALL myoption(Integer icomm[], double com[], Integer *istat,
                                Nag_Comm *comm);

  static void NAG_CALL mymonit(Integer ncv, Integer niter, Integer nconv,
                               const Complex w[], const double rzest[],
                               Integer *istat, Nag_Comm *comm);
#ifdef __cplusplus
}
#endif

int main(void)
{
  /* Scalars */
  double one = 1.0, two = 2.0, three = 3.0;
  double h, rho, s, sigma;
  Integer exit_status = 0;
  Integer fileid, fmode, i, imon, k, maxit, mode;
  Integer n, nconv, ncv, nev, nnz, nx, prtlvl, tdv;

  /* Local Arrays */
  Complex *w = 0;
  double *a = 0, *resid = 0, *v = 0;
  double user[1];
  Integer *icolzp = 0, *irowix = 0;
  Integer iuser[5];
  const char *filename = "f02ekce.monit";

  /* Nag Types */
  Nag_Comm comm;
  NagError fail;

  INIT_FAIL(fail);

  comm.user = user;
  comm.iuser = iuser;
  user[0] = 0.0;
  iuser[0] = 0;

  /* Output preamble */
  printf("nag_eigen_real_gen_sparse_arnoldi (f02ekc) ");
  printf("Example Program Results\n\n");
  fflush(stdout);

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

  /* Read in problem size and parameters */
  scanf("%" NAG_IFMT "%*[^\n]%" NAG_IFMT "%*[^\n]%" NAG_IFMT "", &nx, &nev,
        &ncv);
  scanf("%*[^\n]%lf%*[^\n]%lf%*[^\n]", &rho, &sigma);

  n = nx * nx;
  nnz = 3 * n - 2;
  tdv = n;

  if (!(resid = NAG_ALLOC((ncv), double)) ||
      !(a = NAG_ALLOC((nnz), double)) ||
      !(icolzp = NAG_ALLOC((n + 1), Integer)) ||
      !(irowix = NAG_ALLOC((nnz), Integer)) ||
      !(w = NAG_ALLOC((ncv), Complex)) ||
      !(v = NAG_ALLOC((tdv) * (ncv), double))
         )
  {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  /* Construct  A in compressed column storage (CCS) format where:
   *   A_{i,i}   = 2 + i
   *   A_{i+1,i) = 3
   *   A_{i,i+1} = rho/(2n+2) - 1 
   */
  h = one / (double) (n + 1);
  s = rho * h / two - one;
  a[0] = two + one;
  a[1] = three;
  icolzp[0] = 1;
  irowix[0] = 1;
  irowix[1] = 2;

  k = 3;
  for (i = 2; i <= n - 1; i++) {
    icolzp[i - 1] = k;
    irowix[k - 1] = i - 1;
    irowix[k] = i;
    irowix[k + 1] = i + 1;
    a[k - 1] = s;
    a[k] = two + (double) (i);
    a[k + 1] = three;
    k = k + 3;
  }

  icolzp[n - 1] = k;
  icolzp[n] = k + 2;
  irowix[k - 1] = n - 1;
  irowix[k] = n;
  a[k - 1] = s;
  a[k] = two + (double) (n);

  /* Set some options via iuser array and routine argument OPTION.
   * iuser[0] = print level, iuser[1] = iteration limit,
   * iuser[2]>0 means shifted-invert mode
   * iuser[3]>0 means print monitoring info.
   */
  scanf("%" NAG_IFMT "%*[^\n]%" NAG_IFMT "%*[^\n]", &prtlvl, &maxit);
  scanf("%" NAG_IFMT "%*[^\n]%" NAG_IFMT "%*[^\n]", &mode, &imon);

  if (imon > 0) {
    /* Open the monitoring file for writing using
     * nag_open_file (x04acc):
     *   open unit number for reading, writing or appending, and
     *   associate unit with named file
     */
    /* If prtlvl >=10 internal monitoring information in addition to whatever is
     * written to fileid using mymonit.
     */
    fmode = 1;
    nag_open_file(filename, fmode, &fileid, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_open_file (x04acc) %s\n", fail.message);
      exit_status = 1;
      goto END;
    }
    iuser[4] = fileid;
  }

  iuser[0] = prtlvl;
  iuser[1] = maxit;
  iuser[2] = mode;
  iuser[3] = imon;

  /* Compute eigenvalues and eigenvectors using
   * nag_eigen_real_gen_sparse_arnoldi (f02ekc):
   *   selected eigenvalues of real general matrix driver.
   */
  nag_eigen_real_gen_sparse_arnoldi(n, nnz, a, icolzp, irowix, nev, ncv,
                                    sigma, mymonit, myoption, &nconv, w, v,
                                    tdv, resid, &comm, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_eigen_real_gen_sparse_arnoldi (f02ekc)\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }

  if (imon > 0) {
    /* Close the monitoring file using
     * nag_close_file (x04adc):
     *   close file associated with given unit number
     */
    nag_close_file(fileid, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_close_file (x04adc) %s\n", fail.message);
      exit_status = 1;
      goto END;
    }
  }

  printf(" The %4" NAG_IFMT " ", nconv);
  printf(" Ritz values of closest to %13.5e are \n", sigma);
  for (i = 1; i <= nconv; i++) {
    /* nag_machine_precision (x02ajc) */
    if (resid[i - 1] > (double) (100 * n) * nag_machine_precision) {
      printf("%7" NAG_IFMT "  ( %13.5e, %13.5e) ", i, w[i - 1].re,
             w[i - 1].im);
      printf("%13.5e\n", resid[i - 1]);
    }
    else {
      printf("%8" NAG_IFMT "  ( %13.5e, %13.5e) \n", i, w[i - 1].re,
             w[i - 1].im);
    }
  }

END:

  NAG_FREE(w);
  NAG_FREE(a);
  NAG_FREE(v);
  NAG_FREE(resid);
  NAG_FREE(icolzp);
  NAG_FREE(irowix);
  return exit_status;
}

static void NAG_CALL myoption(Integer icomm[], double com[], Integer *istat,
                              Nag_Comm *comm)
{
  NagError fail1;
  char rec[26];

  INIT_FAIL(fail1);

  if (comm->iuser[0] > 0) {
    sprintf(rec, "Print Level=%5" NAG_IFMT, comm->iuser[0]);
    fail1.code = 1;
    /* Set print level using
     * nag_real_sparse_eigensystem_option (f12adc)
     * Set a single option from a string.
     */
    nag_real_sparse_eigensystem_option(rec, icomm, com, &fail1);
    *istat = MAX(*istat, fail1.code);
  }

  if (comm->iuser[1] > 100) {
    sprintf(rec, "Iteration Limit=%5" NAG_IFMT, comm->iuser[1]);
    fail1.code = 1;
    /* Set iteration limit using
     * nag_real_sparse_eigensystem_option (f12adc)
     * Set a single option from a string.
     */
    nag_real_sparse_eigensystem_option(rec, icomm, com, &fail1);
    *istat = MAX(*istat, fail1.code);
  }

  if (comm->iuser[2] > 0) {
    fail1.code = 1;
    /* Set computational mode to shifted inverse real. */
    nag_real_sparse_eigensystem_option("Shifted Inverse Real", icomm, com,
                                       &fail1);
    *istat = MAX(*istat, fail1.code);
  }

  if (comm->iuser[3] > 0) {
    fail1.code = 1;
    /* Switch monitoring on and use the fileid stored in iuser[4]. */
    sprintf(rec, "Monitoring=%5" NAG_IFMT, comm->iuser[4]);
    nag_real_sparse_eigensystem_option(rec, icomm, com, &fail1);
    *istat = MAX(*istat, fail1.code);
  }
}

static void NAG_CALL mymonit(Integer ncv, Integer niter, Integer nconv,
                             const Complex w[], const double rzest[],
                             Integer *istat, Nag_Comm *comm)
{
  Integer i;
  char line[100];

  if (comm->iuser[3] > 0) {

    /* Write lines to the file we opened for monitoring using
     * nag_write_line (x04bac):
     *   write formatted record to external file.
     */

    if (niter == 1 && comm->iuser[2] > 0) {
      sprintf(line, " Arnoldi basis vectors used: %4" NAG_IFMT "\n", ncv);
      nag_write_line(comm->iuser[4], line);
      sprintf(line, " The following Ritz values (mu) are related to the\n");
      nag_write_line(comm->iuser[4], line);
      sprintf(line, " true eigenvalues (lambda) by lambda = sigma + 1/mu\n");
      nag_write_line(comm->iuser[4], line);
    }
    sprintf(line, "\n Iteration number %4" NAG_IFMT "\n", niter);
    nag_write_line(comm->iuser[4], line);
    sprintf(line,
            " Ritz values converged so far (%4" NAG_IFMT ") and their Ritz "
            "estimates:\n", nconv);
    nag_write_line(comm->iuser[4], line);

    for (i = 1; i <= nconv; i++) {
      sprintf(line,
              "  %4" NAG_IFMT " (%13.5e,%13.5e) %13.5e\n",
              i, w[i - 1].re, w[i - 1].im, rzest[i - 1]);
      nag_write_line(comm->iuser[4], line);
    }
    sprintf(line, " Next (unconverged) Ritz value:\n");
    nag_write_line(comm->iuser[4], line);
    sprintf(line,
            "  %4" NAG_IFMT " (%13.5e,%13.5e)\n",
            nconv + 1, w[nconv].re, w[nconv].im);
    nag_write_line(comm->iuser[4], line);
  }
  *istat = 0;
}