1: *> \brief \b ZTGEVC
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download ZTGEVC + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztgevc.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztgevc.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztgevc.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZTGEVC( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL,
22: * LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO )
23: *
24: * .. Scalar Arguments ..
25: * CHARACTER HOWMNY, SIDE
26: * INTEGER INFO, LDP, LDS, LDVL, LDVR, M, MM, N
27: * ..
28: * .. Array Arguments ..
29: * LOGICAL SELECT( * )
30: * DOUBLE PRECISION RWORK( * )
31: * COMPLEX*16 P( LDP, * ), S( LDS, * ), VL( LDVL, * ),
32: * $ VR( LDVR, * ), WORK( * )
33: * ..
34: *
35: *
36: *
37: *> \par Purpose:
38: * =============
39: *>
40: *> \verbatim
41: *>
42: *> ZTGEVC computes some or all of the right and/or left eigenvectors of
43: *> a pair of complex matrices (S,P), where S and P are upper triangular.
44: *> Matrix pairs of this type are produced by the generalized Schur
45: *> factorization of a complex matrix pair (A,B):
46: *>
47: *> A = Q*S*Z**H, B = Q*P*Z**H
48: *>
49: *> as computed by ZGGHRD + ZHGEQZ.
50: *>
51: *> The right eigenvector x and the left eigenvector y of (S,P)
52: *> corresponding to an eigenvalue w are defined by:
53: *>
54: *> S*x = w*P*x, (y**H)*S = w*(y**H)*P,
55: *>
56: *> where y**H denotes the conjugate tranpose of y.
57: *> The eigenvalues are not input to this routine, but are computed
58: *> directly from the diagonal elements of S and P.
59: *>
60: *> This routine returns the matrices X and/or Y of right and left
61: *> eigenvectors of (S,P), or the products Z*X and/or Q*Y,
62: *> where Z and Q are input matrices.
63: *> If Q and Z are the unitary factors from the generalized Schur
64: *> factorization of a matrix pair (A,B), then Z*X and Q*Y
65: *> are the matrices of right and left eigenvectors of (A,B).
66: *> \endverbatim
67: *
68: * Arguments:
69: * ==========
70: *
71: *> \param[in] SIDE
72: *> \verbatim
73: *> SIDE is CHARACTER*1
74: *> = 'R': compute right eigenvectors only;
75: *> = 'L': compute left eigenvectors only;
76: *> = 'B': compute both right and left eigenvectors.
77: *> \endverbatim
78: *>
79: *> \param[in] HOWMNY
80: *> \verbatim
81: *> HOWMNY is CHARACTER*1
82: *> = 'A': compute all right and/or left eigenvectors;
83: *> = 'B': compute all right and/or left eigenvectors,
84: *> backtransformed by the matrices in VR and/or VL;
85: *> = 'S': compute selected right and/or left eigenvectors,
86: *> specified by the logical array SELECT.
87: *> \endverbatim
88: *>
89: *> \param[in] SELECT
90: *> \verbatim
91: *> SELECT is LOGICAL array, dimension (N)
92: *> If HOWMNY='S', SELECT specifies the eigenvectors to be
93: *> computed. The eigenvector corresponding to the j-th
94: *> eigenvalue is computed if SELECT(j) = .TRUE..
95: *> Not referenced if HOWMNY = 'A' or 'B'.
96: *> \endverbatim
97: *>
98: *> \param[in] N
99: *> \verbatim
100: *> N is INTEGER
101: *> The order of the matrices S and P. N >= 0.
102: *> \endverbatim
103: *>
104: *> \param[in] S
105: *> \verbatim
106: *> S is COMPLEX*16 array, dimension (LDS,N)
107: *> The upper triangular matrix S from a generalized Schur
108: *> factorization, as computed by ZHGEQZ.
109: *> \endverbatim
110: *>
111: *> \param[in] LDS
112: *> \verbatim
113: *> LDS is INTEGER
114: *> The leading dimension of array S. LDS >= max(1,N).
115: *> \endverbatim
116: *>
117: *> \param[in] P
118: *> \verbatim
119: *> P is COMPLEX*16 array, dimension (LDP,N)
120: *> The upper triangular matrix P from a generalized Schur
121: *> factorization, as computed by ZHGEQZ. P must have real
122: *> diagonal elements.
123: *> \endverbatim
124: *>
125: *> \param[in] LDP
126: *> \verbatim
127: *> LDP is INTEGER
128: *> The leading dimension of array P. LDP >= max(1,N).
129: *> \endverbatim
130: *>
131: *> \param[in,out] VL
132: *> \verbatim
133: *> VL is COMPLEX*16 array, dimension (LDVL,MM)
134: *> On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
135: *> contain an N-by-N matrix Q (usually the unitary matrix Q
136: *> of left Schur vectors returned by ZHGEQZ).
137: *> On exit, if SIDE = 'L' or 'B', VL contains:
138: *> if HOWMNY = 'A', the matrix Y of left eigenvectors of (S,P);
139: *> if HOWMNY = 'B', the matrix Q*Y;
140: *> if HOWMNY = 'S', the left eigenvectors of (S,P) specified by
141: *> SELECT, stored consecutively in the columns of
142: *> VL, in the same order as their eigenvalues.
143: *> Not referenced if SIDE = 'R'.
144: *> \endverbatim
145: *>
146: *> \param[in] LDVL
147: *> \verbatim
148: *> LDVL is INTEGER
149: *> The leading dimension of array VL. LDVL >= 1, and if
150: *> SIDE = 'L' or 'l' or 'B' or 'b', LDVL >= N.
151: *> \endverbatim
152: *>
153: *> \param[in,out] VR
154: *> \verbatim
155: *> VR is COMPLEX*16 array, dimension (LDVR,MM)
156: *> On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
157: *> contain an N-by-N matrix Q (usually the unitary matrix Z
158: *> of right Schur vectors returned by ZHGEQZ).
159: *> On exit, if SIDE = 'R' or 'B', VR contains:
160: *> if HOWMNY = 'A', the matrix X of right eigenvectors of (S,P);
161: *> if HOWMNY = 'B', the matrix Z*X;
162: *> if HOWMNY = 'S', the right eigenvectors of (S,P) specified by
163: *> SELECT, stored consecutively in the columns of
164: *> VR, in the same order as their eigenvalues.
165: *> Not referenced if SIDE = 'L'.
166: *> \endverbatim
167: *>
168: *> \param[in] LDVR
169: *> \verbatim
170: *> LDVR is INTEGER
171: *> The leading dimension of the array VR. LDVR >= 1, and if
172: *> SIDE = 'R' or 'B', LDVR >= N.
173: *> \endverbatim
174: *>
175: *> \param[in] MM
176: *> \verbatim
177: *> MM is INTEGER
178: *> The number of columns in the arrays VL and/or VR. MM >= M.
179: *> \endverbatim
180: *>
181: *> \param[out] M
182: *> \verbatim
183: *> M is INTEGER
184: *> The number of columns in the arrays VL and/or VR actually
185: *> used to store the eigenvectors. If HOWMNY = 'A' or 'B', M
186: *> is set to N. Each selected eigenvector occupies one column.
187: *> \endverbatim
188: *>
189: *> \param[out] WORK
190: *> \verbatim
191: *> WORK is COMPLEX*16 array, dimension (2*N)
192: *> \endverbatim
193: *>
194: *> \param[out] RWORK
195: *> \verbatim
196: *> RWORK is DOUBLE PRECISION array, dimension (2*N)
197: *> \endverbatim
198: *>
199: *> \param[out] INFO
200: *> \verbatim
201: *> INFO is INTEGER
202: *> = 0: successful exit.
203: *> < 0: if INFO = -i, the i-th argument had an illegal value.
204: *> \endverbatim
205: *
206: * Authors:
207: * ========
208: *
209: *> \author Univ. of Tennessee
210: *> \author Univ. of California Berkeley
211: *> \author Univ. of Colorado Denver
212: *> \author NAG Ltd.
213: *
214: *> \ingroup complex16GEcomputational
215: *
216: * =====================================================================
217: SUBROUTINE ZTGEVC( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL,
218: $ LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO )
219: *
220: * -- LAPACK computational routine --
221: * -- LAPACK is a software package provided by Univ. of Tennessee, --
222: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
223: *
224: * .. Scalar Arguments ..
225: CHARACTER HOWMNY, SIDE
226: INTEGER INFO, LDP, LDS, LDVL, LDVR, M, MM, N
227: * ..
228: * .. Array Arguments ..
229: LOGICAL SELECT( * )
230: DOUBLE PRECISION RWORK( * )
231: COMPLEX*16 P( LDP, * ), S( LDS, * ), VL( LDVL, * ),
232: $ VR( LDVR, * ), WORK( * )
233: * ..
234: *
235: *
236: * =====================================================================
237: *
238: * .. Parameters ..
239: DOUBLE PRECISION ZERO, ONE
240: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
241: COMPLEX*16 CZERO, CONE
242: PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
243: $ CONE = ( 1.0D+0, 0.0D+0 ) )
244: * ..
245: * .. Local Scalars ..
246: LOGICAL COMPL, COMPR, ILALL, ILBACK, ILBBAD, ILCOMP,
247: $ LSA, LSB
248: INTEGER I, IBEG, IEIG, IEND, IHWMNY, IM, ISIDE, ISRC,
249: $ J, JE, JR
250: DOUBLE PRECISION ACOEFA, ACOEFF, ANORM, ASCALE, BCOEFA, BIG,
251: $ BIGNUM, BNORM, BSCALE, DMIN, SAFMIN, SBETA,
252: $ SCALE, SMALL, TEMP, ULP, XMAX
253: COMPLEX*16 BCOEFF, CA, CB, D, SALPHA, SUM, SUMA, SUMB, X
254: * ..
255: * .. External Functions ..
256: LOGICAL LSAME
257: DOUBLE PRECISION DLAMCH
258: COMPLEX*16 ZLADIV
259: EXTERNAL LSAME, DLAMCH, ZLADIV
260: * ..
261: * .. External Subroutines ..
262: EXTERNAL DLABAD, XERBLA, ZGEMV
263: * ..
264: * .. Intrinsic Functions ..
265: INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX, MIN
266: * ..
267: * .. Statement Functions ..
268: DOUBLE PRECISION ABS1
269: * ..
270: * .. Statement Function definitions ..
271: ABS1( X ) = ABS( DBLE( X ) ) + ABS( DIMAG( X ) )
272: * ..
273: * .. Executable Statements ..
274: *
275: * Decode and Test the input parameters
276: *
277: IF( LSAME( HOWMNY, 'A' ) ) THEN
278: IHWMNY = 1
279: ILALL = .TRUE.
280: ILBACK = .FALSE.
281: ELSE IF( LSAME( HOWMNY, 'S' ) ) THEN
282: IHWMNY = 2
283: ILALL = .FALSE.
284: ILBACK = .FALSE.
285: ELSE IF( LSAME( HOWMNY, 'B' ) ) THEN
286: IHWMNY = 3
287: ILALL = .TRUE.
288: ILBACK = .TRUE.
289: ELSE
290: IHWMNY = -1
291: END IF
292: *
293: IF( LSAME( SIDE, 'R' ) ) THEN
294: ISIDE = 1
295: COMPL = .FALSE.
296: COMPR = .TRUE.
297: ELSE IF( LSAME( SIDE, 'L' ) ) THEN
298: ISIDE = 2
299: COMPL = .TRUE.
300: COMPR = .FALSE.
301: ELSE IF( LSAME( SIDE, 'B' ) ) THEN
302: ISIDE = 3
303: COMPL = .TRUE.
304: COMPR = .TRUE.
305: ELSE
306: ISIDE = -1
307: END IF
308: *
309: INFO = 0
310: IF( ISIDE.LT.0 ) THEN
311: INFO = -1
312: ELSE IF( IHWMNY.LT.0 ) THEN
313: INFO = -2
314: ELSE IF( N.LT.0 ) THEN
315: INFO = -4
316: ELSE IF( LDS.LT.MAX( 1, N ) ) THEN
317: INFO = -6
318: ELSE IF( LDP.LT.MAX( 1, N ) ) THEN
319: INFO = -8
320: END IF
321: IF( INFO.NE.0 ) THEN
322: CALL XERBLA( 'ZTGEVC', -INFO )
323: RETURN
324: END IF
325: *
326: * Count the number of eigenvectors
327: *
328: IF( .NOT.ILALL ) THEN
329: IM = 0
330: DO 10 J = 1, N
331: IF( SELECT( J ) )
332: $ IM = IM + 1
333: 10 CONTINUE
334: ELSE
335: IM = N
336: END IF
337: *
338: * Check diagonal of B
339: *
340: ILBBAD = .FALSE.
341: DO 20 J = 1, N
342: IF( DIMAG( P( J, J ) ).NE.ZERO )
343: $ ILBBAD = .TRUE.
344: 20 CONTINUE
345: *
346: IF( ILBBAD ) THEN
347: INFO = -7
348: ELSE IF( COMPL .AND. LDVL.LT.N .OR. LDVL.LT.1 ) THEN
349: INFO = -10
350: ELSE IF( COMPR .AND. LDVR.LT.N .OR. LDVR.LT.1 ) THEN
351: INFO = -12
352: ELSE IF( MM.LT.IM ) THEN
353: INFO = -13
354: END IF
355: IF( INFO.NE.0 ) THEN
356: CALL XERBLA( 'ZTGEVC', -INFO )
357: RETURN
358: END IF
359: *
360: * Quick return if possible
361: *
362: M = IM
363: IF( N.EQ.0 )
364: $ RETURN
365: *
366: * Machine Constants
367: *
368: SAFMIN = DLAMCH( 'Safe minimum' )
369: BIG = ONE / SAFMIN
370: CALL DLABAD( SAFMIN, BIG )
371: ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' )
372: SMALL = SAFMIN*N / ULP
373: BIG = ONE / SMALL
374: BIGNUM = ONE / ( SAFMIN*N )
375: *
376: * Compute the 1-norm of each column of the strictly upper triangular
377: * part of A and B to check for possible overflow in the triangular
378: * solver.
379: *
380: ANORM = ABS1( S( 1, 1 ) )
381: BNORM = ABS1( P( 1, 1 ) )
382: RWORK( 1 ) = ZERO
383: RWORK( N+1 ) = ZERO
384: DO 40 J = 2, N
385: RWORK( J ) = ZERO
386: RWORK( N+J ) = ZERO
387: DO 30 I = 1, J - 1
388: RWORK( J ) = RWORK( J ) + ABS1( S( I, J ) )
389: RWORK( N+J ) = RWORK( N+J ) + ABS1( P( I, J ) )
390: 30 CONTINUE
391: ANORM = MAX( ANORM, RWORK( J )+ABS1( S( J, J ) ) )
392: BNORM = MAX( BNORM, RWORK( N+J )+ABS1( P( J, J ) ) )
393: 40 CONTINUE
394: *
395: ASCALE = ONE / MAX( ANORM, SAFMIN )
396: BSCALE = ONE / MAX( BNORM, SAFMIN )
397: *
398: * Left eigenvectors
399: *
400: IF( COMPL ) THEN
401: IEIG = 0
402: *
403: * Main loop over eigenvalues
404: *
405: DO 140 JE = 1, N
406: IF( ILALL ) THEN
407: ILCOMP = .TRUE.
408: ELSE
409: ILCOMP = SELECT( JE )
410: END IF
411: IF( ILCOMP ) THEN
412: IEIG = IEIG + 1
413: *
414: IF( ABS1( S( JE, JE ) ).LE.SAFMIN .AND.
415: $ ABS( DBLE( P( JE, JE ) ) ).LE.SAFMIN ) THEN
416: *
417: * Singular matrix pencil -- return unit eigenvector
418: *
419: DO 50 JR = 1, N
420: VL( JR, IEIG ) = CZERO
421: 50 CONTINUE
422: VL( IEIG, IEIG ) = CONE
423: GO TO 140
424: END IF
425: *
426: * Non-singular eigenvalue:
427: * Compute coefficients a and b in
428: * H
429: * y ( a A - b B ) = 0
430: *
431: TEMP = ONE / MAX( ABS1( S( JE, JE ) )*ASCALE,
432: $ ABS( DBLE( P( JE, JE ) ) )*BSCALE, SAFMIN )
433: SALPHA = ( TEMP*S( JE, JE ) )*ASCALE
434: SBETA = ( TEMP*DBLE( P( JE, JE ) ) )*BSCALE
435: ACOEFF = SBETA*ASCALE
436: BCOEFF = SALPHA*BSCALE
437: *
438: * Scale to avoid underflow
439: *
440: LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEFF ).LT.SMALL
441: LSB = ABS1( SALPHA ).GE.SAFMIN .AND. ABS1( BCOEFF ).LT.
442: $ SMALL
443: *
444: SCALE = ONE
445: IF( LSA )
446: $ SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG )
447: IF( LSB )
448: $ SCALE = MAX( SCALE, ( SMALL / ABS1( SALPHA ) )*
449: $ MIN( BNORM, BIG ) )
450: IF( LSA .OR. LSB ) THEN
451: SCALE = MIN( SCALE, ONE /
452: $ ( SAFMIN*MAX( ONE, ABS( ACOEFF ),
453: $ ABS1( BCOEFF ) ) ) )
454: IF( LSA ) THEN
455: ACOEFF = ASCALE*( SCALE*SBETA )
456: ELSE
457: ACOEFF = SCALE*ACOEFF
458: END IF
459: IF( LSB ) THEN
460: BCOEFF = BSCALE*( SCALE*SALPHA )
461: ELSE
462: BCOEFF = SCALE*BCOEFF
463: END IF
464: END IF
465: *
466: ACOEFA = ABS( ACOEFF )
467: BCOEFA = ABS1( BCOEFF )
468: XMAX = ONE
469: DO 60 JR = 1, N
470: WORK( JR ) = CZERO
471: 60 CONTINUE
472: WORK( JE ) = CONE
473: DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN )
474: *
475: * H
476: * Triangular solve of (a A - b B) y = 0
477: *
478: * H
479: * (rowwise in (a A - b B) , or columnwise in a A - b B)
480: *
481: DO 100 J = JE + 1, N
482: *
483: * Compute
484: * j-1
485: * SUM = sum conjg( a*S(k,j) - b*P(k,j) )*x(k)
486: * k=je
487: * (Scale if necessary)
488: *
489: TEMP = ONE / XMAX
490: IF( ACOEFA*RWORK( J )+BCOEFA*RWORK( N+J ).GT.BIGNUM*
491: $ TEMP ) THEN
492: DO 70 JR = JE, J - 1
493: WORK( JR ) = TEMP*WORK( JR )
494: 70 CONTINUE
495: XMAX = ONE
496: END IF
497: SUMA = CZERO
498: SUMB = CZERO
499: *
500: DO 80 JR = JE, J - 1
501: SUMA = SUMA + DCONJG( S( JR, J ) )*WORK( JR )
502: SUMB = SUMB + DCONJG( P( JR, J ) )*WORK( JR )
503: 80 CONTINUE
504: SUM = ACOEFF*SUMA - DCONJG( BCOEFF )*SUMB
505: *
506: * Form x(j) = - SUM / conjg( a*S(j,j) - b*P(j,j) )
507: *
508: * with scaling and perturbation of the denominator
509: *
510: D = DCONJG( ACOEFF*S( J, J )-BCOEFF*P( J, J ) )
511: IF( ABS1( D ).LE.DMIN )
512: $ D = DCMPLX( DMIN )
513: *
514: IF( ABS1( D ).LT.ONE ) THEN
515: IF( ABS1( SUM ).GE.BIGNUM*ABS1( D ) ) THEN
516: TEMP = ONE / ABS1( SUM )
517: DO 90 JR = JE, J - 1
518: WORK( JR ) = TEMP*WORK( JR )
519: 90 CONTINUE
520: XMAX = TEMP*XMAX
521: SUM = TEMP*SUM
522: END IF
523: END IF
524: WORK( J ) = ZLADIV( -SUM, D )
525: XMAX = MAX( XMAX, ABS1( WORK( J ) ) )
526: 100 CONTINUE
527: *
528: * Back transform eigenvector if HOWMNY='B'.
529: *
530: IF( ILBACK ) THEN
531: CALL ZGEMV( 'N', N, N+1-JE, CONE, VL( 1, JE ), LDVL,
532: $ WORK( JE ), 1, CZERO, WORK( N+1 ), 1 )
533: ISRC = 2
534: IBEG = 1
535: ELSE
536: ISRC = 1
537: IBEG = JE
538: END IF
539: *
540: * Copy and scale eigenvector into column of VL
541: *
542: XMAX = ZERO
543: DO 110 JR = IBEG, N
544: XMAX = MAX( XMAX, ABS1( WORK( ( ISRC-1 )*N+JR ) ) )
545: 110 CONTINUE
546: *
547: IF( XMAX.GT.SAFMIN ) THEN
548: TEMP = ONE / XMAX
549: DO 120 JR = IBEG, N
550: VL( JR, IEIG ) = TEMP*WORK( ( ISRC-1 )*N+JR )
551: 120 CONTINUE
552: ELSE
553: IBEG = N + 1
554: END IF
555: *
556: DO 130 JR = 1, IBEG - 1
557: VL( JR, IEIG ) = CZERO
558: 130 CONTINUE
559: *
560: END IF
561: 140 CONTINUE
562: END IF
563: *
564: * Right eigenvectors
565: *
566: IF( COMPR ) THEN
567: IEIG = IM + 1
568: *
569: * Main loop over eigenvalues
570: *
571: DO 250 JE = N, 1, -1
572: IF( ILALL ) THEN
573: ILCOMP = .TRUE.
574: ELSE
575: ILCOMP = SELECT( JE )
576: END IF
577: IF( ILCOMP ) THEN
578: IEIG = IEIG - 1
579: *
580: IF( ABS1( S( JE, JE ) ).LE.SAFMIN .AND.
581: $ ABS( DBLE( P( JE, JE ) ) ).LE.SAFMIN ) THEN
582: *
583: * Singular matrix pencil -- return unit eigenvector
584: *
585: DO 150 JR = 1, N
586: VR( JR, IEIG ) = CZERO
587: 150 CONTINUE
588: VR( IEIG, IEIG ) = CONE
589: GO TO 250
590: END IF
591: *
592: * Non-singular eigenvalue:
593: * Compute coefficients a and b in
594: *
595: * ( a A - b B ) x = 0
596: *
597: TEMP = ONE / MAX( ABS1( S( JE, JE ) )*ASCALE,
598: $ ABS( DBLE( P( JE, JE ) ) )*BSCALE, SAFMIN )
599: SALPHA = ( TEMP*S( JE, JE ) )*ASCALE
600: SBETA = ( TEMP*DBLE( P( JE, JE ) ) )*BSCALE
601: ACOEFF = SBETA*ASCALE
602: BCOEFF = SALPHA*BSCALE
603: *
604: * Scale to avoid underflow
605: *
606: LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEFF ).LT.SMALL
607: LSB = ABS1( SALPHA ).GE.SAFMIN .AND. ABS1( BCOEFF ).LT.
608: $ SMALL
609: *
610: SCALE = ONE
611: IF( LSA )
612: $ SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG )
613: IF( LSB )
614: $ SCALE = MAX( SCALE, ( SMALL / ABS1( SALPHA ) )*
615: $ MIN( BNORM, BIG ) )
616: IF( LSA .OR. LSB ) THEN
617: SCALE = MIN( SCALE, ONE /
618: $ ( SAFMIN*MAX( ONE, ABS( ACOEFF ),
619: $ ABS1( BCOEFF ) ) ) )
620: IF( LSA ) THEN
621: ACOEFF = ASCALE*( SCALE*SBETA )
622: ELSE
623: ACOEFF = SCALE*ACOEFF
624: END IF
625: IF( LSB ) THEN
626: BCOEFF = BSCALE*( SCALE*SALPHA )
627: ELSE
628: BCOEFF = SCALE*BCOEFF
629: END IF
630: END IF
631: *
632: ACOEFA = ABS( ACOEFF )
633: BCOEFA = ABS1( BCOEFF )
634: XMAX = ONE
635: DO 160 JR = 1, N
636: WORK( JR ) = CZERO
637: 160 CONTINUE
638: WORK( JE ) = CONE
639: DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN )
640: *
641: * Triangular solve of (a A - b B) x = 0 (columnwise)
642: *
643: * WORK(1:j-1) contains sums w,
644: * WORK(j+1:JE) contains x
645: *
646: DO 170 JR = 1, JE - 1
647: WORK( JR ) = ACOEFF*S( JR, JE ) - BCOEFF*P( JR, JE )
648: 170 CONTINUE
649: WORK( JE ) = CONE
650: *
651: DO 210 J = JE - 1, 1, -1
652: *
653: * Form x(j) := - w(j) / d
654: * with scaling and perturbation of the denominator
655: *
656: D = ACOEFF*S( J, J ) - BCOEFF*P( J, J )
657: IF( ABS1( D ).LE.DMIN )
658: $ D = DCMPLX( DMIN )
659: *
660: IF( ABS1( D ).LT.ONE ) THEN
661: IF( ABS1( WORK( J ) ).GE.BIGNUM*ABS1( D ) ) THEN
662: TEMP = ONE / ABS1( WORK( J ) )
663: DO 180 JR = 1, JE
664: WORK( JR ) = TEMP*WORK( JR )
665: 180 CONTINUE
666: END IF
667: END IF
668: *
669: WORK( J ) = ZLADIV( -WORK( J ), D )
670: *
671: IF( J.GT.1 ) THEN
672: *
673: * w = w + x(j)*(a S(*,j) - b P(*,j) ) with scaling
674: *
675: IF( ABS1( WORK( J ) ).GT.ONE ) THEN
676: TEMP = ONE / ABS1( WORK( J ) )
677: IF( ACOEFA*RWORK( J )+BCOEFA*RWORK( N+J ).GE.
678: $ BIGNUM*TEMP ) THEN
679: DO 190 JR = 1, JE
680: WORK( JR ) = TEMP*WORK( JR )
681: 190 CONTINUE
682: END IF
683: END IF
684: *
685: CA = ACOEFF*WORK( J )
686: CB = BCOEFF*WORK( J )
687: DO 200 JR = 1, J - 1
688: WORK( JR ) = WORK( JR ) + CA*S( JR, J ) -
689: $ CB*P( JR, J )
690: 200 CONTINUE
691: END IF
692: 210 CONTINUE
693: *
694: * Back transform eigenvector if HOWMNY='B'.
695: *
696: IF( ILBACK ) THEN
697: CALL ZGEMV( 'N', N, JE, CONE, VR, LDVR, WORK, 1,
698: $ CZERO, WORK( N+1 ), 1 )
699: ISRC = 2
700: IEND = N
701: ELSE
702: ISRC = 1
703: IEND = JE
704: END IF
705: *
706: * Copy and scale eigenvector into column of VR
707: *
708: XMAX = ZERO
709: DO 220 JR = 1, IEND
710: XMAX = MAX( XMAX, ABS1( WORK( ( ISRC-1 )*N+JR ) ) )
711: 220 CONTINUE
712: *
713: IF( XMAX.GT.SAFMIN ) THEN
714: TEMP = ONE / XMAX
715: DO 230 JR = 1, IEND
716: VR( JR, IEIG ) = TEMP*WORK( ( ISRC-1 )*N+JR )
717: 230 CONTINUE
718: ELSE
719: IEND = 0
720: END IF
721: *
722: DO 240 JR = IEND + 1, N
723: VR( JR, IEIG ) = CZERO
724: 240 CONTINUE
725: *
726: END IF
727: 250 CONTINUE
728: END IF
729: *
730: RETURN
731: *
732: * End of ZTGEVC
733: *
734: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to ztgevc.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack