Annotation of rpl/lapack/lapack/zgghrd.f, revision 1.10
1.8 bertrand 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 Q;
83: *> = 'I': Q is initialized to the unit matrix, and the
84: *> unitary matrix Q is returned;
85: *> = 'V': Q must contain a unitary matrix Q1 on entry,
86: *> and the product Q1*Q 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: *> \date November 2011
190: *
191: *> \ingroup complex16OTHERcomputational
192: *
193: *> \par Further Details:
194: * =====================
195: *>
196: *> \verbatim
197: *>
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).
201: *> \endverbatim
202: *>
203: * =====================================================================
1.1 bertrand 204: SUBROUTINE ZGGHRD( COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q,
205: $ LDQ, Z, LDZ, INFO )
206: *
1.8 bertrand 207: * -- LAPACK computational routine (version 3.4.0) --
1.1 bertrand 208: * -- LAPACK is a software package provided by Univ. of Tennessee, --
209: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.8 bertrand 210: * November 2011
1.1 bertrand 211: *
212: * .. Scalar Arguments ..
213: CHARACTER COMPQ, COMPZ
214: INTEGER IHI, ILO, INFO, LDA, LDB, LDQ, LDZ, N
215: * ..
216: * .. Array Arguments ..
217: COMPLEX*16 A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
218: $ Z( LDZ, * )
219: * ..
220: *
221: * =====================================================================
222: *
223: * .. Parameters ..
224: COMPLEX*16 CONE, CZERO
225: PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ),
226: $ CZERO = ( 0.0D+0, 0.0D+0 ) )
227: * ..
228: * .. Local Scalars ..
229: LOGICAL ILQ, ILZ
230: INTEGER ICOMPQ, ICOMPZ, JCOL, JROW
231: DOUBLE PRECISION C
232: COMPLEX*16 CTEMP, S
233: * ..
234: * .. External Functions ..
235: LOGICAL LSAME
236: EXTERNAL LSAME
237: * ..
238: * .. External Subroutines ..
239: EXTERNAL XERBLA, ZLARTG, ZLASET, ZROT
240: * ..
241: * .. Intrinsic Functions ..
242: INTRINSIC DCONJG, MAX
243: * ..
244: * .. Executable Statements ..
245: *
246: * Decode COMPQ
247: *
248: IF( LSAME( COMPQ, 'N' ) ) THEN
249: ILQ = .FALSE.
250: ICOMPQ = 1
251: ELSE IF( LSAME( COMPQ, 'V' ) ) THEN
252: ILQ = .TRUE.
253: ICOMPQ = 2
254: ELSE IF( LSAME( COMPQ, 'I' ) ) THEN
255: ILQ = .TRUE.
256: ICOMPQ = 3
257: ELSE
258: ICOMPQ = 0
259: END IF
260: *
261: * Decode COMPZ
262: *
263: IF( LSAME( COMPZ, 'N' ) ) THEN
264: ILZ = .FALSE.
265: ICOMPZ = 1
266: ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
267: ILZ = .TRUE.
268: ICOMPZ = 2
269: ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
270: ILZ = .TRUE.
271: ICOMPZ = 3
272: ELSE
273: ICOMPZ = 0
274: END IF
275: *
276: * Test the input parameters.
277: *
278: INFO = 0
279: IF( ICOMPQ.LE.0 ) THEN
280: INFO = -1
281: ELSE IF( ICOMPZ.LE.0 ) THEN
282: INFO = -2
283: ELSE IF( N.LT.0 ) THEN
284: INFO = -3
285: ELSE IF( ILO.LT.1 ) THEN
286: INFO = -4
287: ELSE IF( IHI.GT.N .OR. IHI.LT.ILO-1 ) THEN
288: INFO = -5
289: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
290: INFO = -7
291: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
292: INFO = -9
293: ELSE IF( ( ILQ .AND. LDQ.LT.N ) .OR. LDQ.LT.1 ) THEN
294: INFO = -11
295: ELSE IF( ( ILZ .AND. LDZ.LT.N ) .OR. LDZ.LT.1 ) THEN
296: INFO = -13
297: END IF
298: IF( INFO.NE.0 ) THEN
299: CALL XERBLA( 'ZGGHRD', -INFO )
300: RETURN
301: END IF
302: *
303: * Initialize Q and Z if desired.
304: *
305: IF( ICOMPQ.EQ.3 )
306: $ CALL ZLASET( 'Full', N, N, CZERO, CONE, Q, LDQ )
307: IF( ICOMPZ.EQ.3 )
308: $ CALL ZLASET( 'Full', N, N, CZERO, CONE, Z, LDZ )
309: *
310: * Quick return if possible
311: *
312: IF( N.LE.1 )
313: $ RETURN
314: *
315: * Zero out lower triangle of B
316: *
317: DO 20 JCOL = 1, N - 1
318: DO 10 JROW = JCOL + 1, N
319: B( JROW, JCOL ) = CZERO
320: 10 CONTINUE
321: 20 CONTINUE
322: *
323: * Reduce A and B
324: *
325: DO 40 JCOL = ILO, IHI - 2
326: *
327: DO 30 JROW = IHI, JCOL + 2, -1
328: *
329: * Step 1: rotate rows JROW-1, JROW to kill A(JROW,JCOL)
330: *
331: CTEMP = A( JROW-1, JCOL )
332: CALL ZLARTG( CTEMP, A( JROW, JCOL ), C, S,
333: $ A( JROW-1, JCOL ) )
334: A( JROW, JCOL ) = CZERO
335: CALL ZROT( N-JCOL, A( JROW-1, JCOL+1 ), LDA,
336: $ A( JROW, JCOL+1 ), LDA, C, S )
337: CALL ZROT( N+2-JROW, B( JROW-1, JROW-1 ), LDB,
338: $ B( JROW, JROW-1 ), LDB, C, S )
339: IF( ILQ )
340: $ CALL ZROT( N, Q( 1, JROW-1 ), 1, Q( 1, JROW ), 1, C,
341: $ DCONJG( S ) )
342: *
343: * Step 2: rotate columns JROW, JROW-1 to kill B(JROW,JROW-1)
344: *
345: CTEMP = B( JROW, JROW )
346: CALL ZLARTG( CTEMP, B( JROW, JROW-1 ), C, S,
347: $ B( JROW, JROW ) )
348: B( JROW, JROW-1 ) = CZERO
349: CALL ZROT( IHI, A( 1, JROW ), 1, A( 1, JROW-1 ), 1, C, S )
350: CALL ZROT( JROW-1, B( 1, JROW ), 1, B( 1, JROW-1 ), 1, C,
351: $ S )
352: IF( ILZ )
353: $ CALL ZROT( N, Z( 1, JROW ), 1, Z( 1, JROW-1 ), 1, C, S )
354: 30 CONTINUE
355: 40 CONTINUE
356: *
357: RETURN
358: *
359: * End of ZGGHRD
360: *
361: 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