LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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 DHGEQZ,
                =N+2: error return from DTGEVC.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
January 2015

Definition at line 218 of file zggev3.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: