/* nag_eigen_real_symm_sparse_arnoldi (f02fkc) Example Program.
 *
 * NAGPRODCODE Version.
 *
 * Copyright 2016 Numerical Algorithms Group.
 *
 * Mark 26, 2016.
 */

#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 double w[], const double rzest[],
                               Integer *istat, Nag_Comm *comm);
#ifdef __cplusplus
}
#endif

int main(void)
{

  /* Scalars */
  double h2, sigma;
  Integer exit_status = 0;
  Integer fileid, fmode, i, imon, j, k, lo, maxit, mode;
  Integer n, nconv, ncv, nev, nnz, nx, prtlvl, tdv;

  /* Local Arrays */
  double *w = 0, *a = 0, *resid = 0, *v = 0;
  double user[1];
  Integer *icol = 0, *irow = 0;
  Integer iuser[5];
  const char *filename = "f02fkce.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_symm_sparse_arnoldi (f02fkc) ");
  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]", &sigma);

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

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

  /* Construct A in sparse (SCS) format where:
   *   A_{i,i}   = 4/(h*h)
   *   A_{i+1,i) = -1/(h*h)
   *   A_{i+nx,i} = -1/(h*h) 
   */
  h2 = 1.0 / (double) ((nx + 1) * (nx + 1));

  /* Main Diagonal of A */
  k = 0;
  for (i = 1; i <= n; i++) {
    irow[k] = i;
    icol[k] = i;
    a[k] = 4.0 / h2;
    k++;
  }

  /* First subdiagonal of A. */
  for (i = 1; i <= nx; i++) {
    lo = (i - 1) * nx;
    for (j = lo + 1; j <= lo + nx - 1; j++) {
      irow[k] = j + 1;
      icol[k] = j;
      a[k] = -1.0 / h2;
      k++;
    }
  }

  /* nx-th subdiagonal of A. */
  for (i = 1; i < nx; i++) {
    lo = (i - 1) * nx;
    for (j = lo + 1; j <= lo + nx; j++) {
      irow[k] = j + nx;
      icol[k] = j;
      a[k] = -1.0 / h2;
      k++;
    }
  }

  /* 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).
     * If prtlvl >=10 internal monitoring information is also written.
     */
    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_symm_sparse_arnoldi (f02fkc).
   *   selected eigenvalues of real general matrix (driver).
   */
  nag_eigen_real_symm_sparse_arnoldi(n, nnz, a, irow, icol, nev, ncv, sigma,
                                     mymonit, myoption, &nconv, w, v, tdv,
                                     resid, &comm, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from  nag_eigen_real_symm_sparse_arnoldi (f02fkc)\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }

  if (imon > 0) {
    /* Close the monitoring file using 
     * nag_close_file (x04adc). 
     */
    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 = 0; i < nconv; i++) {
    /* nag_machine_precision (x02ajc) */
    if (resid[i] > (double) (100 * n) * nag_machine_precision) {
      printf("%7" NAG_IFMT "  %13.5e %13.5e\n", i + 1, w[i], resid[i]);
    }
    else {
      printf("%8" NAG_IFMT "  %13.5e\n", i + 1, w[i]);
    }
  }

END:

  NAG_FREE(w);
  NAG_FREE(a);
  NAG_FREE(v);
  NAG_FREE(resid);
  NAG_FREE(icol);
  NAG_FREE(irow);
  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);

  /* Set options using
   * nag_real_symm_sparse_eigensystem_option (f12fdc).
   */
  if (comm->iuser[0] > 0) {
    sprintf(rec, "Print Level=%5" NAG_IFMT, comm->iuser[0]);
    fail1.code = 1;
    nag_real_symm_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;
    nag_real_symm_sparse_eigensystem_option(rec, icomm, com, &fail1);
    *istat = MAX(*istat, fail1.code);
  }

  if (comm->iuser[2] > 0) {
    fail1.code = 1;
    nag_real_symm_sparse_eigensystem_option("Shifted Inverse", icomm, com,
                                            &fail1);
    *istat = MAX(*istat, fail1.code);
  }

  if (comm->iuser[3] > 0) {
    fail1.code = 1;
    sprintf(rec, "Monitoring=%5" NAG_IFMT, comm->iuser[4]);
    nag_real_symm_sparse_eigensystem_option(rec, icomm, com, &fail1);
    *istat = MAX(*istat, fail1.code);
  }
}

static void NAG_CALL mymonit(Integer ncv, Integer niter, Integer nconv,
                             const double 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).
     */

    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 = 0; i < nconv; i++) {
      sprintf(line, "  %4" NAG_IFMT " %13.5e %13.5e\n", i + 1, w[i],
              rzest[i]);
      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\n", nconv + 1, w[nconv]);
    nag_write_line(comm->iuser[4], line);
  }
  *istat = 0;
}