/* nag_sparseig_real_band_solve (f12agc) Example Program.
 *
 * Copyright 2020 Numerical Algorithms Group.
 *
 * Mark 27.1, 2020.
 */

#include <math.h>
#include <nag.h>
#include <stdio.h>

int main(void) {

  /* Constants */
  double rho = 100.0;

  /* Scalars */
  Complex res, eig;
  double alpha, h, resr, resi, sigmai, sigmar;
  Integer exit_status, i, j, k, kl, ksub, ksup, ku, lcomm, licomm;
  Integer ldab, n, nconv, ncv, nev, nx;
  /* Nag types */
  Nag_Boolean first;
  NagError fail;

  /* Arrays */
  double *ab = 0, *ax = 0, *comm = 0, *eigvr = 0, *eigvi = 0, *eigest = 0;
  double *mb = 0, *mx = 0, *resid = 0, *v = 0;
  Integer *icomm = 0;

  exit_status = 0;
  INIT_FAIL(fail);

  printf("nag_sparseig_real_band_solve (f12agc) Example "
         "Program Results\n");
  /* Skip heading in data file */
  scanf("%*[^\n] ");
  scanf("%" NAG_IFMT "%" NAG_IFMT "%" NAG_IFMT "%lf%lf%*[^\n] ", &nx, &nev,
        &ncv, &sigmar, &sigmai);
  n = nx * nx;
  /* ku, kl are number of superdiagonals and subdiagonals within the */
  /* band of matrices A and M. */
  kl = nx;
  ku = nx;
  /* ldab is the width of the band required to store the */
  /* factorization of the matrix A. */
  ldab = 2 * kl + ku + 1;
  /* Allocate memory. */
  if (!(ab = NAG_ALLOC(ldab * n, double)) || !(ax = NAG_ALLOC(n, double)) ||
      !(eigvr = NAG_ALLOC(ncv, double)) || !(eigvi = NAG_ALLOC(ncv, double)) ||
      !(eigest = NAG_ALLOC(ncv, double)) ||
      !(mb = NAG_ALLOC(ldab * n, double)) || !(mx = NAG_ALLOC(n, double)) ||
      !(resid = NAG_ALLOC(n, double)) || !(v = NAG_ALLOC(ncv * n, double))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }
  /* Initialize communication arrays for problem using
     nag_sparseig_real_band_init (f12afc).
     The first call sets lcomm = licomm = -1 to perform a workspace
     query. */
  lcomm = licomm = -1;
  if (!(comm = NAG_ALLOC(1, double)) || !(icomm = NAG_ALLOC(1, Integer))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }
  nag_sparseig_real_band_init(n, nev, ncv, icomm, licomm, comm, lcomm, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_sparseig_real_band_init "
           "(f12afc).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }
  lcomm = (Integer)comm[0];
  licomm = icomm[0];
  NAG_FREE(comm);
  NAG_FREE(icomm);
  if (!(comm = NAG_ALLOC(lcomm, double)) ||
      !(icomm = NAG_ALLOC(licomm, Integer))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }
  nag_sparseig_real_band_init(n, nev, ncv, icomm, licomm, comm, lcomm, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_sparseig_real_band_init "
           "(f12afc).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }
  /* Select the required spectrum using
     nag_sparseig_real_option (f12adc). */
  nag_sparseig_real_option("SHIFTED IMAGINARY", icomm, comm, &fail);
  /* Select the problem type using
     nag_sparseig_real_option (f12adc). */
  nag_sparseig_real_option("GENERALIZED", icomm, comm, &fail);

  /* Construct the matrix A in banded form and store in AB. */
  /* Zero out matrices first */
  for (j = 0; j < n * ldab; ++j) {
    ab[j] = 0.0;
    mb[j] = 0.0;
  }
  /* Main diagonal of A. */
  k = kl + ku;
  for (j = 0; j < n; ++j) {
    ab[k] = 4.;
    mb[k] = 4.;
    k = k + ldab;
  }
  /* First subdiagonal and superdiagonal of A. */
  ksup = kl + ku - 1;
  ksub = kl + ku + 1;
  h = .5 * rho / (double)(nx + 1);
  for (i = 1; i <= nx; ++i) {
    ksub = ksub + ldab;
    for (j = 1; j <= nx - 1; ++j) {
      ab[ksub] = -1. + h;
      ab[ksup] = -1. - h;
      ksup = ksup + ldab;
      ksub = ksub + ldab;
    }
    ksup = ksup + ldab;
  }

  ksub = kl + ku + 1 + ldab;
  ksup = kl + ku - 1;
  for (j = 1; j <= n - 1; ++j) {
    mb[ksup] = 1.;
    mb[ksub] = 1.;
    ksup = ksup + ldab;
    ksub = ksub + ldab;
  }
  /* KL-th subdiagonal and KU-th superdiagonal. */
  ksup = kl + nx * ldab;
  ksub = 2 * kl + ku;
  for (i = 1; i <= nx - 1; ++i) {
    for (j = 1; j <= nx; ++j) {
      ab[ksup] = -1.;
      ab[ksub] = -1.;
      ksup = ksup + ldab;
      ksub = ksub + ldab;
    }
  }

  /* Find eigenvalues closest in value to SIGMA and corresponding
     eigenvectors using nag_sparseig_real_band_solve (f12agc). */
  nag_sparseig_real_band_solve(kl, ku, ab, mb, sigmar, sigmai, &nconv, eigvr,
                               eigvi, v, resid, v, comm, icomm, &fail);
  if (fail.code == NE_NOERROR) {
    /* Compute the residual norm  ||A*x - lambda*Mx||. */
    first = Nag_TRUE;
    k = 0;
    for (j = 0; j <= nconv - 1; ++j) {
      if (eigvi[j] == 0.) {
        /* Eigenvalue is Real. */
        /* ax <- AV_k, where V_k is k-th Ritz vector. */
        /* Compute matrix-vector multiply using nag_blast_dgbmv
           (f16pbc). */
        nag_blast_dgbmv(Nag_ColMajor, Nag_NoTrans, n, n, kl, ku, 1., &ab[kl],
                        ldab, &v[k], 1, 0., ax, 1, &fail);
        /* mx <- MV_k. */
        nag_blast_dgbmv(Nag_ColMajor, Nag_NoTrans, n, n, kl, ku, 1., &mb[kl],
                        ldab, &v[k], 1, 0., mx, 1, &fail);
        /* ax <- ax - (lambda_j) * mx. */
        alpha = -eigvr[j];
        /* Compute vector update using nag_blast_daxpby (f16ecc). */
        nag_blast_daxpby(n, alpha, mx, 1, 1., ax, 1, &fail);
        /* resr = norm(ax). */
        /* Compute 2-norm of Ritz estimates using nag_blast_dge_norm
           (f16rac). */
        nag_blast_dge_norm(Nag_ColMajor, Nag_FrobeniusNorm, n, 1, ax, n, &resr,
                           &fail);
        /* Scale. */
        eigest[j] = resr / fabs(eigvr[j]);
      } else if (first) {
        /* First or a conjugate pair of eigenvalues. */

        /* Real part of residual Re[Ax-lamdaMx]. */
        /* ax <- AV_k, where V_k is real part of k-th Ritz vector. */
        /* Compute matrix-vector multiply using nag_blast_dgbmv
           (f16pbc). */
        nag_blast_dgbmv(Nag_ColMajor, Nag_NoTrans, n, n, kl, ku, 1., &ab[kl],
                        ldab, &v[k], 1, 0., ax, 1, &fail);
        /* mx <- MV_k */
        nag_blast_dgbmv(Nag_ColMajor, Nag_NoTrans, n, n, kl, ku, 1., &mb[kl],
                        ldab, &v[k], 1, 0., mx, 1, &fail);
        /* ax <- ax - (lambda_j).re * mx. */
        /* Compute vector update using nag_blast_daxpby (f16ecc). */
        alpha = -eigvr[j];
        nag_blast_daxpby(n, alpha, mx, 1, 1., ax, 1, &fail);
        /* mx <- MW_k where W_k is imaginary part k-th Ritz vector. */
        nag_blast_dgbmv(Nag_ColMajor, Nag_NoTrans, n, n, kl, ku, 1., &mb[kl],
                        ldab, &v[k + n], 1, 0., mx, 1, &fail);
        /* ax <- ax - (lambda_j).im * mx. */
        alpha = eigvi[j];
        nag_blast_daxpby(n, alpha, mx, 1, 1., ax, 1, &fail);
        /* resr = norm(ax). */
        /* Compute 2-norm of Ritz estimates using nag_blast_dge_norm
           (f16rac). */
        nag_blast_dge_norm(Nag_ColMajor, Nag_FrobeniusNorm, n, 1, ax, n, &resr,
                           &fail);
        /* Imaginary part of residual Im[Ax-lamdaMx]. */
        /* ax <- AW_k. */
        /* Compute matrix-vector multiply using nag_blast_dgbmv
           (f16pbc). */
        nag_blast_dgbmv(Nag_ColMajor, Nag_NoTrans, n, n, kl, ku, 1., &ab[kl],
                        ldab, &v[k + n], 1, 0., ax, 1, &fail);
        /* ax <- ax - (lambda_j).re * mx. */
        alpha = -eigvr[j];
        /* Compute vector update using nag_blast_daxpby (f16ecc). */
        nag_blast_daxpby(n, alpha, mx, 1, 1., ax, 1, &fail);
        /* mx <- MV_k. */
        nag_blast_dgbmv(Nag_ColMajor, Nag_NoTrans, n, n, kl, ku, 1., &mb[kl],
                        ldab, &v[k], 1, 0., mx, 1, &fail);
        /* ax <- ax - (lambda_j).im * mx. */
        alpha = -eigvi[j];
        nag_blast_daxpby(n, alpha, mx, 1, 1., ax, 1, &fail);
        /* resi = norm(ax). */
        /* Compute 2-norm of Ritz estimates using nag_blast_dge_norm
           (f16rac). */
        nag_blast_dge_norm(Nag_ColMajor, Nag_FrobeniusNorm, n, 1, ax, n, &resi,
                           &fail);
        /* Scale residual using Ritz value. */
        /* Assign to Complex type using nag_complex_create (a02bac) */
        res = nag_complex_create(resr, resi);
        eig = nag_complex_create(eigvr[j], eigvi[j]);
        eigest[j] = nag_complex_abs(res) / nag_complex_abs(eig);
        /* Set residual for second in conjugate pair. */
        eigest[j + 1] = eigest[j];
        first = Nag_FALSE;
      } else {
        /* Second of complex conjugate pair; */
        /* Already got residual from first in pair. */
        first = Nag_TRUE;
      }
      k = k + n;
    }

    /* Print computed eigenvalues. */
    printf("\n The %4" NAG_IFMT " generalized Ritz values", nconv);
    printf(" closest to \n");
    printf(" ( %7.4f, +-%7.4f ) are:\n\n", sigmar, sigmai);
    for (j = 0; j <= nconv - 1; ++j) {
      if (eigest[j] <= 1.0e-10) {
        printf("%8" NAG_IFMT "%5s( %7.4f, %7.4f ).\n", j + 1, "", eigvr[j],
               eigvi[j]);
      } else {
        printf("%8" NAG_IFMT "%5s( %7.4f, %7.4f )%5s%18.16f.\n", j + 1, "",
               eigvr[j], eigvi[j], " *** ", eigest[j]);
      }
    }
  } else {
    printf("Error from "
           "nag_sparseig_real_band_solve (f12agc).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }
END:
  NAG_FREE(ab);
  NAG_FREE(ax);
  NAG_FREE(comm);
  NAG_FREE(eigvr);
  NAG_FREE(eigvi);
  NAG_FREE(eigest);
  NAG_FREE(mb);
  NAG_FREE(mx);
  NAG_FREE(resid);
  NAG_FREE(v);
  NAG_FREE(icomm);

  return exit_status;
}