Annotation of rpl/lapack/lapack/dgehrd.f, revision 1.10
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