|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE ZLARFG( 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*16 ALPHA, TAU 00011 * .. 00012 * .. Array Arguments .. 00013 COMPLEX*16 X( * ) 00014 * .. 00015 * 00016 * Purpose 00017 * ======= 00018 * 00019 * ZLARFG 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*16 00046 * On entry, the value alpha. 00047 * On exit, it is overwritten with the value beta. 00048 * 00049 * X (input/output) COMPLEX*16 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*16 00058 * The value tau. 00059 * 00060 * ===================================================================== 00061 * 00062 * .. Parameters .. 00063 DOUBLE PRECISION ONE, ZERO 00064 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 00065 * .. 00066 * .. Local Scalars .. 00067 INTEGER J, KNT 00068 DOUBLE PRECISION ALPHI, ALPHR, BETA, RSAFMN, SAFMIN, XNORM 00069 * .. 00070 * .. External Functions .. 00071 DOUBLE PRECISION DLAMCH, DLAPY3, DZNRM2 00072 COMPLEX*16 ZLADIV 00073 EXTERNAL DLAMCH, DLAPY3, DZNRM2, ZLADIV 00074 * .. 00075 * .. Intrinsic Functions .. 00076 INTRINSIC ABS, DBLE, DCMPLX, DIMAG, SIGN 00077 * .. 00078 * .. External Subroutines .. 00079 EXTERNAL ZDSCAL, ZSCAL 00080 * .. 00081 * .. Executable Statements .. 00082 * 00083 IF( N.LE.0 ) THEN 00084 TAU = ZERO 00085 RETURN 00086 END IF 00087 * 00088 XNORM = DZNRM2( N-1, X, INCX ) 00089 ALPHR = DBLE( ALPHA ) 00090 ALPHI = DIMAG( 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( DLAPY3( ALPHR, ALPHI, XNORM ), ALPHR ) 00102 SAFMIN = DLAMCH( 'S' ) / DLAMCH( '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 ZDSCAL( 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 = DZNRM2( N-1, X, INCX ) 00122 ALPHA = DCMPLX( ALPHR, ALPHI ) 00123 BETA = -SIGN( DLAPY3( ALPHR, ALPHI, XNORM ), ALPHR ) 00124 END IF 00125 TAU = DCMPLX( ( BETA-ALPHR ) / BETA, -ALPHI / BETA ) 00126 ALPHA = ZLADIV( DCMPLX( ONE ), ALPHA-BETA ) 00127 CALL ZSCAL( 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 ZLARFG 00140 * 00141 END
1.7.3 
