NAG Library Manual, Mark 27.2
```/* nag_eigen_complex_gen_quad (f02jqc) Example Program.
*
* Copyright 2021 Numerical Algorithms Group.
*
* Mark 27.2, 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);

/* 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] ");

* 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) {
fail.message);
exit_status = -1;
goto END;
} else if (iwarn != 0) {
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;
}
```