1 SUBROUTINE ZGGHRD( COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q,
2 $ LDQ, Z, LDZ, INFO )
3 *
4 * -- LAPACK routine (version 3.2) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * November 2006
8 *
9 * .. Scalar Arguments ..
10 CHARACTER COMPQ, COMPZ
11 INTEGER IHI, ILO, INFO, LDA, LDB, LDQ, LDZ, N
12 * ..
13 * .. Array Arguments ..
14 COMPLEX*16 A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
15 $ Z( LDZ, * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * ZGGHRD reduces a pair of complex matrices (A,B) to generalized upper
22 * Hessenberg form using unitary transformations, where A is a
23 * general matrix and B is upper triangular. The form of the
24 * generalized eigenvalue problem is
25 * A*x = lambda*B*x,
26 * and B is typically made upper triangular by computing its QR
27 * factorization and moving the unitary matrix Q to the left side
28 * of the equation.
29 *
30 * This subroutine simultaneously reduces A to a Hessenberg matrix H:
31 * Q**H*A*Z = H
32 * and transforms B to another upper triangular matrix T:
33 * Q**H*B*Z = T
34 * in order to reduce the problem to its standard form
35 * H*y = lambda*T*y
36 * where y = Z**H*x.
37 *
38 * The unitary matrices Q and Z are determined as products of Givens
39 * rotations. They may either be formed explicitly, or they may be
40 * postmultiplied into input matrices Q1 and Z1, so that
41 * Q1 * A * Z1**H = (Q1*Q) * H * (Z1*Z)**H
42 * Q1 * B * Z1**H = (Q1*Q) * T * (Z1*Z)**H
43 * If Q1 is the unitary matrix from the QR factorization of B in the
44 * original equation A*x = lambda*B*x, then ZGGHRD reduces the original
45 * problem to generalized Hessenberg form.
46 *
47 * Arguments
48 * =========
49 *
50 * COMPQ (input) CHARACTER*1
51 * = 'N': do not compute Q;
52 * = 'I': Q is initialized to the unit matrix, and the
53 * unitary matrix Q is returned;
54 * = 'V': Q must contain a unitary matrix Q1 on entry,
55 * and the product Q1*Q is returned.
56 *
57 * COMPZ (input) CHARACTER*1
58 * = 'N': do not compute Q;
59 * = 'I': Q is initialized to the unit matrix, and the
60 * unitary matrix Q is returned;
61 * = 'V': Q must contain a unitary matrix Q1 on entry,
62 * and the product Q1*Q is returned.
63 *
64 * N (input) INTEGER
65 * The order of the matrices A and B. N >= 0.
66 *
67 * ILO (input) INTEGER
68 * IHI (input) INTEGER
69 * ILO and IHI mark the rows and columns of A which are to be
70 * reduced. It is assumed that A is already upper triangular
71 * in rows and columns 1:ILO-1 and IHI+1:N. ILO and IHI are
72 * normally set by a previous call to ZGGBAL; otherwise they
73 * should be set to 1 and N respectively.
74 * 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
75 *
76 * A (input/output) COMPLEX*16 array, dimension (LDA, N)
77 * On entry, the N-by-N general matrix to be reduced.
78 * On exit, the upper triangle and the first subdiagonal of A
79 * are overwritten with the upper Hessenberg matrix H, and the
80 * rest is set to zero.
81 *
82 * LDA (input) INTEGER
83 * The leading dimension of the array A. LDA >= max(1,N).
84 *
85 * B (input/output) COMPLEX*16 array, dimension (LDB, N)
86 * On entry, the N-by-N upper triangular matrix B.
87 * On exit, the upper triangular matrix T = Q**H B Z. The
88 * elements below the diagonal are set to zero.
89 *
90 * LDB (input) INTEGER
91 * The leading dimension of the array B. LDB >= max(1,N).
92 *
93 * Q (input/output) COMPLEX*16 array, dimension (LDQ, N)
94 * On entry, if COMPQ = 'V', the unitary matrix Q1, typically
95 * from the QR factorization of B.
96 * On exit, if COMPQ='I', the unitary matrix Q, and if
97 * COMPQ = 'V', the product Q1*Q.
98 * Not referenced if COMPQ='N'.
99 *
100 * LDQ (input) INTEGER
101 * The leading dimension of the array Q.
102 * LDQ >= N if COMPQ='V' or 'I'; LDQ >= 1 otherwise.
103 *
104 * Z (input/output) COMPLEX*16 array, dimension (LDZ, N)
105 * On entry, if COMPZ = 'V', the unitary matrix Z1.
106 * On exit, if COMPZ='I', the unitary matrix Z, and if
107 * COMPZ = 'V', the product Z1*Z.
108 * Not referenced if COMPZ='N'.
109 *
110 * LDZ (input) INTEGER
111 * The leading dimension of the array Z.
112 * LDZ >= N if COMPZ='V' or 'I'; LDZ >= 1 otherwise.
113 *
114 * INFO (output) INTEGER
115 * = 0: successful exit.
116 * < 0: if INFO = -i, the i-th argument had an illegal value.
117 *
118 * Further Details
119 * ===============
120 *
121 * This routine reduces A to Hessenberg and B to triangular form by
122 * an unblocked reduction, as described in _Matrix_Computations_,
123 * by Golub and van Loan (Johns Hopkins Press).
124 *
125 * =====================================================================
126 *
127 * .. Parameters ..
128 COMPLEX*16 CONE, CZERO
129 PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ),
130 $ CZERO = ( 0.0D+0, 0.0D+0 ) )
131 * ..
132 * .. Local Scalars ..
133 LOGICAL ILQ, ILZ
134 INTEGER ICOMPQ, ICOMPZ, JCOL, JROW
135 DOUBLE PRECISION C
136 COMPLEX*16 CTEMP, S
137 * ..
138 * .. External Functions ..
139 LOGICAL LSAME
140 EXTERNAL LSAME
141 * ..
142 * .. External Subroutines ..
143 EXTERNAL XERBLA, ZLARTG, ZLASET, ZROT
144 * ..
145 * .. Intrinsic Functions ..
146 INTRINSIC DCONJG, MAX
147 * ..
148 * .. Executable Statements ..
149 *
150 * Decode COMPQ
151 *
152 IF( LSAME( COMPQ, 'N' ) ) THEN
153 ILQ = .FALSE.
154 ICOMPQ = 1
155 ELSE IF( LSAME( COMPQ, 'V' ) ) THEN
156 ILQ = .TRUE.
157 ICOMPQ = 2
158 ELSE IF( LSAME( COMPQ, 'I' ) ) THEN
159 ILQ = .TRUE.
160 ICOMPQ = 3
161 ELSE
162 ICOMPQ = 0
163 END IF
164 *
165 * Decode COMPZ
166 *
167 IF( LSAME( COMPZ, 'N' ) ) THEN
168 ILZ = .FALSE.
169 ICOMPZ = 1
170 ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
171 ILZ = .TRUE.
172 ICOMPZ = 2
173 ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
174 ILZ = .TRUE.
175 ICOMPZ = 3
176 ELSE
177 ICOMPZ = 0
178 END IF
179 *
180 * Test the input parameters.
181 *
182 INFO = 0
183 IF( ICOMPQ.LE.0 ) THEN
184 INFO = -1
185 ELSE IF( ICOMPZ.LE.0 ) THEN
186 INFO = -2
187 ELSE IF( N.LT.0 ) THEN
188 INFO = -3
189 ELSE IF( ILO.LT.1 ) THEN
190 INFO = -4
191 ELSE IF( IHI.GT.N .OR. IHI.LT.ILO-1 ) THEN
192 INFO = -5
193 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
194 INFO = -7
195 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
196 INFO = -9
197 ELSE IF( ( ILQ .AND. LDQ.LT.N ) .OR. LDQ.LT.1 ) THEN
198 INFO = -11
199 ELSE IF( ( ILZ .AND. LDZ.LT.N ) .OR. LDZ.LT.1 ) THEN
200 INFO = -13
201 END IF
202 IF( INFO.NE.0 ) THEN
203 CALL XERBLA( 'ZGGHRD', -INFO )
204 RETURN
205 END IF
206 *
207 * Initialize Q and Z if desired.
208 *
209 IF( ICOMPQ.EQ.3 )
210 $ CALL ZLASET( 'Full', N, N, CZERO, CONE, Q, LDQ )
211 IF( ICOMPZ.EQ.3 )
212 $ CALL ZLASET( 'Full', N, N, CZERO, CONE, Z, LDZ )
213 *
214 * Quick return if possible
215 *
216 IF( N.LE.1 )
217 $ RETURN
218 *
219 * Zero out lower triangle of B
220 *
221 DO 20 JCOL = 1, N - 1
222 DO 10 JROW = JCOL + 1, N
223 B( JROW, JCOL ) = CZERO
224 10 CONTINUE
225 20 CONTINUE
226 *
227 * Reduce A and B
228 *
229 DO 40 JCOL = ILO, IHI - 2
230 *
231 DO 30 JROW = IHI, JCOL + 2, -1
232 *
233 * Step 1: rotate rows JROW-1, JROW to kill A(JROW,JCOL)
234 *
235 CTEMP = A( JROW-1, JCOL )
236 CALL ZLARTG( CTEMP, A( JROW, JCOL ), C, S,
237 $ A( JROW-1, JCOL ) )
238 A( JROW, JCOL ) = CZERO
239 CALL ZROT( N-JCOL, A( JROW-1, JCOL+1 ), LDA,
240 $ A( JROW, JCOL+1 ), LDA, C, S )
241 CALL ZROT( N+2-JROW, B( JROW-1, JROW-1 ), LDB,
242 $ B( JROW, JROW-1 ), LDB, C, S )
243 IF( ILQ )
244 $ CALL ZROT( N, Q( 1, JROW-1 ), 1, Q( 1, JROW ), 1, C,
245 $ DCONJG( S ) )
246 *
247 * Step 2: rotate columns JROW, JROW-1 to kill B(JROW,JROW-1)
248 *
249 CTEMP = B( JROW, JROW )
250 CALL ZLARTG( CTEMP, B( JROW, JROW-1 ), C, S,
251 $ B( JROW, JROW ) )
252 B( JROW, JROW-1 ) = CZERO
253 CALL ZROT( IHI, A( 1, JROW ), 1, A( 1, JROW-1 ), 1, C, S )
254 CALL ZROT( JROW-1, B( 1, JROW ), 1, B( 1, JROW-1 ), 1, C,
255 $ S )
256 IF( ILZ )
257 $ CALL ZROT( N, Z( 1, JROW ), 1, Z( 1, JROW-1 ), 1, C, S )
258 30 CONTINUE
259 40 CONTINUE
260 *
261 RETURN
262 *
263 * End of ZGGHRD
264 *
265 END
2 $ LDQ, Z, LDZ, INFO )
3 *
4 * -- LAPACK routine (version 3.2) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * November 2006
8 *
9 * .. Scalar Arguments ..
10 CHARACTER COMPQ, COMPZ
11 INTEGER IHI, ILO, INFO, LDA, LDB, LDQ, LDZ, N
12 * ..
13 * .. Array Arguments ..
14 COMPLEX*16 A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
15 $ Z( LDZ, * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * ZGGHRD reduces a pair of complex matrices (A,B) to generalized upper
22 * Hessenberg form using unitary transformations, where A is a
23 * general matrix and B is upper triangular. The form of the
24 * generalized eigenvalue problem is
25 * A*x = lambda*B*x,
26 * and B is typically made upper triangular by computing its QR
27 * factorization and moving the unitary matrix Q to the left side
28 * of the equation.
29 *
30 * This subroutine simultaneously reduces A to a Hessenberg matrix H:
31 * Q**H*A*Z = H
32 * and transforms B to another upper triangular matrix T:
33 * Q**H*B*Z = T
34 * in order to reduce the problem to its standard form
35 * H*y = lambda*T*y
36 * where y = Z**H*x.
37 *
38 * The unitary matrices Q and Z are determined as products of Givens
39 * rotations. They may either be formed explicitly, or they may be
40 * postmultiplied into input matrices Q1 and Z1, so that
41 * Q1 * A * Z1**H = (Q1*Q) * H * (Z1*Z)**H
42 * Q1 * B * Z1**H = (Q1*Q) * T * (Z1*Z)**H
43 * If Q1 is the unitary matrix from the QR factorization of B in the
44 * original equation A*x = lambda*B*x, then ZGGHRD reduces the original
45 * problem to generalized Hessenberg form.
46 *
47 * Arguments
48 * =========
49 *
50 * COMPQ (input) CHARACTER*1
51 * = 'N': do not compute Q;
52 * = 'I': Q is initialized to the unit matrix, and the
53 * unitary matrix Q is returned;
54 * = 'V': Q must contain a unitary matrix Q1 on entry,
55 * and the product Q1*Q is returned.
56 *
57 * COMPZ (input) CHARACTER*1
58 * = 'N': do not compute Q;
59 * = 'I': Q is initialized to the unit matrix, and the
60 * unitary matrix Q is returned;
61 * = 'V': Q must contain a unitary matrix Q1 on entry,
62 * and the product Q1*Q is returned.
63 *
64 * N (input) INTEGER
65 * The order of the matrices A and B. N >= 0.
66 *
67 * ILO (input) INTEGER
68 * IHI (input) INTEGER
69 * ILO and IHI mark the rows and columns of A which are to be
70 * reduced. It is assumed that A is already upper triangular
71 * in rows and columns 1:ILO-1 and IHI+1:N. ILO and IHI are
72 * normally set by a previous call to ZGGBAL; otherwise they
73 * should be set to 1 and N respectively.
74 * 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
75 *
76 * A (input/output) COMPLEX*16 array, dimension (LDA, N)
77 * On entry, the N-by-N general matrix to be reduced.
78 * On exit, the upper triangle and the first subdiagonal of A
79 * are overwritten with the upper Hessenberg matrix H, and the
80 * rest is set to zero.
81 *
82 * LDA (input) INTEGER
83 * The leading dimension of the array A. LDA >= max(1,N).
84 *
85 * B (input/output) COMPLEX*16 array, dimension (LDB, N)
86 * On entry, the N-by-N upper triangular matrix B.
87 * On exit, the upper triangular matrix T = Q**H B Z. The
88 * elements below the diagonal are set to zero.
89 *
90 * LDB (input) INTEGER
91 * The leading dimension of the array B. LDB >= max(1,N).
92 *
93 * Q (input/output) COMPLEX*16 array, dimension (LDQ, N)
94 * On entry, if COMPQ = 'V', the unitary matrix Q1, typically
95 * from the QR factorization of B.
96 * On exit, if COMPQ='I', the unitary matrix Q, and if
97 * COMPQ = 'V', the product Q1*Q.
98 * Not referenced if COMPQ='N'.
99 *
100 * LDQ (input) INTEGER
101 * The leading dimension of the array Q.
102 * LDQ >= N if COMPQ='V' or 'I'; LDQ >= 1 otherwise.
103 *
104 * Z (input/output) COMPLEX*16 array, dimension (LDZ, N)
105 * On entry, if COMPZ = 'V', the unitary matrix Z1.
106 * On exit, if COMPZ='I', the unitary matrix Z, and if
107 * COMPZ = 'V', the product Z1*Z.
108 * Not referenced if COMPZ='N'.
109 *
110 * LDZ (input) INTEGER
111 * The leading dimension of the array Z.
112 * LDZ >= N if COMPZ='V' or 'I'; LDZ >= 1 otherwise.
113 *
114 * INFO (output) INTEGER
115 * = 0: successful exit.
116 * < 0: if INFO = -i, the i-th argument had an illegal value.
117 *
118 * Further Details
119 * ===============
120 *
121 * This routine reduces A to Hessenberg and B to triangular form by
122 * an unblocked reduction, as described in _Matrix_Computations_,
123 * by Golub and van Loan (Johns Hopkins Press).
124 *
125 * =====================================================================
126 *
127 * .. Parameters ..
128 COMPLEX*16 CONE, CZERO
129 PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ),
130 $ CZERO = ( 0.0D+0, 0.0D+0 ) )
131 * ..
132 * .. Local Scalars ..
133 LOGICAL ILQ, ILZ
134 INTEGER ICOMPQ, ICOMPZ, JCOL, JROW
135 DOUBLE PRECISION C
136 COMPLEX*16 CTEMP, S
137 * ..
138 * .. External Functions ..
139 LOGICAL LSAME
140 EXTERNAL LSAME
141 * ..
142 * .. External Subroutines ..
143 EXTERNAL XERBLA, ZLARTG, ZLASET, ZROT
144 * ..
145 * .. Intrinsic Functions ..
146 INTRINSIC DCONJG, MAX
147 * ..
148 * .. Executable Statements ..
149 *
150 * Decode COMPQ
151 *
152 IF( LSAME( COMPQ, 'N' ) ) THEN
153 ILQ = .FALSE.
154 ICOMPQ = 1
155 ELSE IF( LSAME( COMPQ, 'V' ) ) THEN
156 ILQ = .TRUE.
157 ICOMPQ = 2
158 ELSE IF( LSAME( COMPQ, 'I' ) ) THEN
159 ILQ = .TRUE.
160 ICOMPQ = 3
161 ELSE
162 ICOMPQ = 0
163 END IF
164 *
165 * Decode COMPZ
166 *
167 IF( LSAME( COMPZ, 'N' ) ) THEN
168 ILZ = .FALSE.
169 ICOMPZ = 1
170 ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
171 ILZ = .TRUE.
172 ICOMPZ = 2
173 ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
174 ILZ = .TRUE.
175 ICOMPZ = 3
176 ELSE
177 ICOMPZ = 0
178 END IF
179 *
180 * Test the input parameters.
181 *
182 INFO = 0
183 IF( ICOMPQ.LE.0 ) THEN
184 INFO = -1
185 ELSE IF( ICOMPZ.LE.0 ) THEN
186 INFO = -2
187 ELSE IF( N.LT.0 ) THEN
188 INFO = -3
189 ELSE IF( ILO.LT.1 ) THEN
190 INFO = -4
191 ELSE IF( IHI.GT.N .OR. IHI.LT.ILO-1 ) THEN
192 INFO = -5
193 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
194 INFO = -7
195 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
196 INFO = -9
197 ELSE IF( ( ILQ .AND. LDQ.LT.N ) .OR. LDQ.LT.1 ) THEN
198 INFO = -11
199 ELSE IF( ( ILZ .AND. LDZ.LT.N ) .OR. LDZ.LT.1 ) THEN
200 INFO = -13
201 END IF
202 IF( INFO.NE.0 ) THEN
203 CALL XERBLA( 'ZGGHRD', -INFO )
204 RETURN
205 END IF
206 *
207 * Initialize Q and Z if desired.
208 *
209 IF( ICOMPQ.EQ.3 )
210 $ CALL ZLASET( 'Full', N, N, CZERO, CONE, Q, LDQ )
211 IF( ICOMPZ.EQ.3 )
212 $ CALL ZLASET( 'Full', N, N, CZERO, CONE, Z, LDZ )
213 *
214 * Quick return if possible
215 *
216 IF( N.LE.1 )
217 $ RETURN
218 *
219 * Zero out lower triangle of B
220 *
221 DO 20 JCOL = 1, N - 1
222 DO 10 JROW = JCOL + 1, N
223 B( JROW, JCOL ) = CZERO
224 10 CONTINUE
225 20 CONTINUE
226 *
227 * Reduce A and B
228 *
229 DO 40 JCOL = ILO, IHI - 2
230 *
231 DO 30 JROW = IHI, JCOL + 2, -1
232 *
233 * Step 1: rotate rows JROW-1, JROW to kill A(JROW,JCOL)
234 *
235 CTEMP = A( JROW-1, JCOL )
236 CALL ZLARTG( CTEMP, A( JROW, JCOL ), C, S,
237 $ A( JROW-1, JCOL ) )
238 A( JROW, JCOL ) = CZERO
239 CALL ZROT( N-JCOL, A( JROW-1, JCOL+1 ), LDA,
240 $ A( JROW, JCOL+1 ), LDA, C, S )
241 CALL ZROT( N+2-JROW, B( JROW-1, JROW-1 ), LDB,
242 $ B( JROW, JROW-1 ), LDB, C, S )
243 IF( ILQ )
244 $ CALL ZROT( N, Q( 1, JROW-1 ), 1, Q( 1, JROW ), 1, C,
245 $ DCONJG( S ) )
246 *
247 * Step 2: rotate columns JROW, JROW-1 to kill B(JROW,JROW-1)
248 *
249 CTEMP = B( JROW, JROW )
250 CALL ZLARTG( CTEMP, B( JROW, JROW-1 ), C, S,
251 $ B( JROW, JROW ) )
252 B( JROW, JROW-1 ) = CZERO
253 CALL ZROT( IHI, A( 1, JROW ), 1, A( 1, JROW-1 ), 1, C, S )
254 CALL ZROT( JROW-1, B( 1, JROW ), 1, B( 1, JROW-1 ), 1, C,
255 $ S )
256 IF( ILZ )
257 $ CALL ZROT( N, Z( 1, JROW ), 1, Z( 1, JROW-1 ), 1, C, S )
258 30 CONTINUE
259 40 CONTINUE
260 *
261 RETURN
262 *
263 * End of ZGGHRD
264 *
265 END