![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SLARFGP 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download SLARFGP + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slarfgp.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slarfgp.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slarfgp.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE SLARFGP( N, ALPHA, X, INCX, TAU ) 00022 * 00023 * .. Scalar Arguments .. 00024 * INTEGER INCX, N 00025 * REAL ALPHA, TAU 00026 * .. 00027 * .. Array Arguments .. 00028 * REAL X( * ) 00029 * .. 00030 * 00031 * 00032 *> \par Purpose: 00033 * ============= 00034 *> 00035 *> \verbatim 00036 *> 00037 *> SLARFGP generates a real elementary reflector H of order n, such 00038 *> that 00039 *> 00040 *> H * ( alpha ) = ( beta ), H**T * H = I. 00041 *> ( x ) ( 0 ) 00042 *> 00043 *> where alpha and beta are scalars, beta is non-negative, and x is 00044 *> an (n-1)-element real vector. H is represented in the form 00045 *> 00046 *> H = I - tau * ( 1 ) * ( 1 v**T ) , 00047 *> ( v ) 00048 *> 00049 *> where tau is a real scalar and v is a real (n-1)-element 00050 *> vector. 00051 *> 00052 *> If the elements of x are all zero, then tau = 0 and H is taken to be 00053 *> 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 REAL 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 REAL 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 REAL 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 realOTHERauxiliary 00103 * 00104 * ===================================================================== 00105 SUBROUTINE SLARFGP( 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 REAL ALPHA, TAU 00115 * .. 00116 * .. Array Arguments .. 00117 REAL 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 BETA, BIGNUM, SAVEALPHA, SMLNUM, XNORM 00129 * .. 00130 * .. External Functions .. 00131 REAL SLAMCH, SLAPY2, SNRM2 00132 EXTERNAL SLAMCH, SLAPY2, SNRM2 00133 * .. 00134 * .. Intrinsic Functions .. 00135 INTRINSIC ABS, SIGN 00136 * .. 00137 * .. External Subroutines .. 00138 EXTERNAL SSCAL 00139 * .. 00140 * .. Executable Statements .. 00141 * 00142 IF( N.LE.0 ) THEN 00143 TAU = ZERO 00144 RETURN 00145 END IF 00146 * 00147 XNORM = SNRM2( N-1, X, INCX ) 00148 * 00149 IF( XNORM.EQ.ZERO ) THEN 00150 * 00151 * H = [+/-1, 0; I], sign chosen so ALPHA >= 0. 00152 * 00153 IF( ALPHA.GE.ZERO ) THEN 00154 * When TAU.eq.ZERO, the vector is special-cased to be 00155 * all zeros in the application routines. We do not need 00156 * to clear it. 00157 TAU = ZERO 00158 ELSE 00159 * However, the application routines rely on explicit 00160 * zero checks when TAU.ne.ZERO, and we must clear X. 00161 TAU = TWO 00162 DO J = 1, N-1 00163 X( 1 + (J-1)*INCX ) = 0 00164 END DO 00165 ALPHA = -ALPHA 00166 END IF 00167 ELSE 00168 * 00169 * general case 00170 * 00171 BETA = SIGN( SLAPY2( ALPHA, XNORM ), ALPHA ) 00172 SMLNUM = SLAMCH( 'S' ) / SLAMCH( 'E' ) 00173 KNT = 0 00174 IF( ABS( BETA ).LT.SMLNUM ) THEN 00175 * 00176 * XNORM, BETA may be inaccurate; scale X and recompute them 00177 * 00178 BIGNUM = ONE / SMLNUM 00179 10 CONTINUE 00180 KNT = KNT + 1 00181 CALL SSCAL( N-1, BIGNUM, X, INCX ) 00182 BETA = BETA*BIGNUM 00183 ALPHA = ALPHA*BIGNUM 00184 IF( ABS( BETA ).LT.SMLNUM ) 00185 $ GO TO 10 00186 * 00187 * New BETA is at most 1, at least SMLNUM 00188 * 00189 XNORM = SNRM2( N-1, X, INCX ) 00190 BETA = SIGN( SLAPY2( ALPHA, XNORM ), ALPHA ) 00191 END IF 00192 SAVEALPHA = ALPHA 00193 ALPHA = ALPHA + BETA 00194 IF( BETA.LT.ZERO ) THEN 00195 BETA = -BETA 00196 TAU = -ALPHA / BETA 00197 ELSE 00198 ALPHA = XNORM * (XNORM/ALPHA) 00199 TAU = ALPHA / BETA 00200 ALPHA = -ALPHA 00201 END IF 00202 * 00203 IF ( ABS(TAU).LE.SMLNUM ) THEN 00204 * 00205 * In the case where the computed TAU ends up being a denormalized number, 00206 * it loses relative accuracy. This is a BIG problem. Solution: flush TAU 00207 * to ZERO. This explains the next IF statement. 00208 * 00209 * (Bug report provided by Pat Quillen from MathWorks on Jul 29, 2009.) 00210 * (Thanks Pat. Thanks MathWorks.) 00211 * 00212 IF( SAVEALPHA.GE.ZERO ) THEN 00213 TAU = ZERO 00214 ELSE 00215 TAU = TWO 00216 DO J = 1, N-1 00217 X( 1 + (J-1)*INCX ) = 0 00218 END DO 00219 BETA = -SAVEALPHA 00220 END IF 00221 * 00222 ELSE 00223 * 00224 * This is the general case. 00225 * 00226 CALL SSCAL( N-1, ONE / ALPHA, X, INCX ) 00227 * 00228 END IF 00229 * 00230 * If BETA is subnormal, it may lose relative accuracy 00231 * 00232 DO 20 J = 1, KNT 00233 BETA = BETA*SMLNUM 00234 20 CONTINUE 00235 ALPHA = BETA 00236 END IF 00237 * 00238 RETURN 00239 * 00240 * End of SLARFGP 00241 * 00242 END