/* nag_trans_hessenberg_observer (g13ewc) 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 C(I,J) c[(I)*tdc + J] #define U(I,J) u[(I)*tdu + J] int main(void) { Integer exit_status=0, i, j, n, p, tda, tdc, tdu; NagError fail; Nag_ObserverForm reduceto; double *a=0, *c=0, one=1.0, *u=0, zero=0.0; INIT_FAIL(fail); Vprintf("nag_trans_hessenberg_observer (g13ewc) Example Program Results\n"); /* Skip the heading in the data file and read the data. */ Vscanf("%*[^\n]"); Vscanf("%ld%ld",&n,&p); if (n>=1 || p>=1) { if ( !( a = NAG_ALLOC(n*n, double)) || !( c = NAG_ALLOC(p*n, double)) || !( u = NAG_ALLOC(n*n, double)) ) { Vprintf("Allocation failure\n"); exit_status = -1; goto END; } tda = n; tdc = n; tdu = n; } else { Vprintf("Invalid n or p.\n"); exit_status = 1; return exit_status; } reduceto = Nag_UH_Observer; for (j = 0; j < n; ++j) for (i = 0; i < n; ++i) Vscanf("%lf", &A(i,j)); for (i = 0; i < p; ++i) for (j = 0; j < n; ++j) Vscanf("%lf", &C(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 (A,C) to reduceto observer Hessenberg form. */ /* nag_trans_hessenberg_observer (g13ewc). * Unitary state-space transformation to reduce (AC) to * lower or upper observer Hessenberg form */ nag_trans_hessenberg_observer(n, p, reduceto, a, tda, c, tdc, u, tdu, &fail); if (fail.code != NE_NOERROR) { Vprintf("Error from nag_trans_hessenberg_observer (g13ewc).\n%s\n", fail.message); exit_status = 1; goto END; } Vprintf("\nThe transformed state transition matrix is \n\n"); for (i = 0; i < n; ++i) { for (j = 0; j < n; ++j) Vprintf("%8.4f ",A(i,j)); Vprintf("\n"); } Vprintf("\nThe transformed input matrix is \n\n"); for (i = 0; i < p; ++i) { for (j = 0; j < n; ++j) Vprintf("%8.4f ",C(i,j)); Vprintf("\n"); } if (u) { Vprintf("\nThe transformation matrix that reduces (A,C) " "to observer Hessenberg form is \n\n"); for (i = 0; i < n; ++i) { for (j = 0; j < n; ++j) Vprintf("%8.4f ",U(i,j)); Vprintf("\n"); } } END: if (a) NAG_FREE(a); if (c) NAG_FREE(c); if (u) NAG_FREE(u); return exit_status; }