Annotation of rpl/lapack/lapack/zgghrd.f, revision 1.1.1.1
1.1 bertrand 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
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