/* nag_trans_hessenberg_controller (g13exc) Example Program. * * Copyright 1993 Numerical Algorithms Group * * Mark 3, 1993 * Mark 8 revised, 2004. */ #include #include #include #include #define A(I, J) a[(I) *tda + J] #define B(I, J) b[(I) *tdb + J] #define U(I, J) u[(I) *tdu + J] int main(void) { Integer exit_status = 0, i, j, m, n, tda, tdb, tdu; double *a = 0, *b = 0, one = 1.0, *u = 0, zero = 0.0; Nag_ControllerForm reduceto; NagError fail; INIT_FAIL(fail); printf( "nag_trans_hessenberg_controller (g13exc) Example Program Results\n"); /* Skip the heading in the data file and read the data. */ scanf("%*[^\n]"); scanf("%ld%ld", &n, &m); if (n >= 1 || m >= 1) { if (!(a = NAG_ALLOC(n*n, double)) || !(b = NAG_ALLOC(n*m, double)) || !(u = NAG_ALLOC(n*n, double))) { printf("Allocation failure\n"); exit_status = -1; goto END; } tda = n; tdb = m; tdu = n; } else { printf("Invalid n or m.\n"); exit_status = 1; return exit_status; } reduceto = Nag_UH_Controller; for (j = 0; j < n; ++j) for (i = 0; i < n; ++i) scanf("%lf", &A(i, j)); for (j = 0; j < m; ++j) for (i = 0; i < n; ++i) scanf("%lf", &B(i, j)); if (u) /* Initialise U as the identity matrix. */ for (i = 0; i < n; ++i) { for (j = 0; j < n; ++j) U(i, j) = zero; U(i, i) = one; } /* Reduce the pair (B,A) to reduceto controller Hessenberg form. */ /* nag_trans_hessenberg_controller (g13exc). * Unitary state-space transformation to reduce (BA) to * lower or upper controller Hessenberg form */ nag_trans_hessenberg_controller(n, m, reduceto, a, tda, b, tdb, u, tdu, &fail); if (fail.code != NE_NOERROR) { printf( "Error from nag_trans_hessenberg_controller (g13exc).\n%s\n", fail.message); exit_status = 1; goto END; } printf("\nThe transformed state transition matrix is\n\n"); for (i = 0; i < n; ++i) { for (j = 0; j < n; ++j) printf("%8.4f ", A(i, j)); printf("\n"); } printf("\nThe transformed input matrix is\n\n"); for (i = 0; i < n; ++i) { for (j = 0; j < m; ++j) printf("%8.4f ", B(i, j)); printf("\n"); } if (u) { printf("\nThe matrix that reduces (B,A) to "); printf("controller Hessenberg form is\n\n"); for (i = 0; i < n; ++i) { for (j = 0; j < n; ++j) printf("%8.4f ", U(i, j)); printf("\n"); } } END: if (a) NAG_FREE(a); if (b) NAG_FREE(b); if (u) NAG_FREE(u); return exit_status; }