LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dlarfgp.f
Go to the documentation of this file.
00001 *> \brief \b DLARFGP
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download DLARFGP + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarfgp.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarfgp.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarfgp.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE DLARFGP( N, ALPHA, X, INCX, TAU )
00022 * 
00023 *       .. Scalar Arguments ..
00024 *       INTEGER            INCX, N
00025 *       DOUBLE PRECISION   ALPHA, TAU
00026 *       ..
00027 *       .. Array Arguments ..
00028 *       DOUBLE PRECISION   X( * )
00029 *       ..
00030 *  
00031 *
00032 *> \par Purpose:
00033 *  =============
00034 *>
00035 *> \verbatim
00036 *>
00037 *> DLARFGP generates a real elementary reflector H of order n, such
00038 *> that
00039 *>
00040 *>       H * ( alpha ) = ( beta ),   H**T * H = I.
00041 *>           (   x   )   (   0  )
00042 *>
00043 *> where alpha and beta are scalars, beta is non-negative, and x is
00044 *> an (n-1)-element real vector.  H is represented in the form
00045 *>
00046 *>       H = I - tau * ( 1 ) * ( 1 v**T ) ,
00047 *>                     ( v )
00048 *>
00049 *> where tau is a real scalar and v is a real (n-1)-element
00050 *> vector.
00051 *>
00052 *> If the elements of x are all zero, then tau = 0 and H is taken to be
00053 *> the unit matrix.
00054 *> \endverbatim
00055 *
00056 *  Arguments:
00057 *  ==========
00058 *
00059 *> \param[in] N
00060 *> \verbatim
00061 *>          N is INTEGER
00062 *>          The order of the elementary reflector.
00063 *> \endverbatim
00064 *>
00065 *> \param[in,out] ALPHA
00066 *> \verbatim
00067 *>          ALPHA is DOUBLE PRECISION
00068 *>          On entry, the value alpha.
00069 *>          On exit, it is overwritten with the value beta.
00070 *> \endverbatim
00071 *>
00072 *> \param[in,out] X
00073 *> \verbatim
00074 *>          X is DOUBLE PRECISION array, dimension
00075 *>                         (1+(N-2)*abs(INCX))
00076 *>          On entry, the vector x.
00077 *>          On exit, it is overwritten with the vector v.
00078 *> \endverbatim
00079 *>
00080 *> \param[in] INCX
00081 *> \verbatim
00082 *>          INCX is INTEGER
00083 *>          The increment between elements of X. INCX > 0.
00084 *> \endverbatim
00085 *>
00086 *> \param[out] TAU
00087 *> \verbatim
00088 *>          TAU is DOUBLE PRECISION
00089 *>          The value tau.
00090 *> \endverbatim
00091 *
00092 *  Authors:
00093 *  ========
00094 *
00095 *> \author Univ. of Tennessee 
00096 *> \author Univ. of California Berkeley 
00097 *> \author Univ. of Colorado Denver 
00098 *> \author NAG Ltd. 
00099 *
00100 *> \date November 2011
00101 *
00102 *> \ingroup doubleOTHERauxiliary
00103 *
00104 *  =====================================================================
00105       SUBROUTINE DLARFGP( N, ALPHA, X, INCX, TAU )
00106 *
00107 *  -- LAPACK auxiliary routine (version 3.4.0) --
00108 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00109 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00110 *     November 2011
00111 *
00112 *     .. Scalar Arguments ..
00113       INTEGER            INCX, N
00114       DOUBLE PRECISION   ALPHA, TAU
00115 *     ..
00116 *     .. Array Arguments ..
00117       DOUBLE PRECISION   X( * )
00118 *     ..
00119 *
00120 *  =====================================================================
00121 *
00122 *     .. Parameters ..
00123       DOUBLE PRECISION   TWO, ONE, ZERO
00124       PARAMETER          ( TWO = 2.0D+0, ONE = 1.0D+0, ZERO = 0.0D+0 )
00125 *     ..
00126 *     .. Local Scalars ..
00127       INTEGER            J, KNT
00128       DOUBLE PRECISION   BETA, BIGNUM, SAVEALPHA, SMLNUM, XNORM
00129 *     ..
00130 *     .. External Functions ..
00131       DOUBLE PRECISION   DLAMCH, DLAPY2, DNRM2
00132       EXTERNAL           DLAMCH, DLAPY2, DNRM2
00133 *     ..
00134 *     .. Intrinsic Functions ..
00135       INTRINSIC          ABS, SIGN
00136 *     ..
00137 *     .. External Subroutines ..
00138       EXTERNAL           DSCAL
00139 *     ..
00140 *     .. Executable Statements ..
00141 *
00142       IF( N.LE.0 ) THEN
00143          TAU = ZERO
00144          RETURN
00145       END IF
00146 *
00147       XNORM = DNRM2( N-1, X, INCX )
00148 *
00149       IF( XNORM.EQ.ZERO ) THEN
00150 *
00151 *        H  =  [+/-1, 0; I], sign chosen so ALPHA >= 0
00152 *
00153          IF( ALPHA.GE.ZERO ) THEN
00154 *           When TAU.eq.ZERO, the vector is special-cased to be
00155 *           all zeros in the application routines.  We do not need
00156 *           to clear it.
00157             TAU = ZERO
00158          ELSE
00159 *           However, the application routines rely on explicit
00160 *           zero checks when TAU.ne.ZERO, and we must clear X.
00161             TAU = TWO
00162             DO J = 1, N-1
00163                X( 1 + (J-1)*INCX ) = 0
00164             END DO
00165             ALPHA = -ALPHA
00166          END IF
00167       ELSE
00168 *
00169 *        general case
00170 *
00171          BETA = SIGN( DLAPY2( ALPHA, XNORM ), ALPHA )
00172          SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'E' )
00173          KNT = 0
00174          IF( ABS( BETA ).LT.SMLNUM ) THEN
00175 *
00176 *           XNORM, BETA may be inaccurate; scale X and recompute them
00177 *
00178             BIGNUM = ONE / SMLNUM
00179    10       CONTINUE
00180             KNT = KNT + 1
00181             CALL DSCAL( N-1, BIGNUM, X, INCX )
00182             BETA = BETA*BIGNUM
00183             ALPHA = ALPHA*BIGNUM
00184             IF( ABS( BETA ).LT.SMLNUM )
00185      $         GO TO 10
00186 *
00187 *           New BETA is at most 1, at least SMLNUM
00188 *
00189             XNORM = DNRM2( N-1, X, INCX )
00190             BETA = SIGN( DLAPY2( ALPHA, XNORM ), ALPHA )
00191          END IF
00192          SAVEALPHA = ALPHA
00193          ALPHA = ALPHA + BETA
00194          IF( BETA.LT.ZERO ) THEN
00195             BETA = -BETA
00196             TAU = -ALPHA / BETA
00197          ELSE
00198             ALPHA = XNORM * (XNORM/ALPHA)
00199             TAU = ALPHA / BETA
00200             ALPHA = -ALPHA
00201          END IF
00202 *
00203          IF ( ABS(TAU).LE.SMLNUM ) THEN
00204 *
00205 *           In the case where the computed TAU ends up being a denormalized number,
00206 *           it loses relative accuracy. This is a BIG problem. Solution: flush TAU 
00207 *           to ZERO. This explains the next IF statement.
00208 *
00209 *           (Bug report provided by Pat Quillen from MathWorks on Jul 29, 2009.)
00210 *           (Thanks Pat. Thanks MathWorks.)
00211 *
00212             IF( SAVEALPHA.GE.ZERO ) THEN
00213                TAU = ZERO
00214             ELSE
00215                TAU = TWO
00216                DO J = 1, N-1
00217                   X( 1 + (J-1)*INCX ) = 0
00218                END DO
00219                BETA = -SAVEALPHA
00220             END IF
00221 *
00222          ELSE 
00223 *
00224 *           This is the general case.
00225 *
00226             CALL DSCAL( N-1, ONE / ALPHA, X, INCX )
00227 *
00228          END IF
00229 *
00230 *        If BETA is subnormal, it may lose relative accuracy
00231 *
00232          DO 20 J = 1, KNT
00233             BETA = BETA*SMLNUM
00234  20      CONTINUE
00235          ALPHA = BETA
00236       END IF
00237 *
00238       RETURN
00239 *
00240 *     End of DLARFGP
00241 *
00242       END
 All Files Functions