/* nag_real_symm_banded_sparse_eigensystem_sol (f12fgc) Example Program.
 *
 * Copyright 2005 Numerical Algorithms Group.
 *
 * Mark 8, 2005.
 */

#include <math.h>
#include <nag.h>
#include <nagx04.h>
#include <nag_stdlib.h>
#include <stdio.h>
#include <nagf12.h>
#include <nagf16.h>

int main(int argc, char *argv[])
{
  FILE     *fpin, *fpout;
  /* 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);

  /* Check for command-line IO options */
  fpin = nag_example_file_io(argc, argv, "-data", NULL);
  fpout = nag_example_file_io(argc, argv, "-results", NULL);

  fprintf(fpout, "nag_real_symm_banded_sparse_eigensystem_sol (f12fgc)"
          " Example Program Results\n\n");

  /* Skip heading in data file */
  fscanf(fpin, "%*[^\n] ");
  fscanf(fpin, "%ld%ld%ld%*[^\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)))
    {
      fprintf(fpout, "Allocation failure\n");
      exit_status = -1;
      goto END;
    }
  /* Initialise communication arrays for problem using
     nag_real_symm_banded_sparse_eigensystem_init (f12ffc). */
  nag_real_symm_banded_sparse_eigensystem_init(n, nev, ncv, icomm, licomm,
                                               comm, lcomm, &fail);
  if (fail.code != NE_NOERROR)
    {
      fprintf(fpout, "Error from nag_real_symm_banded_sparse_eigensystem_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 super-diagonal. */
  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_real_symm_banded_sparse_eigensystem_sol (f12fgc).
   */

  nag_real_symm_banded_sparse_eigensystem_sol(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 kth Ritz vector. */
          /* Compute matrix-vector multiply using nag_dgbmv (f16pbc). */
          nag_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_daxpby (f16ecc). */
          nag_daxpby(n, alpha, &v[k], 1, 1.0, ax, 1, &fail);
          /* Compute 2-norm of Ritz estimates using nag_dge_norm (f16rac).*/
          nag_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. */
      fprintf(fpout, "\n The %4ld Ritz values are:\n\n", nconv);
      for (j = 0; j <= nconv-1; ++j)
        {
          if (eigest[j] <= 1.0e-10)
            {
              fprintf(fpout, "%8ld%5s%7.4f\n", j+1, "", eigv[j]);
            }
          else
            {
              fprintf(fpout, "%8ld%5s%7.4f%5s%18.16f\n", j+1, "", eigv[j],
                      " *** ", eigest[j]);
            }
        }
    }
  else
    {
      fprintf(fpout, " Error from "
              "nag_real_symm_banded_sparse_eigensystem_sol (f12fgc).\n%s\n",
              fail.message);
      exit_status = 1;
      goto END;
    }
 END:
  if (fpin != stdin) fclose(fpin);
  if (fpout != stdout) fclose(fpout);
  if (ab) NAG_FREE(ab);
  if (ax) NAG_FREE(ax);
  if (comm) NAG_FREE(comm);
  if (eigv) NAG_FREE(eigv);
  if (eigest) NAG_FREE(eigest);
  if (mb) NAG_FREE(mb);
  if (resid) NAG_FREE(resid);
  if (v) NAG_FREE(v);
  if (icomm) NAG_FREE(icomm);

  return exit_status;
}