Reduction of a square matrix to upper Hessenberg form: Q'*A*Q = H, where Q is an orthogonal matrix, H - Hessenberg matrix. Input parameters: A - matrix A with elements [0..N-1, 0..N-1] N - size of matrix A. Output parameters: A - matrices Q and P in compact form (see below). Tau - array of scalar factors which are used to form matrix Q.
| 2087 | October 31, 1992 |
| 2088 | *************************************************************************/ |
| 2089 | void rmatrixhessenberg(ap::real_2d_array& a, int n, ap::real_1d_array& tau) |
| 2090 | { |
| 2091 | int i; |
| 2092 | double v; |
| 2093 | ap::real_1d_array t; |
| 2094 | ap::real_1d_array work; |
| 2095 | |
| 2096 | ap::ap_error::make_assertion(n>=0, "RMatrixHessenberg: incorrect N!"); |
| 2097 | |
| 2098 | // |
| 2099 | // Quick return if possible |
| 2100 | // |
| 2101 | if( n<=1 ) |
| 2102 | { |
| 2103 | return; |
| 2104 | } |
| 2105 | tau.setbounds(0, n-2); |
| 2106 | t.setbounds(1, n); |
| 2107 | work.setbounds(0, n-1); |
| 2108 | for(i = 0; i <= n-2; i++) |
| 2109 | { |
| 2110 | |
| 2111 | // |
| 2112 | // Compute elementary reflector H(i) to annihilate A(i+2:ihi,i) |
| 2113 | // |
| 2114 | ap::vmove(&t(1), 1, &a(i+1, i), a.getstride(), ap::vlen(1,n-i-1)); |
| 2115 | generatereflection(t, n-i-1, v); |
| 2116 | ap::vmove(&a(i+1, i), a.getstride(), &t(1), 1, ap::vlen(i+1,n-1)); |
| 2117 | tau(i) = v; |
| 2118 | t(1) = 1; |
| 2119 | |
| 2120 | // |
| 2121 | // Apply H(i) to A(1:ihi,i+1:ihi) from the right |
| 2122 | // |
| 2123 | applyreflectionfromtheright(a, v, t, 0, n-1, i+1, n-1, work); |
| 2124 | |
| 2125 | // |
| 2126 | // Apply H(i) to A(i+1:ihi,i+1:n) from the left |
| 2127 | // |
| 2128 | applyreflectionfromtheleft(a, v, t, i+1, n-1, i+1, n-1, work); |
| 2129 | } |
| 2130 | } |
| 2131 | |
| 2132 | |
| 2133 | /************************************************************************* |
nothing calls this directly
no test coverage detected