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 good performance, LWORK should generally be larger.
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*
118*> \author Univ. of Tennessee
119*> \author Univ. of California Berkeley
120*> \author Univ. of Colorado Denver
121*> \author NAG Ltd.
122*
123*> \ingroup complex16GEcomputational
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**H
138*>
139*>  where tau is a complex scalar, and v is a complex 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 ZGEHRD
161*>  subroutine incorporating improvements proposed by Quintana-Orti and
162*>  Van de Geijn (2006). (See ZLAHR2.)
163*> \endverbatim
164*>
165*  =====================================================================
166      SUBROUTINE ZGEHRD( N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO )
167*
168*  -- LAPACK computational routine --
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      COMPLEX*16        A( LDA, * ), TAU( * ), WORK( * )
177*     ..
178*
179*  =====================================================================
180*
181*     .. Parameters ..
182      INTEGER            NBMAX, LDT, TSIZE
183      PARAMETER          ( NBMAX = 64, LDT = NBMAX+1,
184     $                     TSIZE = LDT*NBMAX )
185      COMPLEX*16        ZERO, ONE
186      PARAMETER          ( ZERO = ( 0.0D+0, 0.0D+0 ),
187     $                     ONE = ( 1.0D+0, 0.0D+0 ) )
188*     ..
189*     .. Local Scalars ..
190      LOGICAL            LQUERY
191      INTEGER            I, IB, IINFO, IWT, J, LDWORK, LWKOPT, NB,
192     $                   NBMIN, NH, NX
193      COMPLEX*16        EI
194*     ..
195*     .. External Subroutines ..
196      EXTERNAL           ZAXPY, ZGEHD2, ZGEMM, ZLAHR2, ZLARFB, ZTRMM,
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
223*
224      IF( INFO.EQ.0 ) THEN
225*
226*        Compute the workspace requirements
227*
228         NB = MIN( NBMAX, ILAENV( 1, 'ZGEHRD', ' ', N, ILO, IHI, -1 ) )
229         LWKOPT = N*NB + TSIZE
230         WORK( 1 ) = LWKOPT
231      ENDIF
232*
233      IF( INFO.NE.0 ) THEN
234         CALL XERBLA( 'ZGEHRD', -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, 'ZGEHRD', ' ', 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, 'ZGEHRD', ' ', N, ILO, IHI, -1 ) )
267         IF( NX.LT.NH ) THEN
268*
269*           Determine if workspace is large enough for blocked code
270*
271            IF( LWORK.LT.N*NB+TSIZE ) THEN
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, 'ZGEHRD', ' ', N, ILO, IHI,
278     $                 -1 ) )
279               IF( LWORK.GE.(N*NBMIN + TSIZE) ) THEN
280                  NB = (LWORK-TSIZE) / N
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*
299         IWT = 1 + N*NB
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 ),
308     $                   WORK( IWT ), LDT, 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,
339     $                   WORK( IWT ), LDT, A( I+1, I+IB ), LDA,
340     $                   WORK, LDWORK )
341   40    CONTINUE
342      END IF
343*
344*     Use unblocked code to reduce the rest of the matrix
345*
346      CALL ZGEHD2( N, I, IHI, A, LDA, TAU, WORK, IINFO )
347      WORK( 1 ) = LWKOPT
348*
349      RETURN
350*
351*     End of ZGEHRD
352*
353      END
354