Annotation of rpl/lapack/lapack/dlarfgp.f, revision 1.6
1.6 ! bertrand 1: *> \brief \b DLARFGP
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download DLARFGP + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarfgp.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarfgp.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarfgp.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE DLARFGP( N, ALPHA, X, INCX, TAU )
! 22: *
! 23: * .. Scalar Arguments ..
! 24: * INTEGER INCX, N
! 25: * DOUBLE PRECISION ALPHA, TAU
! 26: * ..
! 27: * .. Array Arguments ..
! 28: * DOUBLE PRECISION X( * )
! 29: * ..
! 30: *
! 31: *
! 32: *> \par Purpose:
! 33: * =============
! 34: *>
! 35: *> \verbatim
! 36: *>
! 37: *> DLARFGP generates a real elementary reflector H of order n, such
! 38: *> that
! 39: *>
! 40: *> H * ( alpha ) = ( beta ), H**T * H = I.
! 41: *> ( x ) ( 0 )
! 42: *>
! 43: *> where alpha and beta are scalars, beta is non-negative, and x is
! 44: *> an (n-1)-element real vector. H is represented in the form
! 45: *>
! 46: *> H = I - tau * ( 1 ) * ( 1 v**T ) ,
! 47: *> ( v )
! 48: *>
! 49: *> where tau is a real scalar and v is a real (n-1)-element
! 50: *> vector.
! 51: *>
! 52: *> If the elements of x are all zero, then tau = 0 and H is taken to be
! 53: *> the unit matrix.
! 54: *> \endverbatim
! 55: *
! 56: * Arguments:
! 57: * ==========
! 58: *
! 59: *> \param[in] N
! 60: *> \verbatim
! 61: *> N is INTEGER
! 62: *> The order of the elementary reflector.
! 63: *> \endverbatim
! 64: *>
! 65: *> \param[in,out] ALPHA
! 66: *> \verbatim
! 67: *> ALPHA is DOUBLE PRECISION
! 68: *> On entry, the value alpha.
! 69: *> On exit, it is overwritten with the value beta.
! 70: *> \endverbatim
! 71: *>
! 72: *> \param[in,out] X
! 73: *> \verbatim
! 74: *> X is DOUBLE PRECISION array, dimension
! 75: *> (1+(N-2)*abs(INCX))
! 76: *> On entry, the vector x.
! 77: *> On exit, it is overwritten with the vector v.
! 78: *> \endverbatim
! 79: *>
! 80: *> \param[in] INCX
! 81: *> \verbatim
! 82: *> INCX is INTEGER
! 83: *> The increment between elements of X. INCX > 0.
! 84: *> \endverbatim
! 85: *>
! 86: *> \param[out] TAU
! 87: *> \verbatim
! 88: *> TAU is DOUBLE PRECISION
! 89: *> The value tau.
! 90: *> \endverbatim
! 91: *
! 92: * Authors:
! 93: * ========
! 94: *
! 95: *> \author Univ. of Tennessee
! 96: *> \author Univ. of California Berkeley
! 97: *> \author Univ. of Colorado Denver
! 98: *> \author NAG Ltd.
! 99: *
! 100: *> \date November 2011
! 101: *
! 102: *> \ingroup doubleOTHERauxiliary
! 103: *
! 104: * =====================================================================
1.1 bertrand 105: SUBROUTINE DLARFGP( N, ALPHA, X, INCX, TAU )
106: *
1.6 ! bertrand 107: * -- LAPACK auxiliary routine (version 3.4.0) --
1.1 bertrand 108: * -- LAPACK is a software package provided by Univ. of Tennessee, --
109: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.6 ! bertrand 110: * November 2011
1.1 bertrand 111: *
112: * .. Scalar Arguments ..
113: INTEGER INCX, N
114: DOUBLE PRECISION ALPHA, TAU
115: * ..
116: * .. Array Arguments ..
117: DOUBLE PRECISION X( * )
118: * ..
119: *
120: * =====================================================================
121: *
122: * .. Parameters ..
123: DOUBLE PRECISION TWO, ONE, ZERO
124: PARAMETER ( TWO = 2.0D+0, ONE = 1.0D+0, ZERO = 0.0D+0 )
125: * ..
126: * .. Local Scalars ..
127: INTEGER J, KNT
128: DOUBLE PRECISION BETA, BIGNUM, SAVEALPHA, SMLNUM, XNORM
129: * ..
130: * .. External Functions ..
131: DOUBLE PRECISION DLAMCH, DLAPY2, DNRM2
132: EXTERNAL DLAMCH, DLAPY2, DNRM2
133: * ..
134: * .. Intrinsic Functions ..
135: INTRINSIC ABS, SIGN
136: * ..
137: * .. External Subroutines ..
138: EXTERNAL DSCAL
139: * ..
140: * .. Executable Statements ..
141: *
142: IF( N.LE.0 ) THEN
143: TAU = ZERO
144: RETURN
145: END IF
146: *
147: XNORM = DNRM2( N-1, X, INCX )
148: *
149: IF( XNORM.EQ.ZERO ) THEN
150: *
151: * H = [+/-1, 0; I], sign chosen so ALPHA >= 0
152: *
153: IF( ALPHA.GE.ZERO ) THEN
154: * When TAU.eq.ZERO, the vector is special-cased to be
155: * all zeros in the application routines. We do not need
156: * to clear it.
157: TAU = ZERO
158: ELSE
159: * However, the application routines rely on explicit
160: * zero checks when TAU.ne.ZERO, and we must clear X.
161: TAU = TWO
162: DO J = 1, N-1
163: X( 1 + (J-1)*INCX ) = 0
164: END DO
165: ALPHA = -ALPHA
166: END IF
167: ELSE
168: *
169: * general case
170: *
171: BETA = SIGN( DLAPY2( ALPHA, XNORM ), ALPHA )
172: SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'E' )
173: KNT = 0
174: IF( ABS( BETA ).LT.SMLNUM ) THEN
175: *
176: * XNORM, BETA may be inaccurate; scale X and recompute them
177: *
178: BIGNUM = ONE / SMLNUM
179: 10 CONTINUE
180: KNT = KNT + 1
181: CALL DSCAL( N-1, BIGNUM, X, INCX )
182: BETA = BETA*BIGNUM
183: ALPHA = ALPHA*BIGNUM
184: IF( ABS( BETA ).LT.SMLNUM )
185: $ GO TO 10
186: *
187: * New BETA is at most 1, at least SMLNUM
188: *
189: XNORM = DNRM2( N-1, X, INCX )
190: BETA = SIGN( DLAPY2( ALPHA, XNORM ), ALPHA )
191: END IF
192: SAVEALPHA = ALPHA
193: ALPHA = ALPHA + BETA
194: IF( BETA.LT.ZERO ) THEN
195: BETA = -BETA
196: TAU = -ALPHA / BETA
197: ELSE
198: ALPHA = XNORM * (XNORM/ALPHA)
199: TAU = ALPHA / BETA
200: ALPHA = -ALPHA
201: END IF
202: *
203: IF ( ABS(TAU).LE.SMLNUM ) THEN
204: *
205: * In the case where the computed TAU ends up being a denormalized number,
206: * it loses relative accuracy. This is a BIG problem. Solution: flush TAU
207: * to ZERO. This explains the next IF statement.
208: *
209: * (Bug report provided by Pat Quillen from MathWorks on Jul 29, 2009.)
210: * (Thanks Pat. Thanks MathWorks.)
211: *
212: IF( SAVEALPHA.GE.ZERO ) THEN
213: TAU = ZERO
214: ELSE
215: TAU = TWO
216: DO J = 1, N-1
217: X( 1 + (J-1)*INCX ) = 0
218: END DO
219: BETA = -SAVEALPHA
220: END IF
221: *
222: ELSE
223: *
224: * This is the general case.
225: *
226: CALL DSCAL( N-1, ONE / ALPHA, X, INCX )
227: *
228: END IF
229: *
230: * If BETA is subnormal, it may lose relative accuracy
231: *
232: DO 20 J = 1, KNT
233: BETA = BETA*SMLNUM
234: 20 CONTINUE
235: ALPHA = BETA
236: END IF
237: *
238: RETURN
239: *
240: * End of DLARFGP
241: *
242: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlarfgp.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack