OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
dggevx.f
Go to the documentation of this file.
1*> \brief <b> DGGEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices</b>
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DGGEVX + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dggevx.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dggevx.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dggevx.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DGGEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, B, LDB,
22* ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR, ILO,
23* IHI, LSCALE, RSCALE, ABNRM, BBNRM, RCONDE,
24* RCONDV, WORK, LWORK, IWORK, BWORK, INFO )
25*
26* .. Scalar Arguments ..
27* CHARACTER BALANC, JOBVL, JOBVR, SENSE
28* INTEGER IHI, ILO, INFO, LDA, LDB, LDVL, LDVR, LWORK, N
29* DOUBLE PRECISION ABNRM, BBNRM
30* ..
31* .. Array Arguments ..
32* LOGICAL BWORK( * )
33* INTEGER IWORK( * )
34* DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
35* $ B( LDB, * ), BETA( * ), LSCALE( * ),
36* $ RCONDE( * ), RCONDV( * ), RSCALE( * ),
37* $ VL( LDVL, * ), VR( LDVR, * ), WORK( * )
38* ..
39*
40*
41*> \par Purpose:
42* =============
43*>
44*> \verbatim
45*>
46*> DGGEVX computes for a pair of N-by-N real nonsymmetric matrices (A,B)
47*> the generalized eigenvalues, and optionally, the left and/or right
48*> generalized eigenvectors.
49*>
50*> Optionally also, it computes a balancing transformation to improve
51*> the conditioning of the eigenvalues and eigenvectors (ILO, IHI,
52*> LSCALE, RSCALE, ABNRM, and BBNRM), reciprocal condition numbers for
53*> the eigenvalues (RCONDE), and reciprocal condition numbers for the
54*> right eigenvectors (RCONDV).
55*>
56*> A generalized eigenvalue for a pair of matrices (A,B) is a scalar
57*> lambda or a ratio alpha/beta = lambda, such that A - lambda*B is
58*> singular. It is usually represented as the pair (alpha,beta), as
59*> there is a reasonable interpretation for beta=0, and even for both
60*> being zero.
61*>
62*> The right eigenvector v(j) corresponding to the eigenvalue lambda(j)
63*> of (A,B) satisfies
64*>
65*> A * v(j) = lambda(j) * B * v(j) .
66*>
67*> The left eigenvector u(j) corresponding to the eigenvalue lambda(j)
68*> of (A,B) satisfies
69*>
70*> u(j)**H * A = lambda(j) * u(j)**H * B.
71*>
72*> where u(j)**H is the conjugate-transpose of u(j).
73*>
74*> \endverbatim
75*
76* Arguments:
77* ==========
78*
79*> \param[in] BALANC
80*> \verbatim
81*> BALANC is CHARACTER*1
82*> Specifies the balance option to be performed.
83*> = 'N': do not diagonally scale or permute;
84*> = 'P': permute only;
85*> = 'S': scale only;
86*> = 'B': both permute and scale.
87*> Computed reciprocal condition numbers will be for the
88*> matrices after permuting and/or balancing. Permuting does
89*> not change condition numbers (in exact arithmetic), but
90*> balancing does.
91*> \endverbatim
92*>
93*> \param[in] JOBVL
94*> \verbatim
95*> JOBVL is CHARACTER*1
96*> = 'N': do not compute the left generalized eigenvectors;
97*> = 'V': compute the left generalized eigenvectors.
98*> \endverbatim
99*>
100*> \param[in] JOBVR
101*> \verbatim
102*> JOBVR is CHARACTER*1
103*> = 'N': do not compute the right generalized eigenvectors;
104*> = 'V': compute the right generalized eigenvectors.
105*> \endverbatim
106*>
107*> \param[in] SENSE
108*> \verbatim
109*> SENSE is CHARACTER*1
110*> Determines which reciprocal condition numbers are computed.
111*> = 'N': none are computed;
112*> = 'E': computed for eigenvalues only;
113*> = 'V': computed for eigenvectors only;
114*> = 'B': computed for eigenvalues and eigenvectors.
115*> \endverbatim
116*>
117*> \param[in] N
118*> \verbatim
119*> N is INTEGER
120*> The order of the matrices A, B, VL, and VR. N >= 0.
121*> \endverbatim
122*>
123*> \param[in,out] A
124*> \verbatim
125*> A is DOUBLE PRECISION array, dimension (LDA, N)
126*> On entry, the matrix A in the pair (A,B).
127*> On exit, A has been overwritten. If JOBVL='V' or JOBVR='V'
128*> or both, then A contains the first part of the real Schur
129*> form of the "balanced" versions of the input A and B.
130*> \endverbatim
131*>
132*> \param[in] LDA
133*> \verbatim
134*> LDA is INTEGER
135*> The leading dimension of A. LDA >= max(1,N).
136*> \endverbatim
137*>
138*> \param[in,out] B
139*> \verbatim
140*> B is DOUBLE PRECISION array, dimension (LDB, N)
141*> On entry, the matrix B in the pair (A,B).
142*> On exit, B has been overwritten. If JOBVL='V' or JOBVR='V'
143*> or both, then B contains the second part of the real Schur
144*> form of the "balanced" versions of the input A and B.
145*> \endverbatim
146*>
147*> \param[in] LDB
148*> \verbatim
149*> LDB is INTEGER
150*> The leading dimension of B. LDB >= max(1,N).
151*> \endverbatim
152*>
153*> \param[out] ALPHAR
154*> \verbatim
155*> ALPHAR is DOUBLE PRECISION array, dimension (N)
156*> \endverbatim
157*>
158*> \param[out] ALPHAI
159*> \verbatim
160*> ALPHAI is DOUBLE PRECISION array, dimension (N)
161*> \endverbatim
162*>
163*> \param[out] BETA
164*> \verbatim
165*> BETA is DOUBLE PRECISION array, dimension (N)
166*> On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will
167*> be the generalized eigenvalues. If ALPHAI(j) is zero, then
168*> the j-th eigenvalue is real; if positive, then the j-th and
169*> (j+1)-st eigenvalues are a complex conjugate pair, with
170*> ALPHAI(j+1) negative.
171*>
172*> Note: the quotients ALPHAR(j)/BETA(j) and ALPHAI(j)/BETA(j)
173*> may easily over- or underflow, and BETA(j) may even be zero.
174*> Thus, the user should avoid naively computing the ratio
175*> ALPHA/BETA. However, ALPHAR and ALPHAI will be always less
176*> than and usually comparable with norm(A) in magnitude, and
177*> BETA always less than and usually comparable with norm(B).
178*> \endverbatim
179*>
180*> \param[out] VL
181*> \verbatim
182*> VL is DOUBLE PRECISION array, dimension (LDVL,N)
183*> If JOBVL = 'V', the left eigenvectors u(j) are stored one
184*> after another in the columns of VL, in the same order as
185*> their eigenvalues. If the j-th eigenvalue is real, then
186*> u(j) = VL(:,j), the j-th column of VL. If the j-th and
187*> (j+1)-th eigenvalues form a complex conjugate pair, then
188*> u(j) = VL(:,j)+i*VL(:,j+1) and u(j+1) = VL(:,j)-i*VL(:,j+1).
189*> Each eigenvector will be scaled so the largest component have
190*> abs(real part) + abs(imag. part) = 1.
191*> Not referenced if JOBVL = 'N'.
192*> \endverbatim
193*>
194*> \param[in] LDVL
195*> \verbatim
196*> LDVL is INTEGER
197*> The leading dimension of the matrix VL. LDVL >= 1, and
198*> if JOBVL = 'V', LDVL >= N.
199*> \endverbatim
200*>
201*> \param[out] VR
202*> \verbatim
203*> VR is DOUBLE PRECISION array, dimension (LDVR,N)
204*> If JOBVR = 'V', the right eigenvectors v(j) are stored one
205*> after another in the columns of VR, in the same order as
206*> their eigenvalues. If the j-th eigenvalue is real, then
207*> v(j) = VR(:,j), the j-th column of VR. If the j-th and
208*> (j+1)-th eigenvalues form a complex conjugate pair, then
209*> v(j) = VR(:,j)+i*VR(:,j+1) and v(j+1) = VR(:,j)-i*VR(:,j+1).
210*> Each eigenvector will be scaled so the largest component have
211*> abs(real part) + abs(imag. part) = 1.
212*> Not referenced if JOBVR = 'N'.
213*> \endverbatim
214*>
215*> \param[in] LDVR
216*> \verbatim
217*> LDVR is INTEGER
218*> The leading dimension of the matrix VR. LDVR >= 1, and
219*> if JOBVR = 'V', LDVR >= N.
220*> \endverbatim
221*>
222*> \param[out] ILO
223*> \verbatim
224*> ILO is INTEGER
225*> \endverbatim
226*>
227*> \param[out] IHI
228*> \verbatim
229*> IHI is INTEGER
230*> ILO and IHI are integer values such that on exit
231*> A(i,j) = 0 and B(i,j) = 0 if i > j and
232*> j = 1,...,ILO-1 or i = IHI+1,...,N.
233*> If BALANC = 'N' or 'S', ILO = 1 and IHI = N.
234*> \endverbatim
235*>
236*> \param[out] LSCALE
237*> \verbatim
238*> LSCALE is DOUBLE PRECISION array, dimension (N)
239*> Details of the permutations and scaling factors applied
240*> to the left side of A and B. If PL(j) is the index of the
241*> row interchanged with row j, and DL(j) is the scaling
242*> factor applied to row j, then
243*> LSCALE(j) = PL(j) for j = 1,...,ILO-1
244*> = DL(j) for j = ILO,...,IHI
245*> = PL(j) for j = IHI+1,...,N.
246*> The order in which the interchanges are made is N to IHI+1,
247*> then 1 to ILO-1.
248*> \endverbatim
249*>
250*> \param[out] RSCALE
251*> \verbatim
252*> RSCALE is DOUBLE PRECISION array, dimension (N)
253*> Details of the permutations and scaling factors applied
254*> to the right side of A and B. If PR(j) is the index of the
255*> column interchanged with column j, and DR(j) is the scaling
256*> factor applied to column j, then
257*> RSCALE(j) = PR(j) for j = 1,...,ILO-1
258*> = DR(j) for j = ILO,...,IHI
259*> = PR(j) for j = IHI+1,...,N
260*> The order in which the interchanges are made is N to IHI+1,
261*> then 1 to ILO-1.
262*> \endverbatim
263*>
264*> \param[out] ABNRM
265*> \verbatim
266*> ABNRM is DOUBLE PRECISION
267*> The one-norm of the balanced matrix A.
268*> \endverbatim
269*>
270*> \param[out] BBNRM
271*> \verbatim
272*> BBNRM is DOUBLE PRECISION
273*> The one-norm of the balanced matrix B.
274*> \endverbatim
275*>
276*> \param[out] RCONDE
277*> \verbatim
278*> RCONDE is DOUBLE PRECISION array, dimension (N)
279*> If SENSE = 'E' or 'B', the reciprocal condition numbers of
280*> the eigenvalues, stored in consecutive elements of the array.
281*> For a complex conjugate pair of eigenvalues two consecutive
282*> elements of RCONDE are set to the same value. Thus RCONDE(j),
283*> RCONDV(j), and the j-th columns of VL and VR all correspond
284*> to the j-th eigenpair.
285*> If SENSE = 'N or 'V', RCONDE is not referenced.
286*> \endverbatim
287*>
288*> \param[out] RCONDV
289*> \verbatim
290*> RCONDV is DOUBLE PRECISION array, dimension (N)
291*> If SENSE = 'V' or 'B', the estimated reciprocal condition
292*> numbers of the eigenvectors, stored in consecutive elements
293*> of the array. For a complex eigenvector two consecutive
294*> elements of RCONDV are set to the same value. If the
295*> eigenvalues cannot be reordered to compute RCONDV(j),
296*> RCONDV(j) is set to 0; this can only occur when the true
297*> value would be very small anyway.
298*> If SENSE = 'N' or 'E', RCONDV is not referenced.
299*> \endverbatim
300*>
301*> \param[out] WORK
302*> \verbatim
303*> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
304*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
305*> \endverbatim
306*>
307*> \param[in] LWORK
308*> \verbatim
309*> LWORK is INTEGER
310*> The dimension of the array WORK. LWORK >= max(1,2*N).
311*> If BALANC = 'S' or 'B', or JOBVL = 'V', or JOBVR = 'V',
312*> LWORK >= max(1,6*N).
313*> If SENSE = 'E' or 'B', LWORK >= max(1,10*N).
314*> If SENSE = 'V' or 'B', LWORK >= 2*N*N+8*N+16.
315*>
316*> If LWORK = -1, then a workspace query is assumed; the routine
317*> only calculates the optimal size of the WORK array, returns
318*> this value as the first entry of the WORK array, and no error
319*> message related to LWORK is issued by XERBLA.
320*> \endverbatim
321*>
322*> \param[out] IWORK
323*> \verbatim
324*> IWORK is INTEGER array, dimension (N+6)
325*> If SENSE = 'E', IWORK is not referenced.
326*> \endverbatim
327*>
328*> \param[out] BWORK
329*> \verbatim
330*> BWORK is LOGICAL array, dimension (N)
331*> If SENSE = 'N', BWORK is not referenced.
332*> \endverbatim
333*>
334*> \param[out] INFO
335*> \verbatim
336*> INFO is INTEGER
337*> = 0: successful exit
338*> < 0: if INFO = -i, the i-th argument had an illegal value.
339*> = 1,...,N:
340*> The QZ iteration failed. No eigenvectors have been
341*> calculated, but ALPHAR(j), ALPHAI(j), and BETA(j)
342*> should be correct for j=INFO+1,...,N.
343*> > N: =N+1: other than QZ iteration failed in DHGEQZ.
344*> =N+2: error return from DTGEVC.
345*> \endverbatim
346*
347* Authors:
348* ========
349*
350*> \author Univ. of Tennessee
351*> \author Univ. of California Berkeley
352*> \author Univ. of Colorado Denver
353*> \author NAG Ltd.
354*
355*> \ingroup doubleGEeigen
356*
357*> \par Further Details:
358* =====================
359*>
360*> \verbatim
361*>
362*> Balancing a matrix pair (A,B) includes, first, permuting rows and
363*> columns to isolate eigenvalues, second, applying diagonal similarity
364*> transformation to the rows and columns to make the rows and columns
365*> as close in norm as possible. The computed reciprocal condition
366*> numbers correspond to the balanced matrix. Permuting rows and columns
367*> will not change the condition numbers (in exact arithmetic) but
368*> diagonal scaling will. For further explanation of balancing, see
369*> section 4.11.1.2 of LAPACK Users' Guide.
370*>
371*> An approximate error bound on the chordal distance between the i-th
372*> computed generalized eigenvalue w and the corresponding exact
373*> eigenvalue lambda is
374*>
375*> chord(w, lambda) <= EPS * norm(ABNRM, BBNRM) / RCONDE(I)
376*>
377*> An approximate error bound for the angle between the i-th computed
378*> eigenvector VL(i) or VR(i) is given by
379*>
380*> EPS * norm(ABNRM, BBNRM) / DIF(i).
381*>
382*> For further explanation of the reciprocal condition numbers RCONDE
383*> and RCONDV, see section 4.11 of LAPACK User's Guide.
384*> \endverbatim
385*>
386* =====================================================================
387 SUBROUTINE dggevx( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, B, LDB,
388 $ ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR, ILO,
389 $ IHI, LSCALE, RSCALE, ABNRM, BBNRM, RCONDE,
390 $ RCONDV, WORK, LWORK, IWORK, BWORK, INFO )
391*
392* -- LAPACK driver routine --
393* -- LAPACK is a software package provided by Univ. of Tennessee, --
394* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
395*
396* .. Scalar Arguments ..
397 CHARACTER BALANC, JOBVL, JOBVR, SENSE
398 INTEGER IHI, ILO, INFO, LDA, LDB, LDVL, LDVR, LWORK, N
399 DOUBLE PRECISION ABNRM, BBNRM
400* ..
401* .. Array Arguments ..
402 LOGICAL BWORK( * )
403 INTEGER IWORK( * )
404 DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
405 $ b( ldb, * ), beta( * ), lscale( * ),
406 $ rconde( * ), rcondv( * ), rscale( * ),
407 $ vl( ldvl, * ), vr( ldvr, * ), work( * )
408* ..
409*
410* =====================================================================
411*
412* .. Parameters ..
413 DOUBLE PRECISION ZERO, ONE
414 PARAMETER ( ZERO = 0.0d+0, one = 1.0d+0 )
415* ..
416* .. Local Scalars ..
417 LOGICAL ILASCL, ILBSCL, ILV, ILVL, ILVR, LQUERY, NOSCL,
418 $ PAIR, WANTSB, WANTSE, WANTSN, WANTSV
419 CHARACTER CHTEMP
420 INTEGER I, ICOLS, IERR, IJOBVL, IJOBVR, IN, IROWS,
421 $ ITAU, IWRK, IWRK1, J, JC, JR, M, MAXWRK,
422 $ minwrk, mm
423 DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
424 $ SMLNUM, TEMP
425* ..
426* .. Local Arrays ..
427 LOGICAL LDUMMA( 1 )
428* ..
429* .. External Subroutines ..
430 EXTERNAL dgeqrf, dggbak, dggbal, dgghrd, dhgeqz, dlabad,
431 $ dlacpy, dlascl, dlaset, dorgqr, dormqr, dtgevc,
432 $ dtgsna, xerbla
433* ..
434* .. External Functions ..
435 LOGICAL LSAME
436 INTEGER ILAENV
437 DOUBLE PRECISION DLAMCH, DLANGE
438 EXTERNAL lsame, ilaenv, dlamch, dlange
439* ..
440* .. Intrinsic Functions ..
441 INTRINSIC abs, max, sqrt
442* ..
443* .. Executable Statements ..
444*
445* Decode the input arguments
446*
447 IF( lsame( jobvl, 'N' ) ) THEN
448 ijobvl = 1
449 ilvl = .false.
450 ELSE IF( lsame( jobvl, 'V' ) ) THEN
451 ijobvl = 2
452 ilvl = .true.
453 ELSE
454 ijobvl = -1
455 ilvl = .false.
456 END IF
457*
458 IF( lsame( jobvr, 'N' ) ) THEN
459 ijobvr = 1
460 ilvr = .false.
461 ELSE IF( lsame( jobvr, 'V' ) ) THEN
462 ijobvr = 2
463 ilvr = .true.
464 ELSE
465 ijobvr = -1
466 ilvr = .false.
467 END IF
468 ilv = ilvl .OR. ilvr
469*
470 noscl = lsame( balanc, 'N' ) .OR. lsame( balanc, 'P' )
471 wantsn = lsame( sense, 'N' )
472 wantse = lsame( sense, 'E' )
473 wantsv = lsame( sense, 'V' )
474 wantsb = lsame( sense, 'B' )
475*
476* Test the input arguments
477*
478 info = 0
479 lquery = ( lwork.EQ.-1 )
480 IF( .NOT.( lsame( balanc, 'N' ) .OR. lsame( balanc,
481 $ 'S' ) .OR. lsame( balanc, 'P' ) .OR. lsame( balanc, 'B' ) ) )
482 $ THEN
483 info = -1
484 ELSE IF( ijobvl.LE.0 ) THEN
485 info = -2
486 ELSE IF( ijobvr.LE.0 ) THEN
487 info = -3
488 ELSE IF( .NOT.( wantsn .OR. wantse .OR. wantsb .OR. wantsv ) )
489 $ THEN
490 info = -4
491 ELSE IF( n.LT.0 ) THEN
492 info = -5
493 ELSE IF( lda.LT.max( 1, n ) ) THEN
494 info = -7
495 ELSE IF( ldb.LT.max( 1, n ) ) THEN
496 info = -9
497 ELSE IF( ldvl.LT.1 .OR. ( ilvl .AND. ldvl.LT.n ) ) THEN
498 info = -14
499 ELSE IF( ldvr.LT.1 .OR. ( ilvr .AND. ldvr.LT.n ) ) THEN
500 info = -16
501 END IF
502*
503* Compute workspace
504* (Note: Comments in the code beginning "Workspace:" describe the
505* minimal amount of workspace needed at that point in the code,
506* as well as the preferred amount for good performance.
507* NB refers to the optimal block size for the immediately
508* following subroutine, as returned by ILAENV. The workspace is
509* computed assuming ILO = 1 and IHI = N, the worst case.)
510*
511 IF( info.EQ.0 ) THEN
512 IF( n.EQ.0 ) THEN
513 minwrk = 1
514 maxwrk = 1
515 ELSE
516 IF( noscl .AND. .NOT.ilv ) THEN
517 minwrk = 2*n
518 ELSE
519 minwrk = 6*n
520 END IF
521 IF( wantse .OR. wantsb ) THEN
522 minwrk = 10*n
523 END IF
524 IF( wantsv .OR. wantsb ) THEN
525 minwrk = max( minwrk, 2*n*( n + 4 ) + 16 )
526 END IF
527 maxwrk = minwrk
528 maxwrk = max( maxwrk,
529 $ n + n*ilaenv( 1, 'DGEQRF', ' ', n, 1, n, 0 ) )
530 maxwrk = max( maxwrk,
531 $ n + n*ilaenv( 1, 'DORMQR', ' ', n, 1, n, 0 ) )
532 IF( ilvl ) THEN
533 maxwrk = max( maxwrk, n +
534 $ n*ilaenv( 1, 'DORGQR', ' ', n, 1, n, 0 ) )
535 END IF
536 END IF
537 work( 1 ) = maxwrk
538*
539 IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
540 info = -26
541 END IF
542 END IF
543*
544 IF( info.NE.0 ) THEN
545 CALL xerbla( 'DGGEVX', -info )
546 RETURN
547 ELSE IF( lquery ) THEN
548 RETURN
549 END IF
550*
551* Quick return if possible
552*
553 IF( n.EQ.0 )
554 $ RETURN
555*
556*
557* Get machine constants
558*
559 eps = dlamch( 'P' )
560 smlnum = dlamch( 'S' )
561 bignum = one / smlnum
562 CALL dlabad( smlnum, bignum )
563 smlnum = sqrt( smlnum ) / eps
564 bignum = one / smlnum
565*
566* Scale A if max element outside range [SMLNUM,BIGNUM]
567*
568 anrm = dlange( 'M', n, n, a, lda, work )
569 ilascl = .false.
570 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
571 anrmto = smlnum
572 ilascl = .true.
573 ELSE IF( anrm.GT.bignum ) THEN
574 anrmto = bignum
575 ilascl = .true.
576 END IF
577 IF( ilascl )
578 $ CALL dlascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
579*
580* Scale B if max element outside range [SMLNUM,BIGNUM]
581*
582 bnrm = dlange( 'M', n, n, b, ldb, work )
583 ilbscl = .false.
584 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
585 bnrmto = smlnum
586 ilbscl = .true.
587 ELSE IF( bnrm.GT.bignum ) THEN
588 bnrmto = bignum
589 ilbscl = .true.
590 END IF
591 IF( ilbscl )
592 $ CALL dlascl( 'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
593*
594* Permute and/or balance the matrix pair (A,B)
595* (Workspace: need 6*N if BALANC = 'S' or 'B', 1 otherwise)
596*
597 CALL dggbal( balanc, n, a, lda, b, ldb, ilo, ihi, lscale, rscale,
598 $ work, ierr )
599*
600* Compute ABNRM and BBNRM
601*
602 abnrm = dlange( '1', n, n, a, lda, work( 1 ) )
603 IF( ilascl ) THEN
604 work( 1 ) = abnrm
605 CALL dlascl( 'G', 0, 0, anrmto, anrm, 1, 1, work( 1 ), 1,
606 $ ierr )
607 abnrm = work( 1 )
608 END IF
609*
610 bbnrm = dlange( '1', n, n, b, ldb, work( 1 ) )
611 IF( ilbscl ) THEN
612 work( 1 ) = bbnrm
613 CALL dlascl( 'G', 0, 0, bnrmto, bnrm, 1, 1, work( 1 ), 1,
614 $ ierr )
615 bbnrm = work( 1 )
616 END IF
617*
618* Reduce B to triangular form (QR decomposition of B)
619* (Workspace: need N, prefer N*NB )
620*
621 irows = ihi + 1 - ilo
622 IF( ilv .OR. .NOT.wantsn ) THEN
623 icols = n + 1 - ilo
624 ELSE
625 icols = irows
626 END IF
627 itau = 1
628 iwrk = itau + irows
629 CALL dgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
630 $ work( iwrk ), lwork+1-iwrk, ierr )
631*
632* Apply the orthogonal transformation to A
633* (Workspace: need N, prefer N*NB)
634*
635 CALL dormqr( 'L', 'T', irows, icols, irows, b( ilo, ilo ), ldb,
636 $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
637 $ lwork+1-iwrk, ierr )
638*
639* Initialize VL and/or VR
640* (Workspace: need N, prefer N*NB)
641*
642 IF( ilvl ) THEN
643 CALL dlaset( 'Full', n, n, zero, one, vl, ldvl )
644 IF( irows.GT.1 ) THEN
645 CALL dlacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
646 $ vl( ilo+1, ilo ), ldvl )
647 END IF
648 CALL dorgqr( irows, irows, irows, vl( ilo, ilo ), ldvl,
649 $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
650 END IF
651*
652 IF( ilvr )
653 $ CALL dlaset( 'Full', n, n, zero, one, vr, ldvr )
654*
655* Reduce to generalized Hessenberg form
656* (Workspace: none needed)
657*
658 IF( ilv .OR. .NOT.wantsn ) THEN
659*
660* Eigenvectors requested -- work on whole matrix.
661*
662 CALL dgghrd( jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb, vl,
663 $ ldvl, vr, ldvr, ierr )
664 ELSE
665 CALL dgghrd( 'N', 'N', irows, 1, irows, a( ilo, ilo ), lda,
666 $ b( ilo, ilo ), ldb, vl, ldvl, vr, ldvr, ierr )
667 END IF
668*
669* Perform QZ algorithm (Compute eigenvalues, and optionally, the
670* Schur forms and Schur vectors)
671* (Workspace: need N)
672*
673 IF( ilv .OR. .NOT.wantsn ) THEN
674 chtemp = 'S'
675 ELSE
676 chtemp = 'E'
677 END IF
678*
679 CALL dhgeqz( chtemp, jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb,
680 $ alphar, alphai, beta, vl, ldvl, vr, ldvr, work,
681 $ lwork, ierr )
682 IF( ierr.NE.0 ) THEN
683 IF( ierr.GT.0 .AND. ierr.LE.n ) THEN
684 info = ierr
685 ELSE IF( ierr.GT.n .AND. ierr.LE.2*n ) THEN
686 info = ierr - n
687 ELSE
688 info = n + 1
689 END IF
690 GO TO 130
691 END IF
692*
693* Compute Eigenvectors and estimate condition numbers if desired
694* (Workspace: DTGEVC: need 6*N
695* DTGSNA: need 2*N*(N+2)+16 if SENSE = 'V' or 'B',
696* need N otherwise )
697*
698 IF( ilv .OR. .NOT.wantsn ) THEN
699 IF( ilv ) THEN
700 IF( ilvl ) THEN
701 IF( ilvr ) THEN
702 chtemp = 'B'
703 ELSE
704 chtemp = 'L'
705 END IF
706 ELSE
707 chtemp = 'R'
708 END IF
709*
710 CALL dtgevc( chtemp, 'B', ldumma, n, a, lda, b, ldb, vl,
711 $ ldvl, vr, ldvr, n, in, work, ierr )
712 IF( ierr.NE.0 ) THEN
713 info = n + 2
714 GO TO 130
715 END IF
716 END IF
717*
718 IF( .NOT.wantsn ) THEN
719*
720* compute eigenvectors (DTGEVC) and estimate condition
721* numbers (DTGSNA). Note that the definition of the condition
722* number is not invariant under transformation (u,v) to
723* (Q*u, Z*v), where (u,v) are eigenvectors of the generalized
724* Schur form (S,T), Q and Z are orthogonal matrices. In order
725* to avoid using extra 2*N*N workspace, we have to recalculate
726* eigenvectors and estimate one condition numbers at a time.
727*
728 pair = .false.
729 DO 20 i = 1, n
730*
731 IF( pair ) THEN
732 pair = .false.
733 GO TO 20
734 END IF
735 mm = 1
736 IF( i.LT.n ) THEN
737 IF( a( i+1, i ).NE.zero ) THEN
738 pair = .true.
739 mm = 2
740 END IF
741 END IF
742*
743 DO 10 j = 1, n
744 bwork( j ) = .false.
745 10 CONTINUE
746 IF( mm.EQ.1 ) THEN
747 bwork( i ) = .true.
748 ELSE IF( mm.EQ.2 ) THEN
749 bwork( i ) = .true.
750 bwork( i+1 ) = .true.
751 END IF
752*
753 iwrk = mm*n + 1
754 iwrk1 = iwrk + mm*n
755*
756* Compute a pair of left and right eigenvectors.
757* (compute workspace: need up to 4*N + 6*N)
758*
759 IF( wantse .OR. wantsb ) THEN
760 CALL dtgevc( 'B', 'S', bwork, n, a, lda, b, ldb,
761 $ work( 1 ), n, work( iwrk ), n, mm, m,
762 $ work( iwrk1 ), ierr )
763 IF( ierr.NE.0 ) THEN
764 info = n + 2
765 GO TO 130
766 END IF
767 END IF
768*
769 CALL dtgsna( sense, 'S', bwork, n, a, lda, b, ldb,
770 $ work( 1 ), n, work( iwrk ), n, rconde( i ),
771 $ rcondv( i ), mm, m, work( iwrk1 ),
772 $ lwork-iwrk1+1, iwork, ierr )
773*
774 20 CONTINUE
775 END IF
776 END IF
777*
778* Undo balancing on VL and VR and normalization
779* (Workspace: none needed)
780*
781 IF( ilvl ) THEN
782 CALL dggbak( balanc, 'L', n, ilo, ihi, lscale, rscale, n, vl,
783 $ ldvl, ierr )
784*
785 DO 70 jc = 1, n
786 IF( alphai( jc ).LT.zero )
787 $ GO TO 70
788 temp = zero
789 IF( alphai( jc ).EQ.zero ) THEN
790 DO 30 jr = 1, n
791 temp = max( temp, abs( vl( jr, jc ) ) )
792 30 CONTINUE
793 ELSE
794 DO 40 jr = 1, n
795 temp = max( temp, abs( vl( jr, jc ) )+
796 $ abs( vl( jr, jc+1 ) ) )
797 40 CONTINUE
798 END IF
799 IF( temp.LT.smlnum )
800 $ GO TO 70
801 temp = one / temp
802 IF( alphai( jc ).EQ.zero ) THEN
803 DO 50 jr = 1, n
804 vl( jr, jc ) = vl( jr, jc )*temp
805 50 CONTINUE
806 ELSE
807 DO 60 jr = 1, n
808 vl( jr, jc ) = vl( jr, jc )*temp
809 vl( jr, jc+1 ) = vl( jr, jc+1 )*temp
810 60 CONTINUE
811 END IF
812 70 CONTINUE
813 END IF
814 IF( ilvr ) THEN
815 CALL dggbak( balanc, 'R', n, ilo, ihi, lscale, rscale, n, vr,
816 $ ldvr, ierr )
817 DO 120 jc = 1, n
818 IF( alphai( jc ).LT.zero )
819 $ GO TO 120
820 temp = zero
821 IF( alphai( jc ).EQ.zero ) THEN
822 DO 80 jr = 1, n
823 temp = max( temp, abs( vr( jr, jc ) ) )
824 80 CONTINUE
825 ELSE
826 DO 90 jr = 1, n
827 temp = max( temp, abs( vr( jr, jc ) )+
828 $ abs( vr( jr, jc+1 ) ) )
829 90 CONTINUE
830 END IF
831 IF( temp.LT.smlnum )
832 $ GO TO 120
833 temp = one / temp
834 IF( alphai( jc ).EQ.zero ) THEN
835 DO 100 jr = 1, n
836 vr( jr, jc ) = vr( jr, jc )*temp
837 100 CONTINUE
838 ELSE
839 DO 110 jr = 1, n
840 vr( jr, jc ) = vr( jr, jc )*temp
841 vr( jr, jc+1 ) = vr( jr, jc+1 )*temp
842 110 CONTINUE
843 END IF
844 120 CONTINUE
845 END IF
846*
847* Undo scaling if necessary
848*
849 130 CONTINUE
850*
851 IF( ilascl ) THEN
852 CALL dlascl( 'G', 0, 0, anrmto, anrm, n, 1, alphar, n, ierr )
853 CALL dlascl( 'G', 0, 0, anrmto, anrm, n, 1, alphai, n, ierr )
854 END IF
855*
856 IF( ilbscl ) THEN
857 CALL dlascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
858 END IF
859*
860 work( 1 ) = maxwrk
861 RETURN
862*
863* End of DGGEVX
864*
865 END
dlabad
subroutine dlabad(small, large)
DLABAD
Definition dlabad.f:74
dlacpy
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
Definition dlacpy.f:103
dlascl
subroutine dlascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition dlascl.f:143
dlaset
subroutine dlaset(uplo, m, n, alpha, beta, a, lda)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition dlaset.f:110
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
dggbal
subroutine dggbal(job, n, a, lda, b, ldb, ilo, ihi, lscale, rscale, work, info)
DGGBAL
Definition dggbal.f:177
dggbak
subroutine dggbak(job, side, n, ilo, ihi, lscale, rscale, m, v, ldv, info)
DGGBAK
Definition dggbak.f:147
dhgeqz
subroutine dhgeqz(job, compq, compz, n, ilo, ihi, h, ldh, t, ldt, alphar, alphai, beta, q, ldq, z, ldz, work, lwork, info)
DHGEQZ
Definition dhgeqz.f:304
dtgevc
subroutine dtgevc(side, howmny, select, n, s, lds, p, ldp, vl, ldvl, vr, ldvr, mm, m, work, info)
DTGEVC
Definition dtgevc.f:295
dgeqrf
subroutine dgeqrf(m, n, a, lda, tau, work, lwork, info)
DGEQRF
Definition dgeqrf.f:146
dggevx
subroutine dggevx(balanc, jobvl, jobvr, sense, n, a, lda, b, ldb, alphar, alphai, beta, vl, ldvl, vr, ldvr, ilo, ihi, lscale, rscale, abnrm, bbnrm, rconde, rcondv, work, lwork, iwork, bwork, info)
DGGEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices
Definition dggevx.f:391
dgghrd
subroutine dgghrd(compq, compz, n, ilo, ihi, a, lda, b, ldb, q, ldq, z, ldz, info)
DGGHRD
Definition dgghrd.f:207
dormqr
subroutine dormqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
DORMQR
Definition dormqr.f:167
dtgsna
subroutine dtgsna(job, howmny, select, n, a, lda, b, ldb, vl, ldvl, vr, ldvr, s, dif, mm, m, work, lwork, iwork, info)
DTGSNA
Definition dtgsna.f:381
dorgqr
subroutine dorgqr(m, n, k, a, lda, tau, work, lwork, info)
DORGQR
Definition dorgqr.f:128
max
#define max(a, b)
Definition macros.h:21