Annotation of rpl/lapack/lapack/dgegv.f, revision 1.17
1.8 bertrand 1: *> \brief <b> DGEEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices</b>
2: *
3: * =========== DOCUMENTATION ===========
4: *
1.14 bertrand 5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
1.8 bertrand 7: *
8: *> \htmlonly
1.14 bertrand 9: *> Download DGEGV + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgegv.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgegv.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgegv.f">
1.8 bertrand 15: *> [TXT]</a>
1.14 bertrand 16: *> \endhtmlonly
1.8 bertrand 17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DGEGV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR, ALPHAI,
22: * BETA, VL, LDVL, VR, LDVR, WORK, LWORK, INFO )
1.14 bertrand 23: *
1.8 bertrand 24: * .. Scalar Arguments ..
25: * CHARACTER JOBVL, JOBVR
26: * INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, N
27: * ..
28: * .. Array Arguments ..
29: * DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
30: * $ B( LDB, * ), BETA( * ), VL( LDVL, * ),
31: * $ VR( LDVR, * ), WORK( * )
32: * ..
1.14 bertrand 33: *
1.8 bertrand 34: *
35: *> \par Purpose:
36: * =============
37: *>
38: *> \verbatim
39: *>
40: *> This routine is deprecated and has been replaced by routine DGGEV.
41: *>
42: *> DGEGV computes the eigenvalues and, optionally, the left and/or right
43: *> eigenvectors of a real matrix pair (A,B).
44: *> Given two square matrices A and B,
45: *> the generalized nonsymmetric eigenvalue problem (GNEP) is to find the
46: *> eigenvalues lambda and corresponding (non-zero) eigenvectors x such
47: *> that
48: *>
49: *> A*x = lambda*B*x.
50: *>
51: *> An alternate form is to find the eigenvalues mu and corresponding
52: *> eigenvectors y such that
53: *>
54: *> mu*A*y = B*y.
55: *>
56: *> These two forms are equivalent with mu = 1/lambda and x = y if
57: *> neither lambda nor mu is zero. In order to deal with the case that
58: *> lambda or mu is zero or small, two values alpha and beta are returned
59: *> for each eigenvalue, such that lambda = alpha/beta and
60: *> mu = beta/alpha.
61: *>
62: *> The vectors x and y in the above equations are right eigenvectors of
63: *> the matrix pair (A,B). Vectors u and v satisfying
64: *>
65: *> u**H*A = lambda*u**H*B or mu*v**H*A = v**H*B
66: *>
67: *> are left eigenvectors of (A,B).
68: *>
69: *> Note: this routine performs "full balancing" on A and B
70: *> \endverbatim
71: *
72: * Arguments:
73: * ==========
74: *
75: *> \param[in] JOBVL
76: *> \verbatim
77: *> JOBVL is CHARACTER*1
78: *> = 'N': do not compute the left generalized eigenvectors;
79: *> = 'V': compute the left generalized eigenvectors (returned
80: *> in VL).
81: *> \endverbatim
82: *>
83: *> \param[in] JOBVR
84: *> \verbatim
85: *> JOBVR is CHARACTER*1
86: *> = 'N': do not compute the right generalized eigenvectors;
87: *> = 'V': compute the right generalized eigenvectors (returned
88: *> in VR).
89: *> \endverbatim
90: *>
91: *> \param[in] N
92: *> \verbatim
93: *> N is INTEGER
94: *> The order of the matrices A, B, VL, and VR. N >= 0.
95: *> \endverbatim
96: *>
97: *> \param[in,out] A
98: *> \verbatim
99: *> A is DOUBLE PRECISION array, dimension (LDA, N)
100: *> On entry, the matrix A.
101: *> If JOBVL = 'V' or JOBVR = 'V', then on exit A
102: *> contains the real Schur form of A from the generalized Schur
103: *> factorization of the pair (A,B) after balancing.
104: *> If no eigenvectors were computed, then only the diagonal
105: *> blocks from the Schur form will be correct. See DGGHRD and
106: *> DHGEQZ for details.
107: *> \endverbatim
108: *>
109: *> \param[in] LDA
110: *> \verbatim
111: *> LDA is INTEGER
112: *> The leading dimension of A. LDA >= max(1,N).
113: *> \endverbatim
114: *>
115: *> \param[in,out] B
116: *> \verbatim
117: *> B is DOUBLE PRECISION array, dimension (LDB, N)
118: *> On entry, the matrix B.
119: *> If JOBVL = 'V' or JOBVR = 'V', then on exit B contains the
120: *> upper triangular matrix obtained from B in the generalized
121: *> Schur factorization of the pair (A,B) after balancing.
122: *> If no eigenvectors were computed, then only those elements of
123: *> B corresponding to the diagonal blocks from the Schur form of
124: *> A will be correct. See DGGHRD and DHGEQZ for details.
125: *> \endverbatim
126: *>
127: *> \param[in] LDB
128: *> \verbatim
129: *> LDB is INTEGER
130: *> The leading dimension of B. LDB >= max(1,N).
131: *> \endverbatim
132: *>
133: *> \param[out] ALPHAR
134: *> \verbatim
135: *> ALPHAR is DOUBLE PRECISION array, dimension (N)
136: *> The real parts of each scalar alpha defining an eigenvalue of
137: *> GNEP.
138: *> \endverbatim
139: *>
140: *> \param[out] ALPHAI
141: *> \verbatim
142: *> ALPHAI is DOUBLE PRECISION array, dimension (N)
143: *> The imaginary parts of each scalar alpha defining an
144: *> eigenvalue of GNEP. If ALPHAI(j) is zero, then the j-th
145: *> eigenvalue is real; if positive, then the j-th and
146: *> (j+1)-st eigenvalues are a complex conjugate pair, with
147: *> ALPHAI(j+1) = -ALPHAI(j).
148: *> \endverbatim
149: *>
150: *> \param[out] BETA
151: *> \verbatim
152: *> BETA is DOUBLE PRECISION array, dimension (N)
153: *> The scalars beta that define the eigenvalues of GNEP.
1.14 bertrand 154: *>
1.8 bertrand 155: *> Together, the quantities alpha = (ALPHAR(j),ALPHAI(j)) and
156: *> beta = BETA(j) represent the j-th eigenvalue of the matrix
157: *> pair (A,B), in one of the forms lambda = alpha/beta or
158: *> mu = beta/alpha. Since either lambda or mu may overflow,
159: *> they should not, in general, be computed.
160: *> \endverbatim
161: *>
162: *> \param[out] VL
163: *> \verbatim
164: *> VL is DOUBLE PRECISION array, dimension (LDVL,N)
165: *> If JOBVL = 'V', the left eigenvectors u(j) are stored
166: *> in the columns of VL, in the same order as their eigenvalues.
167: *> If the j-th eigenvalue is real, then u(j) = VL(:,j).
168: *> If the j-th and (j+1)-st eigenvalues form a complex conjugate
169: *> pair, then
170: *> u(j) = VL(:,j) + i*VL(:,j+1)
171: *> and
172: *> u(j+1) = VL(:,j) - i*VL(:,j+1).
173: *>
174: *> Each eigenvector is scaled so that its largest component has
175: *> abs(real part) + abs(imag. part) = 1, except for eigenvectors
176: *> corresponding to an eigenvalue with alpha = beta = 0, which
177: *> are set to zero.
178: *> Not referenced if JOBVL = 'N'.
179: *> \endverbatim
180: *>
181: *> \param[in] LDVL
182: *> \verbatim
183: *> LDVL is INTEGER
184: *> The leading dimension of the matrix VL. LDVL >= 1, and
185: *> if JOBVL = 'V', LDVL >= N.
186: *> \endverbatim
187: *>
188: *> \param[out] VR
189: *> \verbatim
190: *> VR is DOUBLE PRECISION array, dimension (LDVR,N)
191: *> If JOBVR = 'V', the right eigenvectors x(j) are stored
192: *> in the columns of VR, in the same order as their eigenvalues.
193: *> If the j-th eigenvalue is real, then x(j) = VR(:,j).
194: *> If the j-th and (j+1)-st eigenvalues form a complex conjugate
195: *> pair, then
196: *> x(j) = VR(:,j) + i*VR(:,j+1)
197: *> and
198: *> x(j+1) = VR(:,j) - i*VR(:,j+1).
199: *>
200: *> Each eigenvector is scaled so that its largest component has
201: *> abs(real part) + abs(imag. part) = 1, except for eigenvalues
202: *> corresponding to an eigenvalue with alpha = beta = 0, which
203: *> are set to zero.
204: *> Not referenced if JOBVR = 'N'.
205: *> \endverbatim
206: *>
207: *> \param[in] LDVR
208: *> \verbatim
209: *> LDVR is INTEGER
210: *> The leading dimension of the matrix VR. LDVR >= 1, and
211: *> if JOBVR = 'V', LDVR >= N.
212: *> \endverbatim
213: *>
214: *> \param[out] WORK
215: *> \verbatim
216: *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
217: *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
218: *> \endverbatim
219: *>
220: *> \param[in] LWORK
221: *> \verbatim
222: *> LWORK is INTEGER
223: *> The dimension of the array WORK. LWORK >= max(1,8*N).
224: *> For good performance, LWORK must generally be larger.
225: *> To compute the optimal value of LWORK, call ILAENV to get
226: *> blocksizes (for DGEQRF, DORMQR, and DORGQR.) Then compute:
227: *> NB -- MAX of the blocksizes for DGEQRF, DORMQR, and DORGQR;
228: *> The optimal LWORK is:
229: *> 2*N + MAX( 6*N, N*(NB+1) ).
230: *>
231: *> If LWORK = -1, then a workspace query is assumed; the routine
232: *> only calculates the optimal size of the WORK array, returns
233: *> this value as the first entry of the WORK array, and no error
234: *> message related to LWORK is issued by XERBLA.
235: *> \endverbatim
236: *>
237: *> \param[out] INFO
238: *> \verbatim
239: *> INFO is INTEGER
240: *> = 0: successful exit
241: *> < 0: if INFO = -i, the i-th argument had an illegal value.
242: *> = 1,...,N:
243: *> The QZ iteration failed. No eigenvectors have been
244: *> calculated, but ALPHAR(j), ALPHAI(j), and BETA(j)
245: *> should be correct for j=INFO+1,...,N.
246: *> > N: errors that usually indicate LAPACK problems:
247: *> =N+1: error return from DGGBAL
248: *> =N+2: error return from DGEQRF
249: *> =N+3: error return from DORMQR
250: *> =N+4: error return from DORGQR
251: *> =N+5: error return from DGGHRD
252: *> =N+6: error return from DHGEQZ (other than failed
253: *> iteration)
254: *> =N+7: error return from DTGEVC
255: *> =N+8: error return from DGGBAK (computing VL)
256: *> =N+9: error return from DGGBAK (computing VR)
257: *> =N+10: error return from DLASCL (various calls)
258: *> \endverbatim
259: *
260: * Authors:
261: * ========
262: *
1.14 bertrand 263: *> \author Univ. of Tennessee
264: *> \author Univ. of California Berkeley
265: *> \author Univ. of Colorado Denver
266: *> \author NAG Ltd.
1.8 bertrand 267: *
268: *> \ingroup doubleGEeigen
269: *
270: *> \par Further Details:
271: * =====================
272: *>
273: *> \verbatim
274: *>
275: *> Balancing
276: *> ---------
277: *>
278: *> This driver calls DGGBAL to both permute and scale rows and columns
279: *> of A and B. The permutations PL and PR are chosen so that PL*A*PR
280: *> and PL*B*R will be upper triangular except for the diagonal blocks
281: *> A(i:j,i:j) and B(i:j,i:j), with i and j as close together as
282: *> possible. The diagonal scaling matrices DL and DR are chosen so
283: *> that the pair DL*PL*A*PR*DR, DL*PL*B*PR*DR have elements close to
284: *> one (except for the elements that start out zero.)
285: *>
286: *> After the eigenvalues and eigenvectors of the balanced matrices
287: *> have been computed, DGGBAK transforms the eigenvectors back to what
288: *> they would have been (in perfect arithmetic) if they had not been
289: *> balanced.
290: *>
291: *> Contents of A and B on Exit
292: *> -------- -- - --- - -- ----
293: *>
294: *> If any eigenvectors are computed (either JOBVL='V' or JOBVR='V' or
295: *> both), then on exit the arrays A and B will contain the real Schur
296: *> form[*] of the "balanced" versions of A and B. If no eigenvectors
297: *> are computed, then only the diagonal blocks will be correct.
298: *>
299: *> [*] See DHGEQZ, DGEGS, or read the book "Matrix Computations",
300: *> by Golub & van Loan, pub. by Johns Hopkins U. Press.
301: *> \endverbatim
302: *>
303: * =====================================================================
1.1 bertrand 304: SUBROUTINE DGEGV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR, ALPHAI,
305: $ BETA, VL, LDVL, VR, LDVR, WORK, LWORK, INFO )
306: *
1.17 ! bertrand 307: * -- LAPACK driver routine --
1.1 bertrand 308: * -- LAPACK is a software package provided by Univ. of Tennessee, --
309: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
310: *
311: * .. Scalar Arguments ..
312: CHARACTER JOBVL, JOBVR
313: INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, N
314: * ..
315: * .. Array Arguments ..
316: DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
317: $ B( LDB, * ), BETA( * ), VL( LDVL, * ),
318: $ VR( LDVR, * ), WORK( * )
319: * ..
320: *
321: * =====================================================================
322: *
323: * .. Parameters ..
324: DOUBLE PRECISION ZERO, ONE
325: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
326: * ..
327: * .. Local Scalars ..
328: LOGICAL ILIMIT, ILV, ILVL, ILVR, LQUERY
329: CHARACTER CHTEMP
330: INTEGER ICOLS, IHI, IINFO, IJOBVL, IJOBVR, ILEFT, ILO,
331: $ IN, IRIGHT, IROWS, ITAU, IWORK, JC, JR, LOPT,
332: $ LWKMIN, LWKOPT, NB, NB1, NB2, NB3
333: DOUBLE PRECISION ABSAI, ABSAR, ABSB, ANRM, ANRM1, ANRM2, BNRM,
334: $ BNRM1, BNRM2, EPS, ONEPLS, SAFMAX, SAFMIN,
335: $ SALFAI, SALFAR, SBETA, SCALE, TEMP
336: * ..
337: * .. Local Arrays ..
338: LOGICAL LDUMMA( 1 )
339: * ..
340: * .. External Subroutines ..
341: EXTERNAL DGEQRF, DGGBAK, DGGBAL, DGGHRD, DHGEQZ, DLACPY,
342: $ DLASCL, DLASET, DORGQR, DORMQR, DTGEVC, XERBLA
343: * ..
344: * .. External Functions ..
345: LOGICAL LSAME
346: INTEGER ILAENV
347: DOUBLE PRECISION DLAMCH, DLANGE
348: EXTERNAL LSAME, ILAENV, DLAMCH, DLANGE
349: * ..
350: * .. Intrinsic Functions ..
351: INTRINSIC ABS, INT, MAX
352: * ..
353: * .. Executable Statements ..
354: *
355: * Decode the input arguments
356: *
357: IF( LSAME( JOBVL, 'N' ) ) THEN
358: IJOBVL = 1
359: ILVL = .FALSE.
360: ELSE IF( LSAME( JOBVL, 'V' ) ) THEN
361: IJOBVL = 2
362: ILVL = .TRUE.
363: ELSE
364: IJOBVL = -1
365: ILVL = .FALSE.
366: END IF
367: *
368: IF( LSAME( JOBVR, 'N' ) ) THEN
369: IJOBVR = 1
370: ILVR = .FALSE.
371: ELSE IF( LSAME( JOBVR, 'V' ) ) THEN
372: IJOBVR = 2
373: ILVR = .TRUE.
374: ELSE
375: IJOBVR = -1
376: ILVR = .FALSE.
377: END IF
378: ILV = ILVL .OR. ILVR
379: *
380: * Test the input arguments
381: *
382: LWKMIN = MAX( 8*N, 1 )
383: LWKOPT = LWKMIN
384: WORK( 1 ) = LWKOPT
385: LQUERY = ( LWORK.EQ.-1 )
386: INFO = 0
387: IF( IJOBVL.LE.0 ) THEN
388: INFO = -1
389: ELSE IF( IJOBVR.LE.0 ) THEN
390: INFO = -2
391: ELSE IF( N.LT.0 ) THEN
392: INFO = -3
393: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
394: INFO = -5
395: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
396: INFO = -7
397: ELSE IF( LDVL.LT.1 .OR. ( ILVL .AND. LDVL.LT.N ) ) THEN
398: INFO = -12
399: ELSE IF( LDVR.LT.1 .OR. ( ILVR .AND. LDVR.LT.N ) ) THEN
400: INFO = -14
401: ELSE IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
402: INFO = -16
403: END IF
404: *
405: IF( INFO.EQ.0 ) THEN
406: NB1 = ILAENV( 1, 'DGEQRF', ' ', N, N, -1, -1 )
407: NB2 = ILAENV( 1, 'DORMQR', ' ', N, N, N, -1 )
408: NB3 = ILAENV( 1, 'DORGQR', ' ', N, N, N, -1 )
409: NB = MAX( NB1, NB2, NB3 )
410: LOPT = 2*N + MAX( 6*N, N*( NB+1 ) )
411: WORK( 1 ) = LOPT
412: END IF
413: *
414: IF( INFO.NE.0 ) THEN
415: CALL XERBLA( 'DGEGV ', -INFO )
416: RETURN
417: ELSE IF( LQUERY ) THEN
418: RETURN
419: END IF
420: *
421: * Quick return if possible
422: *
423: IF( N.EQ.0 )
424: $ RETURN
425: *
426: * Get machine constants
427: *
428: EPS = DLAMCH( 'E' )*DLAMCH( 'B' )
429: SAFMIN = DLAMCH( 'S' )
430: SAFMIN = SAFMIN + SAFMIN
431: SAFMAX = ONE / SAFMIN
432: ONEPLS = ONE + ( 4*EPS )
433: *
434: * Scale A
435: *
436: ANRM = DLANGE( 'M', N, N, A, LDA, WORK )
437: ANRM1 = ANRM
438: ANRM2 = ONE
439: IF( ANRM.LT.ONE ) THEN
440: IF( SAFMAX*ANRM.LT.ONE ) THEN
441: ANRM1 = SAFMIN
442: ANRM2 = SAFMAX*ANRM
443: END IF
444: END IF
445: *
446: IF( ANRM.GT.ZERO ) THEN
447: CALL DLASCL( 'G', -1, -1, ANRM, ONE, N, N, A, LDA, IINFO )
448: IF( IINFO.NE.0 ) THEN
449: INFO = N + 10
450: RETURN
451: END IF
452: END IF
453: *
454: * Scale B
455: *
456: BNRM = DLANGE( 'M', N, N, B, LDB, WORK )
457: BNRM1 = BNRM
458: BNRM2 = ONE
459: IF( BNRM.LT.ONE ) THEN
460: IF( SAFMAX*BNRM.LT.ONE ) THEN
461: BNRM1 = SAFMIN
462: BNRM2 = SAFMAX*BNRM
463: END IF
464: END IF
465: *
466: IF( BNRM.GT.ZERO ) THEN
467: CALL DLASCL( 'G', -1, -1, BNRM, ONE, N, N, B, LDB, IINFO )
468: IF( IINFO.NE.0 ) THEN
469: INFO = N + 10
470: RETURN
471: END IF
472: END IF
473: *
474: * Permute the matrix to make it more nearly triangular
475: * Workspace layout: (8*N words -- "work" requires 6*N words)
476: * left_permutation, right_permutation, work...
477: *
478: ILEFT = 1
479: IRIGHT = N + 1
480: IWORK = IRIGHT + N
481: CALL DGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, WORK( ILEFT ),
482: $ WORK( IRIGHT ), WORK( IWORK ), IINFO )
483: IF( IINFO.NE.0 ) THEN
484: INFO = N + 1
485: GO TO 120
486: END IF
487: *
488: * Reduce B to triangular form, and initialize VL and/or VR
489: * Workspace layout: ("work..." must have at least N words)
490: * left_permutation, right_permutation, tau, work...
491: *
492: IROWS = IHI + 1 - ILO
493: IF( ILV ) THEN
494: ICOLS = N + 1 - ILO
495: ELSE
496: ICOLS = IROWS
497: END IF
498: ITAU = IWORK
499: IWORK = ITAU + IROWS
500: CALL DGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
501: $ WORK( IWORK ), LWORK+1-IWORK, IINFO )
502: IF( IINFO.GE.0 )
503: $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
504: IF( IINFO.NE.0 ) THEN
505: INFO = N + 2
506: GO TO 120
507: END IF
508: *
509: CALL DORMQR( 'L', 'T', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
510: $ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWORK ),
511: $ LWORK+1-IWORK, IINFO )
512: IF( IINFO.GE.0 )
513: $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
514: IF( IINFO.NE.0 ) THEN
515: INFO = N + 3
516: GO TO 120
517: END IF
518: *
519: IF( ILVL ) THEN
520: CALL DLASET( 'Full', N, N, ZERO, ONE, VL, LDVL )
521: CALL DLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
522: $ VL( ILO+1, ILO ), LDVL )
523: CALL DORGQR( IROWS, IROWS, IROWS, VL( ILO, ILO ), LDVL,
524: $ WORK( ITAU ), WORK( IWORK ), LWORK+1-IWORK,
525: $ IINFO )
526: IF( IINFO.GE.0 )
527: $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
528: IF( IINFO.NE.0 ) THEN
529: INFO = N + 4
530: GO TO 120
531: END IF
532: END IF
533: *
534: IF( ILVR )
535: $ CALL DLASET( 'Full', N, N, ZERO, ONE, VR, LDVR )
536: *
537: * Reduce to generalized Hessenberg form
538: *
539: IF( ILV ) THEN
540: *
541: * Eigenvectors requested -- work on whole matrix.
542: *
543: CALL DGGHRD( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL,
544: $ LDVL, VR, LDVR, IINFO )
545: ELSE
546: CALL DGGHRD( 'N', 'N', IROWS, 1, IROWS, A( ILO, ILO ), LDA,
547: $ B( ILO, ILO ), LDB, VL, LDVL, VR, LDVR, IINFO )
548: END IF
549: IF( IINFO.NE.0 ) THEN
550: INFO = N + 5
551: GO TO 120
552: END IF
553: *
554: * Perform QZ algorithm
555: * Workspace layout: ("work..." must have at least 1 word)
556: * left_permutation, right_permutation, work...
557: *
558: IWORK = ITAU
559: IF( ILV ) THEN
560: CHTEMP = 'S'
561: ELSE
562: CHTEMP = 'E'
563: END IF
564: CALL DHGEQZ( CHTEMP, JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB,
565: $ ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR,
566: $ WORK( IWORK ), LWORK+1-IWORK, IINFO )
567: IF( IINFO.GE.0 )
568: $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
569: IF( IINFO.NE.0 ) THEN
570: IF( IINFO.GT.0 .AND. IINFO.LE.N ) THEN
571: INFO = IINFO
572: ELSE IF( IINFO.GT.N .AND. IINFO.LE.2*N ) THEN
573: INFO = IINFO - N
574: ELSE
575: INFO = N + 6
576: END IF
577: GO TO 120
578: END IF
579: *
580: IF( ILV ) THEN
581: *
582: * Compute Eigenvectors (DTGEVC requires 6*N words of workspace)
583: *
584: IF( ILVL ) THEN
585: IF( ILVR ) THEN
586: CHTEMP = 'B'
587: ELSE
588: CHTEMP = 'L'
589: END IF
590: ELSE
591: CHTEMP = 'R'
592: END IF
593: *
594: CALL DTGEVC( CHTEMP, 'B', LDUMMA, N, A, LDA, B, LDB, VL, LDVL,
595: $ VR, LDVR, N, IN, WORK( IWORK ), IINFO )
596: IF( IINFO.NE.0 ) THEN
597: INFO = N + 7
598: GO TO 120
599: END IF
600: *
601: * Undo balancing on VL and VR, rescale
602: *
603: IF( ILVL ) THEN
604: CALL DGGBAK( 'P', 'L', N, ILO, IHI, WORK( ILEFT ),
605: $ WORK( IRIGHT ), N, VL, LDVL, IINFO )
606: IF( IINFO.NE.0 ) THEN
607: INFO = N + 8
608: GO TO 120
609: END IF
610: DO 50 JC = 1, N
611: IF( ALPHAI( JC ).LT.ZERO )
612: $ GO TO 50
613: TEMP = ZERO
614: IF( ALPHAI( JC ).EQ.ZERO ) THEN
615: DO 10 JR = 1, N
616: TEMP = MAX( TEMP, ABS( VL( JR, JC ) ) )
617: 10 CONTINUE
618: ELSE
619: DO 20 JR = 1, N
620: TEMP = MAX( TEMP, ABS( VL( JR, JC ) )+
621: $ ABS( VL( JR, JC+1 ) ) )
622: 20 CONTINUE
623: END IF
624: IF( TEMP.LT.SAFMIN )
625: $ GO TO 50
626: TEMP = ONE / TEMP
627: IF( ALPHAI( JC ).EQ.ZERO ) THEN
628: DO 30 JR = 1, N
629: VL( JR, JC ) = VL( JR, JC )*TEMP
630: 30 CONTINUE
631: ELSE
632: DO 40 JR = 1, N
633: VL( JR, JC ) = VL( JR, JC )*TEMP
634: VL( JR, JC+1 ) = VL( JR, JC+1 )*TEMP
635: 40 CONTINUE
636: END IF
637: 50 CONTINUE
638: END IF
639: IF( ILVR ) THEN
640: CALL DGGBAK( 'P', 'R', N, ILO, IHI, WORK( ILEFT ),
641: $ WORK( IRIGHT ), N, VR, LDVR, IINFO )
642: IF( IINFO.NE.0 ) THEN
643: INFO = N + 9
644: GO TO 120
645: END IF
646: DO 100 JC = 1, N
647: IF( ALPHAI( JC ).LT.ZERO )
648: $ GO TO 100
649: TEMP = ZERO
650: IF( ALPHAI( JC ).EQ.ZERO ) THEN
651: DO 60 JR = 1, N
652: TEMP = MAX( TEMP, ABS( VR( JR, JC ) ) )
653: 60 CONTINUE
654: ELSE
655: DO 70 JR = 1, N
656: TEMP = MAX( TEMP, ABS( VR( JR, JC ) )+
657: $ ABS( VR( JR, JC+1 ) ) )
658: 70 CONTINUE
659: END IF
660: IF( TEMP.LT.SAFMIN )
661: $ GO TO 100
662: TEMP = ONE / TEMP
663: IF( ALPHAI( JC ).EQ.ZERO ) THEN
664: DO 80 JR = 1, N
665: VR( JR, JC ) = VR( JR, JC )*TEMP
666: 80 CONTINUE
667: ELSE
668: DO 90 JR = 1, N
669: VR( JR, JC ) = VR( JR, JC )*TEMP
670: VR( JR, JC+1 ) = VR( JR, JC+1 )*TEMP
671: 90 CONTINUE
672: END IF
673: 100 CONTINUE
674: END IF
675: *
676: * End of eigenvector calculation
677: *
678: END IF
679: *
680: * Undo scaling in alpha, beta
681: *
682: * Note: this does not give the alpha and beta for the unscaled
683: * problem.
684: *
685: * Un-scaling is limited to avoid underflow in alpha and beta
686: * if they are significant.
687: *
688: DO 110 JC = 1, N
689: ABSAR = ABS( ALPHAR( JC ) )
690: ABSAI = ABS( ALPHAI( JC ) )
691: ABSB = ABS( BETA( JC ) )
692: SALFAR = ANRM*ALPHAR( JC )
693: SALFAI = ANRM*ALPHAI( JC )
694: SBETA = BNRM*BETA( JC )
695: ILIMIT = .FALSE.
696: SCALE = ONE
697: *
698: * Check for significant underflow in ALPHAI
699: *
700: IF( ABS( SALFAI ).LT.SAFMIN .AND. ABSAI.GE.
701: $ MAX( SAFMIN, EPS*ABSAR, EPS*ABSB ) ) THEN
702: ILIMIT = .TRUE.
703: SCALE = ( ONEPLS*SAFMIN / ANRM1 ) /
704: $ MAX( ONEPLS*SAFMIN, ANRM2*ABSAI )
705: *
706: ELSE IF( SALFAI.EQ.ZERO ) THEN
707: *
708: * If insignificant underflow in ALPHAI, then make the
709: * conjugate eigenvalue real.
710: *
711: IF( ALPHAI( JC ).LT.ZERO .AND. JC.GT.1 ) THEN
712: ALPHAI( JC-1 ) = ZERO
713: ELSE IF( ALPHAI( JC ).GT.ZERO .AND. JC.LT.N ) THEN
714: ALPHAI( JC+1 ) = ZERO
715: END IF
716: END IF
717: *
718: * Check for significant underflow in ALPHAR
719: *
720: IF( ABS( SALFAR ).LT.SAFMIN .AND. ABSAR.GE.
721: $ MAX( SAFMIN, EPS*ABSAI, EPS*ABSB ) ) THEN
722: ILIMIT = .TRUE.
723: SCALE = MAX( SCALE, ( ONEPLS*SAFMIN / ANRM1 ) /
724: $ MAX( ONEPLS*SAFMIN, ANRM2*ABSAR ) )
725: END IF
726: *
727: * Check for significant underflow in BETA
728: *
729: IF( ABS( SBETA ).LT.SAFMIN .AND. ABSB.GE.
730: $ MAX( SAFMIN, EPS*ABSAR, EPS*ABSAI ) ) THEN
731: ILIMIT = .TRUE.
732: SCALE = MAX( SCALE, ( ONEPLS*SAFMIN / BNRM1 ) /
733: $ MAX( ONEPLS*SAFMIN, BNRM2*ABSB ) )
734: END IF
735: *
736: * Check for possible overflow when limiting scaling
737: *
738: IF( ILIMIT ) THEN
739: TEMP = ( SCALE*SAFMIN )*MAX( ABS( SALFAR ), ABS( SALFAI ),
740: $ ABS( SBETA ) )
741: IF( TEMP.GT.ONE )
742: $ SCALE = SCALE / TEMP
743: IF( SCALE.LT.ONE )
744: $ ILIMIT = .FALSE.
745: END IF
746: *
747: * Recompute un-scaled ALPHAR, ALPHAI, BETA if necessary.
748: *
749: IF( ILIMIT ) THEN
750: SALFAR = ( SCALE*ALPHAR( JC ) )*ANRM
751: SALFAI = ( SCALE*ALPHAI( JC ) )*ANRM
752: SBETA = ( SCALE*BETA( JC ) )*BNRM
753: END IF
754: ALPHAR( JC ) = SALFAR
755: ALPHAI( JC ) = SALFAI
756: BETA( JC ) = SBETA
757: 110 CONTINUE
758: *
759: 120 CONTINUE
760: WORK( 1 ) = LWKOPT
761: *
762: RETURN
763: *
764: * End of DGEGV
765: *
766: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dgegv.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack