 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE SLARFGP( 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 REAL ALPHA, TAU 00011 * .. 00012 * .. Array Arguments .. 00013 REAL X( * ) 00014 * .. 00015 * 00016 * Purpose 00017 * ======= 00018 * 00019 * SLARFGP generates a real elementary reflector H of order n, such 00020 * that 00021 * 00022 * H * ( alpha ) = ( beta ), H**T * H = I. 00023 * ( x ) ( 0 ) 00024 * 00025 * where alpha and beta are scalars, beta is non-negative, and x is 00026 * an (n-1)-element real vector. H is represented in the form 00027 * 00028 * H = I - tau * ( 1 ) * ( 1 v**T ) , 00029 * ( v ) 00030 * 00031 * where tau is a real scalar and v is a real (n-1)-element 00032 * vector. 00033 * 00034 * If the elements of x are all zero, then tau = 0 and H is taken to be 00035 * the unit matrix. 00036 * 00037 * Arguments 00038 * ========= 00039 * 00040 * N (input) INTEGER 00041 * The order of the elementary reflector. 00042 * 00043 * ALPHA (input/output) REAL 00044 * On entry, the value alpha. 00045 * On exit, it is overwritten with the value beta. 00046 * 00047 * X (input/output) REAL 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) REAL 00056 * The value tau. 00057 * 00058 * ===================================================================== 00059 * 00060 * .. Parameters .. 00061 REAL TWO, ONE, ZERO 00062 PARAMETER ( TWO = 2.0E+0, ONE = 1.0E+0, ZERO = 0.0E+0 ) 00063 * .. 00064 * .. Local Scalars .. 00065 INTEGER J, KNT 00066 REAL BETA, BIGNUM, SAVEALPHA, SMLNUM, XNORM 00067 * .. 00068 * .. External Functions .. 00069 REAL SLAMCH, SLAPY2, SNRM2 00070 EXTERNAL SLAMCH, SLAPY2, SNRM2 00071 * .. 00072 * .. Intrinsic Functions .. 00073 INTRINSIC ABS, SIGN 00074 * .. 00075 * .. External Subroutines .. 00076 EXTERNAL SSCAL 00077 * .. 00078 * .. Executable Statements .. 00079 * 00080 IF( N.LE.0 ) THEN 00081 TAU = ZERO 00082 RETURN 00083 END IF 00084 * 00085 XNORM = SNRM2( N-1, X, INCX ) 00086 * 00087 IF( XNORM.EQ.ZERO ) THEN 00088 * 00089 * H = [+/-1, 0; I], sign chosen so ALPHA >= 0. 00090 * 00091 IF( ALPHA.GE.ZERO ) THEN 00092 * When TAU.eq.ZERO, the vector is special-cased to be 00093 * all zeros in the application routines. We do not need 00094 * to clear it. 00095 TAU = ZERO 00096 ELSE 00097 * However, the application routines rely on explicit 00098 * zero checks when TAU.ne.ZERO, and we must clear X. 00099 TAU = TWO 00100 DO J = 1, N-1 00101 X( 1 + (J-1)*INCX ) = 0 00102 END DO 00103 ALPHA = -ALPHA 00104 END IF 00105 ELSE 00106 * 00107 * general case 00108 * 00109 BETA = SIGN( SLAPY2( ALPHA, XNORM ), ALPHA ) 00110 SMLNUM = SLAMCH( 'S' ) / SLAMCH( 'E' ) 00111 KNT = 0 00112 IF( ABS( BETA ).LT.SMLNUM ) THEN 00113 * 00114 * XNORM, BETA may be inaccurate; scale X and recompute them 00115 * 00116 BIGNUM = ONE / SMLNUM 00117 10 CONTINUE 00118 KNT = KNT + 1 00119 CALL SSCAL( N-1, BIGNUM, X, INCX ) 00120 BETA = BETA*BIGNUM 00121 ALPHA = ALPHA*BIGNUM 00122 IF( ABS( BETA ).LT.SMLNUM ) 00123 $ GO TO 10 00124 * 00125 * New BETA is at most 1, at least SMLNUM 00126 * 00127 XNORM = SNRM2( N-1, X, INCX ) 00128 BETA = SIGN( SLAPY2( ALPHA, XNORM ), ALPHA ) 00129 END IF 00130 SAVEALPHA = ALPHA 00131 ALPHA = ALPHA + BETA 00132 IF( BETA.LT.ZERO ) THEN 00133 BETA = -BETA 00134 TAU = -ALPHA / BETA 00135 ELSE 00136 ALPHA = XNORM * (XNORM/ALPHA) 00137 TAU = ALPHA / BETA 00138 ALPHA = -ALPHA 00139 END IF 00140 * 00141 IF ( ABS(TAU).LE.SMLNUM ) THEN 00142 * 00143 * In the case where the computed TAU ends up being a denormalized number, 00144 * it loses relative accuracy. This is a BIG problem. Solution: flush TAU 00145 * to ZERO. This explains the next IF statement. 00146 * 00147 * (Bug report provided by Pat Quillen from MathWorks on Jul 29, 2009.) 00148 * (Thanks Pat. Thanks MathWorks.) 00149 * 00150 IF( SAVEALPHA.GE.ZERO ) THEN 00151 TAU = ZERO 00152 ELSE 00153 TAU = TWO 00154 DO J = 1, N-1 00155 X( 1 + (J-1)*INCX ) = 0 00156 END DO 00157 BETA = -SAVEALPHA 00158 END IF 00159 * 00160 ELSE 00161 * 00162 * This is the general case. 00163 * 00164 CALL SSCAL( N-1, ONE / ALPHA, X, INCX ) 00165 * 00166 END IF 00167 * 00168 * If BETA is subnormal, it may lose relative accuracy 00169 * 00170 DO 20 J = 1, KNT 00171 BETA = BETA*SMLNUM 00172 20 CONTINUE 00173 ALPHA = BETA 00174 END IF 00175 * 00176 RETURN 00177 * 00178 * End of SLARFGP 00179 * 00180 END
1.7.3 
