Chimie Théorique » OpenPath

      SUBROUTINE DLARFG( 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

*  =======

*

*  DLARFG 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, 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   ONE, ZERO

      PARAMETER          ( 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.1 ) THEN

         TAU = ZERO

         RETURN

      END IF

*

      XNORM = DNRM2( N-1, X, INCX )

*

      IF( XNORM.EQ.ZERO ) THEN

*

*        H  =  I

*

         TAU = ZERO

      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

         TAU = ( BETA-ALPHA ) / BETA

         CALL DSCAL( N-1, ONE / ( ALPHA-BETA ), X, INCX )

*

*        If ALPHA is subnormal, it may lose relative accuracy

*

         DO 20 J = 1, KNT

            BETA = BETA*SAFMIN

 20      CONTINUE

         ALPHA = BETA

      END IF

*

      RETURN

*

*     End of DLARFG

*

      END