LAPACK 3.3.0

clarfg.f

Go to the documentation of this file.
00001       SUBROUTINE CLARFG( N, ALPHA, X, INCX, TAU )
00002 *
00003 *  -- LAPACK auxiliary routine (version 3.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 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       INTEGER            INCX, N
00010       COMPLEX            ALPHA, TAU
00011 *     ..
00012 *     .. Array Arguments ..
00013       COMPLEX            X( * )
00014 *     ..
00015 *
00016 *  Purpose
00017 *  =======
00018 *
00019 *  CLARFG 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, with beta real, and x is an
00026 *  (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 *  Otherwise  1 <= real(tau) <= 2  and  abs(tau-1) <= 1 .
00038 *
00039 *  Arguments
00040 *  =========
00041 *
00042 *  N       (input) INTEGER
00043 *          The order of the elementary reflector.
00044 *
00045 *  ALPHA   (input/output) COMPLEX
00046 *          On entry, the value alpha.
00047 *          On exit, it is overwritten with the value beta.
00048 *
00049 *  X       (input/output) COMPLEX 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) COMPLEX
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               ALPHI, ALPHR, BETA, RSAFMN, SAFMIN, XNORM
00069 *     ..
00070 *     .. External Functions ..
00071       REAL               SCNRM2, SLAMCH, SLAPY3
00072       COMPLEX            CLADIV
00073       EXTERNAL           SCNRM2, SLAMCH, SLAPY3, CLADIV
00074 *     ..
00075 *     .. Intrinsic Functions ..
00076       INTRINSIC          ABS, AIMAG, CMPLX, REAL, SIGN
00077 *     ..
00078 *     .. External Subroutines ..
00079       EXTERNAL           CSCAL, CSSCAL
00080 *     ..
00081 *     .. Executable Statements ..
00082 *
00083       IF( N.LE.0 ) THEN
00084          TAU = ZERO
00085          RETURN
00086       END IF
00087 *
00088       XNORM = SCNRM2( N-1, X, INCX )
00089       ALPHR = REAL( ALPHA )
00090       ALPHI = AIMAG( ALPHA )
00091 *
00092       IF( XNORM.EQ.ZERO .AND. ALPHI.EQ.ZERO ) THEN
00093 *
00094 *        H  =  I
00095 *
00096          TAU = ZERO
00097       ELSE
00098 *
00099 *        general case
00100 *
00101          BETA = -SIGN( SLAPY3( ALPHR, ALPHI, XNORM ), ALPHR )
00102          SAFMIN = SLAMCH( 'S' ) / SLAMCH( 'E' )
00103          RSAFMN = ONE / SAFMIN
00104 *
00105          KNT = 0
00106          IF( ABS( BETA ).LT.SAFMIN ) THEN
00107 *
00108 *           XNORM, BETA may be inaccurate; scale X and recompute them
00109 *
00110    10       CONTINUE
00111             KNT = KNT + 1
00112             CALL CSSCAL( N-1, RSAFMN, X, INCX )
00113             BETA = BETA*RSAFMN
00114             ALPHI = ALPHI*RSAFMN
00115             ALPHR = ALPHR*RSAFMN
00116             IF( ABS( BETA ).LT.SAFMIN )
00117      $         GO TO 10
00118 *
00119 *           New BETA is at most 1, at least SAFMIN
00120 *
00121             XNORM = SCNRM2( N-1, X, INCX )
00122             ALPHA = CMPLX( ALPHR, ALPHI )
00123             BETA = -SIGN( SLAPY3( ALPHR, ALPHI, XNORM ), ALPHR )
00124          END IF
00125          TAU = CMPLX( ( BETA-ALPHR ) / BETA, -ALPHI / BETA )
00126          ALPHA = CLADIV( CMPLX( ONE ), ALPHA-BETA )
00127          CALL CSCAL( N-1, ALPHA, X, INCX )
00128 *
00129 *        If ALPHA is subnormal, it may lose relative accuracy
00130 *
00131          DO 20 J = 1, KNT
00132             BETA = BETA*SAFMIN
00133  20      CONTINUE
00134          ALPHA = BETA
00135       END IF
00136 *
00137       RETURN
00138 *
00139 *     End of CLARFG
00140 *
00141       END
 All Files Functions