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