/* nag_eigen_real_gen_quad (f02jcc) Example Program.
*
* Copyright 2014 Numerical Algorithms Group.
*
* Mark 24, 2013.
*/
#include <stdio.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagf02.h>
#include <nagf16.h>
#include <nagx02.h>
#include <nagx04.h>
#include <math.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];
/* Initialise 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("%ld%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 = %ld\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(%3ld) = (%11.4e, %11.4e)\n",j+1,er[j],ei[j]);
} else {
printf("Eigenvalue(%3ld) is infinte\n",j+1);
}
}
if (jobvr == Nag_RightVecs) {
printf("\n");
fflush(stdout);
/* x04cac: Print out the right eigenvectors */
nag_gen_real_mat_print (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_gen_real_mat_print (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_gen_real_mat_print (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("%2ld %11.4e\n", j+1, s[j]);
}
if (sense==Nag_BackErrRight || sense==Nag_BackErrBoth ||
sense==Nag_CondBackErrRight || sense==Nag_CondBackErrBoth) {
/* nag_dmax_val (f16jnc).
* Get maximum value (bmax) and location of that value (j) of bevr.
*/
nag_dmax_val(2*n, bevr, 1, &j, &bmax, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_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_dmax_val (f16jnc).
* Get maximum value (bmax) and location of that value (j) of bevl.
*/
nag_dmax_val(2*n, bevl, 1, &j, &bmax, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_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);
}