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