|
1
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 |
SUBROUTINE ZLARFGP( N, ALPHA, X, INCX, TAU )
* * -- LAPACK auxiliary routine (version 3.3.1) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * -- April 2011 -- * * .. Scalar Arguments .. INTEGER INCX, N COMPLEX*16 ALPHA, TAU * .. * .. Array Arguments .. COMPLEX*16 X( * ) * .. * * Purpose * ======= * * ZLARFGP generates a complex elementary reflector H of order n, such * that * * H**H * ( alpha ) = ( beta ), H**H * H = I. * ( x ) ( 0 ) * * where alpha and beta are scalars, beta is real and non-negative, and * x is an (n-1)-element complex vector. H is represented in the form * * H = I - tau * ( 1 ) * ( 1 v**H ) , * ( v ) * * where tau is a complex scalar and v is a complex (n-1)-element * vector. Note that H is not hermitian. * * If the elements of x are all zero and alpha is real, then tau = 0 * and H is taken to be the unit matrix. * * Arguments * ========= * * N (input) INTEGER * The order of the elementary reflector. * * ALPHA (input/output) COMPLEX*16 * On entry, the value alpha. * On exit, it is overwritten with the value beta. * * X (input/output) COMPLEX*16 array, dimension * (1+(N-2)*abs(INCX)) * On entry, the vector x. * On exit, it is overwritten with the vector v. * * INCX (input) INTEGER * The increment between elements of X. INCX > 0. * * TAU (output) COMPLEX*16 * The value tau. * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION TWO, ONE, ZERO PARAMETER ( TWO = 2.0D+0, ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. * .. Local Scalars .. INTEGER J, KNT DOUBLE PRECISION ALPHI, ALPHR, BETA, BIGNUM, SMLNUM, XNORM COMPLEX*16 SAVEALPHA * .. * .. External Functions .. DOUBLE PRECISION DLAMCH, DLAPY3, DLAPY2, DZNRM2 COMPLEX*16 ZLADIV EXTERNAL DLAMCH, DLAPY3, DLAPY2, DZNRM2, ZLADIV * .. * .. Intrinsic Functions .. INTRINSIC ABS, DBLE, DCMPLX, DIMAG, SIGN * .. * .. External Subroutines .. EXTERNAL ZDSCAL, ZSCAL * .. * .. Executable Statements .. * IF( N.LE.0 ) THEN TAU = ZERO RETURN END IF * XNORM = DZNRM2( N-1, X, INCX ) ALPHR = DBLE( ALPHA ) ALPHI = DIMAG( ALPHA ) * IF( XNORM.EQ.ZERO ) THEN * * H = [1-alpha/abs(alpha) 0; 0 I], sign chosen so ALPHA >= 0. * IF( ALPHI.EQ.ZERO ) THEN IF( ALPHR.GE.ZERO ) THEN * When TAU.eq.ZERO, the vector is special-cased to be * all zeros in the application routines. We do not need * to clear it. TAU = ZERO ELSE * However, the application routines rely on explicit * zero checks when TAU.ne.ZERO, and we must clear X. TAU = TWO DO J = 1, N-1 X( 1 + (J-1)*INCX ) = ZERO END DO ALPHA = -ALPHA END IF ELSE * Only "reflecting" the diagonal entry to be real and non-negative. XNORM = DLAPY2( ALPHR, ALPHI ) TAU = DCMPLX( ONE - ALPHR / XNORM, -ALPHI / XNORM ) DO J = 1, N-1 X( 1 + (J-1)*INCX ) = ZERO END DO ALPHA = XNORM END IF ELSE * * general case * BETA = SIGN( DLAPY3( ALPHR, ALPHI, XNORM ), ALPHR ) SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'E' ) BIGNUM = ONE / SMLNUM * KNT = 0 IF( ABS( BETA ).LT.SMLNUM ) THEN * * XNORM, BETA may be inaccurate; scale X and recompute them * 10 CONTINUE KNT = KNT + 1 CALL ZDSCAL( N-1, BIGNUM, X, INCX ) BETA = BETA*BIGNUM ALPHI = ALPHI*BIGNUM ALPHR = ALPHR*BIGNUM IF( ABS( BETA ).LT.SMLNUM ) $ GO TO 10 * * New BETA is at most 1, at least SMLNUM * XNORM = DZNRM2( N-1, X, INCX ) ALPHA = DCMPLX( ALPHR, ALPHI ) BETA = SIGN( DLAPY3( ALPHR, ALPHI, XNORM ), ALPHR ) END IF SAVEALPHA = ALPHA ALPHA = ALPHA + BETA IF( BETA.LT.ZERO ) THEN BETA = -BETA TAU = -ALPHA / BETA ELSE ALPHR = ALPHI * (ALPHI/DBLE( ALPHA )) ALPHR = ALPHR + XNORM * (XNORM/DBLE( ALPHA )) TAU = DCMPLX( ALPHR/BETA, -ALPHI/BETA ) ALPHA = DCMPLX( -ALPHR, ALPHI ) END IF ALPHA = ZLADIV( DCMPLX( ONE ), ALPHA ) * IF ( ABS(TAU).LE.SMLNUM ) THEN * * In the case where the computed TAU ends up being a denormalized number, * it loses relative accuracy. This is a BIG problem. Solution: flush TAU * to ZERO (or TWO or whatever makes a nonnegative real number for BETA). * * (Bug report provided by Pat Quillen from MathWorks on Jul 29, 2009.) * (Thanks Pat. Thanks MathWorks.) * ALPHR = DBLE( SAVEALPHA ) ALPHI = DIMAG( SAVEALPHA ) IF( ALPHI.EQ.ZERO ) THEN IF( ALPHR.GE.ZERO ) THEN TAU = ZERO ELSE TAU = TWO DO J = 1, N-1 X( 1 + (J-1)*INCX ) = ZERO END DO BETA = -SAVEALPHA END IF ELSE XNORM = DLAPY2( ALPHR, ALPHI ) TAU = DCMPLX( ONE - ALPHR / XNORM, -ALPHI / XNORM ) DO J = 1, N-1 X( 1 + (J-1)*INCX ) = ZERO END DO BETA = XNORM END IF * ELSE * * This is the general case. * CALL ZSCAL( N-1, ALPHA, X, INCX ) * END IF * * If BETA is subnormal, it may lose relative accuracy * DO 20 J = 1, KNT BETA = BETA*SMLNUM 20 CONTINUE ALPHA = BETA END IF * RETURN * * End of ZLARFGP * END |