1*> \brief \b ZGEHRD
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
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">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE ZGEHRD( 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*       COMPLEX*16        A( LDA, * ), TAU( * ), WORK( * )
28*       ..
29*
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).
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 complex16GEcomputational
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**H
141*>
142*>  where tau is a complex scalar, and v is a complex 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*  =====================================================================
169      SUBROUTINE ZGEHRD( N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO )
170*
171*  -- LAPACK computational routine (version 3.4.0) --
172*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
173*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
174*     November 2011
175*
176*     .. Scalar Arguments ..
177      INTEGER            IHI, ILO, INFO, LDA, LWORK, N
178*     ..
179*     .. Array Arguments ..
180      COMPLEX*16        A( LDA, * ), TAU( * ), WORK( * )
181*     ..
182*
183*  =====================================================================
184*
185*     .. Parameters ..
186      INTEGER            NBMAX, LDT
187      PARAMETER          ( NBMAX = 64, LDT = NBMAX+1 )
188      COMPLEX*16        ZERO, ONE
189      PARAMETER          ( ZERO = ( 0.0D+0, 0.0D+0 ),
190     $                     ONE = ( 1.0D+0, 0.0D+0 ) )
191*     ..
192*     .. Local Scalars ..
193      LOGICAL            LQUERY
194      INTEGER            I, IB, IINFO, IWS, J, LDWORK, LWKOPT, NB,
195     $                   NBMIN, NH, NX
196      COMPLEX*16        EI
197*     ..
198*     .. Local Arrays ..
199      COMPLEX*16        T( LDT, NBMAX )
200*     ..
201*     .. External Subroutines ..
202      EXTERNAL           ZAXPY, ZGEHD2, ZGEMM, ZLAHR2, ZLARFB, ZTRMM,
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, 'ZGEHRD', ' ', 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( 'ZGEHRD', -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, 'ZGEHRD', ' ', 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, 'ZGEHRD', ' ', 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, 'ZGEHRD', ' ', 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
304*           matrices V and T of the block reflector H = I - V*T*V**H
305*           which performs the reduction, and also the matrix Y = A*V*T
306*
307            CALL ZLAHR2( 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
311*           right, computing  A := A - Y * V**H. V(i+ib,ib-1) must be set
312*           to 1
313*
314            EI = A( I+IB, I+IB-1 )
315            A( I+IB, I+IB-1 ) = ONE
316            CALL ZGEMM( 'No transpose', 'Conjugate 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 ZTRMM( 'Right', 'Lower', 'Conjugate transpose',
326     $                  'Unit', I, IB-1,
327     $                  ONE, A( I+1, I ), LDA, WORK, LDWORK )
328            DO 30 J = 0, IB-2
329               CALL ZAXPY( 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 ZLARFB( 'Left', 'Conjugate 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 ZGEHD2( N, I, IHI, A, LDA, TAU, WORK, IINFO )
346      WORK( 1 ) = IWS
347*
348      RETURN
349*
350*     End of ZGEHRD
351*
352      END
353c $Id$
354