Annotation of rpl/lapack/lapack/zgehrd.f, revision 1.17
1.9 bertrand 1: *> \brief \b ZGEHRD
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 ZGEHRD + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgehrd.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgehrd.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgehrd.f">
1.9 bertrand 15: *> [TXT]</a>
1.16 bertrand 16: *> \endhtmlonly
1.9 bertrand 17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZGEHRD( 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: * COMPLEX*16 A( LDA, * ), TAU( * ), WORK( * )
28: * ..
1.16 bertrand 29: *
1.9 bertrand 30: *
31: *> \par Purpose:
32: * =============
33: *>
34: *> \verbatim
35: *>
36: *> ZGEHRD reduces a complex general matrix A to upper Hessenberg form H by
37: *> an unitary similarity transformation: Q**H * 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 ZGEBAL; 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 COMPLEX*16 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 unitary 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 COMPLEX*16 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 COMPLEX*16 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: *
1.16 bertrand 123: *> \date December 2016
1.9 bertrand 124: *
125: *> \ingroup complex16GEcomputational
126: *
127: *> \par Further Details:
128: * =====================
129: *>
130: *> \verbatim
131: *>
132: *> The matrix Q is represented as a product of (ihi-ilo) elementary
133: *> reflectors
134: *>
135: *> Q = H(ilo) H(ilo+1) . . . H(ihi-1).
136: *>
137: *> Each H(i) has the form
138: *>
139: *> H(i) = I - tau * v * v**H
140: *>
141: *> where tau is a complex scalar, and v is a complex vector with
142: *> v(1:i) = 0, v(i+1) = 1 and v(ihi+1:n) = 0; v(i+2:ihi) is stored on
143: *> exit in A(i+2:ihi,i), and tau in TAU(i).
144: *>
145: *> The contents of A are illustrated by the following example, with
146: *> n = 7, ilo = 2 and ihi = 6:
147: *>
148: *> on entry, on exit,
149: *>
150: *> ( a a a a a a a ) ( a a h h h h a )
151: *> ( a a a a a a ) ( a h h h h a )
152: *> ( a a a a a a ) ( h h h h h h )
153: *> ( a a a a a a ) ( v2 h h h h h )
154: *> ( a a a a a a ) ( v2 v3 h h h h )
155: *> ( a a a a a a ) ( v2 v3 v4 h h h )
156: *> ( a ) ( a )
157: *>
158: *> where a denotes an element of the original matrix A, h denotes a
159: *> modified element of the upper Hessenberg matrix H, and vi denotes an
160: *> element of the vector defining H(i).
161: *>
162: *> This file is a slight modification of LAPACK-3.0's DGEHRD
163: *> subroutine incorporating improvements proposed by Quintana-Orti and
164: *> Van de Geijn (2006). (See DLAHR2.)
165: *> \endverbatim
166: *>
167: * =====================================================================
1.1 bertrand 168: SUBROUTINE ZGEHRD( N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO )
169: *
1.16 bertrand 170: * -- LAPACK computational routine (version 3.7.0) --
1.1 bertrand 171: * -- LAPACK is a software package provided by Univ. of Tennessee, --
172: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.16 bertrand 173: * December 2016
1.1 bertrand 174: *
175: * .. Scalar Arguments ..
176: INTEGER IHI, ILO, INFO, LDA, LWORK, N
177: * ..
178: * .. Array Arguments ..
179: COMPLEX*16 A( LDA, * ), TAU( * ), WORK( * )
180: * ..
181: *
182: * =====================================================================
183: *
184: * .. Parameters ..
1.14 bertrand 185: INTEGER NBMAX, LDT, TSIZE
186: PARAMETER ( NBMAX = 64, LDT = NBMAX+1,
187: $ TSIZE = LDT*NBMAX )
1.1 bertrand 188: COMPLEX*16 ZERO, ONE
1.16 bertrand 189: PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ),
1.1 bertrand 190: $ ONE = ( 1.0D+0, 0.0D+0 ) )
191: * ..
192: * .. Local Scalars ..
193: LOGICAL LQUERY
1.14 bertrand 194: INTEGER I, IB, IINFO, IWT, J, LDWORK, LWKOPT, NB,
1.1 bertrand 195: $ NBMIN, NH, NX
196: COMPLEX*16 EI
197: * ..
198: * .. External Subroutines ..
199: EXTERNAL ZAXPY, ZGEHD2, ZGEMM, ZLAHR2, ZLARFB, ZTRMM,
200: $ XERBLA
201: * ..
202: * .. Intrinsic Functions ..
203: INTRINSIC MAX, MIN
204: * ..
205: * .. External Functions ..
206: INTEGER ILAENV
207: EXTERNAL ILAENV
208: * ..
209: * .. Executable Statements ..
210: *
211: * Test the input parameters
212: *
213: INFO = 0
214: LQUERY = ( LWORK.EQ.-1 )
215: IF( N.LT.0 ) THEN
216: INFO = -1
217: ELSE IF( ILO.LT.1 .OR. ILO.GT.MAX( 1, N ) ) THEN
218: INFO = -2
219: ELSE IF( IHI.LT.MIN( ILO, N ) .OR. IHI.GT.N ) THEN
220: INFO = -3
221: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
222: INFO = -5
223: ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
224: INFO = -8
225: END IF
1.14 bertrand 226: *
227: IF( INFO.EQ.0 ) THEN
228: *
229: * Compute the workspace requirements
230: *
231: NB = MIN( NBMAX, ILAENV( 1, 'ZGEHRD', ' ', N, ILO, IHI, -1 ) )
232: LWKOPT = N*NB + TSIZE
233: WORK( 1 ) = LWKOPT
234: ENDIF
235: *
1.1 bertrand 236: IF( INFO.NE.0 ) THEN
237: CALL XERBLA( 'ZGEHRD', -INFO )
238: RETURN
239: ELSE IF( LQUERY ) THEN
240: RETURN
241: END IF
242: *
243: * Set elements 1:ILO-1 and IHI:N-1 of TAU to zero
244: *
245: DO 10 I = 1, ILO - 1
246: TAU( I ) = ZERO
247: 10 CONTINUE
248: DO 20 I = MAX( 1, IHI ), N - 1
249: TAU( I ) = ZERO
250: 20 CONTINUE
251: *
252: * Quick return if possible
253: *
254: NH = IHI - ILO + 1
255: IF( NH.LE.1 ) THEN
256: WORK( 1 ) = 1
257: RETURN
258: END IF
259: *
260: * Determine the block size
261: *
262: NB = MIN( NBMAX, ILAENV( 1, 'ZGEHRD', ' ', N, ILO, IHI, -1 ) )
263: NBMIN = 2
264: IF( NB.GT.1 .AND. NB.LT.NH ) THEN
265: *
266: * Determine when to cross over from blocked to unblocked code
267: * (last block is always handled by unblocked code)
268: *
269: NX = MAX( NB, ILAENV( 3, 'ZGEHRD', ' ', N, ILO, IHI, -1 ) )
270: IF( NX.LT.NH ) THEN
271: *
272: * Determine if workspace is large enough for blocked code
273: *
1.14 bertrand 274: IF( LWORK.LT.N*NB+TSIZE ) THEN
1.1 bertrand 275: *
276: * Not enough workspace to use optimal NB: determine the
277: * minimum value of NB, and reduce NB or force use of
278: * unblocked code
279: *
280: NBMIN = MAX( 2, ILAENV( 2, 'ZGEHRD', ' ', N, ILO, IHI,
281: $ -1 ) )
1.14 bertrand 282: IF( LWORK.GE.(N*NBMIN + TSIZE) ) THEN
283: NB = (LWORK-TSIZE) / N
1.1 bertrand 284: ELSE
285: NB = 1
286: END IF
287: END IF
288: END IF
289: END IF
290: LDWORK = N
291: *
292: IF( NB.LT.NBMIN .OR. NB.GE.NH ) THEN
293: *
294: * Use unblocked code below
295: *
296: I = ILO
297: *
298: ELSE
299: *
300: * Use blocked code
301: *
1.14 bertrand 302: IWT = 1 + N*NB
1.1 bertrand 303: DO 40 I = ILO, IHI - 1 - NX, NB
304: IB = MIN( NB, IHI-I )
305: *
306: * Reduce columns i:i+ib-1 to Hessenberg form, returning the
1.8 bertrand 307: * matrices V and T of the block reflector H = I - V*T*V**H
1.1 bertrand 308: * which performs the reduction, and also the matrix Y = A*V*T
309: *
1.14 bertrand 310: CALL ZLAHR2( IHI, I, IB, A( 1, I ), LDA, TAU( I ),
311: $ WORK( IWT ), LDT, WORK, LDWORK )
1.1 bertrand 312: *
313: * Apply the block reflector H to A(1:ihi,i+ib:ihi) from the
1.8 bertrand 314: * right, computing A := A - Y * V**H. V(i+ib,ib-1) must be set
1.1 bertrand 315: * to 1
316: *
317: EI = A( I+IB, I+IB-1 )
318: A( I+IB, I+IB-1 ) = ONE
1.16 bertrand 319: CALL ZGEMM( 'No transpose', 'Conjugate transpose',
1.1 bertrand 320: $ IHI, IHI-I-IB+1,
321: $ IB, -ONE, WORK, LDWORK, A( I+IB, I ), LDA, ONE,
322: $ A( 1, I+IB ), LDA )
323: A( I+IB, I+IB-1 ) = EI
324: *
325: * Apply the block reflector H to A(1:i,i+1:i+ib-1) from the
326: * right
327: *
328: CALL ZTRMM( 'Right', 'Lower', 'Conjugate transpose',
329: $ 'Unit', I, IB-1,
330: $ ONE, A( I+1, I ), LDA, WORK, LDWORK )
331: DO 30 J = 0, IB-2
332: CALL ZAXPY( I, -ONE, WORK( LDWORK*J+1 ), 1,
333: $ A( 1, I+J+1 ), 1 )
334: 30 CONTINUE
335: *
336: * Apply the block reflector H to A(i+1:ihi,i+ib:n) from the
337: * left
338: *
339: CALL ZLARFB( 'Left', 'Conjugate transpose', 'Forward',
340: $ 'Columnwise',
1.14 bertrand 341: $ IHI-I, N-I-IB+1, IB, A( I+1, I ), LDA,
342: $ WORK( IWT ), LDT, A( I+1, I+IB ), LDA,
343: $ WORK, LDWORK )
1.1 bertrand 344: 40 CONTINUE
345: END IF
346: *
347: * Use unblocked code to reduce the rest of the matrix
348: *
349: CALL ZGEHD2( N, I, IHI, A, LDA, TAU, WORK, IINFO )
1.14 bertrand 350: WORK( 1 ) = LWKOPT
1.1 bertrand 351: *
352: RETURN
353: *
354: * End of ZGEHRD
355: *
356: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zgehrd.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack