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

NAG CL Interface Introduction
Example description
/* nag_eigen_real_gen_quad (f02jcc) Example Program.
 *
 * Copyright 2022 Numerical Algorithms Group.
 *
 * Mark 28.7, 2022.
 */

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

int main(void) {
  /* Integer scalar and array declarations */
  Integer i, iwarn, j, pda, pdb, pdc, pdvl, pdvr, n;
  Integer exit_status = 0;

  /* Nag Types */
  NagError fail;
  Nag_ScaleType scal;
  Nag_LeftVecsType jobvl;
  Nag_RightVecsType jobvr;
  Nag_CondErrType sense;

  /* Double scalar and array declarations */
  double bmax, inf, tmp;
  double tol = 0.0;
  double *a = 0, *alphai = 0, *alphar = 0, *b = 0, *beta = 0, *bevl = 0;
  double *bevr = 0, *c = 0, *ei = 0, *er = 0, *s = 0, *vl = 0, *vr = 0;

  /* Character scalar declarations */
  char cjobvl[40], cjobvr[40], cscal[40], csense[40];

  /* Initialize the error structure */
  INIT_FAIL(fail);

  printf("nag_eigen_real_gen_quad (f02jcc) Example Program Results\n\n");
  fflush(stdout);

  /* Skip heading in data file */
  scanf("%*[^\n] ");

  /* Read in the problem size, scaling and output required */
  scanf("%" NAG_IFMT "%39s%39s%*[^\n] ", &n, cscal, csense);
  scal = (Nag_ScaleType)nag_enum_name_to_value(cscal);
  sense = (Nag_CondErrType)nag_enum_name_to_value(csense);
  scanf("%39s%39s%*[^\n] ", cjobvl, cjobvr);
  jobvl = (Nag_LeftVecsType)nag_enum_name_to_value(cjobvl);
  jobvr = (Nag_RightVecsType)nag_enum_name_to_value(cjobvr);

  pda = n;
  pdb = n;
  pdc = n;
  pdvl = n;
  pdvr = n;

  if (!(a = NAG_ALLOC(n * pda, double)) || !(b = NAG_ALLOC(n * pdb, double)) ||
      !(c = NAG_ALLOC(n * pdc, double)) ||
      !(alphai = NAG_ALLOC(2 * n, double)) ||
      !(alphar = NAG_ALLOC(2 * n, double)) ||
      !(beta = NAG_ALLOC(2 * n, double)) || !(ei = NAG_ALLOC(2 * n, double)) ||
      !(er = NAG_ALLOC(2 * n, double)) ||
      !(vl = NAG_ALLOC(2 * n * pdvl, double)) ||
      !(vr = NAG_ALLOC(2 * n * pdvr, double)) ||
      !(s = NAG_ALLOC(2 * n, double)) || !(bevr = NAG_ALLOC(2 * n, double)) ||
      !(bevl = NAG_ALLOC(2 * n, double))) {
    printf("Allocation failure\n");
    exit_status = -1;
    goto END;
  }

  /* Read in the matrix A */
  for (i = 0; i < n; i++)
    for (j = 0; j < n; j++)
      scanf("%lf", &a[j * pda + i]);
  scanf("%*[^\n] ");

  /* Read in the matrix B */
  for (i = 0; i < n; i++)
    for (j = 0; j < n; j++)
      scanf("%lf", &b[j * pdb + i]);
  scanf("%*[^\n] ");

  /* Read in the matrix C */
  for (i = 0; i < n; i++)
    for (j = 0; j < n; j++)
      scanf("%lf", &c[j * pdc + i]);
  scanf("%*[^\n] ");

  /* nag_eigen_real_gen_quad (f02jcc): Solve the quadratic eigenvalue problem */
  nag_eigen_real_gen_quad(scal, jobvl, jobvr, sense, tol, n, a, pda, b, pdb, c,
                          pdc, alphar, alphai, beta, vl, pdvl, vr, pdvr, s,
                          bevl, bevr, &iwarn, &fail);
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_eigen_real_gen_quad (f02jcc).\n%s\n", fail.message);
    exit_status = -1;
    goto END;
  } else if (iwarn != 0) {
    printf("Warning from nag_eigen_real_gen_quad (f02jcc).");
    printf("  iwarn = %" NAG_IFMT "\n", iwarn);
  }

  /* Infinity */
  inf = X02ALC;

  /* Display eigenvalues */
  for (j = 0; j < 2 * n; j++) {
    if (beta[j] >= 1.0) {
      er[j] = alphar[j] / beta[j];
      ei[j] = alphai[j] / beta[j];
    } else {
      tmp = inf * beta[j];
      if ((fabs(alphar[j]) < tmp) && (fabs(alphai[j]) < tmp)) {
        er[j] = alphar[j] / beta[j];
        ei[j] = alphai[j] / beta[j];
      } else {
        er[j] = inf;
        ei[j] = 0.0;
      }
    }
    if (er[j] < inf) {
      printf("Eigenvalue(%3" NAG_IFMT ") = (%11.4e, %11.4e)\n", j + 1, er[j],
             ei[j]);
    } else {
      printf("Eigenvalue(%3" NAG_IFMT ") is infinte\n", j + 1);
    }
  }

  if (jobvr == Nag_RightVecs) {
    printf("\n");
    fflush(stdout);
    /* x04cac: Print out the right eigenvectors */
    nag_file_print_matrix_real_gen(
        Nag_ColMajor, Nag_GeneralMatrix, Nag_NonUnitDiag, n, 2 * n, vr, pdvr,
        "Right eigenvectors (matrix VR)", NULL, &fail);
  }

  if (jobvl == Nag_LeftVecs && fail.code == NE_NOERROR) {
    printf("\n");
    fflush(stdout);
    /* x04cac: Print out the left eigenvectors */
    nag_file_print_matrix_real_gen(
        Nag_ColMajor, Nag_GeneralMatrix, Nag_NonUnitDiag, n, 2 * n, vl, pdvl,
        "Left eigenvectors (matrix VL)", NULL, &fail);
  }
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_file_print_matrix_real_gen (x04cac).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }
  /* Display condition numbers and errors, as required */
  if (sense == Nag_CondOnly || sense == Nag_CondBackErrLeft ||
      sense == Nag_CondBackErrRight || sense == Nag_CondBackErrBoth) {
    printf("\n");
    printf("Eigenvalue Condition numbers\n");
    for (j = 0; j < 2 * n; j++)
      printf("%2" NAG_IFMT "  %11.4e\n", j + 1, s[j]);
  }

  if (sense == Nag_BackErrRight || sense == Nag_BackErrBoth ||
      sense == Nag_CondBackErrRight || sense == Nag_CondBackErrBoth) {
    /* nag_blast_dmax_val (f16jnc).
     * Get maximum value (bmax) and location of that value (j) of bevr.
     */
    nag_blast_dmax_val(2 * n, bevr, 1, &j, &bmax, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_blast_dmax_val (f16jnc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }
    printf("\n");
    printf("Backward errors for eigenvalues and right eigenvectors\n");
    if (bmax < 10.0 * X02AJC) {
      printf("  All errors are less than 10 times machine precision\n");
    } else {
      for (j = 0; j < 2 * n; j++)
        printf("%11.4e\n", bevr[j]);
    }
  }

  if (sense == Nag_BackErrLeft || sense == Nag_BackErrBoth ||
      sense == Nag_CondBackErrLeft || sense == Nag_CondBackErrBoth) {
    /* nag_blast_dmax_val (f16jnc).
     * Get maximum value (bmax) and location of that value (j) of bevl.
     */
    nag_blast_dmax_val(2 * n, bevl, 1, &j, &bmax, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_blast_dmax_val (f16jnc).\n%s\n", fail.message);
      exit_status = 2;
      goto END;
    }
    printf("\n");
    printf("Backward errors for eigenvalues and left eigenvectors\n");
    if (bmax < 10.0 * X02AJC) {
      printf("  All errors are less than 10 times machine precision\n");
    } else {
      for (j = 0; j < 2 * n; j++)
        printf("%11.4e\n", bevl[j]);
    }
  }

END:
  NAG_FREE(a);
  NAG_FREE(b);
  NAG_FREE(c);
  NAG_FREE(alphai);
  NAG_FREE(alphar);
  NAG_FREE(beta);
  NAG_FREE(ei);
  NAG_FREE(er);
  NAG_FREE(vl);
  NAG_FREE(vr);
  NAG_FREE(s);
  NAG_FREE(bevr);
  NAG_FREE(bevl);

  return (exit_status);
}