Generation of an elementary complex reflection transformation The subroutine generates elementary complex reflection H of order N, so that, for a given X, the following equality holds true: ( X(1) ) ( Beta ) H' * ( .. ) = ( 0 ), H'*H = I, Beta is a real number ( X(n) ) ( 0 ) where ( V(1) ) H = 1 - Tau * ( .. ) * ( conj(V(1)), ..., conj(V(n)) )
| 64 | September 30, 1994 |
| 65 | *************************************************************************/ |
| 66 | void complexgeneratereflection(ap::complex_1d_array& x, |
| 67 | int n, |
| 68 | ap::complex& tau) |
| 69 | { |
| 70 | int j; |
| 71 | ap::complex alpha; |
| 72 | double alphi; |
| 73 | double alphr; |
| 74 | double beta; |
| 75 | double xnorm; |
| 76 | double mx; |
| 77 | ap::complex t; |
| 78 | double s; |
| 79 | ap::complex v; |
| 80 | |
| 81 | if( n<=0 ) |
| 82 | { |
| 83 | tau = 0; |
| 84 | return; |
| 85 | } |
| 86 | |
| 87 | // |
| 88 | // Scale if needed (to avoid overflow/underflow during intermediate |
| 89 | // calculations). |
| 90 | // |
| 91 | mx = 0; |
| 92 | for(j = 1; j <= n; j++) |
| 93 | { |
| 94 | mx = ap::maxreal(ap::abscomplex(x(j)), mx); |
| 95 | } |
| 96 | s = 1; |
| 97 | if( ap::fp_neq(mx,0) ) |
| 98 | { |
| 99 | if( ap::fp_less(mx,1) ) |
| 100 | { |
| 101 | s = sqrt(ap::minrealnumber); |
| 102 | v = 1/s; |
| 103 | ap::vmul(&x(1), 1, ap::vlen(1,n), v); |
| 104 | } |
| 105 | else |
| 106 | { |
| 107 | s = sqrt(ap::maxrealnumber); |
| 108 | v = 1/s; |
| 109 | ap::vmul(&x(1), 1, ap::vlen(1,n), v); |
| 110 | } |
| 111 | } |
| 112 | |
| 113 | // |
| 114 | // calculate |
| 115 | // |
| 116 | alpha = x(1); |
| 117 | mx = 0; |
| 118 | for(j = 2; j <= n; j++) |
| 119 | { |
| 120 | mx = ap::maxreal(ap::abscomplex(x(j)), mx); |
| 121 | } |
| 122 | xnorm = 0; |
| 123 | if( ap::fp_neq(mx,0) ) |
no test coverage detected