/* nag_eigen_complex_gen_quad (f02jqc) Example Program.
*
* Copyright 2017 Numerical Algorithms Group.
*
* Mark 26.2, 2017.
*/
#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;
}