LAPACK: SRC/zlarfgp.f Source File

Go to the documentation of this file.
00001       SUBROUTINE ZLARFGP( 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       COMPLEX*16         ALPHA, TAU
00011 *     ..
00012 *     .. Array Arguments ..
00013       COMPLEX*16         X( * )
00014 *     ..
00015 *
00016 *  Purpose
00017 *  =======
00018 *
00019 *  ZLARFGP generates a complex 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 real and non-negative, and
00026 *  x is an (n-1)-element complex vector.  H is represented in the form
00027 *
00028 *        H = I - tau * ( 1 ) * ( 1 v' ) ,
00029 *                      ( v )
00030 *
00031 *  where tau is a complex scalar and v is a complex (n-1)-element
00032 *  vector. Note that H is not hermitian.
00033 *
00034 *  If the elements of x are all zero and alpha is real, then tau = 0
00035 *  and H is taken to be the unit matrix.
00036 *
00037 *  Arguments
00038 *  =========
00039 *
00040 *  N       (input) INTEGER
00041 *          The order of the elementary reflector.
00042 *
00043 *  ALPHA   (input/output) COMPLEX*16
00044 *          On entry, the value alpha.
00045 *          On exit, it is overwritten with the value beta.
00046 *
00047 *  X       (input/output) COMPLEX*16 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) COMPLEX*16
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   ALPHI, ALPHR, BETA, BIGNUM, SMLNUM, XNORM
00067       COMPLEX*16         SAVEALPHA
00068 *     ..
00069 *     .. External Functions ..
00070       DOUBLE PRECISION   DLAMCH, DLAPY3, DLAPY2, DZNRM2
00071       COMPLEX*16         ZLADIV
00072       EXTERNAL           DLAMCH, DLAPY3, DLAPY2, DZNRM2, ZLADIV
00073 *     ..
00074 *     .. Intrinsic Functions ..
00075       INTRINSIC          ABS, DBLE, DCMPLX, DIMAG, SIGN
00076 *     ..
00077 *     .. External Subroutines ..
00078       EXTERNAL           ZDSCAL, ZSCAL
00079 *     ..
00080 *     .. Executable Statements ..
00081 *
00082       IF( N.LE.0 ) THEN
00083          TAU = ZERO
00084          RETURN
00085       END IF
00086 *
00087       XNORM = DZNRM2( N-1, X, INCX )
00088       ALPHR = DBLE( ALPHA )
00089       ALPHI = DIMAG( ALPHA )
00090 *
00091       IF( XNORM.EQ.ZERO ) THEN
00092 *
00093 *        H  =  [1-alpha/abs(alpha) 0; 0 I], sign chosen so ALPHA >= 0.
00094 *
00095          IF( ALPHI.EQ.ZERO ) THEN
00096             IF( ALPHR.GE.ZERO ) THEN
00097 *              When TAU.eq.ZERO, the vector is special-cased to be
00098 *              all zeros in the application routines.  We do not need
00099 *              to clear it.
00100                TAU = ZERO
00101             ELSE
00102 *              However, the application routines rely on explicit
00103 *              zero checks when TAU.ne.ZERO, and we must clear X.
00104                TAU = TWO
00105                DO J = 1, N-1
00106                   X( 1 + (J-1)*INCX ) = ZERO
00107                END DO
00108                ALPHA = -ALPHA
00109             END IF
00110          ELSE
00111 *           Only "reflecting" the diagonal entry to be real and non-negative.
00112             XNORM = DLAPY2( ALPHR, ALPHI )
00113             TAU = DCMPLX( ONE - ALPHR / XNORM, -ALPHI / XNORM )
00114             DO J = 1, N-1
00115                X( 1 + (J-1)*INCX ) = ZERO
00116             END DO
00117             ALPHA = XNORM
00118          END IF
00119       ELSE
00120 *
00121 *        general case
00122 *
00123          BETA = SIGN( DLAPY3( ALPHR, ALPHI, XNORM ), ALPHR )
00124          SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'E' )
00125          BIGNUM = ONE / SMLNUM
00126 *
00127          KNT = 0
00128          IF( ABS( BETA ).LT.SMLNUM ) THEN
00129 *
00130 *           XNORM, BETA may be inaccurate; scale X and recompute them
00131 *
00132    10       CONTINUE
00133             KNT = KNT + 1
00134             CALL ZDSCAL( N-1, BIGNUM, X, INCX )
00135             BETA = BETA*BIGNUM
00136             ALPHI = ALPHI*BIGNUM
00137             ALPHR = ALPHR*BIGNUM
00138             IF( ABS( BETA ).LT.SMLNUM )
00139      $         GO TO 10
00140 *
00141 *           New BETA is at most 1, at least SMLNUM
00142 *
00143             XNORM = DZNRM2( N-1, X, INCX )
00144             ALPHA = DCMPLX( ALPHR, ALPHI )
00145             BETA = SIGN( DLAPY3( ALPHR, ALPHI, XNORM ), ALPHR )
00146          END IF
00147          SAVEALPHA = ALPHA
00148          ALPHA = ALPHA + BETA
00149          IF( BETA.LT.ZERO ) THEN
00150             BETA = -BETA
00151             TAU = -ALPHA / BETA
00152          ELSE
00153             ALPHR = ALPHI * (ALPHI/DBLE( ALPHA ))
00154             ALPHR = ALPHR + XNORM * (XNORM/DBLE( ALPHA ))
00155             TAU = DCMPLX( ALPHR/BETA, -ALPHI/BETA )
00156             ALPHA = DCMPLX( -ALPHR, ALPHI )
00157          END IF
00158          ALPHA = ZLADIV( DCMPLX( ONE ), ALPHA )
00159 *
00160          IF ( ABS(TAU).LE.SMLNUM ) THEN
00161 *
00162 *           In the case where the computed TAU ends up being a denormalized number,
00163 *           it loses relative accuracy. This is a BIG problem. Solution: flush TAU 
00164 *           to ZERO (or TWO or whatever makes a nonnegative real number for BETA).
00165 *
00166 *           (Bug report provided by Pat Quillen from MathWorks on Jul 29, 2009.)
00167 *           (Thanks Pat. Thanks MathWorks.)
00168 *
00169             ALPHR = DBLE( SAVEALPHA )
00170             ALPHI = DIMAG( SAVEALPHA )
00171             IF( ALPHI.EQ.ZERO ) THEN
00172                IF( ALPHR.GE.ZERO ) THEN
00173                   TAU = ZERO
00174                ELSE
00175                   TAU = TWO
00176                   DO J = 1, N-1
00177                      X( 1 + (J-1)*INCX ) = ZERO
00178                   END DO
00179                   BETA = -SAVEALPHA
00180                END IF
00181             ELSE
00182                XNORM = DLAPY2( ALPHR, ALPHI )
00183                TAU = DCMPLX( ONE - ALPHR / XNORM, -ALPHI / XNORM )
00184                DO J = 1, N-1
00185                   X( 1 + (J-1)*INCX ) = ZERO
00186                END DO
00187                BETA = XNORM
00188             END IF
00189 *
00190          ELSE 
00191 *
00192 *           This is the general case.
00193 *
00194             CALL ZSCAL( N-1, ALPHA, X, INCX )
00195 *
00196          END IF
00197 *
00198 *        If BETA is subnormal, it may lose relative accuracy
00199 *
00200          DO 20 J = 1, KNT
00201             BETA = BETA*SMLNUM
00202  20      CONTINUE
00203          ALPHA = BETA
00204       END IF
00205 *
00206       RETURN
00207 *
00208 *     End of ZLARFGP
00209 *
00210       END