subroutine dlarfgp	(	integer	N,
		double precision	ALPHA,
		double precision, dimension( * )	X,
		integer	INCX,
		double precision	TAU
	)

DLARFGP generates an elementary reflector (Householder matrix) with non-negative beta.

Download DLARFGP + dependencies [TGZ] [ZIP] [TXT]

Purpose:

 DLARFGP generates a real elementary reflector H of order n, such
 that

       H * ( alpha ) = ( beta ),   H**T * 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**T ) ,
                     ( 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.

Parameters

[in]	N	N is INTEGER The order of the elementary reflector.
[in,out]	ALPHA	ALPHA is DOUBLE PRECISION On entry, the value alpha. On exit, it is overwritten with the value beta.
[in,out]	X	X is DOUBLE PRECISION array, dimension (1+(N-2)*abs(INCX)) On entry, the vector x. On exit, it is overwritten with the vector v.
[in]	INCX	INCX is INTEGER The increment between elements of X. INCX > 0.
[out]	TAU	TAU is DOUBLE PRECISION The value tau.

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Date

November 2015

Definition at line 106 of file dlarfgp.f.

*
*  -- LAPACK auxiliary routine (version 3.6.0) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     November 2015
*
*     .. Scalar Arguments ..
      INTEGER            incx, n
      DOUBLE PRECISION   alpha, tau
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION   x( * )
*     ..
*
*  =====================================================================
*
*     .. 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, bignum, savealpha, smlnum, 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 )
         smlnum = dlamch( 'S' ) / dlamch( 'E' )
         knt = 0
         IF( abs( beta ).LT.smlnum ) THEN
*
*           XNORM, BETA may be inaccurate; scale X and recompute them
*
            bignum = one / smlnum
   10       CONTINUE
            knt = knt + 1
            CALL dscal( n-1, bignum, x, incx )
            beta = beta*bignum
            alpha = alpha*bignum
            IF( abs( beta ).LT.smlnum )
     $         GO TO 10
*
*           New BETA is at most 1, at least SMLNUM
*
            xnorm = dnrm2( n-1, x, incx )
            beta = sign( dlapy2( alpha, xnorm ), alpha )
         END IF
         savealpha = alpha
         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
*
         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. This explains the next IF statement.
*
*           (Bug report provided by Pat Quillen from MathWorks on Jul 29, 2009.)
*           (Thanks Pat. Thanks MathWorks.)
*
            IF( savealpha.GE.zero ) THEN
               tau = zero
            ELSE
               tau = two
               DO j = 1, n-1
                  x( 1 + (j-1)*incx ) = 0
               END DO
               beta = -savealpha
            END IF
*
         ELSE 
*
*           This is the general case.
*
            CALL dscal( n-1, one / 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 DLARFGP
*

Here is the call graph for this function:

Here is the caller graph for this function: