1: *> \brief \b ZGGHRD
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download ZGGHRD + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgghrd.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgghrd.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgghrd.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZGGHRD( 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*16 A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
30: * $ Z( LDZ, * )
31: * ..
32: *
33: *
34: *> \par Purpose:
35: * =============
36: *>
37: *> \verbatim
38: *>
39: *> ZGGHRD 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
42: *> generalized 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 ZGGHRD 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 ZGGBAL; 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*16 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*16 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*16 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*16 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 complex16OTHERcomputational
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 ZGGHRD( 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*16 A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
215: $ Z( LDZ, * )
216: * ..
217: *
218: * =====================================================================
219: *
220: * .. Parameters ..
221: COMPLEX*16 CONE, CZERO
222: PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ),
223: $ CZERO = ( 0.0D+0, 0.0D+0 ) )
224: * ..
225: * .. Local Scalars ..
226: LOGICAL ILQ, ILZ
227: INTEGER ICOMPQ, ICOMPZ, JCOL, JROW
228: DOUBLE PRECISION C
229: COMPLEX*16 CTEMP, S
230: * ..
231: * .. External Functions ..
232: LOGICAL LSAME
233: EXTERNAL LSAME
234: * ..
235: * .. External Subroutines ..
236: EXTERNAL XERBLA, ZLARTG, ZLASET, ZROT
237: * ..
238: * .. Intrinsic Functions ..
239: INTRINSIC DCONJG, 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( 'ZGGHRD', -INFO )
297: RETURN
298: END IF
299: *
300: * Initialize Q and Z if desired.
301: *
302: IF( ICOMPQ.EQ.3 )
303: $ CALL ZLASET( 'Full', N, N, CZERO, CONE, Q, LDQ )
304: IF( ICOMPZ.EQ.3 )
305: $ CALL ZLASET( '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 ZLARTG( CTEMP, A( JROW, JCOL ), C, S,
330: $ A( JROW-1, JCOL ) )
331: A( JROW, JCOL ) = CZERO
332: CALL ZROT( N-JCOL, A( JROW-1, JCOL+1 ), LDA,
333: $ A( JROW, JCOL+1 ), LDA, C, S )
334: CALL ZROT( N+2-JROW, B( JROW-1, JROW-1 ), LDB,
335: $ B( JROW, JROW-1 ), LDB, C, S )
336: IF( ILQ )
337: $ CALL ZROT( N, Q( 1, JROW-1 ), 1, Q( 1, JROW ), 1, C,
338: $ DCONJG( S ) )
339: *
340: * Step 2: rotate columns JROW, JROW-1 to kill B(JROW,JROW-1)
341: *
342: CTEMP = B( JROW, JROW )
343: CALL ZLARTG( CTEMP, B( JROW, JROW-1 ), C, S,
344: $ B( JROW, JROW ) )
345: B( JROW, JROW-1 ) = CZERO
346: CALL ZROT( IHI, A( 1, JROW ), 1, A( 1, JROW-1 ), 1, C, S )
347: CALL ZROT( JROW-1, B( 1, JROW ), 1, B( 1, JROW-1 ), 1, C,
348: $ S )
349: IF( ILZ )
350: $ CALL ZROT( 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 ZGGHRD
357: *
358: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zgghrd.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack