NAG Library Manual, Mark 30
Interfaces:  FL   CL   CPP   AD 

NAG CL Interface Introduction
Example description
/* nag_sparseig_complex_band_solve (f12auc) Example Program.
 *
 * Copyright 2024 Numerical Algorithms Group.
 *
 * Mark 30.0, 2024.
 */
#include <math.h>
#include <nag.h>

#define AB(I, J) ab[(J - 1) * pdab + kl + ku + I - J]
#define MB(I, J) mb[(J - 1) * pdab + kl + ku + I - J]

int main(void) {

  Integer exit_status = 0;
  Complex one = {1.0, 0.0};
  Complex minusone = {-1.0, 0.0};
  Complex zero = {0.0, 0.0};
  Nag_Boolean print_res = Nag_FALSE;
  /* Scalars */
  Complex ch_sub, ch_sup, alpha, sigma;
  double h, rho, anorm;
  Integer j, kl, ku, lcomm, licomm, n, ncol, nconv, ncv, nev, nx;
  Integer pdab, pdmb, pdv;
  /* Arrays */
  Complex commd[1];
  Integer icommd[1];
  Complex *ab = 0, *ax = 0, *comm = 0, *d = 0, *mb = 0;
  Complex *mx = 0, *resid = 0, *v = 0;
  double *d_print = 0;
  Integer *icomm = 0;
  /* NAG types */
  NagError fail;
  Nag_OrderType order = Nag_ColMajor;
  Nag_MatrixType matrix = Nag_GeneralMatrix;
  Nag_DiagType diag = Nag_NonUnitDiag;
  Nag_TransType trans = Nag_NoTrans;
  INIT_FAIL(fail);

  printf("nag_sparseig_complex_band_solve (f12auc) Example Program "
         "Results\n\n");
  fflush(stdout);
  /* Skip heading in data file */
  scanf("%*[^\n] ");
  /* Read mesh size, number of eigenvalues wanted and size of subspace. */
  scanf("%" NAG_IFMT "%" NAG_IFMT "%" NAG_IFMT "%*[^\n] ", &nx, &nev, &ncv);
  /* Read complex value close to which eigenvalues are sought. */
  scanf(" ( %lf , %lf ) %*[^\n] ", &sigma.re, &sigma.im);

  n = nx * nx;
  /* Initialize for the solution of a complex banded eigenproblem using
   * nag_sparseig_complex_band_init (f12atc).
   * But first get the required array lengths using arrays of length 1
   * to store required lengths.
   */
  licomm = -1;
  lcomm = -1;
  nag_sparseig_complex_band_init(n, nev, ncv, icommd, licomm, commd, lcomm,
                                 &fail);
  licomm = icommd[0];
  lcomm = (Integer)(commd[0].re);
  if (!(icomm = NAG_ALLOC((MAX(1, licomm)), Integer)) ||
      !(comm = NAG_ALLOC((MAX(1, lcomm)), Complex))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }
  nag_sparseig_complex_band_init(n, nev, ncv, icomm, licomm, comm, lcomm,
                                 &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_sparseig_complex_init (f12anc).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }

  /* Set options to show that the problem is a generalized eigenproblem and
   * that eigenvalues are to be computed using shifted inverses using
   * nag_sparseig_complex_option (f12arc).
   */
  nag_sparseig_complex_option("SHIFTED INVERSE", icomm, comm, &fail);
  nag_sparseig_complex_option("GENERALIZED", icomm, comm, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_sparseig_complex_option (f12arc).\n%s\n",
           fail.message);
    exit_status = 2;
    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.
   */
  kl = nx;
  ku = nx;
  pdab = 2 * kl + ku + 1;
  pdmb = pdab;
  if (!(ab = NAG_ALLOC(pdab * n, Complex)) ||
      !(mb = NAG_ALLOC(pdmb * n, Complex))) {
    printf("Allocation failure\n");
    exit_status = -2;
    goto END;
  }

  /* Zero out AB and MB. */
  for (j = 0; j < pdab * n; j++) {
    ab[j] = zero;
    mb[j] = zero;
  }

  rho = 100.0;
  h = 1.0 / (double)(nx + 1);
  /* assign complex values using nag_complex_create (a02bac). */
  ch_sub = nag_complex_create(0.5 * h * rho - 1.0, 0.0);
  ch_sup = nag_complex_create(-0.5 * h * rho - 1.0, 0.0);

  /* Set main diagonal of A and B to (4,0);
   * first subdiagonal and superdiagonal of A are (-1 +/- rho/2h,0)..
   * Outermost subdiagonal and superdiagonal are (-1,0).
   */
  for (j = 1; j <= n; j++) {
    AB(j, j) = nag_complex_create(4.0, 0.0);
    MB(j, j) = nag_complex_create(4.0, 0.0);
    if ((j - 1) % nx != 0 && j != n) {
      AB(j + 1, j) = ch_sub;
    }
    if (j % nx != 0 && j != 1) {
      AB(j - 1, j) = ch_sup;
    }
    if (j <= n - kl)
      AB(j + kl, j) = minusone;
    if (j > ku)
      AB(j - ku, j) = minusone;
  }
  for (j = 2; j < n; j++) {
    MB(j + 1, j) = one;
  }
  for (j = 2; j < n; j++) {
    MB(j - 1, j) = one;
  }

  /* Allocate v to hold eigenvectors and d to hold eigenvalues. */
  pdv = n;
  if (!(d = NAG_ALLOC((nev), Complex)) ||
      !(v = NAG_ALLOC((pdv) * (ncv), Complex)) ||
      !(resid = NAG_ALLOC((n), Complex))) {
    printf("Allocation failure\n");
    exit_status = -3;
    goto END;
  }
  /* Find eigenvalues closest in value to sigma and corresponding
   * eigenvectors using nag_sparseig_complex_band_solve (f12auc).
   */
  nag_sparseig_complex_band_solve(kl, ku, ab, mb, sigma, &nconv, d, v, resid, v,
                                  comm, icomm, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_sparseig_complex_band_solve (f12auc).\n%s\n",
           fail.message);
    exit_status = 3;
    goto END;
  }

  /* d_print stores eigenvalues (and possibly residuals) as real array. */
  if (!(d_print = NAG_ALLOC(nconv * 3, double))) {
    printf("Allocation failure\n");
    exit_status = -4;
    goto END;
  }
  for (j = 0; j < nconv; j++) {
    d_print[j] = (d[j]).re;
    d_print[j + nconv] = (d[j]).im;
  }
  /* By default this example will not calculate and print residuals.
   * It will, however if the following line is uncommented.
   */
  /* print_res = Nag_TRUE; */
  if (print_res) {
    ncol = 3;
    /* Compute the residual norm  ||A*x - lambda*B*x||.
     * Matrix vector products Ax and Bx, for complex banded A and B, are formed
     * using nag_blast_zgbmv (f16sbc).
     * The scaled vector subtraction is performed using nag_blast_zaxpby
     * (f16gcc).
     */
    if (!(ax = NAG_ALLOC(n, Complex)) || !(mx = NAG_ALLOC(n, Complex))) {
      printf("Allocation failure\n");
      exit_status = -4;
      goto END;
    }
    for (j = 0; j < nconv; j++) {
      nag_blast_zgbmv(order, trans, n, n, kl, ku, one, &ab[kl], pdab,
                      &v[pdv * j], 1, zero, ax, 1, &fail);
      if (fail.code == NE_NOERROR) {
        nag_blast_zgbmv(order, trans, n, n, kl, ku, one, &mb[kl], pdmb,
                        &v[pdv * j], 1, zero, mx, 1, &fail);
        if (fail.code == NE_NOERROR) {
          /* alpha holds -d[j] using nag_complex_negate (a02cec). */
          alpha = nag_complex_negate(d[j]);
          nag_blast_zaxpby(n, alpha, mx, 1, one, ax, 1, &fail);
        }
      }
      if (fail.code != NE_NOERROR) {
        printf("Error in calculating residual norm.\n%s\n", fail.message);
        exit_status = 4;
        goto END;
      }
      /* Get ||ax|| = ||Ax - lambda Bx|| using nag_blast_zge_norm (f16uac)
       * with the Frobenius norm on a complex general matrix with one column.
       */
      nag_blast_zge_norm(order, Nag_FrobeniusNorm, n, 1, ax, n, &anorm, &fail);
      if (fail.code != NE_NOERROR) {
        printf("Error from nag_blast_zge_norm (f16uac).\n%s\n", fail.message);
        exit_status = 5;
        goto END;
      }
      /* nag_complex_abs (a02dbc) returns the modulus of the eigenvalue */
      d_print[j + 2 * nconv] = anorm / nag_complex_abs(d[j]);
    }
  } else {
    ncol = 2;
  }
  printf("\n");

  /* Print the eigenvalues (and possibly the residuals) using the matrix
   * printer, nag_file_print_matrix_real_gen (x04cac).
   */
  printf("  Number of eigenvalues wanted     = %" NAG_IFMT "\n", nev);
  printf("  Number of eigenvalues converged  = %" NAG_IFMT "\n", nconv);
  printf("  Eigenvalues are closest to Sigma = ( %lf , %lf )\n\n", sigma.re,
         sigma.im);
  fflush(stdout);
  nag_file_print_matrix_real_gen(order, matrix, diag, nconv, ncol, d_print,
                                 nconv, " Ritz values closest to sigma", 0,
                                 &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_file_print_matrix_real_gen (x04cac).\n%s\n",
           fail.message);
    exit_status = 6;
  }
END:
  NAG_FREE(ab);
  NAG_FREE(ax);
  NAG_FREE(comm);
  NAG_FREE(d);
  NAG_FREE(mb);
  NAG_FREE(mx);
  NAG_FREE(resid);
  NAG_FREE(v);
  NAG_FREE(d_print);
  NAG_FREE(icomm);
  return exit_status;
}