SRC/slarfgp.f

   1 *> \brief \b SLARFGP generates an elementary reflector (Householder matrix) with non-negative beta.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download SLARFGP + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slarfgp.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slarfgp.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slarfgp.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE SLARFGP( N, ALPHA, X, INCX, TAU )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       INTEGER            INCX, N
  25 *       REAL               ALPHA, TAU
  26 *       ..
  27 *       .. Array Arguments ..
  28 *       REAL               X( * )
  29 *       ..
  30 *
  31 *
  32 *> \par Purpose:
  33 *  =============
  34 *>
  35 *> \verbatim
  36 *>
  37 *> SLARFGP generates a real elementary reflector H of order n, such
  38 *> that
  39 *>
  40 *>       H * ( alpha ) = ( beta ),   H**T * H = I.
  41 *>           (   x   )   (   0  )
  42 *>
  43 *> where alpha and beta are scalars, beta is non-negative, and x is
  44 *> an (n-1)-element real vector.  H is represented in the form
  45 *>
  46 *>       H = I - tau * ( 1 ) * ( 1 v**T ) ,
  47 *>                     ( v )
  48 *>
  49 *> where tau is a real scalar and v is a real (n-1)-element
  50 *> vector.
  51 *>
  52 *> If the elements of x are all zero, then tau = 0 and H is taken to be
  53 *> the unit matrix.
  54 *> \endverbatim
  55 *
  56 *  Arguments:
  57 *  ==========
  58 *
  59 *> \param[in] N
  60 *> \verbatim
  61 *>          N is INTEGER
  62 *>          The order of the elementary reflector.
  63 *> \endverbatim
  64 *>
  65 *> \param[in,out] ALPHA
  66 *> \verbatim
  67 *>          ALPHA is REAL
  68 *>          On entry, the value alpha.
  69 *>          On exit, it is overwritten with the value beta.
  70 *> \endverbatim
  71 *>
  72 *> \param[in,out] X
  73 *> \verbatim
  74 *>          X is REAL array, dimension
  75 *>                         (1+(N-2)*abs(INCX))
  76 *>          On entry, the vector x.
  77 *>          On exit, it is overwritten with the vector v.
  78 *> \endverbatim
  79 *>
  80 *> \param[in] INCX
  81 *> \verbatim
  82 *>          INCX is INTEGER
  83 *>          The increment between elements of X. INCX > 0.
  84 *> \endverbatim
  85 *>
  86 *> \param[out] TAU
  87 *> \verbatim
  88 *>          TAU is REAL
  89 *>          The value tau.
  90 *> \endverbatim
  91 *
  92 *  Authors:
  93 *  ========
  94 *
  95 *> \author Univ. of Tennessee
  96 *> \author Univ. of California Berkeley
  97 *> \author Univ. of Colorado Denver
  98 *> \author NAG Ltd.
  99 *
 100 *> \date November 2015
 101 *
 102 *> \ingroup realOTHERauxiliary
 103 *
 104 *  =====================================================================
 105       SUBROUTINE SLARFGP( N, ALPHA, X, INCX, TAU )
 106 *
 107 *  -- LAPACK auxiliary routine (version 3.6.0) --
 108 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 109 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 110 *     November 2015
 111 *
 112 *     .. Scalar Arguments ..
 113       INTEGER            INCX, N
 114       REAL               ALPHA, TAU
 115 *     ..
 116 *     .. Array Arguments ..
 117       REAL               X( * )
 118 *     ..
 119 *
 120 *  =====================================================================
 121 *
 122 *     .. Parameters ..
 123       REAL               TWO, ONE, ZERO
 124       PARAMETER          ( TWO = 2.0E+0, ONE = 1.0E+0, ZERO = 0.0E+0 )
 125 *     ..
 126 *     .. Local Scalars ..
 127       INTEGER            J, KNT
 128       REAL               BETA, BIGNUM, SAVEALPHA, SMLNUM, XNORM
 129 *     ..
 130 *     .. External Functions ..
 131       REAL               SLAMCH, SLAPY2, SNRM2
 132       EXTERNAL           SLAMCH, SLAPY2, SNRM2
 133 *     ..
 134 *     .. Intrinsic Functions ..
 135       INTRINSIC          ABS, SIGN
 136 *     ..
 137 *     .. External Subroutines ..
 138       EXTERNAL           SSCAL
 139 *     ..
 140 *     .. Executable Statements ..
 141 *
 142       IF( N.LE.0 ) THEN
 143          TAU = ZERO
 144          RETURN
 145       END IF
 146 *
 147       XNORM = SNRM2( N-1, X, INCX )
 148 *
 149       IF( XNORM.EQ.ZERO ) THEN
 150 *
 151 *        H  =  [+/-1, 0; I], sign chosen so ALPHA >= 0.
 152 *
 153          IF( ALPHA.GE.ZERO ) THEN
 154 *           When TAU.eq.ZERO, the vector is special-cased to be
 155 *           all zeros in the application routines.  We do not need
 156 *           to clear it.
 157             TAU = ZERO
 158          ELSE
 159 *           However, the application routines rely on explicit
 160 *           zero checks when TAU.ne.ZERO, and we must clear X.
 161             TAU = TWO
 162             DO J = 1, N-1
 163                X( 1 + (J-1)*INCX ) = 0
 164             END DO
 165             ALPHA = -ALPHA
 166          END IF
 167       ELSE
 168 *
 169 *        general case
 170 *
 171          BETA = SIGN( SLAPY2( ALPHA, XNORM ), ALPHA )
 172          SMLNUM = SLAMCH( 'S' ) / SLAMCH( 'E' )
 173          KNT = 0
 174          IF( ABS( BETA ).LT.SMLNUM ) THEN
 175 *
 176 *           XNORM, BETA may be inaccurate; scale X and recompute them
 177 *
 178             BIGNUM = ONE / SMLNUM
 179    10       CONTINUE
 180             KNT = KNT + 1
 181             CALL SSCAL( N-1, BIGNUM, X, INCX )
 182             BETA = BETA*BIGNUM
 183             ALPHA = ALPHA*BIGNUM
 184             IF( ABS( BETA ).LT.SMLNUM )
 185      $         GO TO 10
 186 *
 187 *           New BETA is at most 1, at least SMLNUM
 188 *
 189             XNORM = SNRM2( N-1, X, INCX )
 190             BETA = SIGN( SLAPY2( ALPHA, XNORM ), ALPHA )
 191          END IF
 192          SAVEALPHA = ALPHA
 193          ALPHA = ALPHA + BETA
 194          IF( BETA.LT.ZERO ) THEN
 195             BETA = -BETA
 196             TAU = -ALPHA / BETA
 197          ELSE
 198             ALPHA = XNORM * (XNORM/ALPHA)
 199             TAU = ALPHA / BETA
 200             ALPHA = -ALPHA
 201          END IF
 202 *
 203          IF ( ABS(TAU).LE.SMLNUM ) THEN
 204 *
 205 *           In the case where the computed TAU ends up being a denormalized number,
 206 *           it loses relative accuracy. This is a BIG problem. Solution: flush TAU
 207 *           to ZERO. This explains the next IF statement.
 208 *
 209 *           (Bug report provided by Pat Quillen from MathWorks on Jul 29, 2009.)
 210 *           (Thanks Pat. Thanks MathWorks.)
 211 *
 212             IF( SAVEALPHA.GE.ZERO ) THEN
 213                TAU = ZERO
 214             ELSE
 215                TAU = TWO
 216                DO J = 1, N-1
 217                   X( 1 + (J-1)*INCX ) = 0
 218                END DO
 219                BETA = -SAVEALPHA
 220             END IF
 221 *
 222          ELSE
 223 *
 224 *           This is the general case.
 225 *
 226             CALL SSCAL( N-1, ONE / ALPHA, X, INCX )
 227 *
 228          END IF
 229 *
 230 *        If BETA is subnormal, it may lose relative accuracy
 231 *
 232          DO 20 J = 1, KNT
 233             BETA = BETA*SMLNUM
 234  20      CONTINUE
 235          ALPHA = BETA
 236       END IF
 237 *
 238       RETURN
 239 *
 240 *     End of SLARFGP
 241 *
 242       END