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