/* nag_eigen_complex_gen_quad (f02jqc) Example Program. * * Copyright 2013 Numerical Algorithms Group. * * Mark 24, 2013. */ #include #include #include #include #include #include #include #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 structures and 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)