/* nag_sparseig_real_symm_band_solve (f12fgc) Example Program.
 *
 * Copyright 2022 Numerical Algorithms Group.
 *
 * Mark 28.6, 2022.
 */

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

int main(void) {
  /* Constants */
  Integer licomm = 140;

  /* Scalars */
  double alpha, h2, resr, sigma = 0.0;
  Integer exit_status, i, j, k, kl, ku, ksub, ksup, lcomm;
  Integer ldab, n, nconv, ncv, nev, nx;
  /* Nag types */
  NagError fail;

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

  exit_status = 0;
  INIT_FAIL(fail);

  printf("nag_sparseig_real_symm_band_solve (f12fgc)"
         " Example Program Results\n\n");

  /* Skip heading in data file */
  scanf("%*[^\n] ");
  scanf("%" NAG_IFMT "%" NAG_IFMT "%" NAG_IFMT "%*[^\n] ", &nx, &nev, &ncv);
  n = nx * nx;
  kl = nx;
  ku = nx;
  ldab = 2 * kl + ku + 1;
  lcomm = 3 * n + ncv * ncv + 8 * ncv + 60;
  /* Allocate memory */
  if (!(ab = NAG_ALLOC(ldab * n, double)) || !(ax = NAG_ALLOC(n, double)) ||
      !(comm = NAG_ALLOC(lcomm, double)) || !(eigv = NAG_ALLOC(ncv, double)) ||
      !(eigest = NAG_ALLOC(ncv, double)) || !(mb = NAG_ALLOC(1, double)) ||
      !(resid = NAG_ALLOC(n, double)) || !(v = NAG_ALLOC(n * ncv, double)) ||
      !(icomm = NAG_ALLOC(licomm, Integer))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }
  /* Initialize communication arrays for problem using
     nag_sparseig_real_symm_band_init (f12ffc). */
  nag_sparseig_real_symm_band_init(n, nev, ncv, icomm, licomm, comm, lcomm,
                                   &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_sparseig_real_symm_band_init"
           " (f12ffc).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }
  /* Construct the matrix A in banded form and store in AB */
  /* ku, kl are number of superdiagonals and subdiagonals within the
     band of matrices A and M. */
  for (j = 0; j < n * ldab; ++j) {
    ab[j] = 0.0;
  }
  /* Main diagonal of A. */
  h2 = 1.0 / ((nx + 1) * (nx + 1));
  k = kl + ku;
  for (j = 0; j < n; ++j) {
    ab[k] = 4.0 / h2;
    k = k + ldab;
  }
  /* First subdiagonal and superdiagonal of A. */
  ksup = kl + ku - 1;
  ksub = kl + ku + 1;
  for (i = 0; i < nx; ++i) {
    ksup = ksup + ldab;
    for (j = 0; j < nx - 1; ++j) {
      ab[ksub] = -1.0 / h2;
      ab[ksup] = -1.0 / h2;
      ksup = ksup + ldab;
      ksub = ksub + ldab;
    }
    ksub = ksub + ldab;
  }
  /* kl-th subdiagonal and ku-th superdiagonal. */
  ksup = kl + nx * ldab;
  ksub = 2 * kl + ku;
  for (i = 0; i < nx - 1; ++i) {
    for (j = 0; j < nx; ++j) {
      ab[ksup] = -1.0 / h2;
      ab[ksub] = -1.0 / h2;
      ksup = ksup + ldab;
      ksub = ksub + ldab;
    }
  }

  /* Find eigenvalues of largest magnitude and the corresponding
   * eigenvectors using nag_sparseig_real_symm_band_solve (f12fgc).
   */

  nag_sparseig_real_symm_band_solve(kl, ku, ab, mb, sigma, &nconv, eigv, v,
                                    resid, v, comm, icomm, &fail);
  if (fail.code == NE_NOERROR) {
    /* Compute the residual norm  ||A*x - lambda*x||. */
    k = 0;
    for (j = 0; j <= nconv - 1; ++j) {
      /* 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.0, &ab[kl],
                      ldab, &v[k], 1, 0.0, ax, 1, &fail);
      /* ax <- ax - (lambda_j) * V_k. */
      alpha = -eigv[j];
      /* Compute vector update using nag_blast_daxpby (f16ecc). */
      nag_blast_daxpby(n, alpha, &v[k], 1, 1.0, ax, 1, &fail);
      /* 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(eigv[j]);
      k = k + n;
    }

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

  return exit_status;
}