/* nag_eigen_complex_gen_quad (f02jqc) Example Program.
 *
 * NAGPRODCODE Version.
 *
 * Copyright 2016 Numerical Algorithms Group.
 *
 * Mark 26, 2016.
 */

#include <stdio.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagf02.h>
#include <nagx02.h>
#include <nagx04.h>
#include <nagm01.h>
#include <math.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_rank_sort((Pointer) e, 2*n, (ptrdiff_t) (sizeof(Complex)),
                  compare, Nag_Descending, indices, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_rank_sort (m01dsc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }
    /* nag_make_indices (m01zac).
     * Inverts a permutation converting a rank vector to an
     * index vector or vice versa
     */
    nag_make_indices(indices, 2*n, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_make_indices (m01zac).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }
    /* nag_reorder_vector (m01esc).
     * Reorders set of values of arbitrary data type into the
     * order specified by a set of indices
     */
    nag_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_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_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_reorder_vector (m01esc).\n%s\n", fail.message);
      exit_status = 1;
      goto END;
    }

    /* Print eigenvalues as 1 by 2n matrix using 
     * nag_gen_real_mat_print_comp (x04dbc).
     */
    fflush(stdout);
    nag_gen_complx_mat_print_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_gen_complx_mat_print_comp (x04dbc).
     */
    fflush(stdout);
    nag_gen_complx_mat_print_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_gen_complx_mat_print_comp (x04dbc).
     */
    fflush(stdout);
    nag_gen_complx_mat_print_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_gen_complx_mat_print_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_gen_real_mat_print_comp (x04cbc).
     */
    fflush(stdout);
    nag_gen_real_mat_print_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_gen_real_mat_print_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_gen_real_mat_print_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_gen_real_mat_print_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;
}