![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b CLARFGP 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download CLARFGP + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clarfgp.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clarfgp.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clarfgp.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE CLARFGP( N, ALPHA, X, INCX, TAU ) 00022 * 00023 * .. Scalar Arguments .. 00024 * INTEGER INCX, N 00025 * COMPLEX ALPHA, TAU 00026 * .. 00027 * .. Array Arguments .. 00028 * COMPLEX X( * ) 00029 * .. 00030 * 00031 * 00032 *> \par Purpose: 00033 * ============= 00034 *> 00035 *> \verbatim 00036 *> 00037 *> CLARFGP generates a complex elementary reflector H of order n, such 00038 *> that 00039 *> 00040 *> H**H * ( alpha ) = ( beta ), H**H * H = I. 00041 *> ( x ) ( 0 ) 00042 *> 00043 *> where alpha and beta are scalars, beta is real and non-negative, and 00044 *> x is an (n-1)-element complex vector. H is represented in the form 00045 *> 00046 *> H = I - tau * ( 1 ) * ( 1 v**H ) , 00047 *> ( v ) 00048 *> 00049 *> where tau is a complex scalar and v is a complex (n-1)-element 00050 *> vector. Note that H is not hermitian. 00051 *> 00052 *> If the elements of x are all zero and alpha is real, then tau = 0 00053 *> and H is taken to be the unit matrix. 00054 *> \endverbatim 00055 * 00056 * Arguments: 00057 * ========== 00058 * 00059 *> \param[in] N 00060 *> \verbatim 00061 *> N is INTEGER 00062 *> The order of the elementary reflector. 00063 *> \endverbatim 00064 *> 00065 *> \param[in,out] ALPHA 00066 *> \verbatim 00067 *> ALPHA is COMPLEX 00068 *> On entry, the value alpha. 00069 *> On exit, it is overwritten with the value beta. 00070 *> \endverbatim 00071 *> 00072 *> \param[in,out] X 00073 *> \verbatim 00074 *> X is COMPLEX array, dimension 00075 *> (1+(N-2)*abs(INCX)) 00076 *> On entry, the vector x. 00077 *> On exit, it is overwritten with the vector v. 00078 *> \endverbatim 00079 *> 00080 *> \param[in] INCX 00081 *> \verbatim 00082 *> INCX is INTEGER 00083 *> The increment between elements of X. INCX > 0. 00084 *> \endverbatim 00085 *> 00086 *> \param[out] TAU 00087 *> \verbatim 00088 *> TAU is COMPLEX 00089 *> The value tau. 00090 *> \endverbatim 00091 * 00092 * Authors: 00093 * ======== 00094 * 00095 *> \author Univ. of Tennessee 00096 *> \author Univ. of California Berkeley 00097 *> \author Univ. of Colorado Denver 00098 *> \author NAG Ltd. 00099 * 00100 *> \date November 2011 00101 * 00102 *> \ingroup complexOTHERauxiliary 00103 * 00104 * ===================================================================== 00105 SUBROUTINE CLARFGP( N, ALPHA, X, INCX, TAU ) 00106 * 00107 * -- LAPACK auxiliary routine (version 3.4.0) -- 00108 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00109 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00110 * November 2011 00111 * 00112 * .. Scalar Arguments .. 00113 INTEGER INCX, N 00114 COMPLEX ALPHA, TAU 00115 * .. 00116 * .. Array Arguments .. 00117 COMPLEX X( * ) 00118 * .. 00119 * 00120 * ===================================================================== 00121 * 00122 * .. Parameters .. 00123 REAL TWO, ONE, ZERO 00124 PARAMETER ( TWO = 2.0E+0, ONE = 1.0E+0, ZERO = 0.0E+0 ) 00125 * .. 00126 * .. Local Scalars .. 00127 INTEGER J, KNT 00128 REAL ALPHI, ALPHR, BETA, BIGNUM, SMLNUM, XNORM 00129 COMPLEX SAVEALPHA 00130 * .. 00131 * .. External Functions .. 00132 REAL SCNRM2, SLAMCH, SLAPY3, SLAPY2 00133 COMPLEX CLADIV 00134 EXTERNAL SCNRM2, SLAMCH, SLAPY3, SLAPY2, CLADIV 00135 * .. 00136 * .. Intrinsic Functions .. 00137 INTRINSIC ABS, AIMAG, CMPLX, REAL, SIGN 00138 * .. 00139 * .. External Subroutines .. 00140 EXTERNAL CSCAL, CSSCAL 00141 * .. 00142 * .. Executable Statements .. 00143 * 00144 IF( N.LE.0 ) THEN 00145 TAU = ZERO 00146 RETURN 00147 END IF 00148 * 00149 XNORM = SCNRM2( N-1, X, INCX ) 00150 ALPHR = REAL( ALPHA ) 00151 ALPHI = AIMAG( ALPHA ) 00152 * 00153 IF( XNORM.EQ.ZERO ) THEN 00154 * 00155 * H = [1-alpha/abs(alpha) 0; 0 I], sign chosen so ALPHA >= 0. 00156 * 00157 IF( ALPHI.EQ.ZERO ) THEN 00158 IF( ALPHR.GE.ZERO ) THEN 00159 * When TAU.eq.ZERO, the vector is special-cased to be 00160 * all zeros in the application routines. We do not need 00161 * to clear it. 00162 TAU = ZERO 00163 ELSE 00164 * However, the application routines rely on explicit 00165 * zero checks when TAU.ne.ZERO, and we must clear X. 00166 TAU = TWO 00167 DO J = 1, N-1 00168 X( 1 + (J-1)*INCX ) = ZERO 00169 END DO 00170 ALPHA = -ALPHA 00171 END IF 00172 ELSE 00173 * Only "reflecting" the diagonal entry to be real and non-negative. 00174 XNORM = SLAPY2( ALPHR, ALPHI ) 00175 TAU = CMPLX( ONE - ALPHR / XNORM, -ALPHI / XNORM ) 00176 DO J = 1, N-1 00177 X( 1 + (J-1)*INCX ) = ZERO 00178 END DO 00179 ALPHA = XNORM 00180 END IF 00181 ELSE 00182 * 00183 * general case 00184 * 00185 BETA = SIGN( SLAPY3( ALPHR, ALPHI, XNORM ), ALPHR ) 00186 SMLNUM = SLAMCH( 'S' ) / SLAMCH( 'E' ) 00187 BIGNUM = ONE / SMLNUM 00188 * 00189 KNT = 0 00190 IF( ABS( BETA ).LT.SMLNUM ) THEN 00191 * 00192 * XNORM, BETA may be inaccurate; scale X and recompute them 00193 * 00194 10 CONTINUE 00195 KNT = KNT + 1 00196 CALL CSSCAL( N-1, BIGNUM, X, INCX ) 00197 BETA = BETA*BIGNUM 00198 ALPHI = ALPHI*BIGNUM 00199 ALPHR = ALPHR*BIGNUM 00200 IF( ABS( BETA ).LT.SMLNUM ) 00201 $ GO TO 10 00202 * 00203 * New BETA is at most 1, at least SMLNUM 00204 * 00205 XNORM = SCNRM2( N-1, X, INCX ) 00206 ALPHA = CMPLX( ALPHR, ALPHI ) 00207 BETA = SIGN( SLAPY3( ALPHR, ALPHI, XNORM ), ALPHR ) 00208 END IF 00209 SAVEALPHA = ALPHA 00210 ALPHA = ALPHA + BETA 00211 IF( BETA.LT.ZERO ) THEN 00212 BETA = -BETA 00213 TAU = -ALPHA / BETA 00214 ELSE 00215 ALPHR = ALPHI * (ALPHI/REAL( ALPHA )) 00216 ALPHR = ALPHR + XNORM * (XNORM/REAL( ALPHA )) 00217 TAU = CMPLX( ALPHR/BETA, -ALPHI/BETA ) 00218 ALPHA = CMPLX( -ALPHR, ALPHI ) 00219 END IF 00220 ALPHA = CLADIV( CMPLX( ONE ), ALPHA ) 00221 * 00222 IF ( ABS(TAU).LE.SMLNUM ) THEN 00223 * 00224 * In the case where the computed TAU ends up being a denormalized number, 00225 * it loses relative accuracy. This is a BIG problem. Solution: flush TAU 00226 * to ZERO (or TWO or whatever makes a nonnegative real number for BETA). 00227 * 00228 * (Bug report provided by Pat Quillen from MathWorks on Jul 29, 2009.) 00229 * (Thanks Pat. Thanks MathWorks.) 00230 * 00231 ALPHR = REAL( SAVEALPHA ) 00232 ALPHI = AIMAG( SAVEALPHA ) 00233 IF( ALPHI.EQ.ZERO ) THEN 00234 IF( ALPHR.GE.ZERO ) THEN 00235 TAU = ZERO 00236 ELSE 00237 TAU = TWO 00238 DO J = 1, N-1 00239 X( 1 + (J-1)*INCX ) = ZERO 00240 END DO 00241 BETA = -SAVEALPHA 00242 END IF 00243 ELSE 00244 XNORM = SLAPY2( ALPHR, ALPHI ) 00245 TAU = CMPLX( ONE - ALPHR / XNORM, -ALPHI / XNORM ) 00246 DO J = 1, N-1 00247 X( 1 + (J-1)*INCX ) = ZERO 00248 END DO 00249 BETA = XNORM 00250 END IF 00251 * 00252 ELSE 00253 * 00254 * This is the general case. 00255 * 00256 CALL CSCAL( N-1, ALPHA, X, INCX ) 00257 * 00258 END IF 00259 * 00260 * If BETA is subnormal, it may lose relative accuracy 00261 * 00262 DO 20 J = 1, KNT 00263 BETA = BETA*SMLNUM 00264 20 CONTINUE 00265 ALPHA = BETA 00266 END IF 00267 * 00268 RETURN 00269 * 00270 * End of CLARFGP 00271 * 00272 END