/* nag_real_general_eigensystem (f02bjc) Example Program. * * Copyright 1991 Numerical Algorithms Group. * * Mark 2, 1991. * Mark 8 revised, 2004. */ #include #include #include #include #include #include #define A(I, J) a[(I) *tda + J] #define B(I, J) b[(I) *tdb + J] #define V(I, J) v[(I) *tdv + J] int main(void) { Nag_Boolean wantv; Complex *alfa = 0; Integer exit_status = 0, i, ip, *iter = 0, j, k, n, tda, tdb, tdv; double *a = 0, *b = 0, *beta = 0, tol, *v = 0; NagError fail; INIT_FAIL(fail); printf( "nag_real_general_eigensystem (f02bjc) Example Program Results\n"); scanf("%*[^\n]"); /* Skip heading in data file */ scanf("%ld", &n); if (n >= 1) { if (!(beta = NAG_ALLOC(n, double)) || !(a = NAG_ALLOC(n*n, double)) || !(b = NAG_ALLOC(n*n, double)) || !(v = NAG_ALLOC(n*n, double)) || !(iter = NAG_ALLOC(n, Integer)) || !(alfa = NAG_ALLOC(n, Complex))) { printf("Allocation failure\n"); exit_status = -1; goto END; } tda = n; tdb = n; tdv = n; } else { printf("Invalid n.\n"); exit_status = 1; return exit_status; } for (i = 0; i < n; ++i) for (j = 0; j < n; ++j) scanf("%lf", &A(i, j)); for (i = 0; i < n; ++i) for (j = 0; j < n; ++j) scanf("%lf", &B(i, j)); wantv = Nag_TRUE; /* nag_machine_precision (x02ajc). * The machine precision */ tol = nag_machine_precision; /* nag_real_general_eigensystem (f02bjc). * All eigenvalues and optionally eigenvectors of real * generalized eigenproblem, by QZ algorithm */ nag_real_general_eigensystem(n, a, tda, b, tdb, tol, alfa, beta, wantv, v, tdv, iter, &fail); if (fail.code != NE_NOERROR) { printf("Error from nag_real_general_eigensystem (f02bjc).\n%s\n", fail.message); exit_status = 1; goto END; } ip = 0; for (i = 0; i < n; ++i) { printf("Eigensolution %4ld\n", i+1); printf("alfa[%ld].re %7.3f", i, alfa[i].re); printf(" alfa[%ld].im %7.3f", i, alfa[i].im); printf(" beta[%ld] %7.3f\n", i, beta[i]); if (beta[i] == 0.0) printf("lambda is infinite"); else if (alfa[i].im == 0.0) { printf("lambda %7.3f\n", alfa[i].re/beta[i]); printf("Eigenvector\n"); for (j = 0; j < n; ++j) printf("%7.3f\n", V(j, i)); } else { printf("lambda %7.3f %7.3f\n", alfa[i].re/beta[i], alfa[i].im/beta[i]); printf("Eigenvector\n"); k = (Integer) pow((double) -1, (double)(ip+2)); for (j = 0; j < n; ++j) { printf("%7.3f", V(j, i-ip)); printf("%7.3f\n", k*V(j, i-ip+1)); } ip = 1-ip; } } printf("Number of iterations (machine-dependent)\n"); for (i = 0; i < n; ++i) printf("%2ld", iter[i]); printf("\n"); END: NAG_FREE(beta); NAG_FREE(a); NAG_FREE(b); NAG_FREE(v); NAG_FREE(iter); NAG_FREE(alfa); return exit_status; }