1*> \brief \b CGGHRD
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download CGGHRD + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgghrd.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgghrd.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgghrd.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE CGGHRD( COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q,
22*                          LDQ, Z, LDZ, INFO )
23*
24*       .. Scalar Arguments ..
25*       CHARACTER          COMPQ, COMPZ
26*       INTEGER            IHI, ILO, INFO, LDA, LDB, LDQ, LDZ, N
27*       ..
28*       .. Array Arguments ..
29*       COMPLEX            A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
30*      $                   Z( LDZ, * )
31*       ..
32*
33*
34*> \par Purpose:
35*  =============
36*>
37*> \verbatim
38*>
39*> CGGHRD reduces a pair of complex matrices (A,B) to generalized upper
40*> Hessenberg form using unitary transformations, where A is a
41*> general matrix and B is upper triangular.  The form of the generalized
42*> eigenvalue problem is
43*>    A*x = lambda*B*x,
44*> and B is typically made upper triangular by computing its QR
45*> factorization and moving the unitary matrix Q to the left side
46*> of the equation.
47*>
48*> This subroutine simultaneously reduces A to a Hessenberg matrix H:
49*>    Q**H*A*Z = H
50*> and transforms B to another upper triangular matrix T:
51*>    Q**H*B*Z = T
52*> in order to reduce the problem to its standard form
53*>    H*y = lambda*T*y
54*> where y = Z**H*x.
55*>
56*> The unitary matrices Q and Z are determined as products of Givens
57*> rotations.  They may either be formed explicitly, or they may be
58*> postmultiplied into input matrices Q1 and Z1, so that
59*>      Q1 * A * Z1**H = (Q1*Q) * H * (Z1*Z)**H
60*>      Q1 * B * Z1**H = (Q1*Q) * T * (Z1*Z)**H
61*> If Q1 is the unitary matrix from the QR factorization of B in the
62*> original equation A*x = lambda*B*x, then CGGHRD reduces the original
63*> problem to generalized Hessenberg form.
64*> \endverbatim
65*
66*  Arguments:
67*  ==========
68*
69*> \param[in] COMPQ
70*> \verbatim
71*>          COMPQ is CHARACTER*1
72*>          = 'N': do not compute Q;
73*>          = 'I': Q is initialized to the unit matrix, and the
74*>                 unitary matrix Q is returned;
75*>          = 'V': Q must contain a unitary matrix Q1 on entry,
76*>                 and the product Q1*Q is returned.
77*> \endverbatim
78*>
79*> \param[in] COMPZ
80*> \verbatim
81*>          COMPZ is CHARACTER*1
82*>          = 'N': do not compute Z;
83*>          = 'I': Z is initialized to the unit matrix, and the
84*>                 unitary matrix Z is returned;
85*>          = 'V': Z must contain a unitary matrix Z1 on entry,
86*>                 and the product Z1*Z is returned.
87*> \endverbatim
88*>
89*> \param[in] N
90*> \verbatim
91*>          N is INTEGER
92*>          The order of the matrices A and B.  N >= 0.
93*> \endverbatim
94*>
95*> \param[in] ILO
96*> \verbatim
97*>          ILO is INTEGER
98*> \endverbatim
99*>
100*> \param[in] IHI
101*> \verbatim
102*>          IHI is INTEGER
103*>
104*>          ILO and IHI mark the rows and columns of A which are to be
105*>          reduced.  It is assumed that A is already upper triangular
106*>          in rows and columns 1:ILO-1 and IHI+1:N.  ILO and IHI are
107*>          normally set by a previous call to CGGBAL; otherwise they
108*>          should be set to 1 and N respectively.
109*>          1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
110*> \endverbatim
111*>
112*> \param[in,out] A
113*> \verbatim
114*>          A is COMPLEX array, dimension (LDA, N)
115*>          On entry, the N-by-N general matrix to be reduced.
116*>          On exit, the upper triangle and the first subdiagonal of A
117*>          are overwritten with the upper Hessenberg matrix H, and the
118*>          rest is set to zero.
119*> \endverbatim
120*>
121*> \param[in] LDA
122*> \verbatim
123*>          LDA is INTEGER
124*>          The leading dimension of the array A.  LDA >= max(1,N).
125*> \endverbatim
126*>
127*> \param[in,out] B
128*> \verbatim
129*>          B is COMPLEX array, dimension (LDB, N)
130*>          On entry, the N-by-N upper triangular matrix B.
131*>          On exit, the upper triangular matrix T = Q**H B Z.  The
132*>          elements below the diagonal are set to zero.
133*> \endverbatim
134*>
135*> \param[in] LDB
136*> \verbatim
137*>          LDB is INTEGER
138*>          The leading dimension of the array B.  LDB >= max(1,N).
139*> \endverbatim
140*>
141*> \param[in,out] Q
142*> \verbatim
143*>          Q is COMPLEX array, dimension (LDQ, N)
144*>          On entry, if COMPQ = 'V', the unitary matrix Q1, typically
145*>          from the QR factorization of B.
146*>          On exit, if COMPQ='I', the unitary matrix Q, and if
147*>          COMPQ = 'V', the product Q1*Q.
148*>          Not referenced if COMPQ='N'.
149*> \endverbatim
150*>
151*> \param[in] LDQ
152*> \verbatim
153*>          LDQ is INTEGER
154*>          The leading dimension of the array Q.
155*>          LDQ >= N if COMPQ='V' or 'I'; LDQ >= 1 otherwise.
156*> \endverbatim
157*>
158*> \param[in,out] Z
159*> \verbatim
160*>          Z is COMPLEX array, dimension (LDZ, N)
161*>          On entry, if COMPZ = 'V', the unitary matrix Z1.
162*>          On exit, if COMPZ='I', the unitary matrix Z, and if
163*>          COMPZ = 'V', the product Z1*Z.
164*>          Not referenced if COMPZ='N'.
165*> \endverbatim
166*>
167*> \param[in] LDZ
168*> \verbatim
169*>          LDZ is INTEGER
170*>          The leading dimension of the array Z.
171*>          LDZ >= N if COMPZ='V' or 'I'; LDZ >= 1 otherwise.
172*> \endverbatim
173*>
174*> \param[out] INFO
175*> \verbatim
176*>          INFO is INTEGER
177*>          = 0:  successful exit.
178*>          < 0:  if INFO = -i, the i-th argument had an illegal value.
179*> \endverbatim
180*
181*  Authors:
182*  ========
183*
184*> \author Univ. of Tennessee
185*> \author Univ. of California Berkeley
186*> \author Univ. of Colorado Denver
187*> \author NAG Ltd.
188*
189*> \ingroup complexOTHERcomputational
190*
191*> \par Further Details:
192*  =====================
193*>
194*> \verbatim
195*>
196*>  This routine reduces A to Hessenberg and B to triangular form by
197*>  an unblocked reduction, as described in _Matrix_Computations_,
198*>  by Golub and van Loan (Johns Hopkins Press).
199*> \endverbatim
200*>
201*  =====================================================================
202      SUBROUTINE CGGHRD( COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q,
203     $                   LDQ, Z, LDZ, INFO )
204*
205*  -- LAPACK computational routine --
206*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
207*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
208*
209*     .. Scalar Arguments ..
210      CHARACTER          COMPQ, COMPZ
211      INTEGER            IHI, ILO, INFO, LDA, LDB, LDQ, LDZ, N
212*     ..
213*     .. Array Arguments ..
214      COMPLEX            A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
215     $                   Z( LDZ, * )
216*     ..
217*
218*  =====================================================================
219*
220*     .. Parameters ..
221      COMPLEX            CONE, CZERO
222      PARAMETER          ( CONE = ( 1.0E+0, 0.0E+0 ),
223     $                   CZERO = ( 0.0E+0, 0.0E+0 ) )
224*     ..
225*     .. Local Scalars ..
226      LOGICAL            ILQ, ILZ
227      INTEGER            ICOMPQ, ICOMPZ, JCOL, JROW
228      REAL               C
229      COMPLEX            CTEMP, S
230*     ..
231*     .. External Functions ..
232      LOGICAL            LSAME
233      EXTERNAL           LSAME
234*     ..
235*     .. External Subroutines ..
236      EXTERNAL           CLARTG, CLASET, CROT, XERBLA
237*     ..
238*     .. Intrinsic Functions ..
239      INTRINSIC          CONJG, MAX
240*     ..
241*     .. Executable Statements ..
242*
243*     Decode COMPQ
244*
245      IF( LSAME( COMPQ, 'N' ) ) THEN
246         ILQ = .FALSE.
247         ICOMPQ = 1
248      ELSE IF( LSAME( COMPQ, 'V' ) ) THEN
249         ILQ = .TRUE.
250         ICOMPQ = 2
251      ELSE IF( LSAME( COMPQ, 'I' ) ) THEN
252         ILQ = .TRUE.
253         ICOMPQ = 3
254      ELSE
255         ICOMPQ = 0
256      END IF
257*
258*     Decode COMPZ
259*
260      IF( LSAME( COMPZ, 'N' ) ) THEN
261         ILZ = .FALSE.
262         ICOMPZ = 1
263      ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
264         ILZ = .TRUE.
265         ICOMPZ = 2
266      ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
267         ILZ = .TRUE.
268         ICOMPZ = 3
269      ELSE
270         ICOMPZ = 0
271      END IF
272*
273*     Test the input parameters.
274*
275      INFO = 0
276      IF( ICOMPQ.LE.0 ) THEN
277         INFO = -1
278      ELSE IF( ICOMPZ.LE.0 ) THEN
279         INFO = -2
280      ELSE IF( N.LT.0 ) THEN
281         INFO = -3
282      ELSE IF( ILO.LT.1 ) THEN
283         INFO = -4
284      ELSE IF( IHI.GT.N .OR. IHI.LT.ILO-1 ) THEN
285         INFO = -5
286      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
287         INFO = -7
288      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
289         INFO = -9
290      ELSE IF( ( ILQ .AND. LDQ.LT.N ) .OR. LDQ.LT.1 ) THEN
291         INFO = -11
292      ELSE IF( ( ILZ .AND. LDZ.LT.N ) .OR. LDZ.LT.1 ) THEN
293         INFO = -13
294      END IF
295      IF( INFO.NE.0 ) THEN
296         CALL XERBLA( 'CGGHRD', -INFO )
297         RETURN
298      END IF
299*
300*     Initialize Q and Z if desired.
301*
302      IF( ICOMPQ.EQ.3 )
303     $   CALL CLASET( 'Full', N, N, CZERO, CONE, Q, LDQ )
304      IF( ICOMPZ.EQ.3 )
305     $   CALL CLASET( 'Full', N, N, CZERO, CONE, Z, LDZ )
306*
307*     Quick return if possible
308*
309      IF( N.LE.1 )
310     $   RETURN
311*
312*     Zero out lower triangle of B
313*
314      DO 20 JCOL = 1, N - 1
315         DO 10 JROW = JCOL + 1, N
316            B( JROW, JCOL ) = CZERO
317   10    CONTINUE
318   20 CONTINUE
319*
320*     Reduce A and B
321*
322      DO 40 JCOL = ILO, IHI - 2
323*
324         DO 30 JROW = IHI, JCOL + 2, -1
325*
326*           Step 1: rotate rows JROW-1, JROW to kill A(JROW,JCOL)
327*
328            CTEMP = A( JROW-1, JCOL )
329            CALL CLARTG( CTEMP, A( JROW, JCOL ), C, S,
330     $                   A( JROW-1, JCOL ) )
331            A( JROW, JCOL ) = CZERO
332            CALL CROT( N-JCOL, A( JROW-1, JCOL+1 ), LDA,
333     $                 A( JROW, JCOL+1 ), LDA, C, S )
334            CALL CROT( N+2-JROW, B( JROW-1, JROW-1 ), LDB,
335     $                 B( JROW, JROW-1 ), LDB, C, S )
336            IF( ILQ )
337     $         CALL CROT( N, Q( 1, JROW-1 ), 1, Q( 1, JROW ), 1, C,
338     $                    CONJG( S ) )
339*
340*           Step 2: rotate columns JROW, JROW-1 to kill B(JROW,JROW-1)
341*
342            CTEMP = B( JROW, JROW )
343            CALL CLARTG( CTEMP, B( JROW, JROW-1 ), C, S,
344     $                   B( JROW, JROW ) )
345            B( JROW, JROW-1 ) = CZERO
346            CALL CROT( IHI, A( 1, JROW ), 1, A( 1, JROW-1 ), 1, C, S )
347            CALL CROT( JROW-1, B( 1, JROW ), 1, B( 1, JROW-1 ), 1, C,
348     $                 S )
349            IF( ILZ )
350     $         CALL CROT( N, Z( 1, JROW ), 1, Z( 1, JROW-1 ), 1, C, S )
351   30    CONTINUE
352   40 CONTINUE
353*
354      RETURN
355*
356*     End of CGGHRD
357*
358      END
359