|
2480
|
1 SUBROUTINE DGGHRD( COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, |
|
|
2 $ LDQ, Z, LDZ, INFO ) |
|
|
3 * |
|
7034
|
4 * -- LAPACK routine (version 3.1) -- |
|
|
5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. |
|
|
6 * November 2006 |
|
2480
|
7 * |
|
|
8 * .. Scalar Arguments .. |
|
|
9 CHARACTER COMPQ, COMPZ |
|
|
10 INTEGER IHI, ILO, INFO, LDA, LDB, LDQ, LDZ, N |
|
|
11 * .. |
|
|
12 * .. Array Arguments .. |
|
|
13 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), Q( LDQ, * ), |
|
|
14 $ Z( LDZ, * ) |
|
|
15 * .. |
|
|
16 * |
|
|
17 * Purpose |
|
|
18 * ======= |
|
|
19 * |
|
|
20 * DGGHRD reduces a pair of real matrices (A,B) to generalized upper |
|
|
21 * Hessenberg form using orthogonal transformations, where A is a |
|
7034
|
22 * general matrix and B is upper triangular. The form of the |
|
|
23 * generalized eigenvalue problem is |
|
|
24 * A*x = lambda*B*x, |
|
|
25 * and B is typically made upper triangular by computing its QR |
|
|
26 * factorization and moving the orthogonal matrix Q to the left side |
|
|
27 * of the equation. |
|
|
28 * |
|
|
29 * This subroutine simultaneously reduces A to a Hessenberg matrix H: |
|
|
30 * Q**T*A*Z = H |
|
|
31 * and transforms B to another upper triangular matrix T: |
|
|
32 * Q**T*B*Z = T |
|
|
33 * in order to reduce the problem to its standard form |
|
|
34 * H*y = lambda*T*y |
|
|
35 * where y = Z**T*x. |
|
2480
|
36 * |
|
|
37 * The orthogonal matrices Q and Z are determined as products of Givens |
|
|
38 * rotations. They may either be formed explicitly, or they may be |
|
|
39 * postmultiplied into input matrices Q1 and Z1, so that |
|
|
40 * |
|
7034
|
41 * Q1 * A * Z1**T = (Q1*Q) * H * (Z1*Z)**T |
|
|
42 * |
|
|
43 * Q1 * B * Z1**T = (Q1*Q) * T * (Z1*Z)**T |
|
|
44 * |
|
|
45 * If Q1 is the orthogonal matrix from the QR factorization of B in the |
|
|
46 * original equation A*x = lambda*B*x, then DGGHRD reduces the original |
|
|
47 * problem to generalized Hessenberg form. |
|
2480
|
48 * |
|
|
49 * Arguments |
|
|
50 * ========= |
|
|
51 * |
|
|
52 * COMPQ (input) CHARACTER*1 |
|
|
53 * = 'N': do not compute Q; |
|
|
54 * = 'I': Q is initialized to the unit matrix, and the |
|
|
55 * orthogonal matrix Q is returned; |
|
|
56 * = 'V': Q must contain an orthogonal matrix Q1 on entry, |
|
|
57 * and the product Q1*Q is returned. |
|
|
58 * |
|
|
59 * COMPZ (input) CHARACTER*1 |
|
|
60 * = 'N': do not compute Z; |
|
|
61 * = 'I': Z is initialized to the unit matrix, and the |
|
|
62 * orthogonal matrix Z is returned; |
|
|
63 * = 'V': Z must contain an orthogonal matrix Z1 on entry, |
|
|
64 * and the product Z1*Z is returned. |
|
|
65 * |
|
|
66 * N (input) INTEGER |
|
|
67 * The order of the matrices A and B. N >= 0. |
|
|
68 * |
|
|
69 * ILO (input) INTEGER |
|
|
70 * IHI (input) INTEGER |
|
7034
|
71 * ILO and IHI mark the rows and columns of A which are to be |
|
|
72 * reduced. It is assumed that A is already upper triangular |
|
|
73 * in rows and columns 1:ILO-1 and IHI+1:N. ILO and IHI are |
|
|
74 * normally set by a previous call to SGGBAL; otherwise they |
|
|
75 * should be set to 1 and N respectively. |
|
2480
|
76 * 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0. |
|
|
77 * |
|
|
78 * A (input/output) DOUBLE PRECISION array, dimension (LDA, N) |
|
|
79 * On entry, the N-by-N general matrix to be reduced. |
|
|
80 * On exit, the upper triangle and the first subdiagonal of A |
|
|
81 * are overwritten with the upper Hessenberg matrix H, and the |
|
|
82 * rest is set to zero. |
|
|
83 * |
|
|
84 * LDA (input) INTEGER |
|
|
85 * The leading dimension of the array A. LDA >= max(1,N). |
|
|
86 * |
|
|
87 * B (input/output) DOUBLE PRECISION array, dimension (LDB, N) |
|
|
88 * On entry, the N-by-N upper triangular matrix B. |
|
7034
|
89 * On exit, the upper triangular matrix T = Q**T B Z. The |
|
2480
|
90 * elements below the diagonal are set to zero. |
|
|
91 * |
|
|
92 * LDB (input) INTEGER |
|
|
93 * The leading dimension of the array B. LDB >= max(1,N). |
|
|
94 * |
|
|
95 * Q (input/output) DOUBLE PRECISION array, dimension (LDQ, N) |
|
7034
|
96 * On entry, if COMPQ = 'V', the orthogonal matrix Q1, |
|
|
97 * typically from the QR factorization of B. |
|
|
98 * On exit, if COMPQ='I', the orthogonal matrix Q, and if |
|
|
99 * COMPQ = 'V', the product Q1*Q. |
|
|
100 * Not referenced if COMPQ='N'. |
|
2480
|
101 * |
|
|
102 * LDQ (input) INTEGER |
|
|
103 * The leading dimension of the array Q. |
|
|
104 * LDQ >= N if COMPQ='V' or 'I'; LDQ >= 1 otherwise. |
|
|
105 * |
|
|
106 * Z (input/output) DOUBLE PRECISION array, dimension (LDZ, N) |
|
7034
|
107 * On entry, if COMPZ = 'V', the orthogonal matrix Z1. |
|
|
108 * On exit, if COMPZ='I', the orthogonal matrix Z, and if |
|
|
109 * COMPZ = 'V', the product Z1*Z. |
|
|
110 * Not referenced if COMPZ='N'. |
|
2480
|
111 * |
|
|
112 * LDZ (input) INTEGER |
|
|
113 * The leading dimension of the array Z. |
|
|
114 * LDZ >= N if COMPZ='V' or 'I'; LDZ >= 1 otherwise. |
|
|
115 * |
|
|
116 * INFO (output) INTEGER |
|
|
117 * = 0: successful exit. |
|
|
118 * < 0: if INFO = -i, the i-th argument had an illegal value. |
|
|
119 * |
|
|
120 * Further Details |
|
|
121 * =============== |
|
|
122 * |
|
|
123 * This routine reduces A to Hessenberg and B to triangular form by |
|
|
124 * an unblocked reduction, as described in _Matrix_Computations_, |
|
|
125 * by Golub and Van Loan (Johns Hopkins Press.) |
|
|
126 * |
|
|
127 * ===================================================================== |
|
|
128 * |
|
|
129 * .. Parameters .. |
|
|
130 DOUBLE PRECISION ONE, ZERO |
|
|
131 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) |
|
|
132 * .. |
|
|
133 * .. Local Scalars .. |
|
|
134 LOGICAL ILQ, ILZ |
|
|
135 INTEGER ICOMPQ, ICOMPZ, JCOL, JROW |
|
|
136 DOUBLE PRECISION C, S, TEMP |
|
|
137 * .. |
|
|
138 * .. External Functions .. |
|
|
139 LOGICAL LSAME |
|
|
140 EXTERNAL LSAME |
|
|
141 * .. |
|
|
142 * .. External Subroutines .. |
|
|
143 EXTERNAL DLARTG, DLASET, DROT, XERBLA |
|
|
144 * .. |
|
|
145 * .. Intrinsic Functions .. |
|
|
146 INTRINSIC 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( 'DGGHRD', -INFO ) |
|
|
204 RETURN |
|
|
205 END IF |
|
|
206 * |
|
|
207 * Initialize Q and Z if desired. |
|
|
208 * |
|
|
209 IF( ICOMPQ.EQ.3 ) |
|
|
210 $ CALL DLASET( 'Full', N, N, ZERO, ONE, Q, LDQ ) |
|
|
211 IF( ICOMPZ.EQ.3 ) |
|
|
212 $ CALL DLASET( 'Full', N, N, ZERO, ONE, 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 ) = ZERO |
|
|
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 TEMP = A( JROW-1, JCOL ) |
|
|
236 CALL DLARTG( TEMP, A( JROW, JCOL ), C, S, |
|
|
237 $ A( JROW-1, JCOL ) ) |
|
|
238 A( JROW, JCOL ) = ZERO |
|
|
239 CALL DROT( N-JCOL, A( JROW-1, JCOL+1 ), LDA, |
|
|
240 $ A( JROW, JCOL+1 ), LDA, C, S ) |
|
|
241 CALL DROT( N+2-JROW, B( JROW-1, JROW-1 ), LDB, |
|
|
242 $ B( JROW, JROW-1 ), LDB, C, S ) |
|
|
243 IF( ILQ ) |
|
|
244 $ CALL DROT( N, Q( 1, JROW-1 ), 1, Q( 1, JROW ), 1, C, S ) |
|
|
245 * |
|
|
246 * Step 2: rotate columns JROW, JROW-1 to kill B(JROW,JROW-1) |
|
|
247 * |
|
|
248 TEMP = B( JROW, JROW ) |
|
|
249 CALL DLARTG( TEMP, B( JROW, JROW-1 ), C, S, |
|
|
250 $ B( JROW, JROW ) ) |
|
|
251 B( JROW, JROW-1 ) = ZERO |
|
|
252 CALL DROT( IHI, A( 1, JROW ), 1, A( 1, JROW-1 ), 1, C, S ) |
|
|
253 CALL DROT( JROW-1, B( 1, JROW ), 1, B( 1, JROW-1 ), 1, C, |
|
|
254 $ S ) |
|
|
255 IF( ILZ ) |
|
|
256 $ CALL DROT( N, Z( 1, JROW ), 1, Z( 1, JROW-1 ), 1, C, S ) |
|
|
257 30 CONTINUE |
|
|
258 40 CONTINUE |
|
|
259 * |
|
|
260 RETURN |
|
|
261 * |
|
|
262 * End of DGGHRD |
|
|
263 * |
|
|
264 END |
