 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE ZLARFGP( N, ALPHA, X, INCX, TAU ) 00002 * 00003 * -- LAPACK auxiliary routine (version 3.3.1) -- 00004 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00005 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00006 * -- April 2011 -- 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**H * ( alpha ) = ( beta ), H**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**H ) , 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
1.7.3 
