/* nag_eigen_complex_gen_quad (f02jqc) Example Program.
*
* Copyright 2014 Numerical Algorithms Group.
*
* Mark 24, 2013.
*/
#include <stdio.h>
#include <nag.h>
#include <nag_stdlib.h>
#include <nagf02.h>
#include <nagx02.h>
#include <nagx04.h>
#include <math.h>
#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;
/* Nag Types */
NagError fail;
Nag_ScaleType scal;
Nag_LeftVecsType jobvl;
Nag_RightVecsType jobvr;
Nag_CondErrType sense;
/* Double scalar and array declarations */
double inf, tmp, 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];
/* Initialise 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("%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, 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)))
{
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 = %ld\n", iwarn);
}
/* Infinity */
inf = X02ALC;
/* Display eigenvalues */
for (j = 0; j < 2*n; j++)
{
rbetaj = beta[j].re;
if (rbetaj >= 1.0)
{
e[j].re = alpha[j].re/rbetaj;
e[j].im = alpha[j].im/rbetaj;
}
else
{
tmp = inf*rbetaj;
if ((fabs(alpha[j].re)<tmp) && (fabs(alpha[j].im)<tmp))
{
e[j].re = alpha[j].re/rbetaj;
e[j].im = alpha[j].im/rbetaj;
}
else
{
e[j].re = inf;
e[j].im = 0.0;
}
}
if (fabs(e[j].re)<inf)
{
printf("Eigenvalue(%3ld) = (%11.4e, %11.4e)\n",j+1,
COMPLEX(e[j]));
}
else
{
printf("Eigenvalue(%3ld) is infinte\n",j+1);
}
}
if (jobvr == Nag_RightVecs)
{
printf("\n");
printf("Right Eigenvectors\n");
for (j = 0; j < 2*n; j += 3)
{
printf(" Eigenvector %ld",j + 1);
if (j < 2*n-1)
printf(" Eigenvector %ld",j + 2);
if (j < 2*n-2)
printf(" Eigenvector %ld",j + 3);
printf("\n");
for (i = 0; i < n; i++)
{
printf(" (%11.4e,%11.4e)", COMPLEX(vr[j*pdvr+i]));
if (j < 2*n-1)
printf(" (%11.4e,%11.4e)", COMPLEX(vr[(j+1)*pdvr+i]));
if (j < 2*n-2)
printf(" (%11.4e,%11.4e)", COMPLEX(vr[(j+2)*pdvr+i]));
printf("\n");
}
}
}
if (jobvl == Nag_LeftVecs)
{
printf("\n");
printf("Left Eigenvectors\n");
for (j = 0; j < 2*n; j += 3)
{
printf(" Eigenvector %ld",j + 1);
if (j < 2*n-1)
printf(" Eigenvector %ld",j + 2);
if (j < 2*n-2)
printf(" Eigenvector %ld",j + 3);
printf("\n");
for (i = 0; i < n; i++)
{
printf(" (%11.4e,%11.4e)", COMPLEX(vl[j*pdvl+i]));
if (j < 2*n-1)
printf(" (%11.4e,%11.4e)", COMPLEX(vl[(j+1)*pdvl+i]));
if (j < 2*n-2)
printf(" (%11.4e,%11.4e)", COMPLEX(vl[(j+2)*pdvl+i]));
printf("\n");
}
}
}
/* 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("%11.4e\n", s[j]);
}
if (sense==Nag_BackErrRight || sense==Nag_BackErrBoth ||
sense==Nag_CondBackErrRight || sense==Nag_CondBackErrBoth)
{
printf("\n");
printf("Backward errors for eigenvalues and right eigenvectors\n");
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)
{
printf("\n");
printf("Backward errors for eigenvalues and left eigenvectors\n");
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(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);
return(exit_status);
}