|
3274
|
1 SUBROUTINE DTGEVC( SIDE, HOWMNY, SELECT, N, A, LDA, B, LDB, VL, |
|
|
2 $ LDVL, VR, LDVR, MM, M, WORK, INFO ) |
|
|
3 * |
|
3333
|
4 * -- LAPACK routine (version 3.0) -- |
|
3274
|
5 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., |
|
|
6 * Courant Institute, Argonne National Lab, and Rice University |
|
3333
|
7 * June 30, 1999 |
|
3274
|
8 * |
|
|
9 * .. Scalar Arguments .. |
|
|
10 CHARACTER HOWMNY, SIDE |
|
|
11 INTEGER INFO, LDA, LDB, LDVL, LDVR, M, MM, N |
|
|
12 * .. |
|
|
13 * .. Array Arguments .. |
|
|
14 LOGICAL SELECT( * ) |
|
|
15 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), VL( LDVL, * ), |
|
|
16 $ VR( LDVR, * ), WORK( * ) |
|
|
17 * .. |
|
|
18 * |
|
|
19 * |
|
|
20 * Purpose |
|
|
21 * ======= |
|
|
22 * |
|
|
23 * DTGEVC computes some or all of the right and/or left generalized |
|
|
24 * eigenvectors of a pair of real upper triangular matrices (A,B). |
|
|
25 * |
|
|
26 * The right generalized eigenvector x and the left generalized |
|
|
27 * eigenvector y of (A,B) corresponding to a generalized eigenvalue |
|
|
28 * w are defined by: |
|
|
29 * |
|
|
30 * (A - wB) * x = 0 and y**H * (A - wB) = 0 |
|
|
31 * |
|
|
32 * where y**H denotes the conjugate tranpose of y. |
|
|
33 * |
|
|
34 * If an eigenvalue w is determined by zero diagonal elements of both A |
|
|
35 * and B, a unit vector is returned as the corresponding eigenvector. |
|
|
36 * |
|
|
37 * If all eigenvectors are requested, the routine may either return |
|
|
38 * the matrices X and/or Y of right or left eigenvectors of (A,B), or |
|
|
39 * the products Z*X and/or Q*Y, where Z and Q are input orthogonal |
|
|
40 * matrices. If (A,B) was obtained from the generalized real-Schur |
|
|
41 * factorization of an original pair of matrices |
|
|
42 * (A0,B0) = (Q*A*Z**H,Q*B*Z**H), |
|
|
43 * then Z*X and Q*Y are the matrices of right or left eigenvectors of |
|
|
44 * A. |
|
|
45 * |
|
|
46 * A must be block upper triangular, with 1-by-1 and 2-by-2 diagonal |
|
|
47 * blocks. Corresponding to each 2-by-2 diagonal block is a complex |
|
|
48 * conjugate pair of eigenvalues and eigenvectors; only one |
|
|
49 * eigenvector of the pair is computed, namely the one corresponding |
|
|
50 * to the eigenvalue with positive imaginary part. |
|
|
51 * |
|
|
52 * Arguments |
|
|
53 * ========= |
|
|
54 * |
|
|
55 * SIDE (input) CHARACTER*1 |
|
|
56 * = 'R': compute right eigenvectors only; |
|
|
57 * = 'L': compute left eigenvectors only; |
|
|
58 * = 'B': compute both right and left eigenvectors. |
|
|
59 * |
|
|
60 * HOWMNY (input) CHARACTER*1 |
|
|
61 * = 'A': compute all right and/or left eigenvectors; |
|
|
62 * = 'B': compute all right and/or left eigenvectors, and |
|
|
63 * backtransform them using the input matrices supplied |
|
|
64 * in VR and/or VL; |
|
|
65 * = 'S': compute selected right and/or left eigenvectors, |
|
|
66 * specified by the logical array SELECT. |
|
|
67 * |
|
|
68 * SELECT (input) LOGICAL array, dimension (N) |
|
|
69 * If HOWMNY='S', SELECT specifies the eigenvectors to be |
|
|
70 * computed. |
|
|
71 * If HOWMNY='A' or 'B', SELECT is not referenced. |
|
|
72 * To select the real eigenvector corresponding to the real |
|
|
73 * eigenvalue w(j), SELECT(j) must be set to .TRUE. To select |
|
|
74 * the complex eigenvector corresponding to a complex conjugate |
|
|
75 * pair w(j) and w(j+1), either SELECT(j) or SELECT(j+1) must |
|
|
76 * be set to .TRUE.. |
|
|
77 * |
|
|
78 * N (input) INTEGER |
|
|
79 * The order of the matrices A and B. N >= 0. |
|
|
80 * |
|
|
81 * A (input) DOUBLE PRECISION array, dimension (LDA,N) |
|
|
82 * The upper quasi-triangular matrix A. |
|
|
83 * |
|
|
84 * LDA (input) INTEGER |
|
|
85 * The leading dimension of array A. LDA >= max(1, N). |
|
|
86 * |
|
|
87 * B (input) DOUBLE PRECISION array, dimension (LDB,N) |
|
|
88 * The upper triangular matrix B. If A has a 2-by-2 diagonal |
|
|
89 * block, then the corresponding 2-by-2 block of B must be |
|
|
90 * diagonal with positive elements. |
|
|
91 * |
|
|
92 * LDB (input) INTEGER |
|
|
93 * The leading dimension of array B. LDB >= max(1,N). |
|
|
94 * |
|
|
95 * VL (input/output) DOUBLE PRECISION array, dimension (LDVL,MM) |
|
|
96 * On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must |
|
|
97 * contain an N-by-N matrix Q (usually the orthogonal matrix Q |
|
|
98 * of left Schur vectors returned by DHGEQZ). |
|
|
99 * On exit, if SIDE = 'L' or 'B', VL contains: |
|
|
100 * if HOWMNY = 'A', the matrix Y of left eigenvectors of (A,B); |
|
|
101 * if HOWMNY = 'B', the matrix Q*Y; |
|
|
102 * if HOWMNY = 'S', the left eigenvectors of (A,B) specified by |
|
|
103 * SELECT, stored consecutively in the columns of |
|
|
104 * VL, in the same order as their eigenvalues. |
|
|
105 * If SIDE = 'R', VL is not referenced. |
|
|
106 * |
|
|
107 * A complex eigenvector corresponding to a complex eigenvalue |
|
|
108 * is stored in two consecutive columns, the first holding the |
|
|
109 * real part, and the second the imaginary part. |
|
|
110 * |
|
|
111 * LDVL (input) INTEGER |
|
|
112 * The leading dimension of array VL. |
|
|
113 * LDVL >= max(1,N) if SIDE = 'L' or 'B'; LDVL >= 1 otherwise. |
|
|
114 * |
|
3333
|
115 * VR (input/output) DOUBLE PRECISION array, dimension (LDVR,MM) |
|
3274
|
116 * On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must |
|
|
117 * contain an N-by-N matrix Q (usually the orthogonal matrix Z |
|
|
118 * of right Schur vectors returned by DHGEQZ). |
|
|
119 * On exit, if SIDE = 'R' or 'B', VR contains: |
|
|
120 * if HOWMNY = 'A', the matrix X of right eigenvectors of (A,B); |
|
|
121 * if HOWMNY = 'B', the matrix Z*X; |
|
|
122 * if HOWMNY = 'S', the right eigenvectors of (A,B) specified by |
|
|
123 * SELECT, stored consecutively in the columns of |
|
|
124 * VR, in the same order as their eigenvalues. |
|
|
125 * If SIDE = 'L', VR is not referenced. |
|
|
126 * |
|
|
127 * A complex eigenvector corresponding to a complex eigenvalue |
|
|
128 * is stored in two consecutive columns, the first holding the |
|
|
129 * real part and the second the imaginary part. |
|
|
130 * |
|
|
131 * LDVR (input) INTEGER |
|
|
132 * The leading dimension of the array VR. |
|
|
133 * LDVR >= max(1,N) if SIDE = 'R' or 'B'; LDVR >= 1 otherwise. |
|
|
134 * |
|
|
135 * MM (input) INTEGER |
|
|
136 * The number of columns in the arrays VL and/or VR. MM >= M. |
|
|
137 * |
|
|
138 * M (output) INTEGER |
|
|
139 * The number of columns in the arrays VL and/or VR actually |
|
|
140 * used to store the eigenvectors. If HOWMNY = 'A' or 'B', M |
|
|
141 * is set to N. Each selected real eigenvector occupies one |
|
|
142 * column and each selected complex eigenvector occupies two |
|
|
143 * columns. |
|
|
144 * |
|
|
145 * WORK (workspace) DOUBLE PRECISION array, dimension (6*N) |
|
|
146 * |
|
|
147 * INFO (output) INTEGER |
|
|
148 * = 0: successful exit. |
|
|
149 * < 0: if INFO = -i, the i-th argument had an illegal value. |
|
|
150 * > 0: the 2-by-2 block (INFO:INFO+1) does not have a complex |
|
|
151 * eigenvalue. |
|
|
152 * |
|
|
153 * Further Details |
|
|
154 * =============== |
|
|
155 * |
|
|
156 * Allocation of workspace: |
|
|
157 * ---------- -- --------- |
|
|
158 * |
|
|
159 * WORK( j ) = 1-norm of j-th column of A, above the diagonal |
|
|
160 * WORK( N+j ) = 1-norm of j-th column of B, above the diagonal |
|
|
161 * WORK( 2*N+1:3*N ) = real part of eigenvector |
|
|
162 * WORK( 3*N+1:4*N ) = imaginary part of eigenvector |
|
|
163 * WORK( 4*N+1:5*N ) = real part of back-transformed eigenvector |
|
|
164 * WORK( 5*N+1:6*N ) = imaginary part of back-transformed eigenvector |
|
|
165 * |
|
|
166 * Rowwise vs. columnwise solution methods: |
|
|
167 * ------- -- ---------- -------- ------- |
|
|
168 * |
|
|
169 * Finding a generalized eigenvector consists basically of solving the |
|
|
170 * singular triangular system |
|
|
171 * |
|
|
172 * (A - w B) x = 0 (for right) or: (A - w B)**H y = 0 (for left) |
|
|
173 * |
|
|
174 * Consider finding the i-th right eigenvector (assume all eigenvalues |
|
|
175 * are real). The equation to be solved is: |
|
|
176 * n i |
|
|
177 * 0 = sum C(j,k) v(k) = sum C(j,k) v(k) for j = i,. . .,1 |
|
|
178 * k=j k=j |
|
|
179 * |
|
|
180 * where C = (A - w B) (The components v(i+1:n) are 0.) |
|
|
181 * |
|
|
182 * The "rowwise" method is: |
|
|
183 * |
|
|
184 * (1) v(i) := 1 |
|
|
185 * for j = i-1,. . .,1: |
|
|
186 * i |
|
|
187 * (2) compute s = - sum C(j,k) v(k) and |
|
|
188 * k=j+1 |
|
|
189 * |
|
|
190 * (3) v(j) := s / C(j,j) |
|
|
191 * |
|
|
192 * Step 2 is sometimes called the "dot product" step, since it is an |
|
|
193 * inner product between the j-th row and the portion of the eigenvector |
|
|
194 * that has been computed so far. |
|
|
195 * |
|
|
196 * The "columnwise" method consists basically in doing the sums |
|
|
197 * for all the rows in parallel. As each v(j) is computed, the |
|
|
198 * contribution of v(j) times the j-th column of C is added to the |
|
|
199 * partial sums. Since FORTRAN arrays are stored columnwise, this has |
|
|
200 * the advantage that at each step, the elements of C that are accessed |
|
|
201 * are adjacent to one another, whereas with the rowwise method, the |
|
|
202 * elements accessed at a step are spaced LDA (and LDB) words apart. |
|
|
203 * |
|
|
204 * When finding left eigenvectors, the matrix in question is the |
|
|
205 * transpose of the one in storage, so the rowwise method then |
|
|
206 * actually accesses columns of A and B at each step, and so is the |
|
|
207 * preferred method. |
|
|
208 * |
|
|
209 * ===================================================================== |
|
|
210 * |
|
|
211 * .. Parameters .. |
|
|
212 DOUBLE PRECISION ZERO, ONE, SAFETY |
|
|
213 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, |
|
|
214 $ SAFETY = 1.0D+2 ) |
|
|
215 * .. |
|
|
216 * .. Local Scalars .. |
|
|
217 LOGICAL COMPL, COMPR, IL2BY2, ILABAD, ILALL, ILBACK, |
|
|
218 $ ILBBAD, ILCOMP, ILCPLX, LSA, LSB |
|
|
219 INTEGER I, IBEG, IEIG, IEND, IHWMNY, IINFO, IM, ISIDE, |
|
|
220 $ J, JA, JC, JE, JR, JW, NA, NW |
|
|
221 DOUBLE PRECISION ACOEF, ACOEFA, ANORM, ASCALE, BCOEFA, BCOEFI, |
|
|
222 $ BCOEFR, BIG, BIGNUM, BNORM, BSCALE, CIM2A, |
|
|
223 $ CIM2B, CIMAGA, CIMAGB, CRE2A, CRE2B, CREALA, |
|
|
224 $ CREALB, DMIN, SAFMIN, SALFAR, SBETA, SCALE, |
|
|
225 $ SMALL, TEMP, TEMP2, TEMP2I, TEMP2R, ULP, XMAX, |
|
|
226 $ XSCALE |
|
|
227 * .. |
|
|
228 * .. Local Arrays .. |
|
|
229 DOUBLE PRECISION BDIAG( 2 ), SUM( 2, 2 ), SUMA( 2, 2 ), |
|
|
230 $ SUMB( 2, 2 ) |
|
|
231 * .. |
|
|
232 * .. External Functions .. |
|
|
233 LOGICAL LSAME |
|
|
234 DOUBLE PRECISION DLAMCH |
|
|
235 EXTERNAL LSAME, DLAMCH |
|
|
236 * .. |
|
|
237 * .. External Subroutines .. |
|
3333
|
238 EXTERNAL DGEMV, DLACPY, DLAG2, DLALN2, XERBLA |
|
3274
|
239 * .. |
|
|
240 * .. Intrinsic Functions .. |
|
|
241 INTRINSIC ABS, MAX, MIN |
|
|
242 * .. |
|
|
243 * .. Executable Statements .. |
|
|
244 * |
|
|
245 * Decode and Test the input parameters |
|
|
246 * |
|
|
247 IF( LSAME( HOWMNY, 'A' ) ) THEN |
|
|
248 IHWMNY = 1 |
|
|
249 ILALL = .TRUE. |
|
|
250 ILBACK = .FALSE. |
|
|
251 ELSE IF( LSAME( HOWMNY, 'S' ) ) THEN |
|
|
252 IHWMNY = 2 |
|
|
253 ILALL = .FALSE. |
|
|
254 ILBACK = .FALSE. |
|
|
255 ELSE IF( LSAME( HOWMNY, 'B' ) .OR. LSAME( HOWMNY, 'T' ) ) THEN |
|
|
256 IHWMNY = 3 |
|
|
257 ILALL = .TRUE. |
|
|
258 ILBACK = .TRUE. |
|
|
259 ELSE |
|
|
260 IHWMNY = -1 |
|
|
261 ILALL = .TRUE. |
|
|
262 END IF |
|
|
263 * |
|
|
264 IF( LSAME( SIDE, 'R' ) ) THEN |
|
|
265 ISIDE = 1 |
|
|
266 COMPL = .FALSE. |
|
|
267 COMPR = .TRUE. |
|
|
268 ELSE IF( LSAME( SIDE, 'L' ) ) THEN |
|
|
269 ISIDE = 2 |
|
|
270 COMPL = .TRUE. |
|
|
271 COMPR = .FALSE. |
|
|
272 ELSE IF( LSAME( SIDE, 'B' ) ) THEN |
|
|
273 ISIDE = 3 |
|
|
274 COMPL = .TRUE. |
|
|
275 COMPR = .TRUE. |
|
|
276 ELSE |
|
|
277 ISIDE = -1 |
|
|
278 END IF |
|
|
279 * |
|
3333
|
280 INFO = 0 |
|
|
281 IF( ISIDE.LT.0 ) THEN |
|
|
282 INFO = -1 |
|
|
283 ELSE IF( IHWMNY.LT.0 ) THEN |
|
|
284 INFO = -2 |
|
|
285 ELSE IF( N.LT.0 ) THEN |
|
|
286 INFO = -4 |
|
|
287 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN |
|
|
288 INFO = -6 |
|
|
289 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN |
|
|
290 INFO = -8 |
|
|
291 END IF |
|
|
292 IF( INFO.NE.0 ) THEN |
|
|
293 CALL XERBLA( 'DTGEVC', -INFO ) |
|
|
294 RETURN |
|
|
295 END IF |
|
|
296 * |
|
3274
|
297 * Count the number of eigenvectors to be computed |
|
|
298 * |
|
|
299 IF( .NOT.ILALL ) THEN |
|
|
300 IM = 0 |
|
|
301 ILCPLX = .FALSE. |
|
|
302 DO 10 J = 1, N |
|
|
303 IF( ILCPLX ) THEN |
|
|
304 ILCPLX = .FALSE. |
|
|
305 GO TO 10 |
|
|
306 END IF |
|
|
307 IF( J.LT.N ) THEN |
|
|
308 IF( A( J+1, J ).NE.ZERO ) |
|
|
309 $ ILCPLX = .TRUE. |
|
|
310 END IF |
|
|
311 IF( ILCPLX ) THEN |
|
|
312 IF( SELECT( J ) .OR. SELECT( J+1 ) ) |
|
|
313 $ IM = IM + 2 |
|
|
314 ELSE |
|
|
315 IF( SELECT( J ) ) |
|
|
316 $ IM = IM + 1 |
|
|
317 END IF |
|
|
318 10 CONTINUE |
|
|
319 ELSE |
|
|
320 IM = N |
|
|
321 END IF |
|
|
322 * |
|
|
323 * Check 2-by-2 diagonal blocks of A, B |
|
|
324 * |
|
|
325 ILABAD = .FALSE. |
|
|
326 ILBBAD = .FALSE. |
|
|
327 DO 20 J = 1, N - 1 |
|
|
328 IF( A( J+1, J ).NE.ZERO ) THEN |
|
|
329 IF( B( J, J ).EQ.ZERO .OR. B( J+1, J+1 ).EQ.ZERO .OR. |
|
|
330 $ B( J, J+1 ).NE.ZERO )ILBBAD = .TRUE. |
|
|
331 IF( J.LT.N-1 ) THEN |
|
|
332 IF( A( J+2, J+1 ).NE.ZERO ) |
|
|
333 $ ILABAD = .TRUE. |
|
|
334 END IF |
|
|
335 END IF |
|
|
336 20 CONTINUE |
|
|
337 * |
|
3333
|
338 IF( ILABAD ) THEN |
|
3274
|
339 INFO = -5 |
|
|
340 ELSE IF( ILBBAD ) THEN |
|
|
341 INFO = -7 |
|
|
342 ELSE IF( COMPL .AND. LDVL.LT.N .OR. LDVL.LT.1 ) THEN |
|
|
343 INFO = -10 |
|
|
344 ELSE IF( COMPR .AND. LDVR.LT.N .OR. LDVR.LT.1 ) THEN |
|
|
345 INFO = -12 |
|
|
346 ELSE IF( MM.LT.IM ) THEN |
|
|
347 INFO = -13 |
|
|
348 END IF |
|
|
349 IF( INFO.NE.0 ) THEN |
|
|
350 CALL XERBLA( 'DTGEVC', -INFO ) |
|
|
351 RETURN |
|
|
352 END IF |
|
|
353 * |
|
|
354 * Quick return if possible |
|
|
355 * |
|
|
356 M = IM |
|
|
357 IF( N.EQ.0 ) |
|
|
358 $ RETURN |
|
|
359 * |
|
|
360 * Machine Constants |
|
|
361 * |
|
|
362 SAFMIN = DLAMCH( 'Safe minimum' ) |
|
|
363 BIG = ONE / SAFMIN |
|
|
364 CALL DLABAD( SAFMIN, BIG ) |
|
|
365 ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' ) |
|
|
366 SMALL = SAFMIN*N / ULP |
|
|
367 BIG = ONE / SMALL |
|
|
368 BIGNUM = ONE / ( SAFMIN*N ) |
|
|
369 * |
|
|
370 * Compute the 1-norm of each column of the strictly upper triangular |
|
|
371 * part (i.e., excluding all elements belonging to the diagonal |
|
|
372 * blocks) of A and B to check for possible overflow in the |
|
|
373 * triangular solver. |
|
|
374 * |
|
|
375 ANORM = ABS( A( 1, 1 ) ) |
|
|
376 IF( N.GT.1 ) |
|
|
377 $ ANORM = ANORM + ABS( A( 2, 1 ) ) |
|
|
378 BNORM = ABS( B( 1, 1 ) ) |
|
|
379 WORK( 1 ) = ZERO |
|
|
380 WORK( N+1 ) = ZERO |
|
|
381 * |
|
|
382 DO 50 J = 2, N |
|
|
383 TEMP = ZERO |
|
|
384 TEMP2 = ZERO |
|
|
385 IF( A( J, J-1 ).EQ.ZERO ) THEN |
|
|
386 IEND = J - 1 |
|
|
387 ELSE |
|
|
388 IEND = J - 2 |
|
|
389 END IF |
|
|
390 DO 30 I = 1, IEND |
|
|
391 TEMP = TEMP + ABS( A( I, J ) ) |
|
|
392 TEMP2 = TEMP2 + ABS( B( I, J ) ) |
|
|
393 30 CONTINUE |
|
|
394 WORK( J ) = TEMP |
|
|
395 WORK( N+J ) = TEMP2 |
|
|
396 DO 40 I = IEND + 1, MIN( J+1, N ) |
|
|
397 TEMP = TEMP + ABS( A( I, J ) ) |
|
|
398 TEMP2 = TEMP2 + ABS( B( I, J ) ) |
|
|
399 40 CONTINUE |
|
|
400 ANORM = MAX( ANORM, TEMP ) |
|
|
401 BNORM = MAX( BNORM, TEMP2 ) |
|
|
402 50 CONTINUE |
|
|
403 * |
|
|
404 ASCALE = ONE / MAX( ANORM, SAFMIN ) |
|
|
405 BSCALE = ONE / MAX( BNORM, SAFMIN ) |
|
|
406 * |
|
|
407 * Left eigenvectors |
|
|
408 * |
|
|
409 IF( COMPL ) THEN |
|
|
410 IEIG = 0 |
|
|
411 * |
|
|
412 * Main loop over eigenvalues |
|
|
413 * |
|
|
414 ILCPLX = .FALSE. |
|
|
415 DO 220 JE = 1, N |
|
|
416 * |
|
|
417 * Skip this iteration if (a) HOWMNY='S' and SELECT=.FALSE., or |
|
|
418 * (b) this would be the second of a complex pair. |
|
|
419 * Check for complex eigenvalue, so as to be sure of which |
|
|
420 * entry(-ies) of SELECT to look at. |
|
|
421 * |
|
|
422 IF( ILCPLX ) THEN |
|
|
423 ILCPLX = .FALSE. |
|
|
424 GO TO 220 |
|
|
425 END IF |
|
|
426 NW = 1 |
|
|
427 IF( JE.LT.N ) THEN |
|
|
428 IF( A( JE+1, JE ).NE.ZERO ) THEN |
|
|
429 ILCPLX = .TRUE. |
|
|
430 NW = 2 |
|
|
431 END IF |
|
|
432 END IF |
|
|
433 IF( ILALL ) THEN |
|
|
434 ILCOMP = .TRUE. |
|
|
435 ELSE IF( ILCPLX ) THEN |
|
|
436 ILCOMP = SELECT( JE ) .OR. SELECT( JE+1 ) |
|
|
437 ELSE |
|
|
438 ILCOMP = SELECT( JE ) |
|
|
439 END IF |
|
|
440 IF( .NOT.ILCOMP ) |
|
|
441 $ GO TO 220 |
|
|
442 * |
|
|
443 * Decide if (a) singular pencil, (b) real eigenvalue, or |
|
|
444 * (c) complex eigenvalue. |
|
|
445 * |
|
|
446 IF( .NOT.ILCPLX ) THEN |
|
|
447 IF( ABS( A( JE, JE ) ).LE.SAFMIN .AND. |
|
|
448 $ ABS( B( JE, JE ) ).LE.SAFMIN ) THEN |
|
|
449 * |
|
|
450 * Singular matrix pencil -- return unit eigenvector |
|
|
451 * |
|
|
452 IEIG = IEIG + 1 |
|
|
453 DO 60 JR = 1, N |
|
|
454 VL( JR, IEIG ) = ZERO |
|
|
455 60 CONTINUE |
|
|
456 VL( IEIG, IEIG ) = ONE |
|
|
457 GO TO 220 |
|
|
458 END IF |
|
|
459 END IF |
|
|
460 * |
|
|
461 * Clear vector |
|
|
462 * |
|
|
463 DO 70 JR = 1, NW*N |
|
|
464 WORK( 2*N+JR ) = ZERO |
|
|
465 70 CONTINUE |
|
|
466 * T |
|
|
467 * Compute coefficients in ( a A - b B ) y = 0 |
|
|
468 * a is ACOEF |
|
|
469 * b is BCOEFR + i*BCOEFI |
|
|
470 * |
|
|
471 IF( .NOT.ILCPLX ) THEN |
|
|
472 * |
|
|
473 * Real eigenvalue |
|
|
474 * |
|
|
475 TEMP = ONE / MAX( ABS( A( JE, JE ) )*ASCALE, |
|
|
476 $ ABS( B( JE, JE ) )*BSCALE, SAFMIN ) |
|
|
477 SALFAR = ( TEMP*A( JE, JE ) )*ASCALE |
|
|
478 SBETA = ( TEMP*B( JE, JE ) )*BSCALE |
|
|
479 ACOEF = SBETA*ASCALE |
|
|
480 BCOEFR = SALFAR*BSCALE |
|
|
481 BCOEFI = ZERO |
|
|
482 * |
|
|
483 * Scale to avoid underflow |
|
|
484 * |
|
|
485 SCALE = ONE |
|
|
486 LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEF ).LT.SMALL |
|
|
487 LSB = ABS( SALFAR ).GE.SAFMIN .AND. ABS( BCOEFR ).LT. |
|
|
488 $ SMALL |
|
|
489 IF( LSA ) |
|
|
490 $ SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG ) |
|
|
491 IF( LSB ) |
|
|
492 $ SCALE = MAX( SCALE, ( SMALL / ABS( SALFAR ) )* |
|
|
493 $ MIN( BNORM, BIG ) ) |
|
|
494 IF( LSA .OR. LSB ) THEN |
|
|
495 SCALE = MIN( SCALE, ONE / |
|
|
496 $ ( SAFMIN*MAX( ONE, ABS( ACOEF ), |
|
|
497 $ ABS( BCOEFR ) ) ) ) |
|
|
498 IF( LSA ) THEN |
|
|
499 ACOEF = ASCALE*( SCALE*SBETA ) |
|
|
500 ELSE |
|
|
501 ACOEF = SCALE*ACOEF |
|
|
502 END IF |
|
|
503 IF( LSB ) THEN |
|
|
504 BCOEFR = BSCALE*( SCALE*SALFAR ) |
|
|
505 ELSE |
|
|
506 BCOEFR = SCALE*BCOEFR |
|
|
507 END IF |
|
|
508 END IF |
|
|
509 ACOEFA = ABS( ACOEF ) |
|
|
510 BCOEFA = ABS( BCOEFR ) |
|
|
511 * |
|
|
512 * First component is 1 |
|
|
513 * |
|
|
514 WORK( 2*N+JE ) = ONE |
|
|
515 XMAX = ONE |
|
|
516 ELSE |
|
|
517 * |
|
|
518 * Complex eigenvalue |
|
|
519 * |
|
|
520 CALL DLAG2( A( JE, JE ), LDA, B( JE, JE ), LDB, |
|
|
521 $ SAFMIN*SAFETY, ACOEF, TEMP, BCOEFR, TEMP2, |
|
|
522 $ BCOEFI ) |
|
|
523 BCOEFI = -BCOEFI |
|
|
524 IF( BCOEFI.EQ.ZERO ) THEN |
|
|
525 INFO = JE |
|
|
526 RETURN |
|
|
527 END IF |
|
|
528 * |
|
|
529 * Scale to avoid over/underflow |
|
|
530 * |
|
|
531 ACOEFA = ABS( ACOEF ) |
|
|
532 BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI ) |
|
|
533 SCALE = ONE |
|
|
534 IF( ACOEFA*ULP.LT.SAFMIN .AND. ACOEFA.GE.SAFMIN ) |
|
|
535 $ SCALE = ( SAFMIN / ULP ) / ACOEFA |
|
|
536 IF( BCOEFA*ULP.LT.SAFMIN .AND. BCOEFA.GE.SAFMIN ) |
|
|
537 $ SCALE = MAX( SCALE, ( SAFMIN / ULP ) / BCOEFA ) |
|
|
538 IF( SAFMIN*ACOEFA.GT.ASCALE ) |
|
|
539 $ SCALE = ASCALE / ( SAFMIN*ACOEFA ) |
|
|
540 IF( SAFMIN*BCOEFA.GT.BSCALE ) |
|
|
541 $ SCALE = MIN( SCALE, BSCALE / ( SAFMIN*BCOEFA ) ) |
|
|
542 IF( SCALE.NE.ONE ) THEN |
|
|
543 ACOEF = SCALE*ACOEF |
|
|
544 ACOEFA = ABS( ACOEF ) |
|
|
545 BCOEFR = SCALE*BCOEFR |
|
|
546 BCOEFI = SCALE*BCOEFI |
|
|
547 BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI ) |
|
|
548 END IF |
|
|
549 * |
|
|
550 * Compute first two components of eigenvector |
|
|
551 * |
|
|
552 TEMP = ACOEF*A( JE+1, JE ) |
|
|
553 TEMP2R = ACOEF*A( JE, JE ) - BCOEFR*B( JE, JE ) |
|
|
554 TEMP2I = -BCOEFI*B( JE, JE ) |
|
|
555 IF( ABS( TEMP ).GT.ABS( TEMP2R )+ABS( TEMP2I ) ) THEN |
|
|
556 WORK( 2*N+JE ) = ONE |
|
|
557 WORK( 3*N+JE ) = ZERO |
|
|
558 WORK( 2*N+JE+1 ) = -TEMP2R / TEMP |
|
|
559 WORK( 3*N+JE+1 ) = -TEMP2I / TEMP |
|
|
560 ELSE |
|
|
561 WORK( 2*N+JE+1 ) = ONE |
|
|
562 WORK( 3*N+JE+1 ) = ZERO |
|
|
563 TEMP = ACOEF*A( JE, JE+1 ) |
|
|
564 WORK( 2*N+JE ) = ( BCOEFR*B( JE+1, JE+1 )-ACOEF* |
|
|
565 $ A( JE+1, JE+1 ) ) / TEMP |
|
|
566 WORK( 3*N+JE ) = BCOEFI*B( JE+1, JE+1 ) / TEMP |
|
|
567 END IF |
|
|
568 XMAX = MAX( ABS( WORK( 2*N+JE ) )+ABS( WORK( 3*N+JE ) ), |
|
|
569 $ ABS( WORK( 2*N+JE+1 ) )+ABS( WORK( 3*N+JE+1 ) ) ) |
|
|
570 END IF |
|
|
571 * |
|
|
572 DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN ) |
|
|
573 * |
|
|
574 * T |
|
|
575 * Triangular solve of (a A - b B) y = 0 |
|
|
576 * |
|
|
577 * T |
|
|
578 * (rowwise in (a A - b B) , or columnwise in (a A - b B) ) |
|
|
579 * |
|
|
580 IL2BY2 = .FALSE. |
|
|
581 * |
|
|
582 DO 160 J = JE + NW, N |
|
|
583 IF( IL2BY2 ) THEN |
|
|
584 IL2BY2 = .FALSE. |
|
|
585 GO TO 160 |
|
|
586 END IF |
|
|
587 * |
|
|
588 NA = 1 |
|
|
589 BDIAG( 1 ) = B( J, J ) |
|
|
590 IF( J.LT.N ) THEN |
|
|
591 IF( A( J+1, J ).NE.ZERO ) THEN |
|
|
592 IL2BY2 = .TRUE. |
|
|
593 BDIAG( 2 ) = B( J+1, J+1 ) |
|
|
594 NA = 2 |
|
|
595 END IF |
|
|
596 END IF |
|
|
597 * |
|
|
598 * Check whether scaling is necessary for dot products |
|
|
599 * |
|
|
600 XSCALE = ONE / MAX( ONE, XMAX ) |
|
|
601 TEMP = MAX( WORK( J ), WORK( N+J ), |
|
|
602 $ ACOEFA*WORK( J )+BCOEFA*WORK( N+J ) ) |
|
|
603 IF( IL2BY2 ) |
|
|
604 $ TEMP = MAX( TEMP, WORK( J+1 ), WORK( N+J+1 ), |
|
|
605 $ ACOEFA*WORK( J+1 )+BCOEFA*WORK( N+J+1 ) ) |
|
|
606 IF( TEMP.GT.BIGNUM*XSCALE ) THEN |
|
|
607 DO 90 JW = 0, NW - 1 |
|
|
608 DO 80 JR = JE, J - 1 |
|
|
609 WORK( ( JW+2 )*N+JR ) = XSCALE* |
|
|
610 $ WORK( ( JW+2 )*N+JR ) |
|
|
611 80 CONTINUE |
|
|
612 90 CONTINUE |
|
|
613 XMAX = XMAX*XSCALE |
|
|
614 END IF |
|
|
615 * |
|
|
616 * Compute dot products |
|
|
617 * |
|
|
618 * j-1 |
|
|
619 * SUM = sum conjg( a*A(k,j) - b*B(k,j) )*x(k) |
|
|
620 * k=je |
|
|
621 * |
|
|
622 * To reduce the op count, this is done as |
|
|
623 * |
|
|
624 * _ j-1 _ j-1 |
|
|
625 * a*conjg( sum A(k,j)*x(k) ) - b*conjg( sum B(k,j)*x(k) ) |
|
|
626 * k=je k=je |
|
|
627 * |
|
|
628 * which may cause underflow problems if A or B are close |
|
|
629 * to underflow. (E.g., less than SMALL.) |
|
|
630 * |
|
|
631 * |
|
|
632 * A series of compiler directives to defeat vectorization |
|
|
633 * for the next loop |
|
|
634 * |
|
|
635 *$PL$ CMCHAR=' ' |
|
|
636 CDIR$ NEXTSCALAR |
|
|
637 C$DIR SCALAR |
|
|
638 CDIR$ NEXT SCALAR |
|
|
639 CVD$L NOVECTOR |
|
|
640 CDEC$ NOVECTOR |
|
|
641 CVD$ NOVECTOR |
|
|
642 *VDIR NOVECTOR |
|
|
643 *VOCL LOOP,SCALAR |
|
|
644 CIBM PREFER SCALAR |
|
|
645 *$PL$ CMCHAR='*' |
|
|
646 * |
|
|
647 DO 120 JW = 1, NW |
|
|
648 * |
|
|
649 *$PL$ CMCHAR=' ' |
|
|
650 CDIR$ NEXTSCALAR |
|
|
651 C$DIR SCALAR |
|
|
652 CDIR$ NEXT SCALAR |
|
|
653 CVD$L NOVECTOR |
|
|
654 CDEC$ NOVECTOR |
|
|
655 CVD$ NOVECTOR |
|
|
656 *VDIR NOVECTOR |
|
|
657 *VOCL LOOP,SCALAR |
|
|
658 CIBM PREFER SCALAR |
|
|
659 *$PL$ CMCHAR='*' |
|
|
660 * |
|
|
661 DO 110 JA = 1, NA |
|
|
662 SUMA( JA, JW ) = ZERO |
|
|
663 SUMB( JA, JW ) = ZERO |
|
|
664 * |
|
|
665 DO 100 JR = JE, J - 1 |
|
|
666 SUMA( JA, JW ) = SUMA( JA, JW ) + |
|
|
667 $ A( JR, J+JA-1 )* |
|
|
668 $ WORK( ( JW+1 )*N+JR ) |
|
|
669 SUMB( JA, JW ) = SUMB( JA, JW ) + |
|
|
670 $ B( JR, J+JA-1 )* |
|
|
671 $ WORK( ( JW+1 )*N+JR ) |
|
|
672 100 CONTINUE |
|
|
673 110 CONTINUE |
|
|
674 120 CONTINUE |
|
|
675 * |
|
|
676 *$PL$ CMCHAR=' ' |
|
|
677 CDIR$ NEXTSCALAR |
|
|
678 C$DIR SCALAR |
|
|
679 CDIR$ NEXT SCALAR |
|
|
680 CVD$L NOVECTOR |
|
|
681 CDEC$ NOVECTOR |
|
|
682 CVD$ NOVECTOR |
|
|
683 *VDIR NOVECTOR |
|
|
684 *VOCL LOOP,SCALAR |
|
|
685 CIBM PREFER SCALAR |
|
|
686 *$PL$ CMCHAR='*' |
|
|
687 * |
|
|
688 DO 130 JA = 1, NA |
|
|
689 IF( ILCPLX ) THEN |
|
|
690 SUM( JA, 1 ) = -ACOEF*SUMA( JA, 1 ) + |
|
|
691 $ BCOEFR*SUMB( JA, 1 ) - |
|
|
692 $ BCOEFI*SUMB( JA, 2 ) |
|
|
693 SUM( JA, 2 ) = -ACOEF*SUMA( JA, 2 ) + |
|
|
694 $ BCOEFR*SUMB( JA, 2 ) + |
|
|
695 $ BCOEFI*SUMB( JA, 1 ) |
|
|
696 ELSE |
|
|
697 SUM( JA, 1 ) = -ACOEF*SUMA( JA, 1 ) + |
|
|
698 $ BCOEFR*SUMB( JA, 1 ) |
|
|
699 END IF |
|
|
700 130 CONTINUE |
|
|
701 * |
|
|
702 * T |
|
|
703 * Solve ( a A - b B ) y = SUM(,) |
|
|
704 * with scaling and perturbation of the denominator |
|
|
705 * |
|
|
706 CALL DLALN2( .TRUE., NA, NW, DMIN, ACOEF, A( J, J ), LDA, |
|
|
707 $ BDIAG( 1 ), BDIAG( 2 ), SUM, 2, BCOEFR, |
|
|
708 $ BCOEFI, WORK( 2*N+J ), N, SCALE, TEMP, |
|
|
709 $ IINFO ) |
|
|
710 IF( SCALE.LT.ONE ) THEN |
|
|
711 DO 150 JW = 0, NW - 1 |
|
|
712 DO 140 JR = JE, J - 1 |
|
|
713 WORK( ( JW+2 )*N+JR ) = SCALE* |
|
|
714 $ WORK( ( JW+2 )*N+JR ) |
|
|
715 140 CONTINUE |
|
|
716 150 CONTINUE |
|
|
717 XMAX = SCALE*XMAX |
|
|
718 END IF |
|
|
719 XMAX = MAX( XMAX, TEMP ) |
|
|
720 160 CONTINUE |
|
|
721 * |
|
|
722 * Copy eigenvector to VL, back transforming if |
|
|
723 * HOWMNY='B'. |
|
|
724 * |
|
|
725 IEIG = IEIG + 1 |
|
|
726 IF( ILBACK ) THEN |
|
|
727 DO 170 JW = 0, NW - 1 |
|
|
728 CALL DGEMV( 'N', N, N+1-JE, ONE, VL( 1, JE ), LDVL, |
|
|
729 $ WORK( ( JW+2 )*N+JE ), 1, ZERO, |
|
|
730 $ WORK( ( JW+4 )*N+1 ), 1 ) |
|
|
731 170 CONTINUE |
|
|
732 CALL DLACPY( ' ', N, NW, WORK( 4*N+1 ), N, VL( 1, JE ), |
|
|
733 $ LDVL ) |
|
|
734 IBEG = 1 |
|
|
735 ELSE |
|
|
736 CALL DLACPY( ' ', N, NW, WORK( 2*N+1 ), N, VL( 1, IEIG ), |
|
|
737 $ LDVL ) |
|
|
738 IBEG = JE |
|
|
739 END IF |
|
|
740 * |
|
|
741 * Scale eigenvector |
|
|
742 * |
|
|
743 XMAX = ZERO |
|
|
744 IF( ILCPLX ) THEN |
|
|
745 DO 180 J = IBEG, N |
|
|
746 XMAX = MAX( XMAX, ABS( VL( J, IEIG ) )+ |
|
|
747 $ ABS( VL( J, IEIG+1 ) ) ) |
|
|
748 180 CONTINUE |
|
|
749 ELSE |
|
|
750 DO 190 J = IBEG, N |
|
|
751 XMAX = MAX( XMAX, ABS( VL( J, IEIG ) ) ) |
|
|
752 190 CONTINUE |
|
|
753 END IF |
|
|
754 * |
|
|
755 IF( XMAX.GT.SAFMIN ) THEN |
|
|
756 XSCALE = ONE / XMAX |
|
|
757 * |
|
|
758 DO 210 JW = 0, NW - 1 |
|
|
759 DO 200 JR = IBEG, N |
|
|
760 VL( JR, IEIG+JW ) = XSCALE*VL( JR, IEIG+JW ) |
|
|
761 200 CONTINUE |
|
|
762 210 CONTINUE |
|
|
763 END IF |
|
|
764 IEIG = IEIG + NW - 1 |
|
|
765 * |
|
|
766 220 CONTINUE |
|
|
767 END IF |
|
|
768 * |
|
|
769 * Right eigenvectors |
|
|
770 * |
|
|
771 IF( COMPR ) THEN |
|
|
772 IEIG = IM + 1 |
|
|
773 * |
|
|
774 * Main loop over eigenvalues |
|
|
775 * |
|
|
776 ILCPLX = .FALSE. |
|
|
777 DO 500 JE = N, 1, -1 |
|
|
778 * |
|
|
779 * Skip this iteration if (a) HOWMNY='S' and SELECT=.FALSE., or |
|
|
780 * (b) this would be the second of a complex pair. |
|
|
781 * Check for complex eigenvalue, so as to be sure of which |
|
|
782 * entry(-ies) of SELECT to look at -- if complex, SELECT(JE) |
|
|
783 * or SELECT(JE-1). |
|
|
784 * If this is a complex pair, the 2-by-2 diagonal block |
|
|
785 * corresponding to the eigenvalue is in rows/columns JE-1:JE |
|
|
786 * |
|
|
787 IF( ILCPLX ) THEN |
|
|
788 ILCPLX = .FALSE. |
|
|
789 GO TO 500 |
|
|
790 END IF |
|
|
791 NW = 1 |
|
|
792 IF( JE.GT.1 ) THEN |
|
|
793 IF( A( JE, JE-1 ).NE.ZERO ) THEN |
|
|
794 ILCPLX = .TRUE. |
|
|
795 NW = 2 |
|
|
796 END IF |
|
|
797 END IF |
|
|
798 IF( ILALL ) THEN |
|
|
799 ILCOMP = .TRUE. |
|
|
800 ELSE IF( ILCPLX ) THEN |
|
|
801 ILCOMP = SELECT( JE ) .OR. SELECT( JE-1 ) |
|
|
802 ELSE |
|
|
803 ILCOMP = SELECT( JE ) |
|
|
804 END IF |
|
|
805 IF( .NOT.ILCOMP ) |
|
|
806 $ GO TO 500 |
|
|
807 * |
|
|
808 * Decide if (a) singular pencil, (b) real eigenvalue, or |
|
|
809 * (c) complex eigenvalue. |
|
|
810 * |
|
|
811 IF( .NOT.ILCPLX ) THEN |
|
|
812 IF( ABS( A( JE, JE ) ).LE.SAFMIN .AND. |
|
|
813 $ ABS( B( JE, JE ) ).LE.SAFMIN ) THEN |
|
|
814 * |
|
|
815 * Singular matrix pencil -- unit eigenvector |
|
|
816 * |
|
|
817 IEIG = IEIG - 1 |
|
|
818 DO 230 JR = 1, N |
|
|
819 VR( JR, IEIG ) = ZERO |
|
|
820 230 CONTINUE |
|
|
821 VR( IEIG, IEIG ) = ONE |
|
|
822 GO TO 500 |
|
|
823 END IF |
|
|
824 END IF |
|
|
825 * |
|
|
826 * Clear vector |
|
|
827 * |
|
|
828 DO 250 JW = 0, NW - 1 |
|
|
829 DO 240 JR = 1, N |
|
|
830 WORK( ( JW+2 )*N+JR ) = ZERO |
|
|
831 240 CONTINUE |
|
|
832 250 CONTINUE |
|
|
833 * |
|
|
834 * Compute coefficients in ( a A - b B ) x = 0 |
|
|
835 * a is ACOEF |
|
|
836 * b is BCOEFR + i*BCOEFI |
|
|
837 * |
|
|
838 IF( .NOT.ILCPLX ) THEN |
|
|
839 * |
|
|
840 * Real eigenvalue |
|
|
841 * |
|
|
842 TEMP = ONE / MAX( ABS( A( JE, JE ) )*ASCALE, |
|
|
843 $ ABS( B( JE, JE ) )*BSCALE, SAFMIN ) |
|
|
844 SALFAR = ( TEMP*A( JE, JE ) )*ASCALE |
|
|
845 SBETA = ( TEMP*B( JE, JE ) )*BSCALE |
|
|
846 ACOEF = SBETA*ASCALE |
|
|
847 BCOEFR = SALFAR*BSCALE |
|
|
848 BCOEFI = ZERO |
|
|
849 * |
|
|
850 * Scale to avoid underflow |
|
|
851 * |
|
|
852 SCALE = ONE |
|
|
853 LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEF ).LT.SMALL |
|
|
854 LSB = ABS( SALFAR ).GE.SAFMIN .AND. ABS( BCOEFR ).LT. |
|
|
855 $ SMALL |
|
|
856 IF( LSA ) |
|
|
857 $ SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG ) |
|
|
858 IF( LSB ) |
|
|
859 $ SCALE = MAX( SCALE, ( SMALL / ABS( SALFAR ) )* |
|
|
860 $ MIN( BNORM, BIG ) ) |
|
|
861 IF( LSA .OR. LSB ) THEN |
|
|
862 SCALE = MIN( SCALE, ONE / |
|
|
863 $ ( SAFMIN*MAX( ONE, ABS( ACOEF ), |
|
|
864 $ ABS( BCOEFR ) ) ) ) |
|
|
865 IF( LSA ) THEN |
|
|
866 ACOEF = ASCALE*( SCALE*SBETA ) |
|
|
867 ELSE |
|
|
868 ACOEF = SCALE*ACOEF |
|
|
869 END IF |
|
|
870 IF( LSB ) THEN |
|
|
871 BCOEFR = BSCALE*( SCALE*SALFAR ) |
|
|
872 ELSE |
|
|
873 BCOEFR = SCALE*BCOEFR |
|
|
874 END IF |
|
|
875 END IF |
|
|
876 ACOEFA = ABS( ACOEF ) |
|
|
877 BCOEFA = ABS( BCOEFR ) |
|
|
878 * |
|
|
879 * First component is 1 |
|
|
880 * |
|
|
881 WORK( 2*N+JE ) = ONE |
|
|
882 XMAX = ONE |
|
|
883 * |
|
|
884 * Compute contribution from column JE of A and B to sum |
|
|
885 * (See "Further Details", above.) |
|
|
886 * |
|
|
887 DO 260 JR = 1, JE - 1 |
|
|
888 WORK( 2*N+JR ) = BCOEFR*B( JR, JE ) - |
|
|
889 $ ACOEF*A( JR, JE ) |
|
|
890 260 CONTINUE |
|
|
891 ELSE |
|
|
892 * |
|
|
893 * Complex eigenvalue |
|
|
894 * |
|
|
895 CALL DLAG2( A( JE-1, JE-1 ), LDA, B( JE-1, JE-1 ), LDB, |
|
|
896 $ SAFMIN*SAFETY, ACOEF, TEMP, BCOEFR, TEMP2, |
|
|
897 $ BCOEFI ) |
|
|
898 IF( BCOEFI.EQ.ZERO ) THEN |
|
|
899 INFO = JE - 1 |
|
|
900 RETURN |
|
|
901 END IF |
|
|
902 * |
|
|
903 * Scale to avoid over/underflow |
|
|
904 * |
|
|
905 ACOEFA = ABS( ACOEF ) |
|
|
906 BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI ) |
|
|
907 SCALE = ONE |
|
|
908 IF( ACOEFA*ULP.LT.SAFMIN .AND. ACOEFA.GE.SAFMIN ) |
|
|
909 $ SCALE = ( SAFMIN / ULP ) / ACOEFA |
|
|
910 IF( BCOEFA*ULP.LT.SAFMIN .AND. BCOEFA.GE.SAFMIN ) |
|
|
911 $ SCALE = MAX( SCALE, ( SAFMIN / ULP ) / BCOEFA ) |
|
|
912 IF( SAFMIN*ACOEFA.GT.ASCALE ) |
|
|
913 $ SCALE = ASCALE / ( SAFMIN*ACOEFA ) |
|
|
914 IF( SAFMIN*BCOEFA.GT.BSCALE ) |
|
|
915 $ SCALE = MIN( SCALE, BSCALE / ( SAFMIN*BCOEFA ) ) |
|
|
916 IF( SCALE.NE.ONE ) THEN |
|
|
917 ACOEF = SCALE*ACOEF |
|
|
918 ACOEFA = ABS( ACOEF ) |
|
|
919 BCOEFR = SCALE*BCOEFR |
|
|
920 BCOEFI = SCALE*BCOEFI |
|
|
921 BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI ) |
|
|
922 END IF |
|
|
923 * |
|
|
924 * Compute first two components of eigenvector |
|
|
925 * and contribution to sums |
|
|
926 * |
|
|
927 TEMP = ACOEF*A( JE, JE-1 ) |
|
|
928 TEMP2R = ACOEF*A( JE, JE ) - BCOEFR*B( JE, JE ) |
|
|
929 TEMP2I = -BCOEFI*B( JE, JE ) |
|
|
930 IF( ABS( TEMP ).GE.ABS( TEMP2R )+ABS( TEMP2I ) ) THEN |
|
|
931 WORK( 2*N+JE ) = ONE |
|
|
932 WORK( 3*N+JE ) = ZERO |
|
|
933 WORK( 2*N+JE-1 ) = -TEMP2R / TEMP |
|
|
934 WORK( 3*N+JE-1 ) = -TEMP2I / TEMP |
|
|
935 ELSE |
|
|
936 WORK( 2*N+JE-1 ) = ONE |
|
|
937 WORK( 3*N+JE-1 ) = ZERO |
|
|
938 TEMP = ACOEF*A( JE-1, JE ) |
|
|
939 WORK( 2*N+JE ) = ( BCOEFR*B( JE-1, JE-1 )-ACOEF* |
|
|
940 $ A( JE-1, JE-1 ) ) / TEMP |
|
|
941 WORK( 3*N+JE ) = BCOEFI*B( JE-1, JE-1 ) / TEMP |
|
|
942 END IF |
|
|
943 * |
|
|
944 XMAX = MAX( ABS( WORK( 2*N+JE ) )+ABS( WORK( 3*N+JE ) ), |
|
|
945 $ ABS( WORK( 2*N+JE-1 ) )+ABS( WORK( 3*N+JE-1 ) ) ) |
|
|
946 * |
|
|
947 * Compute contribution from columns JE and JE-1 |
|
|
948 * of A and B to the sums. |
|
|
949 * |
|
|
950 CREALA = ACOEF*WORK( 2*N+JE-1 ) |
|
|
951 CIMAGA = ACOEF*WORK( 3*N+JE-1 ) |
|
|
952 CREALB = BCOEFR*WORK( 2*N+JE-1 ) - |
|
|
953 $ BCOEFI*WORK( 3*N+JE-1 ) |
|
|
954 CIMAGB = BCOEFI*WORK( 2*N+JE-1 ) + |
|
|
955 $ BCOEFR*WORK( 3*N+JE-1 ) |
|
|
956 CRE2A = ACOEF*WORK( 2*N+JE ) |
|
|
957 CIM2A = ACOEF*WORK( 3*N+JE ) |
|
|
958 CRE2B = BCOEFR*WORK( 2*N+JE ) - BCOEFI*WORK( 3*N+JE ) |
|
|
959 CIM2B = BCOEFI*WORK( 2*N+JE ) + BCOEFR*WORK( 3*N+JE ) |
|
|
960 DO 270 JR = 1, JE - 2 |
|
|
961 WORK( 2*N+JR ) = -CREALA*A( JR, JE-1 ) + |
|
|
962 $ CREALB*B( JR, JE-1 ) - |
|
|
963 $ CRE2A*A( JR, JE ) + CRE2B*B( JR, JE ) |
|
|
964 WORK( 3*N+JR ) = -CIMAGA*A( JR, JE-1 ) + |
|
|
965 $ CIMAGB*B( JR, JE-1 ) - |
|
|
966 $ CIM2A*A( JR, JE ) + CIM2B*B( JR, JE ) |
|
|
967 270 CONTINUE |
|
|
968 END IF |
|
|
969 * |
|
|
970 DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN ) |
|
|
971 * |
|
|
972 * Columnwise triangular solve of (a A - b B) x = 0 |
|
|
973 * |
|
|
974 IL2BY2 = .FALSE. |
|
|
975 DO 370 J = JE - NW, 1, -1 |
|
|
976 * |
|
|
977 * If a 2-by-2 block, is in position j-1:j, wait until |
|
|
978 * next iteration to process it (when it will be j:j+1) |
|
|
979 * |
|
|
980 IF( .NOT.IL2BY2 .AND. J.GT.1 ) THEN |
|
|
981 IF( A( J, J-1 ).NE.ZERO ) THEN |
|
|
982 IL2BY2 = .TRUE. |
|
|
983 GO TO 370 |
|
|
984 END IF |
|
|
985 END IF |
|
|
986 BDIAG( 1 ) = B( J, J ) |
|
|
987 IF( IL2BY2 ) THEN |
|
|
988 NA = 2 |
|
|
989 BDIAG( 2 ) = B( J+1, J+1 ) |
|
|
990 ELSE |
|
|
991 NA = 1 |
|
|
992 END IF |
|
|
993 * |
|
|
994 * Compute x(j) (and x(j+1), if 2-by-2 block) |
|
|
995 * |
|
|
996 CALL DLALN2( .FALSE., NA, NW, DMIN, ACOEF, A( J, J ), |
|
|
997 $ LDA, BDIAG( 1 ), BDIAG( 2 ), WORK( 2*N+J ), |
|
|
998 $ N, BCOEFR, BCOEFI, SUM, 2, SCALE, TEMP, |
|
|
999 $ IINFO ) |
|
|
1000 IF( SCALE.LT.ONE ) THEN |
|
|
1001 * |
|
|
1002 DO 290 JW = 0, NW - 1 |
|
|
1003 DO 280 JR = 1, JE |
|
|
1004 WORK( ( JW+2 )*N+JR ) = SCALE* |
|
|
1005 $ WORK( ( JW+2 )*N+JR ) |
|
|
1006 280 CONTINUE |
|
|
1007 290 CONTINUE |
|
|
1008 END IF |
|
|
1009 XMAX = MAX( SCALE*XMAX, TEMP ) |
|
|
1010 * |
|
|
1011 DO 310 JW = 1, NW |
|
|
1012 DO 300 JA = 1, NA |
|
|
1013 WORK( ( JW+1 )*N+J+JA-1 ) = SUM( JA, JW ) |
|
|
1014 300 CONTINUE |
|
|
1015 310 CONTINUE |
|
|
1016 * |
|
|
1017 * w = w + x(j)*(a A(*,j) - b B(*,j) ) with scaling |
|
|
1018 * |
|
|
1019 IF( J.GT.1 ) THEN |
|
|
1020 * |
|
|
1021 * Check whether scaling is necessary for sum. |
|
|
1022 * |
|
|
1023 XSCALE = ONE / MAX( ONE, XMAX ) |
|
|
1024 TEMP = ACOEFA*WORK( J ) + BCOEFA*WORK( N+J ) |
|
|
1025 IF( IL2BY2 ) |
|
|
1026 $ TEMP = MAX( TEMP, ACOEFA*WORK( J+1 )+BCOEFA* |
|
|
1027 $ WORK( N+J+1 ) ) |
|
|
1028 TEMP = MAX( TEMP, ACOEFA, BCOEFA ) |
|
|
1029 IF( TEMP.GT.BIGNUM*XSCALE ) THEN |
|
|
1030 * |
|
|
1031 DO 330 JW = 0, NW - 1 |
|
|
1032 DO 320 JR = 1, JE |
|
|
1033 WORK( ( JW+2 )*N+JR ) = XSCALE* |
|
|
1034 $ WORK( ( JW+2 )*N+JR ) |
|
|
1035 320 CONTINUE |
|
|
1036 330 CONTINUE |
|
|
1037 XMAX = XMAX*XSCALE |
|
|
1038 END IF |
|
|
1039 * |
|
|
1040 * Compute the contributions of the off-diagonals of |
|
|
1041 * column j (and j+1, if 2-by-2 block) of A and B to the |
|
|
1042 * sums. |
|
|
1043 * |
|
|
1044 * |
|
|
1045 DO 360 JA = 1, NA |
|
|
1046 IF( ILCPLX ) THEN |
|
|
1047 CREALA = ACOEF*WORK( 2*N+J+JA-1 ) |
|
|
1048 CIMAGA = ACOEF*WORK( 3*N+J+JA-1 ) |
|
|
1049 CREALB = BCOEFR*WORK( 2*N+J+JA-1 ) - |
|
|
1050 $ BCOEFI*WORK( 3*N+J+JA-1 ) |
|
|
1051 CIMAGB = BCOEFI*WORK( 2*N+J+JA-1 ) + |
|
|
1052 $ BCOEFR*WORK( 3*N+J+JA-1 ) |
|
|
1053 DO 340 JR = 1, J - 1 |
|
|
1054 WORK( 2*N+JR ) = WORK( 2*N+JR ) - |
|
|
1055 $ CREALA*A( JR, J+JA-1 ) + |
|
|
1056 $ CREALB*B( JR, J+JA-1 ) |
|
|
1057 WORK( 3*N+JR ) = WORK( 3*N+JR ) - |
|
|
1058 $ CIMAGA*A( JR, J+JA-1 ) + |
|
|
1059 $ CIMAGB*B( JR, J+JA-1 ) |
|
|
1060 340 CONTINUE |
|
|
1061 ELSE |
|
|
1062 CREALA = ACOEF*WORK( 2*N+J+JA-1 ) |
|
|
1063 CREALB = BCOEFR*WORK( 2*N+J+JA-1 ) |
|
|
1064 DO 350 JR = 1, J - 1 |
|
|
1065 WORK( 2*N+JR ) = WORK( 2*N+JR ) - |
|
|
1066 $ CREALA*A( JR, J+JA-1 ) + |
|
|
1067 $ CREALB*B( JR, J+JA-1 ) |
|
|
1068 350 CONTINUE |
|
|
1069 END IF |
|
|
1070 360 CONTINUE |
|
|
1071 END IF |
|
|
1072 * |
|
|
1073 IL2BY2 = .FALSE. |
|
|
1074 370 CONTINUE |
|
|
1075 * |
|
|
1076 * Copy eigenvector to VR, back transforming if |
|
|
1077 * HOWMNY='B'. |
|
|
1078 * |
|
|
1079 IEIG = IEIG - NW |
|
|
1080 IF( ILBACK ) THEN |
|
|
1081 * |
|
|
1082 DO 410 JW = 0, NW - 1 |
|
|
1083 DO 380 JR = 1, N |
|
|
1084 WORK( ( JW+4 )*N+JR ) = WORK( ( JW+2 )*N+1 )* |
|
|
1085 $ VR( JR, 1 ) |
|
|
1086 380 CONTINUE |
|
|
1087 * |
|
|
1088 * A series of compiler directives to defeat |
|
|
1089 * vectorization for the next loop |
|
|
1090 * |
|
|
1091 * |
|
|
1092 DO 400 JC = 2, JE |
|
|
1093 DO 390 JR = 1, N |
|
|
1094 WORK( ( JW+4 )*N+JR ) = WORK( ( JW+4 )*N+JR ) + |
|
|
1095 $ WORK( ( JW+2 )*N+JC )*VR( JR, JC ) |
|
|
1096 390 CONTINUE |
|
|
1097 400 CONTINUE |
|
|
1098 410 CONTINUE |
|
|
1099 * |
|
|
1100 DO 430 JW = 0, NW - 1 |
|
|
1101 DO 420 JR = 1, N |
|
|
1102 VR( JR, IEIG+JW ) = WORK( ( JW+4 )*N+JR ) |
|
|
1103 420 CONTINUE |
|
|
1104 430 CONTINUE |
|
|
1105 * |
|
|
1106 IEND = N |
|
|
1107 ELSE |
|
|
1108 DO 450 JW = 0, NW - 1 |
|
|
1109 DO 440 JR = 1, N |
|
|
1110 VR( JR, IEIG+JW ) = WORK( ( JW+2 )*N+JR ) |
|
|
1111 440 CONTINUE |
|
|
1112 450 CONTINUE |
|
|
1113 * |
|
|
1114 IEND = JE |
|
|
1115 END IF |
|
|
1116 * |
|
|
1117 * Scale eigenvector |
|
|
1118 * |
|
|
1119 XMAX = ZERO |
|
|
1120 IF( ILCPLX ) THEN |
|
|
1121 DO 460 J = 1, IEND |
|
|
1122 XMAX = MAX( XMAX, ABS( VR( J, IEIG ) )+ |
|
|
1123 $ ABS( VR( J, IEIG+1 ) ) ) |
|
|
1124 460 CONTINUE |
|
|
1125 ELSE |
|
|
1126 DO 470 J = 1, IEND |
|
|
1127 XMAX = MAX( XMAX, ABS( VR( J, IEIG ) ) ) |
|
|
1128 470 CONTINUE |
|
|
1129 END IF |
|
|
1130 * |
|
|
1131 IF( XMAX.GT.SAFMIN ) THEN |
|
|
1132 XSCALE = ONE / XMAX |
|
|
1133 DO 490 JW = 0, NW - 1 |
|
|
1134 DO 480 JR = 1, IEND |
|
|
1135 VR( JR, IEIG+JW ) = XSCALE*VR( JR, IEIG+JW ) |
|
|
1136 480 CONTINUE |
|
|
1137 490 CONTINUE |
|
|
1138 END IF |
|
|
1139 500 CONTINUE |
|
|
1140 END IF |
|
|
1141 * |
|
|
1142 RETURN |
|
|
1143 * |
|
|
1144 * End of DTGEVC |
|
|
1145 * |
|
|
1146 END |
