Annotation of rpl/lapack/lapack/dgehrd.f, revision 1.19
1.9 bertrand 1: *> \brief \b DGEHRD
2: *
3: * =========== DOCUMENTATION ===========
4: *
1.16 bertrand 5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
1.9 bertrand 7: *
8: *> \htmlonly
1.16 bertrand 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">
1.9 bertrand 15: *> [TXT]</a>
1.16 bertrand 16: *> \endhtmlonly
1.9 bertrand 17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DGEHRD( N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO )
1.16 bertrand 22: *
1.9 bertrand 23: * .. Scalar Arguments ..
24: * INTEGER IHI, ILO, INFO, LDA, LWORK, N
25: * ..
26: * .. Array Arguments ..
27: * DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * )
28: * ..
1.16 bertrand 29: *
1.9 bertrand 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).
1.14 bertrand 100: *> For good performance, LWORK should generally be larger.
1.9 bertrand 101: *>
102: *> If LWORK = -1, then a workspace query is assumed; the routine
103: *> only calculates the optimal size of the WORK array, returns
104: *> this value as the first entry of the WORK array, and no error
105: *> message related to LWORK is issued by XERBLA.
106: *> \endverbatim
107: *>
108: *> \param[out] INFO
109: *> \verbatim
110: *> INFO is INTEGER
111: *> = 0: successful exit
112: *> < 0: if INFO = -i, the i-th argument had an illegal value.
113: *> \endverbatim
114: *
115: * Authors:
116: * ========
117: *
1.16 bertrand 118: *> \author Univ. of Tennessee
119: *> \author Univ. of California Berkeley
120: *> \author Univ. of Colorado Denver
121: *> \author NAG Ltd.
1.9 bertrand 122: *
123: *> \ingroup doubleGEcomputational
124: *
125: *> \par Further Details:
126: * =====================
127: *>
128: *> \verbatim
129: *>
130: *> The matrix Q is represented as a product of (ihi-ilo) elementary
131: *> reflectors
132: *>
133: *> Q = H(ilo) H(ilo+1) . . . H(ihi-1).
134: *>
135: *> Each H(i) has the form
136: *>
137: *> H(i) = I - tau * v * v**T
138: *>
139: *> where tau is a real scalar, and v is a real vector with
140: *> v(1:i) = 0, v(i+1) = 1 and v(ihi+1:n) = 0; v(i+2:ihi) is stored on
141: *> exit in A(i+2:ihi,i), and tau in TAU(i).
142: *>
143: *> The contents of A are illustrated by the following example, with
144: *> n = 7, ilo = 2 and ihi = 6:
145: *>
146: *> on entry, on exit,
147: *>
148: *> ( a a a a a a a ) ( a a h h h h a )
149: *> ( a a a a a a ) ( a h h h h a )
150: *> ( a a a a a a ) ( h h h h h h )
151: *> ( a a a a a a ) ( v2 h h h h h )
152: *> ( a a a a a a ) ( v2 v3 h h h h )
153: *> ( a a a a a a ) ( v2 v3 v4 h h h )
154: *> ( a ) ( a )
155: *>
156: *> where a denotes an element of the original matrix A, h denotes a
157: *> modified element of the upper Hessenberg matrix H, and vi denotes an
158: *> element of the vector defining H(i).
159: *>
160: *> This file is a slight modification of LAPACK-3.0's DGEHRD
161: *> subroutine incorporating improvements proposed by Quintana-Orti and
162: *> Van de Geijn (2006). (See DLAHR2.)
163: *> \endverbatim
164: *>
165: * =====================================================================
1.1 bertrand 166: SUBROUTINE DGEHRD( N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO )
167: *
1.19 ! bertrand 168: * -- LAPACK computational routine --
1.1 bertrand 169: * -- LAPACK is a software package provided by Univ. of Tennessee, --
170: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
171: *
172: * .. Scalar Arguments ..
173: INTEGER IHI, ILO, INFO, LDA, LWORK, N
174: * ..
175: * .. Array Arguments ..
176: DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * )
177: * ..
178: *
179: * =====================================================================
180: *
181: * .. Parameters ..
1.14 bertrand 182: INTEGER NBMAX, LDT, TSIZE
183: PARAMETER ( NBMAX = 64, LDT = NBMAX+1,
184: $ TSIZE = LDT*NBMAX )
1.1 bertrand 185: DOUBLE PRECISION ZERO, ONE
1.16 bertrand 186: PARAMETER ( ZERO = 0.0D+0,
1.1 bertrand 187: $ ONE = 1.0D+0 )
188: * ..
189: * .. Local Scalars ..
190: LOGICAL LQUERY
1.14 bertrand 191: INTEGER I, IB, IINFO, IWT, J, LDWORK, LWKOPT, NB,
1.1 bertrand 192: $ NBMIN, NH, NX
193: DOUBLE PRECISION EI
194: * ..
195: * .. External Subroutines ..
196: EXTERNAL DAXPY, DGEHD2, DGEMM, DLAHR2, DLARFB, DTRMM,
197: $ XERBLA
198: * ..
199: * .. Intrinsic Functions ..
200: INTRINSIC MAX, MIN
201: * ..
202: * .. External Functions ..
203: INTEGER ILAENV
204: EXTERNAL ILAENV
205: * ..
206: * .. Executable Statements ..
207: *
208: * Test the input parameters
209: *
210: INFO = 0
211: LQUERY = ( LWORK.EQ.-1 )
212: IF( N.LT.0 ) THEN
213: INFO = -1
214: ELSE IF( ILO.LT.1 .OR. ILO.GT.MAX( 1, N ) ) THEN
215: INFO = -2
216: ELSE IF( IHI.LT.MIN( ILO, N ) .OR. IHI.GT.N ) THEN
217: INFO = -3
218: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
219: INFO = -5
220: ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
221: INFO = -8
222: END IF
1.14 bertrand 223: *
224: IF( INFO.EQ.0 ) THEN
225: *
226: * Compute the workspace requirements
227: *
228: NB = MIN( NBMAX, ILAENV( 1, 'DGEHRD', ' ', N, ILO, IHI, -1 ) )
229: LWKOPT = N*NB + TSIZE
230: WORK( 1 ) = LWKOPT
231: END IF
1.16 bertrand 232: *
1.1 bertrand 233: IF( INFO.NE.0 ) THEN
234: CALL XERBLA( 'DGEHRD', -INFO )
235: RETURN
236: ELSE IF( LQUERY ) THEN
237: RETURN
238: END IF
239: *
240: * Set elements 1:ILO-1 and IHI:N-1 of TAU to zero
241: *
242: DO 10 I = 1, ILO - 1
243: TAU( I ) = ZERO
244: 10 CONTINUE
245: DO 20 I = MAX( 1, IHI ), N - 1
246: TAU( I ) = ZERO
247: 20 CONTINUE
248: *
249: * Quick return if possible
250: *
251: NH = IHI - ILO + 1
252: IF( NH.LE.1 ) THEN
253: WORK( 1 ) = 1
254: RETURN
255: END IF
256: *
257: * Determine the block size
258: *
259: NB = MIN( NBMAX, ILAENV( 1, 'DGEHRD', ' ', N, ILO, IHI, -1 ) )
260: NBMIN = 2
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: *
1.14 bertrand 271: IF( LWORK.LT.N*NB+TSIZE ) THEN
1.1 bertrand 272: *
273: * Not enough workspace to use optimal NB: determine the
274: * minimum value of NB, and reduce NB or force use of
275: * unblocked code
276: *
277: NBMIN = MAX( 2, ILAENV( 2, 'DGEHRD', ' ', N, ILO, IHI,
278: $ -1 ) )
1.14 bertrand 279: IF( LWORK.GE.(N*NBMIN + TSIZE) ) THEN
280: NB = (LWORK-TSIZE) / N
1.1 bertrand 281: ELSE
282: NB = 1
283: END IF
284: END IF
285: END IF
286: END IF
287: LDWORK = N
288: *
289: IF( NB.LT.NBMIN .OR. NB.GE.NH ) THEN
290: *
291: * Use unblocked code below
292: *
293: I = ILO
294: *
295: ELSE
296: *
297: * Use blocked code
298: *
1.14 bertrand 299: IWT = 1 + N*NB
1.1 bertrand 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: *
1.14 bertrand 307: CALL DLAHR2( IHI, I, IB, A( 1, I ), LDA, TAU( I ),
308: $ WORK( IWT ), LDT, WORK, LDWORK )
1.1 bertrand 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
1.16 bertrand 316: CALL DGEMM( 'No transpose', 'Transpose',
1.1 bertrand 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',
1.14 bertrand 338: $ IHI-I, N-I-IB+1, IB, A( I+1, I ), LDA,
339: $ WORK( IWT ), LDT, A( I+1, I+IB ), LDA,
340: $ WORK, LDWORK )
1.1 bertrand 341: 40 CONTINUE
342: END IF
343: *
344: * Use unblocked code to reduce the rest of the matrix
345: *
346: CALL DGEHD2( N, I, IHI, A, LDA, TAU, WORK, IINFO )
1.14 bertrand 347: WORK( 1 ) = LWKOPT
1.1 bertrand 348: *
349: RETURN
350: *
351: * End of DGEHRD
352: *
353: 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