/* 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 <nag_stdlib.h>
#include <stdio.h>
#include <nagf12.h>
#include <nagf16.h>

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

  /* Scalars */
  double alpha, h2, resr, sigma;
  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);

  Vprintf("nag_real_symm_banded_sparse_eigensystem_sol (f12fgc) Example "
          "Program Results\n\n");

  /* Skip heading in data file */
  Vscanf("%*[^\n] ");
  Vscanf("%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)) )
    {
      Vprintf("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);
  /* 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. */
  /* Main diagonal of A. */
  h2 = 1.0 / ((nx + 1) * (nx + 1));
  k = kl + ku;
  for (j = 0; j <= n-1; ++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 = 1; i <= nx; ++i)
    {
      ksup =  ksup + ldab;
      for (j = 1; 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 = 1; i <= nx - 1; ++i)
    {
      for (j = 1; 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. */
      Vprintf("\n The %4ld Ritz values are:\n\n", nconv);
      for (j = 0; j <= nconv-1; ++j)
        {
          if (eigest[j] <= 1.0e-10)
            {
              Vprintf("%8ld%5s%7.4f\n", j+1, "", eigv[j]);
            }
          else
            {
              Vprintf("%8ld%5s%7.4f%5s%18.16f\n", j+1, "", eigv[j],
                      " *** ",eigest[j]);
            }
        }
    }
  else
    {
      Vprintf(" Error from nag_real_symm_banded_sparse_eigensystem_sol "
              "(f12fgc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }
 END:
  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;
}