/* nag_ztgsna (f08yyc) Example Program. * * Copyright 2011 Numerical Algorithms Group. * * Mark 23, 2011. */ #include #include #include #include #include #include #include int main(void) { /* Scalars */ double eps, snorm, stnrm, tnorm, tol; Integer i, j, m, n, pds, pdt, pdvl, pdvr; Integer exit_status = 0; /* Arrays */ Complex *s = 0, *t = 0, *vl = 0, *vr = 0; double *dif = 0, *scon = 0; /* Nag Types */ NagError fail; Nag_OrderType order; #ifdef NAG_COLUMN_MAJOR #define S(I, J) s[(J-1)*pds + I - 1] #define T(I, J) t[(J-1)*pdt + I - 1] order = Nag_ColMajor; #else #define S(I, J) s[(I-1)*pds + J - 1] #define T(I, J) t[(I-1)*pdt + J - 1] order = Nag_RowMajor; #endif INIT_FAIL(fail); printf("nag_ztgsna (f08yyc) Example Program Results\n\n"); /* Skip heading in data file */ scanf("%*[^\n]"); scanf("%ld%*[^\n]", &n); if (n < 0) { printf("Invalid n\n"); exit_status = 1; goto END; } m = n; pds = n; pdt = n; pdvl = n; pdvr = n; /* Allocate memory */ if (!(dif = NAG_ALLOC(n, double)) || !(scon = NAG_ALLOC(n, double)) || !(s = NAG_ALLOC(n*n, Complex)) || !(t = NAG_ALLOC(n*n, Complex)) || !(vl = NAG_ALLOC(n*m, Complex)) || !(vr = NAG_ALLOC(n*m, Complex))) { printf("Allocation failure\n"); exit_status = -1; goto END; } /* Read S and T from data file */ for (i = 1; i <= n; ++i) for (j = 1; j <= n; ++j) scanf(" ( %lf , %lf )", &S(i, j).re, &S(i, j).im); scanf("%*[^\n]"); for (i = 1; i <= n; ++i) for (j = 1; j <= n; ++j) scanf(" ( %lf , %lf )", &T(i, j).re, &T(i, j).im); scanf("%*[^\n]"); /* Calculate the left and right generalized eigenvectors of the * matrix pair (S,T) using nag_ztgevc (f08yxc). * NULL may be passed here in place of the select array since all * eigenvectors are requested. */ nag_ztgevc(order, Nag_BothSides, Nag_ComputeAll, NULL, n, s, pds, t, pdt, vl, pdvl, vr, pdvr, n, &m, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_ztgevc (f08yxc).\n%s\n", fail.message); exit_status = 1; goto END; } /* Estimate condition numbers for all the generalized eigenvalues and right * eigenvectors of the pair (S,T) using nag_ztgsna (f08yyc). * NULL may be passed here in place of the select array since all * eigenvectors are requested. */ nag_ztgsna(order, Nag_DoBoth, Nag_ComputeAll, NULL, n, s, pds, t, pdt, vl, pdvl, vr, pdvr, scon, dif, n, &m, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_ztgsna (f08yyc).\n%s\n", fail.message); exit_status = 1; goto END; } /* Print condition numbers of eigenvalues and right eigenvectors */ printf("Condition numbers of eigenvalues (scon) and right eigenvectors " "(diff),\n"); printf("scon: "); for (i = 0; i < m; ++i) printf(" %10.1e%s", scon[i], i%7 == 6?"\n ":""); printf("\ndif: "); for (i = 0; i < m; ++i) printf(" %10.1e%s", dif[i], i%7 == 6?"\n ":""); /* Compute the norm of (S,T) using nag_zge_norm (f16uac). */ eps = nag_machine_precision; nag_zge_norm(order, Nag_OneNorm, n, n, s, pds, &snorm, &fail); nag_zge_norm(order, Nag_OneNorm, n, n, t, pdt, &tnorm, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_zge_norm (f16uac).\n%s\n", fail.message); exit_status = 1; goto END; } if (snorm == 0.0) stnrm = ABS(tnorm); else if (tnorm == 0.0) stnrm = ABS(snorm); else if (ABS(snorm) >= ABS(tnorm)) stnrm = ABS(snorm)*sqrt(1.0+(tnorm/snorm)*(tnorm/snorm)); else stnrm = ABS(tnorm)*sqrt(1.0+(snorm/tnorm)*(snorm/tnorm)); printf("\n"); printf("Approximate error estimates for eigenvalues of (S,T)\n"); /* Calculate approximate error estimates */ tol = eps*stnrm; printf("\n\nError estimates for eigenvalues (errval) and right eigenvectors" " (errvec),\n"); printf("errval: "); for (i = 0; i < m; ++i) printf(" %10.1e%s", tol/scon[i], i%7 == 6?"\n ":""); printf("\nerrvec: "); for (i = 0; i < m; ++i) printf(" %10.1e%s", tol/dif[i], i%7 == 6?"\n ":""); END: if (dif) NAG_FREE(dif); if (scon) NAG_FREE(scon); if (s) NAG_FREE(s); if (t) NAG_FREE(t); if (vl) NAG_FREE(vl); if (vr) NAG_FREE(vr); return exit_status; }