diff options
| author | jason <jason@8a072113-8704-0410-8d35-dd094bca7971> | 2008-10-28 01:38:50 +0000 |
|---|---|---|
| committer | jason <jason@8a072113-8704-0410-8d35-dd094bca7971> | 2008-10-28 01:38:50 +0000 |
| commit | baba851215b44ac3b60b9248eb02bcce7eb76247 (patch) | |
| tree | 8c0f5c006875532a30d4409f5e94b0f310ff00a7 /SRC/clarfg.f | |
| download | lapack-baba851215b44ac3b60b9248eb02bcce7eb76247.tar.gz lapack-baba851215b44ac3b60b9248eb02bcce7eb76247.tar.bz2 lapack-baba851215b44ac3b60b9248eb02bcce7eb76247.zip | |
Move LAPACK trunk into position.
Diffstat (limited to 'SRC/clarfg.f')
| -rw-r--r-- | SRC/clarfg.f | 140 |
1 files changed, 140 insertions, 0 deletions
diff --git a/SRC/clarfg.f b/SRC/clarfg.f new file mode 100644 index 00000000..8867f54b --- /dev/null +++ b/SRC/clarfg.f @@ -0,0 +1,140 @@ + SUBROUTINE CLARFG( N, ALPHA, X, INCX, TAU ) +* +* -- LAPACK auxiliary routine (version 3.1) -- +* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. +* November 2006 +* +* .. Scalar Arguments .. + INTEGER INCX, N + COMPLEX ALPHA, TAU +* .. +* .. Array Arguments .. + COMPLEX X( * ) +* .. +* +* Purpose +* ======= +* +* CLARFG generates a complex elementary reflector H of order n, such +* that +* +* H' * ( alpha ) = ( beta ), H' * H = I. +* ( x ) ( 0 ) +* +* where alpha and beta are scalars, with beta real, and x is an +* (n-1)-element complex vector. H is represented in the form +* +* H = I - tau * ( 1 ) * ( 1 v' ) , +* ( v ) +* +* where tau is a complex scalar and v is a complex (n-1)-element +* vector. Note that H is not hermitian. +* +* If the elements of x are all zero and alpha is real, then tau = 0 +* and H is taken to be the unit matrix. +* +* Otherwise 1 <= real(tau) <= 2 and abs(tau-1) <= 1 . +* +* Arguments +* ========= +* +* N (input) INTEGER +* The order of the elementary reflector. +* +* ALPHA (input/output) COMPLEX +* On entry, the value alpha. +* On exit, it is overwritten with the value beta. +* +* X (input/output) COMPLEX array, dimension +* (1+(N-2)*abs(INCX)) +* On entry, the vector x. +* On exit, it is overwritten with the vector v. +* +* INCX (input) INTEGER +* The increment between elements of X. INCX > 0. +* +* TAU (output) COMPLEX +* The value tau. +* +* ===================================================================== +* +* .. Parameters .. + REAL ONE, ZERO + PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) +* .. +* .. Local Scalars .. + INTEGER J, KNT + REAL ALPHI, ALPHR, BETA, RSAFMN, SAFMIN, XNORM +* .. +* .. External Functions .. + REAL SCNRM2, SLAMCH, SLAPY3 + COMPLEX CLADIV + EXTERNAL SCNRM2, SLAMCH, SLAPY3, CLADIV +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS, AIMAG, CMPLX, REAL, SIGN +* .. +* .. External Subroutines .. + EXTERNAL CSCAL, CSSCAL +* .. +* .. Executable Statements .. +* + IF( N.LE.0 ) THEN + TAU = ZERO + RETURN + END IF +* + XNORM = SCNRM2( N-1, X, INCX ) + ALPHR = REAL( ALPHA ) + ALPHI = AIMAG( ALPHA ) +* + IF( XNORM.EQ.ZERO .AND. ALPHI.EQ.ZERO ) THEN +* +* H = I +* + TAU = ZERO + ELSE +* +* general case +* + BETA = -SIGN( SLAPY3( ALPHR, ALPHI, XNORM ), ALPHR ) + SAFMIN = SLAMCH( 'S' ) / SLAMCH( 'E' ) + RSAFMN = ONE / SAFMIN +* + KNT = 0 + IF( ABS( BETA ).LT.SAFMIN ) THEN +* +* XNORM, BETA may be inaccurate; scale X and recompute them +* + 10 CONTINUE + KNT = KNT + 1 + CALL CSSCAL( N-1, RSAFMN, X, INCX ) + BETA = BETA*RSAFMN + ALPHI = ALPHI*RSAFMN + ALPHR = ALPHR*RSAFMN + IF( ABS( BETA ).LT.SAFMIN ) + $ GO TO 10 +* +* New BETA is at most 1, at least SAFMIN +* + XNORM = SCNRM2( N-1, X, INCX ) + ALPHA = CMPLX( ALPHR, ALPHI ) + BETA = -SIGN( SLAPY3( ALPHR, ALPHI, XNORM ), ALPHR ) + END IF + TAU = CMPLX( ( BETA-ALPHR ) / BETA, -ALPHI / BETA ) + ALPHA = CLADIV( CMPLX( ONE ), ALPHA-BETA ) + CALL CSCAL( N-1, ALPHA, X, INCX ) +* +* If ALPHA is subnormal, it may lose relative accuracy +* + DO 20 J = 1, KNT + BETA = BETA*SAFMIN + 20 CONTINUE + ALPHA = BETA + END IF +* + RETURN +* +* End of CLARFG +* + END |