/* nag_eigen_complex_gen_quad (f02jqc) Example Program.
 *
 * Copyright 2021 Numerical Algorithms Group.
 *
 * Mark 27.3, 2021.
 */

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

#ifdef __cplusplus
extern "C" {
#endif
static Integer NAG_CALL compare(const Nag_Pointer a, const Nag_Pointer b);
#ifdef __cplusplus
}
#endif

#define COMPLEX(A) A.re, A.im

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

  size_t *indices = 0;

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

  /* Double scalar and array declarations */
  double rbetaj;
  double tol = 0.0;
  double *bevl = 0, *bevr = 0, *s = 0;

  /* Complex scalar and array declarations */
  Complex *a = 0, *alpha = 0, *b = 0, *beta = 0, *c = 0, *e = 0, *vl = 0,
          *vr = 0, *cvr = 0;

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

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

  printf("nag_eigen_complex_gen_quad (f02jqc) Example Program Results\n\n");

  /* 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, Complex)) ||
      !(b = NAG_ALLOC(n * pdb, Complex)) ||
      !(c = NAG_ALLOC(n * pdc, Complex)) ||
      !(alpha = NAG_ALLOC(2 * n, Complex)) ||
      !(beta = NAG_ALLOC(2 * n, Complex)) || !(e = NAG_ALLOC(2 * n, Complex)) ||
      !(vl = NAG_ALLOC(2 * n * pdvl, Complex)) ||
      !(vr = NAG_ALLOC(2 * n * pdvr, Complex)) ||
      !(s = NAG_ALLOC(2 * n, double)) || !(bevr = NAG_ALLOC(2 * n, double)) ||
      !(bevl = NAG_ALLOC(2 * n, double)) || !(cvr = NAG_ALLOC(n, Complex)) ||
      !(indices = NAG_ALLOC(2 * n, size_t))) {
    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%lf", COMPLEX(&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%lf", COMPLEX(&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%lf", COMPLEX(&c[j * pdc + i]));
  scanf("%*[^\n] ");

  /* nag_eigen_complex_gen_quad (f02jqc):
   * Solve the quadratic eigenvalue problem */
  nag_eigen_complex_gen_quad(scal, jobvl, jobvr, sense, tol, n, a, pda, b, pdb,
                             c, pdc, alpha, beta, vl, pdvl, vr, pdvr, s, bevl,
                             bevr, &iwarn, &fail);

  if (fail.code != NE_NOERROR) {
    printf("Error from nag_eigen_complex_gen_quad (f02jqc).\n%s\n",
           fail.message);
    exit_status = -1;
    goto END;
  } else if (iwarn != 0) {
    printf("Warning from nag_eigen_complex_gen_quad (f02jqc).");
    printf("  iwarn = %" NAG_IFMT "\n", iwarn);
  }

  /* Display eigenvalues */
  for (j = 0; j < 2 * n; j++) {
    rbetaj = beta[j].re;
    if (rbetaj > 0.0) {
      e[j].re = alpha[j].re / rbetaj;
      e[j].im = alpha[j].im / rbetaj;
    } else {
      isinf = j + 1;
    }
  }
  if (isinf) {
    printf("Eigenvalue(%3" NAG_IFMT ") is infinite\n", isinf);
  } else {

    /* Sort eigenvalues by decscending absolute value and then by
     * descending real part. Sort requested eigenvectors correspondingly.
     */
    nag_sort_rank_sort((Pointer)e, 2 * n, (ptrdiff_t)(sizeof(Complex)), compare,
                       Nag_Descending, indices, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_sort_rank_sort (m01dsc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }
    /* nag_sort_permute_invert (m01zac).
     * Inverts a permutation converting a rank vector to an
     * index vector or vice versa
     */
    nag_sort_permute_invert(indices, 2 * n, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_sort_permute_invert (m01zac).\n%s\n",
             fail.message);
      exit_status = 1;
      goto END;
    }
    /* nag_sort_reorder_vector (m01esc).
     * Reorders set of values of arbitrary data type into the
     * order specified by a set of indices
     */
    nag_sort_reorder_vector((Pointer)e, 2 * n, sizeof(Complex),
                            (ptrdiff_t)(sizeof(Complex)), indices, &fail);
    if (fail.code == NE_NOERROR && jobvl == Nag_LeftVecs) {
      for (i = 0; i < n; i++) {
        nag_sort_reorder_vector((Pointer)&vl[i], 2 * n, sizeof(Complex),
                                (ptrdiff_t)(n * sizeof(Complex)), indices,
                                &fail);
      }
    }
    if (fail.code == NE_NOERROR && jobvr == Nag_RightVecs) {
      for (i = 0; i < n; i++) {
        nag_sort_reorder_vector((Pointer)&vr[i], 2 * n, sizeof(Complex),
                                (ptrdiff_t)(n * sizeof(Complex)), indices,
                                &fail);
      }
    }
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_sort_reorder_vector (m01esc).\n%s\n",
             fail.message);
      exit_status = 1;
      goto END;
    }

    /* Print eigenvalues as 1 by 2n matrix using
     * nag_file_print_matrix_real_gen_comp (x04dbc).
     */
    fflush(stdout);
    nag_file_print_matrix_complex_gen_comp(
        Nag_ColMajor, Nag_GeneralMatrix, Nag_NonUnitDiag, 1, 2 * n, e, 1,
        Nag_BracketForm, "%7.4f", "Eigenvalues:", Nag_NoLabels, 0,
        Nag_IntegerLabels, 0, 60, 0, 0, &fail);
  }
  printf("\n");

  if (fail.code == NE_NOERROR && jobvr == Nag_RightVecs) {
    /* Print right eigenvectors using
     * nag_file_print_matrix_complex_gen_comp (x04dbc).
     */
    fflush(stdout);
    nag_file_print_matrix_complex_gen_comp(
        Nag_ColMajor, Nag_GeneralMatrix, Nag_NonUnitDiag, n, 2 * n, vr, pdvr,
        Nag_BracketForm, "%7.4f", "Right Eigenvectors:", Nag_NoLabels, 0,
        Nag_IntegerLabels, 0, 60, 0, 0, &fail);
    printf("\n");
  }

  if (fail.code == NE_NOERROR && jobvl == Nag_LeftVecs) {
    /* Print left eigenvectors using
     * nag_file_print_matrix_complex_gen_comp (x04dbc).
     */
    fflush(stdout);
    nag_file_print_matrix_complex_gen_comp(
        Nag_ColMajor, Nag_GeneralMatrix, Nag_NonUnitDiag, n, 2 * n, vl, pdvl,
        Nag_BracketForm, "%7.4f", "Left Eigenvectors:", Nag_NoLabels, 0,
        Nag_IntegerLabels, 0, 60, 0, 0, &fail);
    printf("\n");
  }

  if (fail.code != NE_NOERROR) {
    printf("Error from nag_file_print_matrix_complex_gen_comp (x04dbc).\n%s\n",
           fail.message);
    exit_status = 1;
    goto END;
  }

  if (!cond_errors)
    goto END;

  /* Display condition numbers and errors, as required */
  if (sense == Nag_CondOnly || sense == Nag_CondBackErrLeft ||
      sense == Nag_CondBackErrRight || sense == Nag_CondBackErrBoth) {
    /* Print eigenvalue condition numbers as 1 by 2n matrix using
     * nag_file_print_matrix_real_gen_comp (x04cbc).
     */
    fflush(stdout);
    nag_file_print_matrix_real_gen_comp(
        Nag_ColMajor, Nag_GeneralMatrix, Nag_NonUnitDiag, 1, 2 * n, s, 1,
        "%17.2f", "Eigenvalue Condition numbers:", Nag_NoLabels, 0,
        Nag_IntegerLabels, 0, 60, 0, 0, &fail);
    printf("\n");
  }

  if (sense == Nag_BackErrRight || sense == Nag_BackErrBoth ||
      sense == Nag_CondBackErrRight || sense == Nag_CondBackErrBoth) {
    if (fail.code == NE_NOERROR) {
      fflush(stdout);
      nag_file_print_matrix_real_gen_comp(
          Nag_ColMajor, Nag_GeneralMatrix, Nag_NonUnitDiag, 1, 2 * n, bevr, 1,
          "%17.1e",
          "Backward errors for eigenvalues and "
          "right eigenvectors:",
          Nag_NoLabels, 0, Nag_IntegerLabels, 0, 60, 0, 0, &fail);
      printf("\n");
    }
  }

  if (sense == Nag_BackErrLeft || sense == Nag_BackErrBoth ||
      sense == Nag_CondBackErrLeft || sense == Nag_CondBackErrBoth) {
    if (fail.code == NE_NOERROR) {
      fflush(stdout);
      nag_file_print_matrix_real_gen_comp(
          Nag_ColMajor, Nag_GeneralMatrix, Nag_NonUnitDiag, 1, 2 * n, bevl, 1,
          "%17.1e",
          "Backward errors for eigenvalues and "
          "left eigenvectors:",
          Nag_NoLabels, 0, Nag_IntegerLabels, 0, 60, 0, 0, &fail);
      printf("\n");
    }
  }
  if (fail.code != NE_NOERROR) {
    printf("Error from nag_file_print_matrix_real_gen_comp (x04cbc).\n%s\n",
           fail.message);
    exit_status = 2;
  }

END:
  NAG_FREE(a);
  NAG_FREE(b);
  NAG_FREE(c);
  NAG_FREE(alpha);
  NAG_FREE(beta);
  NAG_FREE(e);
  NAG_FREE(vl);
  NAG_FREE(vr);
  NAG_FREE(s);
  NAG_FREE(bevr);
  NAG_FREE(bevl);
  NAG_FREE(cvr);
  NAG_FREE(indices);

  return (exit_status);
}
static Integer NAG_CALL compare(const Nag_Pointer a, const Nag_Pointer b) {
  double a_re = (*((const Complex *)a)).re;
  double a_im = (*((const Complex *)a)).im;
  double b_re = (*((const Complex *)b)).re;
  double b_im = (*((const Complex *)b)).im;
  double diff;
  Integer x;

  diff = (a_re * a_re + a_im * a_im) - (b_re * b_re + b_im * b_im);
  if (fabs(diff) < 1000.0 * x02ajc()) {
    if (a_re < b_re) {
      x = -1;
    } else if (a_re > b_re) {
      x = 1;
    } else {
      x = 0;
    }
  } else if (diff < 0.0) {
    x = -1;
  } else {
    x = 1;
  }
  return x;
}