Annotation of rpl/lapack/lapack/zgegv.f, revision 1.7
1.1 bertrand 1: SUBROUTINE ZGEGV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHA, BETA,
2: $ VL, LDVL, VR, LDVR, WORK, LWORK, RWORK, INFO )
3: *
4: * -- LAPACK driver 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 JOBVL, JOBVR
11: INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, N
12: * ..
13: * .. Array Arguments ..
14: DOUBLE PRECISION RWORK( * )
15: COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDB, * ),
16: $ BETA( * ), VL( LDVL, * ), VR( LDVR, * ),
17: $ WORK( * )
18: * ..
19: *
20: * Purpose
21: * =======
22: *
23: * This routine is deprecated and has been replaced by routine ZGGEV.
24: *
25: * ZGEGV computes the eigenvalues and, optionally, the left and/or right
26: * eigenvectors of a complex matrix pair (A,B).
27: * Given two square matrices A and B,
28: * the generalized nonsymmetric eigenvalue problem (GNEP) is to find the
29: * eigenvalues lambda and corresponding (non-zero) eigenvectors x such
30: * that
31: * A*x = lambda*B*x.
32: *
33: * An alternate form is to find the eigenvalues mu and corresponding
34: * eigenvectors y such that
35: * mu*A*y = B*y.
36: *
37: * These two forms are equivalent with mu = 1/lambda and x = y if
38: * neither lambda nor mu is zero. In order to deal with the case that
39: * lambda or mu is zero or small, two values alpha and beta are returned
40: * for each eigenvalue, such that lambda = alpha/beta and
41: * mu = beta/alpha.
42: *
43: * The vectors x and y in the above equations are right eigenvectors of
44: * the matrix pair (A,B). Vectors u and v satisfying
45: * u**H*A = lambda*u**H*B or mu*v**H*A = v**H*B
46: * are left eigenvectors of (A,B).
47: *
48: * Note: this routine performs "full balancing" on A and B -- see
49: * "Further Details", below.
50: *
51: * Arguments
52: * =========
53: *
54: * JOBVL (input) CHARACTER*1
55: * = 'N': do not compute the left generalized eigenvectors;
56: * = 'V': compute the left generalized eigenvectors (returned
57: * in VL).
58: *
59: * JOBVR (input) CHARACTER*1
60: * = 'N': do not compute the right generalized eigenvectors;
61: * = 'V': compute the right generalized eigenvectors (returned
62: * in VR).
63: *
64: * N (input) INTEGER
65: * The order of the matrices A, B, VL, and VR. N >= 0.
66: *
67: * A (input/output) COMPLEX*16 array, dimension (LDA, N)
68: * On entry, the matrix A.
69: * If JOBVL = 'V' or JOBVR = 'V', then on exit A
70: * contains the Schur form of A from the generalized Schur
71: * factorization of the pair (A,B) after balancing. If no
72: * eigenvectors were computed, then only the diagonal elements
73: * of the Schur form will be correct. See ZGGHRD and ZHGEQZ
74: * for details.
75: *
76: * LDA (input) INTEGER
77: * The leading dimension of A. LDA >= max(1,N).
78: *
79: * B (input/output) COMPLEX*16 array, dimension (LDB, N)
80: * On entry, the matrix B.
81: * If JOBVL = 'V' or JOBVR = 'V', then on exit B contains the
82: * upper triangular matrix obtained from B in the generalized
83: * Schur factorization of the pair (A,B) after balancing.
84: * If no eigenvectors were computed, then only the diagonal
85: * elements of B will be correct. See ZGGHRD and ZHGEQZ for
86: * details.
87: *
88: * LDB (input) INTEGER
89: * The leading dimension of B. LDB >= max(1,N).
90: *
91: * ALPHA (output) COMPLEX*16 array, dimension (N)
92: * The complex scalars alpha that define the eigenvalues of
93: * GNEP.
94: *
95: * BETA (output) COMPLEX*16 array, dimension (N)
96: * The complex scalars beta that define the eigenvalues of GNEP.
97: *
98: * Together, the quantities alpha = ALPHA(j) and beta = BETA(j)
99: * represent the j-th eigenvalue of the matrix pair (A,B), in
100: * one of the forms lambda = alpha/beta or mu = beta/alpha.
101: * Since either lambda or mu may overflow, they should not,
102: * in general, be computed.
103: *
104: * VL (output) COMPLEX*16 array, dimension (LDVL,N)
105: * If JOBVL = 'V', the left eigenvectors u(j) are stored
106: * in the columns of VL, in the same order as their eigenvalues.
107: * Each eigenvector is scaled so that its largest component has
108: * abs(real part) + abs(imag. part) = 1, except for eigenvectors
109: * corresponding to an eigenvalue with alpha = beta = 0, which
110: * are set to zero.
111: * Not referenced if JOBVL = 'N'.
112: *
113: * LDVL (input) INTEGER
114: * The leading dimension of the matrix VL. LDVL >= 1, and
115: * if JOBVL = 'V', LDVL >= N.
116: *
117: * VR (output) COMPLEX*16 array, dimension (LDVR,N)
118: * If JOBVR = 'V', the right eigenvectors x(j) are stored
119: * in the columns of VR, in the same order as their eigenvalues.
120: * Each eigenvector is scaled so that its largest component has
121: * abs(real part) + abs(imag. part) = 1, except for eigenvectors
122: * corresponding to an eigenvalue with alpha = beta = 0, which
123: * are set to zero.
124: * Not referenced if JOBVR = 'N'.
125: *
126: * LDVR (input) INTEGER
127: * The leading dimension of the matrix VR. LDVR >= 1, and
128: * if JOBVR = 'V', LDVR >= N.
129: *
130: * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
131: * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
132: *
133: * LWORK (input) INTEGER
134: * The dimension of the array WORK. LWORK >= max(1,2*N).
135: * For good performance, LWORK must generally be larger.
136: * To compute the optimal value of LWORK, call ILAENV to get
137: * blocksizes (for ZGEQRF, ZUNMQR, and ZUNGQR.) Then compute:
138: * NB -- MAX of the blocksizes for ZGEQRF, ZUNMQR, and ZUNGQR;
139: * The optimal LWORK is MAX( 2*N, N*(NB+1) ).
140: *
141: * If LWORK = -1, then a workspace query is assumed; the routine
142: * only calculates the optimal size of the WORK array, returns
143: * this value as the first entry of the WORK array, and no error
144: * message related to LWORK is issued by XERBLA.
145: *
146: * RWORK (workspace/output) DOUBLE PRECISION array, dimension (8*N)
147: *
148: * INFO (output) INTEGER
149: * = 0: successful exit
150: * < 0: if INFO = -i, the i-th argument had an illegal value.
151: * =1,...,N:
152: * The QZ iteration failed. No eigenvectors have been
153: * calculated, but ALPHA(j) and BETA(j) should be
154: * correct for j=INFO+1,...,N.
155: * > N: errors that usually indicate LAPACK problems:
156: * =N+1: error return from ZGGBAL
157: * =N+2: error return from ZGEQRF
158: * =N+3: error return from ZUNMQR
159: * =N+4: error return from ZUNGQR
160: * =N+5: error return from ZGGHRD
161: * =N+6: error return from ZHGEQZ (other than failed
162: * iteration)
163: * =N+7: error return from ZTGEVC
164: * =N+8: error return from ZGGBAK (computing VL)
165: * =N+9: error return from ZGGBAK (computing VR)
166: * =N+10: error return from ZLASCL (various calls)
167: *
168: * Further Details
169: * ===============
170: *
171: * Balancing
172: * ---------
173: *
174: * This driver calls ZGGBAL to both permute and scale rows and columns
175: * of A and B. The permutations PL and PR are chosen so that PL*A*PR
176: * and PL*B*R will be upper triangular except for the diagonal blocks
177: * A(i:j,i:j) and B(i:j,i:j), with i and j as close together as
178: * possible. The diagonal scaling matrices DL and DR are chosen so
179: * that the pair DL*PL*A*PR*DR, DL*PL*B*PR*DR have elements close to
180: * one (except for the elements that start out zero.)
181: *
182: * After the eigenvalues and eigenvectors of the balanced matrices
183: * have been computed, ZGGBAK transforms the eigenvectors back to what
184: * they would have been (in perfect arithmetic) if they had not been
185: * balanced.
186: *
187: * Contents of A and B on Exit
188: * -------- -- - --- - -- ----
189: *
190: * If any eigenvectors are computed (either JOBVL='V' or JOBVR='V' or
191: * both), then on exit the arrays A and B will contain the complex Schur
192: * form[*] of the "balanced" versions of A and B. If no eigenvectors
193: * are computed, then only the diagonal blocks will be correct.
194: *
195: * [*] In other words, upper triangular form.
196: *
197: * =====================================================================
198: *
199: * .. Parameters ..
200: DOUBLE PRECISION ZERO, ONE
201: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
202: COMPLEX*16 CZERO, CONE
203: PARAMETER ( CZERO = ( 0.0D0, 0.0D0 ),
204: $ CONE = ( 1.0D0, 0.0D0 ) )
205: * ..
206: * .. Local Scalars ..
207: LOGICAL ILIMIT, ILV, ILVL, ILVR, LQUERY
208: CHARACTER CHTEMP
209: INTEGER ICOLS, IHI, IINFO, IJOBVL, IJOBVR, ILEFT, ILO,
210: $ IN, IRIGHT, IROWS, IRWORK, ITAU, IWORK, JC, JR,
211: $ LOPT, LWKMIN, LWKOPT, NB, NB1, NB2, NB3
212: DOUBLE PRECISION ABSAI, ABSAR, ABSB, ANRM, ANRM1, ANRM2, BNRM,
213: $ BNRM1, BNRM2, EPS, SAFMAX, SAFMIN, SALFAI,
214: $ SALFAR, SBETA, SCALE, TEMP
215: COMPLEX*16 X
216: * ..
217: * .. Local Arrays ..
218: LOGICAL LDUMMA( 1 )
219: * ..
220: * .. External Subroutines ..
221: EXTERNAL XERBLA, ZGEQRF, ZGGBAK, ZGGBAL, ZGGHRD, ZHGEQZ,
222: $ ZLACPY, ZLASCL, ZLASET, ZTGEVC, ZUNGQR, ZUNMQR
223: * ..
224: * .. External Functions ..
225: LOGICAL LSAME
226: INTEGER ILAENV
227: DOUBLE PRECISION DLAMCH, ZLANGE
228: EXTERNAL LSAME, ILAENV, DLAMCH, ZLANGE
229: * ..
230: * .. Intrinsic Functions ..
231: INTRINSIC ABS, DBLE, DCMPLX, DIMAG, INT, MAX
232: * ..
233: * .. Statement Functions ..
234: DOUBLE PRECISION ABS1
235: * ..
236: * .. Statement Function definitions ..
237: ABS1( X ) = ABS( DBLE( X ) ) + ABS( DIMAG( X ) )
238: * ..
239: * .. Executable Statements ..
240: *
241: * Decode the input arguments
242: *
243: IF( LSAME( JOBVL, 'N' ) ) THEN
244: IJOBVL = 1
245: ILVL = .FALSE.
246: ELSE IF( LSAME( JOBVL, 'V' ) ) THEN
247: IJOBVL = 2
248: ILVL = .TRUE.
249: ELSE
250: IJOBVL = -1
251: ILVL = .FALSE.
252: END IF
253: *
254: IF( LSAME( JOBVR, 'N' ) ) THEN
255: IJOBVR = 1
256: ILVR = .FALSE.
257: ELSE IF( LSAME( JOBVR, 'V' ) ) THEN
258: IJOBVR = 2
259: ILVR = .TRUE.
260: ELSE
261: IJOBVR = -1
262: ILVR = .FALSE.
263: END IF
264: ILV = ILVL .OR. ILVR
265: *
266: * Test the input arguments
267: *
268: LWKMIN = MAX( 2*N, 1 )
269: LWKOPT = LWKMIN
270: WORK( 1 ) = LWKOPT
271: LQUERY = ( LWORK.EQ.-1 )
272: INFO = 0
273: IF( IJOBVL.LE.0 ) THEN
274: INFO = -1
275: ELSE IF( IJOBVR.LE.0 ) THEN
276: INFO = -2
277: ELSE IF( N.LT.0 ) THEN
278: INFO = -3
279: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
280: INFO = -5
281: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
282: INFO = -7
283: ELSE IF( LDVL.LT.1 .OR. ( ILVL .AND. LDVL.LT.N ) ) THEN
284: INFO = -11
285: ELSE IF( LDVR.LT.1 .OR. ( ILVR .AND. LDVR.LT.N ) ) THEN
286: INFO = -13
287: ELSE IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
288: INFO = -15
289: END IF
290: *
291: IF( INFO.EQ.0 ) THEN
292: NB1 = ILAENV( 1, 'ZGEQRF', ' ', N, N, -1, -1 )
293: NB2 = ILAENV( 1, 'ZUNMQR', ' ', N, N, N, -1 )
294: NB3 = ILAENV( 1, 'ZUNGQR', ' ', N, N, N, -1 )
295: NB = MAX( NB1, NB2, NB3 )
296: LOPT = MAX( 2*N, N*( NB+1 ) )
297: WORK( 1 ) = LOPT
298: END IF
299: *
300: IF( INFO.NE.0 ) THEN
301: CALL XERBLA( 'ZGEGV ', -INFO )
302: RETURN
303: ELSE IF( LQUERY ) THEN
304: RETURN
305: END IF
306: *
307: * Quick return if possible
308: *
309: IF( N.EQ.0 )
310: $ RETURN
311: *
312: * Get machine constants
313: *
314: EPS = DLAMCH( 'E' )*DLAMCH( 'B' )
315: SAFMIN = DLAMCH( 'S' )
316: SAFMIN = SAFMIN + SAFMIN
317: SAFMAX = ONE / SAFMIN
318: *
319: * Scale A
320: *
321: ANRM = ZLANGE( 'M', N, N, A, LDA, RWORK )
322: ANRM1 = ANRM
323: ANRM2 = ONE
324: IF( ANRM.LT.ONE ) THEN
325: IF( SAFMAX*ANRM.LT.ONE ) THEN
326: ANRM1 = SAFMIN
327: ANRM2 = SAFMAX*ANRM
328: END IF
329: END IF
330: *
331: IF( ANRM.GT.ZERO ) THEN
332: CALL ZLASCL( 'G', -1, -1, ANRM, ONE, N, N, A, LDA, IINFO )
333: IF( IINFO.NE.0 ) THEN
334: INFO = N + 10
335: RETURN
336: END IF
337: END IF
338: *
339: * Scale B
340: *
341: BNRM = ZLANGE( 'M', N, N, B, LDB, RWORK )
342: BNRM1 = BNRM
343: BNRM2 = ONE
344: IF( BNRM.LT.ONE ) THEN
345: IF( SAFMAX*BNRM.LT.ONE ) THEN
346: BNRM1 = SAFMIN
347: BNRM2 = SAFMAX*BNRM
348: END IF
349: END IF
350: *
351: IF( BNRM.GT.ZERO ) THEN
352: CALL ZLASCL( 'G', -1, -1, BNRM, ONE, N, N, B, LDB, IINFO )
353: IF( IINFO.NE.0 ) THEN
354: INFO = N + 10
355: RETURN
356: END IF
357: END IF
358: *
359: * Permute the matrix to make it more nearly triangular
360: * Also "balance" the matrix.
361: *
362: ILEFT = 1
363: IRIGHT = N + 1
364: IRWORK = IRIGHT + N
365: CALL ZGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, RWORK( ILEFT ),
366: $ RWORK( IRIGHT ), RWORK( IRWORK ), IINFO )
367: IF( IINFO.NE.0 ) THEN
368: INFO = N + 1
369: GO TO 80
370: END IF
371: *
372: * Reduce B to triangular form, and initialize VL and/or VR
373: *
374: IROWS = IHI + 1 - ILO
375: IF( ILV ) THEN
376: ICOLS = N + 1 - ILO
377: ELSE
378: ICOLS = IROWS
379: END IF
380: ITAU = 1
381: IWORK = ITAU + IROWS
382: CALL ZGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
383: $ WORK( IWORK ), LWORK+1-IWORK, IINFO )
384: IF( IINFO.GE.0 )
385: $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
386: IF( IINFO.NE.0 ) THEN
387: INFO = N + 2
388: GO TO 80
389: END IF
390: *
391: CALL ZUNMQR( 'L', 'C', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
392: $ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWORK ),
393: $ LWORK+1-IWORK, IINFO )
394: IF( IINFO.GE.0 )
395: $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
396: IF( IINFO.NE.0 ) THEN
397: INFO = N + 3
398: GO TO 80
399: END IF
400: *
401: IF( ILVL ) THEN
402: CALL ZLASET( 'Full', N, N, CZERO, CONE, VL, LDVL )
403: CALL ZLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
404: $ VL( ILO+1, ILO ), LDVL )
405: CALL ZUNGQR( IROWS, IROWS, IROWS, VL( ILO, ILO ), LDVL,
406: $ WORK( ITAU ), WORK( IWORK ), LWORK+1-IWORK,
407: $ IINFO )
408: IF( IINFO.GE.0 )
409: $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
410: IF( IINFO.NE.0 ) THEN
411: INFO = N + 4
412: GO TO 80
413: END IF
414: END IF
415: *
416: IF( ILVR )
417: $ CALL ZLASET( 'Full', N, N, CZERO, CONE, VR, LDVR )
418: *
419: * Reduce to generalized Hessenberg form
420: *
421: IF( ILV ) THEN
422: *
423: * Eigenvectors requested -- work on whole matrix.
424: *
425: CALL ZGGHRD( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL,
426: $ LDVL, VR, LDVR, IINFO )
427: ELSE
428: CALL ZGGHRD( 'N', 'N', IROWS, 1, IROWS, A( ILO, ILO ), LDA,
429: $ B( ILO, ILO ), LDB, VL, LDVL, VR, LDVR, IINFO )
430: END IF
431: IF( IINFO.NE.0 ) THEN
432: INFO = N + 5
433: GO TO 80
434: END IF
435: *
436: * Perform QZ algorithm
437: *
438: IWORK = ITAU
439: IF( ILV ) THEN
440: CHTEMP = 'S'
441: ELSE
442: CHTEMP = 'E'
443: END IF
444: CALL ZHGEQZ( CHTEMP, JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB,
445: $ ALPHA, BETA, VL, LDVL, VR, LDVR, WORK( IWORK ),
446: $ LWORK+1-IWORK, RWORK( IRWORK ), IINFO )
447: IF( IINFO.GE.0 )
448: $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
449: IF( IINFO.NE.0 ) THEN
450: IF( IINFO.GT.0 .AND. IINFO.LE.N ) THEN
451: INFO = IINFO
452: ELSE IF( IINFO.GT.N .AND. IINFO.LE.2*N ) THEN
453: INFO = IINFO - N
454: ELSE
455: INFO = N + 6
456: END IF
457: GO TO 80
458: END IF
459: *
460: IF( ILV ) THEN
461: *
462: * Compute Eigenvectors
463: *
464: IF( ILVL ) THEN
465: IF( ILVR ) THEN
466: CHTEMP = 'B'
467: ELSE
468: CHTEMP = 'L'
469: END IF
470: ELSE
471: CHTEMP = 'R'
472: END IF
473: *
474: CALL ZTGEVC( CHTEMP, 'B', LDUMMA, N, A, LDA, B, LDB, VL, LDVL,
475: $ VR, LDVR, N, IN, WORK( IWORK ), RWORK( IRWORK ),
476: $ IINFO )
477: IF( IINFO.NE.0 ) THEN
478: INFO = N + 7
479: GO TO 80
480: END IF
481: *
482: * Undo balancing on VL and VR, rescale
483: *
484: IF( ILVL ) THEN
485: CALL ZGGBAK( 'P', 'L', N, ILO, IHI, RWORK( ILEFT ),
486: $ RWORK( IRIGHT ), N, VL, LDVL, IINFO )
487: IF( IINFO.NE.0 ) THEN
488: INFO = N + 8
489: GO TO 80
490: END IF
491: DO 30 JC = 1, N
492: TEMP = ZERO
493: DO 10 JR = 1, N
494: TEMP = MAX( TEMP, ABS1( VL( JR, JC ) ) )
495: 10 CONTINUE
496: IF( TEMP.LT.SAFMIN )
497: $ GO TO 30
498: TEMP = ONE / TEMP
499: DO 20 JR = 1, N
500: VL( JR, JC ) = VL( JR, JC )*TEMP
501: 20 CONTINUE
502: 30 CONTINUE
503: END IF
504: IF( ILVR ) THEN
505: CALL ZGGBAK( 'P', 'R', N, ILO, IHI, RWORK( ILEFT ),
506: $ RWORK( IRIGHT ), N, VR, LDVR, IINFO )
507: IF( IINFO.NE.0 ) THEN
508: INFO = N + 9
509: GO TO 80
510: END IF
511: DO 60 JC = 1, N
512: TEMP = ZERO
513: DO 40 JR = 1, N
514: TEMP = MAX( TEMP, ABS1( VR( JR, JC ) ) )
515: 40 CONTINUE
516: IF( TEMP.LT.SAFMIN )
517: $ GO TO 60
518: TEMP = ONE / TEMP
519: DO 50 JR = 1, N
520: VR( JR, JC ) = VR( JR, JC )*TEMP
521: 50 CONTINUE
522: 60 CONTINUE
523: END IF
524: *
525: * End of eigenvector calculation
526: *
527: END IF
528: *
529: * Undo scaling in alpha, beta
530: *
531: * Note: this does not give the alpha and beta for the unscaled
532: * problem.
533: *
534: * Un-scaling is limited to avoid underflow in alpha and beta
535: * if they are significant.
536: *
537: DO 70 JC = 1, N
538: ABSAR = ABS( DBLE( ALPHA( JC ) ) )
539: ABSAI = ABS( DIMAG( ALPHA( JC ) ) )
540: ABSB = ABS( DBLE( BETA( JC ) ) )
541: SALFAR = ANRM*DBLE( ALPHA( JC ) )
542: SALFAI = ANRM*DIMAG( ALPHA( JC ) )
543: SBETA = BNRM*DBLE( BETA( JC ) )
544: ILIMIT = .FALSE.
545: SCALE = ONE
546: *
547: * Check for significant underflow in imaginary part of ALPHA
548: *
549: IF( ABS( SALFAI ).LT.SAFMIN .AND. ABSAI.GE.
550: $ MAX( SAFMIN, EPS*ABSAR, EPS*ABSB ) ) THEN
551: ILIMIT = .TRUE.
552: SCALE = ( SAFMIN / ANRM1 ) / MAX( SAFMIN, ANRM2*ABSAI )
553: END IF
554: *
555: * Check for significant underflow in real part of ALPHA
556: *
557: IF( ABS( SALFAR ).LT.SAFMIN .AND. ABSAR.GE.
558: $ MAX( SAFMIN, EPS*ABSAI, EPS*ABSB ) ) THEN
559: ILIMIT = .TRUE.
560: SCALE = MAX( SCALE, ( SAFMIN / ANRM1 ) /
561: $ MAX( SAFMIN, ANRM2*ABSAR ) )
562: END IF
563: *
564: * Check for significant underflow in BETA
565: *
566: IF( ABS( SBETA ).LT.SAFMIN .AND. ABSB.GE.
567: $ MAX( SAFMIN, EPS*ABSAR, EPS*ABSAI ) ) THEN
568: ILIMIT = .TRUE.
569: SCALE = MAX( SCALE, ( SAFMIN / BNRM1 ) /
570: $ MAX( SAFMIN, BNRM2*ABSB ) )
571: END IF
572: *
573: * Check for possible overflow when limiting scaling
574: *
575: IF( ILIMIT ) THEN
576: TEMP = ( SCALE*SAFMIN )*MAX( ABS( SALFAR ), ABS( SALFAI ),
577: $ ABS( SBETA ) )
578: IF( TEMP.GT.ONE )
579: $ SCALE = SCALE / TEMP
580: IF( SCALE.LT.ONE )
581: $ ILIMIT = .FALSE.
582: END IF
583: *
584: * Recompute un-scaled ALPHA, BETA if necessary.
585: *
586: IF( ILIMIT ) THEN
587: SALFAR = ( SCALE*DBLE( ALPHA( JC ) ) )*ANRM
588: SALFAI = ( SCALE*DIMAG( ALPHA( JC ) ) )*ANRM
589: SBETA = ( SCALE*BETA( JC ) )*BNRM
590: END IF
591: ALPHA( JC ) = DCMPLX( SALFAR, SALFAI )
592: BETA( JC ) = SBETA
593: 70 CONTINUE
594: *
595: 80 CONTINUE
596: WORK( 1 ) = LWKOPT
597: *
598: RETURN
599: *
600: * End of ZGEGV
601: *
602: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zgegv.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack