Annotation of rpl/lapack/lapack/dgehrd.f, revision 1.9
1.9 ! bertrand 1: *> \brief \b DGEHRD
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download DGEHRD + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgehrd.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgehrd.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgehrd.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE DGEHRD( N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO )
! 22: *
! 23: * .. Scalar Arguments ..
! 24: * INTEGER IHI, ILO, INFO, LDA, LWORK, N
! 25: * ..
! 26: * .. Array Arguments ..
! 27: * DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * )
! 28: * ..
! 29: *
! 30: *
! 31: *> \par Purpose:
! 32: * =============
! 33: *>
! 34: *> \verbatim
! 35: *>
! 36: *> DGEHRD reduces a real general matrix A to upper Hessenberg form H by
! 37: *> an orthogonal similarity transformation: Q**T * A * Q = H .
! 38: *> \endverbatim
! 39: *
! 40: * Arguments:
! 41: * ==========
! 42: *
! 43: *> \param[in] N
! 44: *> \verbatim
! 45: *> N is INTEGER
! 46: *> The order of the matrix A. N >= 0.
! 47: *> \endverbatim
! 48: *>
! 49: *> \param[in] ILO
! 50: *> \verbatim
! 51: *> ILO is INTEGER
! 52: *> \endverbatim
! 53: *>
! 54: *> \param[in] IHI
! 55: *> \verbatim
! 56: *> IHI is INTEGER
! 57: *>
! 58: *> It is assumed that A is already upper triangular in rows
! 59: *> and columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally
! 60: *> set by a previous call to DGEBAL; otherwise they should be
! 61: *> set to 1 and N respectively. See Further Details.
! 62: *> 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
! 63: *> \endverbatim
! 64: *>
! 65: *> \param[in,out] A
! 66: *> \verbatim
! 67: *> A is DOUBLE PRECISION array, dimension (LDA,N)
! 68: *> On entry, the N-by-N general matrix to be reduced.
! 69: *> On exit, the upper triangle and the first subdiagonal of A
! 70: *> are overwritten with the upper Hessenberg matrix H, and the
! 71: *> elements below the first subdiagonal, with the array TAU,
! 72: *> represent the orthogonal matrix Q as a product of elementary
! 73: *> reflectors. See Further Details.
! 74: *> \endverbatim
! 75: *>
! 76: *> \param[in] LDA
! 77: *> \verbatim
! 78: *> LDA is INTEGER
! 79: *> The leading dimension of the array A. LDA >= max(1,N).
! 80: *> \endverbatim
! 81: *>
! 82: *> \param[out] TAU
! 83: *> \verbatim
! 84: *> TAU is DOUBLE PRECISION array, dimension (N-1)
! 85: *> The scalar factors of the elementary reflectors (see Further
! 86: *> Details). Elements 1:ILO-1 and IHI:N-1 of TAU are set to
! 87: *> zero.
! 88: *> \endverbatim
! 89: *>
! 90: *> \param[out] WORK
! 91: *> \verbatim
! 92: *> WORK is DOUBLE PRECISION array, dimension (LWORK)
! 93: *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
! 94: *> \endverbatim
! 95: *>
! 96: *> \param[in] LWORK
! 97: *> \verbatim
! 98: *> LWORK is INTEGER
! 99: *> The length of the array WORK. LWORK >= max(1,N).
! 100: *> For optimum performance LWORK >= N*NB, where NB is the
! 101: *> optimal blocksize.
! 102: *>
! 103: *> If LWORK = -1, then a workspace query is assumed; the routine
! 104: *> only calculates the optimal size of the WORK array, returns
! 105: *> this value as the first entry of the WORK array, and no error
! 106: *> message related to LWORK is issued by XERBLA.
! 107: *> \endverbatim
! 108: *>
! 109: *> \param[out] INFO
! 110: *> \verbatim
! 111: *> INFO is INTEGER
! 112: *> = 0: successful exit
! 113: *> < 0: if INFO = -i, the i-th argument had an illegal value.
! 114: *> \endverbatim
! 115: *
! 116: * Authors:
! 117: * ========
! 118: *
! 119: *> \author Univ. of Tennessee
! 120: *> \author Univ. of California Berkeley
! 121: *> \author Univ. of Colorado Denver
! 122: *> \author NAG Ltd.
! 123: *
! 124: *> \date November 2011
! 125: *
! 126: *> \ingroup doubleGEcomputational
! 127: *
! 128: *> \par Further Details:
! 129: * =====================
! 130: *>
! 131: *> \verbatim
! 132: *>
! 133: *> The matrix Q is represented as a product of (ihi-ilo) elementary
! 134: *> reflectors
! 135: *>
! 136: *> Q = H(ilo) H(ilo+1) . . . H(ihi-1).
! 137: *>
! 138: *> Each H(i) has the form
! 139: *>
! 140: *> H(i) = I - tau * v * v**T
! 141: *>
! 142: *> where tau is a real scalar, and v is a real vector with
! 143: *> v(1:i) = 0, v(i+1) = 1 and v(ihi+1:n) = 0; v(i+2:ihi) is stored on
! 144: *> exit in A(i+2:ihi,i), and tau in TAU(i).
! 145: *>
! 146: *> The contents of A are illustrated by the following example, with
! 147: *> n = 7, ilo = 2 and ihi = 6:
! 148: *>
! 149: *> on entry, on exit,
! 150: *>
! 151: *> ( a a a a a a a ) ( a a h h h h a )
! 152: *> ( a a a a a a ) ( a h h h h a )
! 153: *> ( a a a a a a ) ( h h h h h h )
! 154: *> ( a a a a a a ) ( v2 h h h h h )
! 155: *> ( a a a a a a ) ( v2 v3 h h h h )
! 156: *> ( a a a a a a ) ( v2 v3 v4 h h h )
! 157: *> ( a ) ( a )
! 158: *>
! 159: *> where a denotes an element of the original matrix A, h denotes a
! 160: *> modified element of the upper Hessenberg matrix H, and vi denotes an
! 161: *> element of the vector defining H(i).
! 162: *>
! 163: *> This file is a slight modification of LAPACK-3.0's DGEHRD
! 164: *> subroutine incorporating improvements proposed by Quintana-Orti and
! 165: *> Van de Geijn (2006). (See DLAHR2.)
! 166: *> \endverbatim
! 167: *>
! 168: * =====================================================================
1.1 bertrand 169: SUBROUTINE DGEHRD( N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO )
170: *
1.9 ! bertrand 171: * -- LAPACK computational routine (version 3.4.0) --
1.1 bertrand 172: * -- LAPACK is a software package provided by Univ. of Tennessee, --
173: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.9 ! bertrand 174: * November 2011
1.1 bertrand 175: *
176: * .. Scalar Arguments ..
177: INTEGER IHI, ILO, INFO, LDA, LWORK, N
178: * ..
179: * .. Array Arguments ..
180: DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * )
181: * ..
182: *
183: * =====================================================================
184: *
185: * .. Parameters ..
186: INTEGER NBMAX, LDT
187: PARAMETER ( NBMAX = 64, LDT = NBMAX+1 )
188: DOUBLE PRECISION ZERO, ONE
189: PARAMETER ( ZERO = 0.0D+0,
190: $ ONE = 1.0D+0 )
191: * ..
192: * .. Local Scalars ..
193: LOGICAL LQUERY
194: INTEGER I, IB, IINFO, IWS, J, LDWORK, LWKOPT, NB,
195: $ NBMIN, NH, NX
196: DOUBLE PRECISION EI
197: * ..
198: * .. Local Arrays ..
199: DOUBLE PRECISION T( LDT, NBMAX )
200: * ..
201: * .. External Subroutines ..
202: EXTERNAL DAXPY, DGEHD2, DGEMM, DLAHR2, DLARFB, DTRMM,
203: $ XERBLA
204: * ..
205: * .. Intrinsic Functions ..
206: INTRINSIC MAX, MIN
207: * ..
208: * .. External Functions ..
209: INTEGER ILAENV
210: EXTERNAL ILAENV
211: * ..
212: * .. Executable Statements ..
213: *
214: * Test the input parameters
215: *
216: INFO = 0
217: NB = MIN( NBMAX, ILAENV( 1, 'DGEHRD', ' ', N, ILO, IHI, -1 ) )
218: LWKOPT = N*NB
219: WORK( 1 ) = LWKOPT
220: LQUERY = ( LWORK.EQ.-1 )
221: IF( N.LT.0 ) THEN
222: INFO = -1
223: ELSE IF( ILO.LT.1 .OR. ILO.GT.MAX( 1, N ) ) THEN
224: INFO = -2
225: ELSE IF( IHI.LT.MIN( ILO, N ) .OR. IHI.GT.N ) THEN
226: INFO = -3
227: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
228: INFO = -5
229: ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
230: INFO = -8
231: END IF
232: IF( INFO.NE.0 ) THEN
233: CALL XERBLA( 'DGEHRD', -INFO )
234: RETURN
235: ELSE IF( LQUERY ) THEN
236: RETURN
237: END IF
238: *
239: * Set elements 1:ILO-1 and IHI:N-1 of TAU to zero
240: *
241: DO 10 I = 1, ILO - 1
242: TAU( I ) = ZERO
243: 10 CONTINUE
244: DO 20 I = MAX( 1, IHI ), N - 1
245: TAU( I ) = ZERO
246: 20 CONTINUE
247: *
248: * Quick return if possible
249: *
250: NH = IHI - ILO + 1
251: IF( NH.LE.1 ) THEN
252: WORK( 1 ) = 1
253: RETURN
254: END IF
255: *
256: * Determine the block size
257: *
258: NB = MIN( NBMAX, ILAENV( 1, 'DGEHRD', ' ', N, ILO, IHI, -1 ) )
259: NBMIN = 2
260: IWS = 1
261: IF( NB.GT.1 .AND. NB.LT.NH ) THEN
262: *
263: * Determine when to cross over from blocked to unblocked code
264: * (last block is always handled by unblocked code)
265: *
266: NX = MAX( NB, ILAENV( 3, 'DGEHRD', ' ', N, ILO, IHI, -1 ) )
267: IF( NX.LT.NH ) THEN
268: *
269: * Determine if workspace is large enough for blocked code
270: *
271: IWS = N*NB
272: IF( LWORK.LT.IWS ) THEN
273: *
274: * Not enough workspace to use optimal NB: determine the
275: * minimum value of NB, and reduce NB or force use of
276: * unblocked code
277: *
278: NBMIN = MAX( 2, ILAENV( 2, 'DGEHRD', ' ', N, ILO, IHI,
279: $ -1 ) )
280: IF( LWORK.GE.N*NBMIN ) THEN
281: NB = LWORK / N
282: ELSE
283: NB = 1
284: END IF
285: END IF
286: END IF
287: END IF
288: LDWORK = N
289: *
290: IF( NB.LT.NBMIN .OR. NB.GE.NH ) THEN
291: *
292: * Use unblocked code below
293: *
294: I = ILO
295: *
296: ELSE
297: *
298: * Use blocked code
299: *
300: DO 40 I = ILO, IHI - 1 - NX, NB
301: IB = MIN( NB, IHI-I )
302: *
303: * Reduce columns i:i+ib-1 to Hessenberg form, returning the
1.8 bertrand 304: * matrices V and T of the block reflector H = I - V*T*V**T
1.1 bertrand 305: * which performs the reduction, and also the matrix Y = A*V*T
306: *
307: CALL DLAHR2( IHI, I, IB, A( 1, I ), LDA, TAU( I ), T, LDT,
308: $ WORK, LDWORK )
309: *
310: * Apply the block reflector H to A(1:ihi,i+ib:ihi) from the
1.8 bertrand 311: * right, computing A := A - Y * V**T. V(i+ib,ib-1) must be set
1.1 bertrand 312: * to 1
313: *
314: EI = A( I+IB, I+IB-1 )
315: A( I+IB, I+IB-1 ) = ONE
316: CALL DGEMM( 'No transpose', 'Transpose',
317: $ IHI, IHI-I-IB+1,
318: $ IB, -ONE, WORK, LDWORK, A( I+IB, I ), LDA, ONE,
319: $ A( 1, I+IB ), LDA )
320: A( I+IB, I+IB-1 ) = EI
321: *
322: * Apply the block reflector H to A(1:i,i+1:i+ib-1) from the
323: * right
324: *
325: CALL DTRMM( 'Right', 'Lower', 'Transpose',
326: $ 'Unit', I, IB-1,
327: $ ONE, A( I+1, I ), LDA, WORK, LDWORK )
328: DO 30 J = 0, IB-2
329: CALL DAXPY( I, -ONE, WORK( LDWORK*J+1 ), 1,
330: $ A( 1, I+J+1 ), 1 )
331: 30 CONTINUE
332: *
333: * Apply the block reflector H to A(i+1:ihi,i+ib:n) from the
334: * left
335: *
336: CALL DLARFB( 'Left', 'Transpose', 'Forward',
337: $ 'Columnwise',
338: $ IHI-I, N-I-IB+1, IB, A( I+1, I ), LDA, T, LDT,
339: $ A( I+1, I+IB ), LDA, WORK, LDWORK )
340: 40 CONTINUE
341: END IF
342: *
343: * Use unblocked code to reduce the rest of the matrix
344: *
345: CALL DGEHD2( N, I, IHI, A, LDA, TAU, WORK, IINFO )
346: WORK( 1 ) = IWS
347: *
348: RETURN
349: *
350: * End of DGEHRD
351: *
352: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dgehrd.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack