LAPACK 3.3.1 Linear Algebra PACKage

# slarfg.f

Go to the documentation of this file.
00001       SUBROUTINE SLARFG( N, ALPHA, X, INCX, TAU )
00002 *
00003 *  -- LAPACK auxiliary routine (version 3.3.1) --
00004 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00005 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00006 *  -- April 2011                                                      --
00007 *
00008 *     .. Scalar Arguments ..
00009       INTEGER            INCX, N
00010       REAL               ALPHA, TAU
00011 *     ..
00012 *     .. Array Arguments ..
00013       REAL               X( * )
00014 *     ..
00015 *
00016 *  Purpose
00017 *  =======
00018 *
00019 *  SLARFG generates a real elementary reflector H of order n, such
00020 *  that
00021 *
00022 *        H * ( alpha ) = ( beta ),   H**T * H = I.
00023 *            (   x   )   (   0  )
00024 *
00025 *  where alpha and beta are scalars, and x is an (n-1)-element real
00026 *  vector. H is represented in the form
00027 *
00028 *        H = I - tau * ( 1 ) * ( 1 v**T ) ,
00029 *                      ( v )
00030 *
00031 *  where tau is a real scalar and v is a real (n-1)-element
00032 *  vector.
00033 *
00034 *  If the elements of x are all zero, then tau = 0 and H is taken to be
00035 *  the unit matrix.
00036 *
00037 *  Otherwise  1 <= tau <= 2.
00038 *
00039 *  Arguments
00040 *  =========
00041 *
00042 *  N       (input) INTEGER
00043 *          The order of the elementary reflector.
00044 *
00045 *  ALPHA   (input/output) REAL
00046 *          On entry, the value alpha.
00047 *          On exit, it is overwritten with the value beta.
00048 *
00049 *  X       (input/output) REAL array, dimension
00050 *                         (1+(N-2)*abs(INCX))
00051 *          On entry, the vector x.
00052 *          On exit, it is overwritten with the vector v.
00053 *
00054 *  INCX    (input) INTEGER
00055 *          The increment between elements of X. INCX > 0.
00056 *
00057 *  TAU     (output) REAL
00058 *          The value tau.
00059 *
00060 *  =====================================================================
00061 *
00062 *     .. Parameters ..
00063       REAL               ONE, ZERO
00064       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
00065 *     ..
00066 *     .. Local Scalars ..
00067       INTEGER            J, KNT
00068       REAL               BETA, RSAFMN, SAFMIN, XNORM
00069 *     ..
00070 *     .. External Functions ..
00071       REAL               SLAMCH, SLAPY2, SNRM2
00072       EXTERNAL           SLAMCH, SLAPY2, SNRM2
00073 *     ..
00074 *     .. Intrinsic Functions ..
00075       INTRINSIC          ABS, SIGN
00076 *     ..
00077 *     .. External Subroutines ..
00078       EXTERNAL           SSCAL
00079 *     ..
00080 *     .. Executable Statements ..
00081 *
00082       IF( N.LE.1 ) THEN
00083          TAU = ZERO
00084          RETURN
00085       END IF
00086 *
00087       XNORM = SNRM2( N-1, X, INCX )
00088 *
00089       IF( XNORM.EQ.ZERO ) THEN
00090 *
00091 *        H  =  I
00092 *
00093          TAU = ZERO
00094       ELSE
00095 *
00096 *        general case
00097 *
00098          BETA = -SIGN( SLAPY2( ALPHA, XNORM ), ALPHA )
00099          SAFMIN = SLAMCH( 'S' ) / SLAMCH( 'E' )
00100          KNT = 0
00101          IF( ABS( BETA ).LT.SAFMIN ) THEN
00102 *
00103 *           XNORM, BETA may be inaccurate; scale X and recompute them
00104 *
00105             RSAFMN = ONE / SAFMIN
00106    10       CONTINUE
00107             KNT = KNT + 1
00108             CALL SSCAL( N-1, RSAFMN, X, INCX )
00109             BETA = BETA*RSAFMN
00110             ALPHA = ALPHA*RSAFMN
00111             IF( ABS( BETA ).LT.SAFMIN )
00112      \$         GO TO 10
00113 *
00114 *           New BETA is at most 1, at least SAFMIN
00115 *
00116             XNORM = SNRM2( N-1, X, INCX )
00117             BETA = -SIGN( SLAPY2( ALPHA, XNORM ), ALPHA )
00118          END IF
00119          TAU = ( BETA-ALPHA ) / BETA
00120          CALL SSCAL( N-1, ONE / ( ALPHA-BETA ), X, INCX )
00121 *
00122 *        If ALPHA is subnormal, it may lose relative accuracy
00123 *
00124          DO 20 J = 1, KNT
00125             BETA = BETA*SAFMIN
00126  20      CONTINUE
00127          ALPHA = BETA
00128       END IF
00129 *
00130       RETURN
00131 *
00132 *     End of SLARFG
00133 *
00134       END