SUBROUTINE DLARFP( N, ALPHA, X, INCX, TAU ) * * -- LAPACK auxiliary routine (version 3.2) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * November 2006 * * .. Scalar Arguments .. INTEGER INCX, N DOUBLE PRECISION ALPHA, TAU * .. * .. Array Arguments .. DOUBLE PRECISION X( * ) * .. * * Purpose * ======= * * DLARFP generates a real elementary reflector H of order n, such * that * * H * ( alpha ) = ( beta ), H' * H = I. * ( x ) ( 0 ) * * where alpha and beta are scalars, beta is non-negative, and x is * an (n-1)-element real vector. H is represented in the form * * H = I - tau * ( 1 ) * ( 1 v' ) , * ( v ) * * where tau is a real scalar and v is a real (n-1)-element * vector. * * If the elements of x are all zero, then tau = 0 and H is taken to be * the unit matrix. * * Otherwise 1 <= tau <= 2. * * Arguments * ========= * * N (input) INTEGER * The order of the elementary reflector. * * ALPHA (input/output) DOUBLE PRECISION * On entry, the value alpha. * On exit, it is overwritten with the value beta. * * X (input/output) DOUBLE PRECISION 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) DOUBLE PRECISION * 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 BETA, RSAFMN, SAFMIN, XNORM * .. * .. External Functions .. DOUBLE PRECISION DLAMCH, DLAPY2, DNRM2 EXTERNAL DLAMCH, DLAPY2, DNRM2 * .. * .. Intrinsic Functions .. INTRINSIC ABS, SIGN * .. * .. External Subroutines .. EXTERNAL DSCAL * .. * .. Executable Statements .. * IF( N.LE.0 ) THEN TAU = ZERO RETURN END IF * XNORM = DNRM2( N-1, X, INCX ) * IF( XNORM.EQ.ZERO ) THEN * * H = [+/-1, 0; I], sign chosen so ALPHA >= 0 * IF( ALPHA.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 ) = 0 END DO ALPHA = -ALPHA END IF ELSE * * general case * BETA = SIGN( DLAPY2( ALPHA, XNORM ), ALPHA ) SAFMIN = DLAMCH( 'S' ) / DLAMCH( 'E' ) KNT = 0 IF( ABS( BETA ).LT.SAFMIN ) THEN * * XNORM, BETA may be inaccurate; scale X and recompute them * RSAFMN = ONE / SAFMIN 10 CONTINUE KNT = KNT + 1 CALL DSCAL( N-1, RSAFMN, X, INCX ) BETA = BETA*RSAFMN ALPHA = ALPHA*RSAFMN IF( ABS( BETA ).LT.SAFMIN ) \$ GO TO 10 * * New BETA is at most 1, at least SAFMIN * XNORM = DNRM2( N-1, X, INCX ) BETA = SIGN( DLAPY2( ALPHA, XNORM ), ALPHA ) END IF ALPHA = ALPHA + BETA IF( BETA.LT.ZERO ) THEN BETA = -BETA TAU = -ALPHA / BETA ELSE ALPHA = XNORM * (XNORM/ALPHA) TAU = ALPHA / BETA ALPHA = -ALPHA END IF CALL DSCAL( N-1, ONE / ALPHA, X, INCX ) * * If BETA is subnormal, it may lose relative accuracy * DO 20 J = 1, KNT BETA = BETA*SAFMIN 20 CONTINUE ALPHA = BETA END IF * RETURN * * End of DLARFP * END