/* nag_eigen_real_gen_quad (f02jcc) Example Program.
*
* Copyright 2022 Numerical Algorithms Group.
*
* Mark 28.3, 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);
}