LAPACK  3.10.1
LAPACK: Linear Algebra PACKage

◆ zggev3()

subroutine zggev3 ( character  JOBVL,
character  JOBVR,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( ldb, * )  B,
integer  LDB,
complex*16, dimension( * )  ALPHA,
complex*16, dimension( * )  BETA,
complex*16, dimension( ldvl, * )  VL,
integer  LDVL,
complex*16, dimension( ldvr, * )  VR,
integer  LDVR,
complex*16, dimension( * )  WORK,
integer  LWORK,
double precision, dimension( * )  RWORK,
integer  INFO 
)

ZGGEV3 computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices (blocked algorithm)

Download ZGGEV3 + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 ZGGEV3 computes for a pair of N-by-N complex nonsymmetric matrices
 (A,B), the generalized eigenvalues, and optionally, the left and/or
 right generalized eigenvectors.

 A generalized eigenvalue for a pair of matrices (A,B) is a scalar
 lambda or a ratio alpha/beta = lambda, such that A - lambda*B is
 singular. It is usually represented as the pair (alpha,beta), as
 there is a reasonable interpretation for beta=0, and even for both
 being zero.

 The right generalized eigenvector v(j) corresponding to the
 generalized eigenvalue lambda(j) of (A,B) satisfies

              A * v(j) = lambda(j) * B * v(j).

 The left generalized eigenvector u(j) corresponding to the
 generalized eigenvalues lambda(j) of (A,B) satisfies

              u(j)**H * A = lambda(j) * u(j)**H * B

 where u(j)**H is the conjugate-transpose of u(j).
Parameters
[in]JOBVL
          JOBVL is CHARACTER*1
          = 'N':  do not compute the left generalized eigenvectors;
          = 'V':  compute the left generalized eigenvectors.
[in]JOBVR
          JOBVR is CHARACTER*1
          = 'N':  do not compute the right generalized eigenvectors;
          = 'V':  compute the right generalized eigenvectors.
[in]N
          N is INTEGER
          The order of the matrices A, B, VL, and VR.  N >= 0.
[in,out]A
          A is COMPLEX*16 array, dimension (LDA, N)
          On entry, the matrix A in the pair (A,B).
          On exit, A has been overwritten.
[in]LDA
          LDA is INTEGER
          The leading dimension of A.  LDA >= max(1,N).
[in,out]B
          B is COMPLEX*16 array, dimension (LDB, N)
          On entry, the matrix B in the pair (A,B).
          On exit, B has been overwritten.
[in]LDB
          LDB is INTEGER
          The leading dimension of B.  LDB >= max(1,N).
[out]ALPHA
          ALPHA is COMPLEX*16 array, dimension (N)
[out]BETA
          BETA is COMPLEX*16 array, dimension (N)
          On exit, ALPHA(j)/BETA(j), j=1,...,N, will be the
          generalized eigenvalues.

          Note: the quotients ALPHA(j)/BETA(j) may easily over- or
          underflow, and BETA(j) may even be zero.  Thus, the user
          should avoid naively computing the ratio alpha/beta.
          However, ALPHA will be always less than and usually
          comparable with norm(A) in magnitude, and BETA always less
          than and usually comparable with norm(B).
[out]VL
          VL is COMPLEX*16 array, dimension (LDVL,N)
          If JOBVL = 'V', the left generalized eigenvectors u(j) are
          stored one after another in the columns of VL, in the same
          order as their eigenvalues.
          Each eigenvector is scaled so the largest component has
          abs(real part) + abs(imag. part) = 1.
          Not referenced if JOBVL = 'N'.
[in]LDVL
          LDVL is INTEGER
          The leading dimension of the matrix VL. LDVL >= 1, and
          if JOBVL = 'V', LDVL >= N.
[out]VR
          VR is COMPLEX*16 array, dimension (LDVR,N)
          If JOBVR = 'V', the right generalized eigenvectors v(j) are
          stored one after another in the columns of VR, in the same
          order as their eigenvalues.
          Each eigenvector is scaled so the largest component has
          abs(real part) + abs(imag. part) = 1.
          Not referenced if JOBVR = 'N'.
[in]LDVR
          LDVR is INTEGER
          The leading dimension of the matrix VR. LDVR >= 1, and
          if JOBVR = 'V', LDVR >= N.
[out]WORK
          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.
[out]RWORK
          RWORK is DOUBLE PRECISION array, dimension (8*N)
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value.
          =1,...,N:
                The QZ iteration failed.  No eigenvectors have been
                calculated, but ALPHA(j) and BETA(j) should be
                correct for j=INFO+1,...,N.
          > N:  =N+1: other then QZ iteration failed in ZHGEQZ,
                =N+2: error return from ZTGEVC.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 214 of file zggev3.f.

216 *
217 * -- LAPACK driver routine --
218 * -- LAPACK is a software package provided by Univ. of Tennessee, --
219 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
220 *
221 * .. Scalar Arguments ..
222  CHARACTER JOBVL, JOBVR
223  INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, N
224 * ..
225 * .. Array Arguments ..
226  DOUBLE PRECISION RWORK( * )
227  COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDB, * ),
228  $ BETA( * ), VL( LDVL, * ), VR( LDVR, * ),
229  $ WORK( * )
230 * ..
231 *
232 * =====================================================================
233 *
234 * .. Parameters ..
235  DOUBLE PRECISION ZERO, ONE
236  parameter( zero = 0.0d0, one = 1.0d0 )
237  COMPLEX*16 CZERO, CONE
238  parameter( czero = ( 0.0d0, 0.0d0 ),
239  $ cone = ( 1.0d0, 0.0d0 ) )
240 * ..
241 * .. Local Scalars ..
242  LOGICAL ILASCL, ILBSCL, ILV, ILVL, ILVR, LQUERY
243  CHARACTER CHTEMP
244  INTEGER ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT, ILO,
245  $ IN, IRIGHT, IROWS, IRWRK, ITAU, IWRK, JC, JR,
246  $ LWKOPT
247  DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
248  $ SMLNUM, TEMP
249  COMPLEX*16 X
250 * ..
251 * .. Local Arrays ..
252  LOGICAL LDUMMA( 1 )
253 * ..
254 * .. External Subroutines ..
255  EXTERNAL dlabad, xerbla, zgeqrf, zggbak, zggbal, zgghd3,
257  $ zunmqr
258 * ..
259 * .. External Functions ..
260  LOGICAL LSAME
261  DOUBLE PRECISION DLAMCH, ZLANGE
262  EXTERNAL lsame, dlamch, zlange
263 * ..
264 * .. Intrinsic Functions ..
265  INTRINSIC abs, dble, dimag, max, sqrt
266 * ..
267 * .. Statement Functions ..
268  DOUBLE PRECISION ABS1
269 * ..
270 * .. Statement Function definitions ..
271  abs1( x ) = abs( dble( x ) ) + abs( dimag( x ) )
272 * ..
273 * .. Executable Statements ..
274 *
275 * Decode the input arguments
276 *
277  IF( lsame( jobvl, 'N' ) ) THEN
278  ijobvl = 1
279  ilvl = .false.
280  ELSE IF( lsame( jobvl, 'V' ) ) THEN
281  ijobvl = 2
282  ilvl = .true.
283  ELSE
284  ijobvl = -1
285  ilvl = .false.
286  END IF
287 *
288  IF( lsame( jobvr, 'N' ) ) THEN
289  ijobvr = 1
290  ilvr = .false.
291  ELSE IF( lsame( jobvr, 'V' ) ) THEN
292  ijobvr = 2
293  ilvr = .true.
294  ELSE
295  ijobvr = -1
296  ilvr = .false.
297  END IF
298  ilv = ilvl .OR. ilvr
299 *
300 * Test the input arguments
301 *
302  info = 0
303  lquery = ( lwork.EQ.-1 )
304  IF( ijobvl.LE.0 ) THEN
305  info = -1
306  ELSE IF( ijobvr.LE.0 ) THEN
307  info = -2
308  ELSE IF( n.LT.0 ) THEN
309  info = -3
310  ELSE IF( lda.LT.max( 1, n ) ) THEN
311  info = -5
312  ELSE IF( ldb.LT.max( 1, n ) ) THEN
313  info = -7
314  ELSE IF( ldvl.LT.1 .OR. ( ilvl .AND. ldvl.LT.n ) ) THEN
315  info = -11
316  ELSE IF( ldvr.LT.1 .OR. ( ilvr .AND. ldvr.LT.n ) ) THEN
317  info = -13
318  ELSE IF( lwork.LT.max( 1, 2*n ) .AND. .NOT.lquery ) THEN
319  info = -15
320  END IF
321 *
322 * Compute workspace
323 *
324  IF( info.EQ.0 ) THEN
325  CALL zgeqrf( n, n, b, ldb, work, work, -1, ierr )
326  lwkopt = max( 1, n+int( work( 1 ) ) )
327  CALL zunmqr( 'L', 'C', n, n, n, b, ldb, work, a, lda, work,
328  $ -1, ierr )
329  lwkopt = max( lwkopt, n+int( work( 1 ) ) )
330  IF( ilvl ) THEN
331  CALL zungqr( n, n, n, vl, ldvl, work, work, -1, ierr )
332  lwkopt = max( lwkopt, n+int( work( 1 ) ) )
333  END IF
334  IF( ilv ) THEN
335  CALL zgghd3( jobvl, jobvr, n, 1, n, a, lda, b, ldb, vl,
336  $ ldvl, vr, ldvr, work, -1, ierr )
337  lwkopt = max( lwkopt, n+int( work( 1 ) ) )
338  CALL zlaqz0( 'S', jobvl, jobvr, n, 1, n, a, lda, b, ldb,
339  $ alpha, beta, vl, ldvl, vr, ldvr, work, -1,
340  $ rwork, 0, ierr )
341  lwkopt = max( lwkopt, n+int( work( 1 ) ) )
342  ELSE
343  CALL zgghd3( jobvl, jobvr, n, 1, n, a, lda, b, ldb, vl,
344  $ ldvl, vr, ldvr, work, -1, ierr )
345  lwkopt = max( lwkopt, n+int( work( 1 ) ) )
346  CALL zlaqz0( 'E', jobvl, jobvr, n, 1, n, a, lda, b, ldb,
347  $ alpha, beta, vl, ldvl, vr, ldvr, work, -1,
348  $ rwork, 0, ierr )
349  lwkopt = max( lwkopt, n+int( work( 1 ) ) )
350  END IF
351  work( 1 ) = dcmplx( lwkopt )
352  END IF
353 *
354  IF( info.NE.0 ) THEN
355  CALL xerbla( 'ZGGEV3 ', -info )
356  RETURN
357  ELSE IF( lquery ) THEN
358  RETURN
359  END IF
360 *
361 * Quick return if possible
362 *
363  IF( n.EQ.0 )
364  $ RETURN
365 *
366 * Get machine constants
367 *
368  eps = dlamch( 'E' )*dlamch( 'B' )
369  smlnum = dlamch( 'S' )
370  bignum = one / smlnum
371  CALL dlabad( smlnum, bignum )
372  smlnum = sqrt( smlnum ) / eps
373  bignum = one / smlnum
374 *
375 * Scale A if max element outside range [SMLNUM,BIGNUM]
376 *
377  anrm = zlange( 'M', n, n, a, lda, rwork )
378  ilascl = .false.
379  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
380  anrmto = smlnum
381  ilascl = .true.
382  ELSE IF( anrm.GT.bignum ) THEN
383  anrmto = bignum
384  ilascl = .true.
385  END IF
386  IF( ilascl )
387  $ CALL zlascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
388 *
389 * Scale B if max element outside range [SMLNUM,BIGNUM]
390 *
391  bnrm = zlange( 'M', n, n, b, ldb, rwork )
392  ilbscl = .false.
393  IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
394  bnrmto = smlnum
395  ilbscl = .true.
396  ELSE IF( bnrm.GT.bignum ) THEN
397  bnrmto = bignum
398  ilbscl = .true.
399  END IF
400  IF( ilbscl )
401  $ CALL zlascl( 'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
402 *
403 * Permute the matrices A, B to isolate eigenvalues if possible
404 *
405  ileft = 1
406  iright = n + 1
407  irwrk = iright + n
408  CALL zggbal( 'P', n, a, lda, b, ldb, ilo, ihi, rwork( ileft ),
409  $ rwork( iright ), rwork( irwrk ), ierr )
410 *
411 * Reduce B to triangular form (QR decomposition of B)
412 *
413  irows = ihi + 1 - ilo
414  IF( ilv ) THEN
415  icols = n + 1 - ilo
416  ELSE
417  icols = irows
418  END IF
419  itau = 1
420  iwrk = itau + irows
421  CALL zgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
422  $ work( iwrk ), lwork+1-iwrk, ierr )
423 *
424 * Apply the orthogonal transformation to matrix A
425 *
426  CALL zunmqr( 'L', 'C', irows, icols, irows, b( ilo, ilo ), ldb,
427  $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
428  $ lwork+1-iwrk, ierr )
429 *
430 * Initialize VL
431 *
432  IF( ilvl ) THEN
433  CALL zlaset( 'Full', n, n, czero, cone, vl, ldvl )
434  IF( irows.GT.1 ) THEN
435  CALL zlacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
436  $ vl( ilo+1, ilo ), ldvl )
437  END IF
438  CALL zungqr( irows, irows, irows, vl( ilo, ilo ), ldvl,
439  $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
440  END IF
441 *
442 * Initialize VR
443 *
444  IF( ilvr )
445  $ CALL zlaset( 'Full', n, n, czero, cone, vr, ldvr )
446 *
447 * Reduce to generalized Hessenberg form
448 *
449  IF( ilv ) THEN
450 *
451 * Eigenvectors requested -- work on whole matrix.
452 *
453  CALL zgghd3( jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb, vl,
454  $ ldvl, vr, ldvr, work( iwrk ), lwork+1-iwrk, ierr )
455  ELSE
456  CALL zgghd3( 'N', 'N', irows, 1, irows, a( ilo, ilo ), lda,
457  $ b( ilo, ilo ), ldb, vl, ldvl, vr, ldvr,
458  $ work( iwrk ), lwork+1-iwrk, ierr )
459  END IF
460 *
461 * Perform QZ algorithm (Compute eigenvalues, and optionally, the
462 * Schur form and Schur vectors)
463 *
464  iwrk = itau
465  IF( ilv ) THEN
466  chtemp = 'S'
467  ELSE
468  chtemp = 'E'
469  END IF
470  CALL zlaqz0( chtemp, jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb,
471  $ alpha, beta, vl, ldvl, vr, ldvr, work( iwrk ),
472  $ lwork+1-iwrk, rwork( irwrk ), 0, ierr )
473  IF( ierr.NE.0 ) THEN
474  IF( ierr.GT.0 .AND. ierr.LE.n ) THEN
475  info = ierr
476  ELSE IF( ierr.GT.n .AND. ierr.LE.2*n ) THEN
477  info = ierr - n
478  ELSE
479  info = n + 1
480  END IF
481  GO TO 70
482  END IF
483 *
484 * Compute Eigenvectors
485 *
486  IF( ilv ) THEN
487  IF( ilvl ) THEN
488  IF( ilvr ) THEN
489  chtemp = 'B'
490  ELSE
491  chtemp = 'L'
492  END IF
493  ELSE
494  chtemp = 'R'
495  END IF
496 *
497  CALL ztgevc( chtemp, 'B', ldumma, n, a, lda, b, ldb, vl, ldvl,
498  $ vr, ldvr, n, in, work( iwrk ), rwork( irwrk ),
499  $ ierr )
500  IF( ierr.NE.0 ) THEN
501  info = n + 2
502  GO TO 70
503  END IF
504 *
505 * Undo balancing on VL and VR and normalization
506 *
507  IF( ilvl ) THEN
508  CALL zggbak( 'P', 'L', n, ilo, ihi, rwork( ileft ),
509  $ rwork( iright ), n, vl, ldvl, ierr )
510  DO 30 jc = 1, n
511  temp = zero
512  DO 10 jr = 1, n
513  temp = max( temp, abs1( vl( jr, jc ) ) )
514  10 CONTINUE
515  IF( temp.LT.smlnum )
516  $ GO TO 30
517  temp = one / temp
518  DO 20 jr = 1, n
519  vl( jr, jc ) = vl( jr, jc )*temp
520  20 CONTINUE
521  30 CONTINUE
522  END IF
523  IF( ilvr ) THEN
524  CALL zggbak( 'P', 'R', n, ilo, ihi, rwork( ileft ),
525  $ rwork( iright ), n, vr, ldvr, ierr )
526  DO 60 jc = 1, n
527  temp = zero
528  DO 40 jr = 1, n
529  temp = max( temp, abs1( vr( jr, jc ) ) )
530  40 CONTINUE
531  IF( temp.LT.smlnum )
532  $ GO TO 60
533  temp = one / temp
534  DO 50 jr = 1, n
535  vr( jr, jc ) = vr( jr, jc )*temp
536  50 CONTINUE
537  60 CONTINUE
538  END IF
539  END IF
540 *
541 * Undo scaling if necessary
542 *
543  70 CONTINUE
544 *
545  IF( ilascl )
546  $ CALL zlascl( 'G', 0, 0, anrmto, anrm, n, 1, alpha, n, ierr )
547 *
548  IF( ilbscl )
549  $ CALL zlascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
550 *
551  work( 1 ) = dcmplx( lwkopt )
552  RETURN
553 *
554 * End of ZGGEV3
555 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:74
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine zggbal(JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE, WORK, INFO)
ZGGBAL
Definition: zggbal.f:177
subroutine zggbak(JOB, SIDE, N, ILO, IHI, LSCALE, RSCALE, M, V, LDV, INFO)
ZGGBAK
Definition: zggbak.f:148
double precision function zlange(NORM, M, N, A, LDA, WORK)
ZLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: zlange.f:115
recursive subroutine zlaqz0(WANTS, WANTQ, WANTZ, N, ILO, IHI, A, LDA, B, LDB, ALPHA, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, RWORK, REC, INFO)
ZLAQZ0
Definition: zlaqz0.f:284
subroutine ztgevc(SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL, LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO)
ZTGEVC
Definition: ztgevc.f:219
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:143
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:103
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:106
subroutine zungqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGQR
Definition: zungqr.f:128
subroutine zgghd3(COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, LDQ, Z, LDZ, WORK, LWORK, INFO)
ZGGHD3
Definition: zgghd3.f:227
subroutine zunmqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMQR
Definition: zunmqr.f:167
subroutine zgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
ZGEQRF VARIANT: left-looking Level 3 BLAS of the algorithm.
Definition: zgeqrf.f:151
Here is the call graph for this function:
Here is the caller graph for this function: