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