3 * =========== DOCUMENTATION ===========
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
9 *> Download CGGHRD + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgghrd.f">
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgghrd.f">
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgghrd.f">
21 * SUBROUTINE CGGHRD( COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q,
24 * .. Scalar Arguments ..
25 * CHARACTER COMPQ, COMPZ
26 * INTEGER IHI, ILO, INFO, LDA, LDB, LDQ, LDZ, N
28 * .. Array Arguments ..
29 * COMPLEX A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
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
44 *> and B is typically made upper triangular by computing its QR
45 *> factorization and moving the unitary matrix Q to the left side
48 *> This subroutine simultaneously reduces A to a Hessenberg matrix H:
50 *> and transforms B to another upper triangular matrix T:
52 *> in order to reduce the problem to its standard form
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.
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.
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.
92 *> The order of the matrices A and B. N >= 0.
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.
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.
124 *> The leading dimension of the array A. LDA >= max(1,N).
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.
138 *> The leading dimension of the array B. LDB >= max(1,N).
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'.
154 *> The leading dimension of the array Q.
155 *> LDQ >= N if COMPQ='V' or 'I'; LDQ >= 1 otherwise.
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'.
170 *> The leading dimension of the array Z.
171 *> LDZ >= N if COMPZ='V' or 'I'; LDZ >= 1 otherwise.
177 *> = 0: successful exit.
178 *> < 0: if INFO = -i, the i-th argument had an illegal value.
184 *> \author Univ. of Tennessee
185 *> \author Univ. of California Berkeley
186 *> \author Univ. of Colorado Denver
189 *> \date November 2015
191 *> \ingroup complexOTHERcomputational
193 *> \par Further Details:
194 * =====================
198 *> This routine reduces A to Hessenberg and B to triangular form by
199 *> an unblocked reduction, as described in _Matrix_Computations_,
200 *> by Golub and van Loan (Johns Hopkins Press).
203 * =====================================================================
204 SUBROUTINE CGGHRD( COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q,
205 $ LDQ, Z, LDZ, INFO )
207 * -- LAPACK computational routine (version 3.6.0) --
208 * -- LAPACK is a software package provided by Univ. of Tennessee, --
209 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
212 * .. Scalar Arguments ..
213 CHARACTER COMPQ, COMPZ
214 INTEGER IHI, ILO, INFO, LDA, LDB, LDQ, LDZ, N
216 * .. Array Arguments ..
217 COMPLEX A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
221 * =====================================================================
225 PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ),
226 $ CZERO = ( 0.0E+0, 0.0E+0 ) )
228 * .. Local Scalars ..
230 INTEGER ICOMPQ, ICOMPZ, JCOL, JROW
234 * .. External Functions ..
238 * .. External Subroutines ..
239 EXTERNAL CLARTG, CLASET, CROT, XERBLA
241 * .. Intrinsic Functions ..
244 * .. Executable Statements ..
248 IF( LSAME( COMPQ, 'N' ) ) THEN
251 ELSE IF( LSAME( COMPQ, 'V' ) ) THEN
254 ELSE IF( LSAME( COMPQ, 'I' ) ) THEN
263 IF( LSAME( COMPZ, 'N' ) ) THEN
266 ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
269 ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
276 * Test the input parameters.
279 IF( ICOMPQ.LE.0 ) THEN
281 ELSE IF( ICOMPZ.LE.0 ) THEN
283 ELSE IF( N.LT.0 ) THEN
285 ELSE IF( ILO.LT.1 ) THEN
287 ELSE IF( IHI.GT.N .OR. IHI.LT.ILO-1 ) THEN
289 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
291 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
293 ELSE IF( ( ILQ .AND. LDQ.LT.N ) .OR. LDQ.LT.1 ) THEN
295 ELSE IF( ( ILZ .AND. LDZ.LT.N ) .OR. LDZ.LT.1 ) THEN
299 CALL XERBLA( 'CGGHRD', -INFO )
303 * Initialize Q and Z if desired.
306 $ CALL CLASET( 'Full', N, N, CZERO, CONE, Q, LDQ )
308 $ CALL CLASET( 'Full', N, N, CZERO, CONE, Z, LDZ )
310 * Quick return if possible
315 * Zero out lower triangle of B
317 DO 20 JCOL = 1, N - 1
318 DO 10 JROW = JCOL + 1, N
319 B( JROW, JCOL ) = CZERO
325 DO 40 JCOL = ILO, IHI - 2
327 DO 30 JROW = IHI, JCOL + 2, -1
329 * Step 1: rotate rows JROW-1, JROW to kill A(JROW,JCOL)
331 CTEMP = A( JROW-1, JCOL )
332 CALL CLARTG( CTEMP, A( JROW, JCOL ), C, S,
333 $ A( JROW-1, JCOL ) )
334 A( JROW, JCOL ) = CZERO
335 CALL CROT( N-JCOL, A( JROW-1, JCOL+1 ), LDA,
336 $ A( JROW, JCOL+1 ), LDA, C, S )
337 CALL CROT( N+2-JROW, B( JROW-1, JROW-1 ), LDB,
338 $ B( JROW, JROW-1 ), LDB, C, S )
340 $ CALL CROT( N, Q( 1, JROW-1 ), 1, Q( 1, JROW ), 1, C,
343 * Step 2: rotate columns JROW, JROW-1 to kill B(JROW,JROW-1)
345 CTEMP = B( JROW, JROW )
346 CALL CLARTG( CTEMP, B( JROW, JROW-1 ), C, S,
348 B( JROW, JROW-1 ) = CZERO
349 CALL CROT( IHI, A( 1, JROW ), 1, A( 1, JROW-1 ), 1, C, S )
350 CALL CROT( JROW-1, B( 1, JROW ), 1, B( 1, JROW-1 ), 1, C,
353 $ CALL CROT( N, Z( 1, JROW ), 1, Z( 1, JROW-1 ), 1, C, S )