/* nag_real_banded_sparse_eigensystem_sol (f12agc) Example Program.
 *
 * Copyright 2005 Numerical Algorithms Group.
 *
 * Mark 8, 2005.
 */

#include <math.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <stdio.h>
#include <naga02.h>
#include <nagf12.h>
#include <nagf16.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_real_banded_sparse_eigensystem_sol (f12agc) Example "
          "Program Results\n");
  /* Skip heading in data file */
  scanf("%*[^\n] ");

  scanf("%ld%ld%ld%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 */
  /* factorisationof 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;
    }
  /* Initialise communication arrays for problem using
     nag_real_banded_sparse_eigensystem_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_real_banded_sparse_eigensystem_init(n, nev, ncv, icomm, licomm,
                                          comm, lcomm, &fail);
  if (fail.code != NE_NOERROR)
    {
      printf("Error from nag_real_banded_sparse_eigensystem_init "
              "(f12afc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }
  lcomm = (Integer)comm[0];
  licomm = icomm[0];
  if (comm) NAG_FREE(comm);
  if (icomm) NAG_FREE(icomm);
  if (!(comm = NAG_ALLOC(lcomm, double)) ||
      !(icomm = NAG_ALLOC(licomm, Integer)))
    {
      printf("Allocation failure\n");
      exit_status = -1;
      goto END;
    }
  nag_real_banded_sparse_eigensystem_init(n, nev, ncv, icomm, licomm,
                                          comm, lcomm, &fail);
  if (fail.code != NE_NOERROR)
    {
      printf("Error from nag_real_banded_sparse_eigensystem_init "
              "(f12afc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }
  /* Select the required spectrum using
     nag_real_sparse_eigensystem_option (f12adc). */
  nag_real_sparse_eigensystem_option("SHIFTED IMAGINARY", icomm, comm,
                                     &fail);
  /* Select the problem type using
     nag_real_sparse_eigensystem_option (f12adc). */
  nag_real_sparse_eigensystem_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 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.;
          ab[ksub] = -1.;
          ksup = ksup + ldab;
          ksub = ksub + ldab;
        }
    }

  /* Find eigenvalues closest in value to SIGMA and corresponding
     eigenvectors using nag_real_banded_sparse_eigensystem_sol (f12agc). */
  nag_real_banded_sparse_eigensystem_sol(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 kth Ritz vector. */
              /* Compute matrix-vector multiply using nag_dgbmv
                 (f16pbc). */
              nag_dgbmv(Nag_ColMajor, Nag_NoTrans, n, n, kl, ku, 1.,
                        &ab[kl], ldab, &v[k], 1, 0., ax, 1, &fail);
              /* mx <- MV_k. */
              nag_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_daxpby (f16ecc). */
              nag_daxpby(n, alpha, mx, 1, 1., ax, 1, &fail);
              /* resr = norm(ax). */
              /* 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(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 kth Ritz vector. */
              /* Compute matrix-vector multiply using nag_dgbmv
                 (f16pbc). */
              nag_dgbmv(Nag_ColMajor, Nag_NoTrans, n, n, kl, ku, 1.,
                        &ab[kl], ldab, &v[k], 1, 0., ax, 1, &fail);
              /* mx <- MV_k */
              nag_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_daxpby (f16ecc). */
              alpha = -eigvr[j];
              nag_daxpby(n, alpha, mx, 1, 1., ax, 1, &fail);
              /* mx <- MW_k where W_k is imaginary part kth Ritz vector. */
              nag_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_daxpby(n, alpha, mx, 1, 1., ax, 1, &fail);
              /* resr = norm(ax). */
              /* Compute 2-norm of Ritz estimates using nag_dge_norm
                 (f16rac).*/
              nag_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_dgbmv
                 (f16pbc). */
              nag_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_daxpby (f16ecc). */
              nag_daxpby(n, alpha, mx, 1, 1., ax, 1, &fail);
              /* mx <- MV_k. */
              nag_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_daxpby(n, alpha, mx, 1, 1., ax, 1, &fail);
              /* resi = norm(ax). */
              /* Compute 2-norm of Ritz estimates using nag_dge_norm
                 (f16rac).*/
              nag_dge_norm(Nag_ColMajor, Nag_FrobeniusNorm, n, 1, ax,
                           n, &resi, &fail);
              /* Scale residual using Ritz value. */
              /* Assign to Complex type using nag_complex (a02bac) */
              res = nag_complex(resr, resi);
              eig = nag_complex(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 %4ld 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("%8ld%5s( %7.4f, %7.4f ).\n", j+1,
                      "", eigvr[j], eigvi[j]);
            }
          else
            {
              printf("%8ld%5s( %7.4f, %7.4f )%5s%18.16f.\n", j+1,
                      "", eigvr[j], eigvi[j], " *** ", eigest[j]);
            }
        }
    }
  else
    {
      printf("Error from "
             "nag_real_banded_sparse_eigensystem_sol (f12agc).\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 (eigvr) NAG_FREE(eigvr);
  if (eigvi) NAG_FREE(eigvi);
  if (eigest) NAG_FREE(eigest);
  if (mb) NAG_FREE(mb);
  if (mx) NAG_FREE(mx);
  if (resid) NAG_FREE(resid);
  if (v) NAG_FREE(v);
  if (icomm) NAG_FREE(icomm);


  return exit_status;
}