Generate block reflector (complex): * fill unused parts of reflectors matrix by zeros * fill diagonal of reflectors matrix by ones * generate triangular factor T -- ALGLIB routine -- 17.02.2010 Bochkanov Sergey *************************************************************************/
| 3200 | Bochkanov Sergey |
| 3201 | *************************************************************************/ |
| 3202 | static void cmatrixblockreflector(ap::complex_2d_array& a, |
| 3203 | ap::complex_1d_array& tau, |
| 3204 | bool columnwisea, |
| 3205 | int lengtha, |
| 3206 | int blocksize, |
| 3207 | ap::complex_2d_array& t, |
| 3208 | ap::complex_1d_array& work) |
| 3209 | { |
| 3210 | int i; |
| 3211 | int k; |
| 3212 | ap::complex v; |
| 3213 | |
| 3214 | |
| 3215 | // |
| 3216 | // Prepare Y (stored in TmpA) and T (stored in TmpT) |
| 3217 | // |
| 3218 | for(k = 0; k <= blocksize-1; k++) |
| 3219 | { |
| 3220 | |
| 3221 | // |
| 3222 | // fill beginning of new column with zeros, |
| 3223 | // load 1.0 in the first non-zero element |
| 3224 | // |
| 3225 | if( columnwisea ) |
| 3226 | { |
| 3227 | for(i = 0; i <= k-1; i++) |
| 3228 | { |
| 3229 | a(i,k) = 0; |
| 3230 | } |
| 3231 | } |
| 3232 | else |
| 3233 | { |
| 3234 | for(i = 0; i <= k-1; i++) |
| 3235 | { |
| 3236 | a(k,i) = 0; |
| 3237 | } |
| 3238 | } |
| 3239 | a(k,k) = 1; |
| 3240 | |
| 3241 | // |
| 3242 | // fill non-zero part of T, |
| 3243 | // |
| 3244 | for(i = 0; i <= k-1; i++) |
| 3245 | { |
| 3246 | if( columnwisea ) |
| 3247 | { |
| 3248 | v = ap::vdotproduct(&a(k, i), a.getstride(), "Conj", &a(k, k), a.getstride(), "N", ap::vlen(k,lengtha-1)); |
| 3249 | } |
| 3250 | else |
| 3251 | { |
| 3252 | v = ap::vdotproduct(&a(i, k), 1, "N", &a(k, k), 1, "Conj", ap::vlen(k,lengtha-1)); |
| 3253 | } |
| 3254 | work(i) = v; |
| 3255 | } |
| 3256 | for(i = 0; i <= k-1; i++) |
| 3257 | { |
| 3258 | v = ap::vdotproduct(&t(i, i), 1, "N", &work(i), 1, "N", ap::vlen(i,k-1)); |
| 3259 | t(i,k) = -tau(k)*v; |
no test coverage detected