LAPACK: SRC/dlarfgp.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DLARFGP( N, ALPHA, X, INCX, TAU )
00002 *
00003 *  -- LAPACK auxiliary routine (version 3.2.2) --
00004 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00005 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00006 *     June 2010
00007 *
00008 *     .. Scalar Arguments ..
00009       INTEGER            INCX, N
00010       DOUBLE PRECISION   ALPHA, TAU
00011 *     ..
00012 *     .. Array Arguments ..
00013       DOUBLE PRECISION   X( * )
00014 *     ..
00015 *
00016 *  Purpose
00017 *  =======
00018 *
00019 *  DLARFGP generates a real elementary reflector H of order n, such
00020 *  that
00021 *
00022 *        H * ( alpha ) = ( beta ),   H' * H = I.
00023 *            (   x   )   (   0  )
00024 *
00025 *  where alpha and beta are scalars, beta is non-negative, and x is
00026 *  an (n-1)-element real vector.  H is represented in the form
00027 *
00028 *        H = I - tau * ( 1 ) * ( 1 v' ) ,
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 *  Arguments
00038 *  =========
00039 *
00040 *  N       (input) INTEGER
00041 *          The order of the elementary reflector.
00042 *
00043 *  ALPHA   (input/output) DOUBLE PRECISION
00044 *          On entry, the value alpha.
00045 *          On exit, it is overwritten with the value beta.
00046 *
00047 *  X       (input/output) DOUBLE PRECISION array, dimension
00048 *                         (1+(N-2)*abs(INCX))
00049 *          On entry, the vector x.
00050 *          On exit, it is overwritten with the vector v.
00051 *
00052 *  INCX    (input) INTEGER
00053 *          The increment between elements of X. INCX > 0.
00054 *
00055 *  TAU     (output) DOUBLE PRECISION
00056 *          The value tau.
00057 *
00058 *  =====================================================================
00059 *
00060 *     .. Parameters ..
00061       DOUBLE PRECISION   TWO, ONE, ZERO
00062       PARAMETER          ( TWO = 2.0D+0, ONE = 1.0D+0, ZERO = 0.0D+0 )
00063 *     ..
00064 *     .. Local Scalars ..
00065       INTEGER            J, KNT
00066       DOUBLE PRECISION   BETA, BIGNUM, SAVEALPHA, SMLNUM, XNORM
00067 *     ..
00068 *     .. External Functions ..
00069       DOUBLE PRECISION   DLAMCH, DLAPY2, DNRM2
00070       EXTERNAL           DLAMCH, DLAPY2, DNRM2
00071 *     ..
00072 *     .. Intrinsic Functions ..
00073       INTRINSIC          ABS, SIGN
00074 *     ..
00075 *     .. External Subroutines ..
00076       EXTERNAL           DSCAL
00077 *     ..
00078 *     .. Executable Statements ..
00079 *
00080       IF( N.LE.0 ) THEN
00081          TAU = ZERO
00082          RETURN
00083       END IF
00084 *
00085       XNORM = DNRM2( N-1, X, INCX )
00086 *
00087       IF( XNORM.EQ.ZERO ) THEN
00088 *
00089 *        H  =  [+/-1, 0; I], sign chosen so ALPHA >= 0
00090 *
00091          IF( ALPHA.GE.ZERO ) THEN
00092 *           When TAU.eq.ZERO, the vector is special-cased to be
00093 *           all zeros in the application routines.  We do not need
00094 *           to clear it.
00095             TAU = ZERO
00096          ELSE
00097 *           However, the application routines rely on explicit
00098 *           zero checks when TAU.ne.ZERO, and we must clear X.
00099             TAU = TWO
00100             DO J = 1, N-1
00101                X( 1 + (J-1)*INCX ) = 0
00102             END DO
00103             ALPHA = -ALPHA
00104          END IF
00105       ELSE
00106 *
00107 *        general case
00108 *
00109          BETA = SIGN( DLAPY2( ALPHA, XNORM ), ALPHA )
00110          SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'E' )
00111          KNT = 0
00112          IF( ABS( BETA ).LT.SMLNUM ) THEN
00113 *
00114 *           XNORM, BETA may be inaccurate; scale X and recompute them
00115 *
00116             BIGNUM = ONE / SMLNUM
00117    10       CONTINUE
00118             KNT = KNT + 1
00119             CALL DSCAL( N-1, BIGNUM, X, INCX )
00120             BETA = BETA*BIGNUM
00121             ALPHA = ALPHA*BIGNUM
00122             IF( ABS( BETA ).LT.SMLNUM )
00123      $         GO TO 10
00124 *
00125 *           New BETA is at most 1, at least SMLNUM
00126 *
00127             XNORM = DNRM2( N-1, X, INCX )
00128             BETA = SIGN( DLAPY2( ALPHA, XNORM ), ALPHA )
00129          END IF
00130          SAVEALPHA = ALPHA
00131          ALPHA = ALPHA + BETA
00132          IF( BETA.LT.ZERO ) THEN
00133             BETA = -BETA
00134             TAU = -ALPHA / BETA
00135          ELSE
00136             ALPHA = XNORM * (XNORM/ALPHA)
00137             TAU = ALPHA / BETA
00138             ALPHA = -ALPHA
00139          END IF
00140 *
00141          IF ( ABS(TAU).LE.SMLNUM ) THEN
00142 *
00143 *           In the case where the computed TAU ends up being a denormalized number,
00144 *           it loses relative accuracy. This is a BIG problem. Solution: flush TAU 
00145 *           to ZERO. This explains the next IF statement.
00146 *
00147 *           (Bug report provided by Pat Quillen from MathWorks on Jul 29, 2009.)
00148 *           (Thanks Pat. Thanks MathWorks.)
00149 *
00150             IF( SAVEALPHA.GE.ZERO ) THEN
00151                TAU = ZERO
00152             ELSE
00153                TAU = TWO
00154                DO J = 1, N-1
00155                   X( 1 + (J-1)*INCX ) = 0
00156                END DO
00157                BETA = -SAVEALPHA
00158             END IF
00159 *
00160          ELSE 
00161 *
00162 *           This is the general case.
00163 *
00164             CALL DSCAL( N-1, ONE / ALPHA, X, INCX )
00165 *
00166          END IF
00167 *
00168 *        If BETA is subnormal, it may lose relative accuracy
00169 *
00170          DO 20 J = 1, KNT
00171             BETA = BETA*SMLNUM
00172  20      CONTINUE
00173          ALPHA = BETA
00174       END IF
00175 *
00176       RETURN
00177 *
00178 *     End of DLARFGP
00179 *
00180       END