LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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. LWORK >= MAX(1,2*N).
!>          For good performance, LWORK must generally be larger.
!>
!>          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 213 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 $ LWKMIN, 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 xerbla, zgeqrf, zggbak, zggbal, zgghd3,
256 $ zlaqz0,
257 $ zlacpy, zlascl, zlaset, ztgevc, zungqr, 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 lwkmin = max( 1, 2*n )
305 IF( ijobvl.LE.0 ) THEN
306 info = -1
307 ELSE IF( ijobvr.LE.0 ) THEN
308 info = -2
309 ELSE IF( n.LT.0 ) THEN
310 info = -3
311 ELSE IF( lda.LT.max( 1, n ) ) THEN
312 info = -5
313 ELSE IF( ldb.LT.max( 1, n ) ) THEN
314 info = -7
315 ELSE IF( ldvl.LT.1 .OR. ( ilvl .AND. ldvl.LT.n ) ) THEN
316 info = -11
317 ELSE IF( ldvr.LT.1 .OR. ( ilvr .AND. ldvr.LT.n ) ) THEN
318 info = -13
319 ELSE IF( lwork.LT.lwkmin .AND. .NOT.lquery ) THEN
320 info = -15
321 END IF
322*
323* Compute workspace
324*
325 IF( info.EQ.0 ) THEN
326 CALL zgeqrf( n, n, b, ldb, work, work, -1, ierr )
327 lwkopt = max( lwkmin, n+int( work( 1 ) ) )
328 CALL zunmqr( 'L', 'C', n, n, n, b, ldb, work, a, lda, work,
329 $ -1, ierr )
330 lwkopt = max( lwkopt, n+int( work( 1 ) ) )
331 IF( ilvl ) THEN
332 CALL zungqr( n, n, n, vl, ldvl, work, work, -1, ierr )
333 lwkopt = max( lwkopt, n+int( work( 1 ) ) )
334 END IF
335 IF( ilv ) THEN
336 CALL zgghd3( jobvl, jobvr, n, 1, n, a, lda, b, ldb, vl,
337 $ ldvl, vr, ldvr, work, -1, ierr )
338 lwkopt = max( lwkopt, n+int( work( 1 ) ) )
339 CALL zlaqz0( 'S', jobvl, jobvr, n, 1, n, a, lda, b, ldb,
340 $ alpha, beta, vl, ldvl, vr, ldvr, work, -1,
341 $ rwork, 0, ierr )
342 lwkopt = max( lwkopt, n+int( work( 1 ) ) )
343 ELSE
344 CALL zgghd3( jobvl, jobvr, n, 1, n, a, lda, b, ldb, vl,
345 $ ldvl, vr, ldvr, work, -1, ierr )
346 lwkopt = max( lwkopt, n+int( work( 1 ) ) )
347 CALL zlaqz0( 'E', jobvl, jobvr, n, 1, n, a, lda, b, ldb,
348 $ alpha, beta, vl, ldvl, vr, ldvr, work, -1,
349 $ rwork, 0, ierr )
350 lwkopt = max( lwkopt, n+int( work( 1 ) ) )
351 END IF
352 IF( n.EQ.0 ) THEN
353 work( 1 ) = 1
354 ELSE
355 work( 1 ) = dcmplx( lwkopt )
356 END IF
357 END IF
358*
359 IF( info.NE.0 ) THEN
360 CALL xerbla( 'ZGGEV3 ', -info )
361 RETURN
362 ELSE IF( lquery ) THEN
363 RETURN
364 END IF
365*
366* Quick return if possible
367*
368 IF( n.EQ.0 )
369 $ RETURN
370*
371* Get machine constants
372*
373 eps = dlamch( 'E' )*dlamch( 'B' )
374 smlnum = dlamch( 'S' )
375 bignum = one / smlnum
376 smlnum = sqrt( smlnum ) / eps
377 bignum = one / smlnum
378*
379* Scale A if max element outside range [SMLNUM,BIGNUM]
380*
381 anrm = zlange( 'M', n, n, a, lda, rwork )
382 ilascl = .false.
383 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
384 anrmto = smlnum
385 ilascl = .true.
386 ELSE IF( anrm.GT.bignum ) THEN
387 anrmto = bignum
388 ilascl = .true.
389 END IF
390 IF( ilascl )
391 $ CALL zlascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
392*
393* Scale B if max element outside range [SMLNUM,BIGNUM]
394*
395 bnrm = zlange( 'M', n, n, b, ldb, rwork )
396 ilbscl = .false.
397 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
398 bnrmto = smlnum
399 ilbscl = .true.
400 ELSE IF( bnrm.GT.bignum ) THEN
401 bnrmto = bignum
402 ilbscl = .true.
403 END IF
404 IF( ilbscl )
405 $ CALL zlascl( 'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
406*
407* Permute the matrices A, B to isolate eigenvalues if possible
408*
409 ileft = 1
410 iright = n + 1
411 irwrk = iright + n
412 CALL zggbal( 'P', n, a, lda, b, ldb, ilo, ihi, rwork( ileft ),
413 $ rwork( iright ), rwork( irwrk ), ierr )
414*
415* Reduce B to triangular form (QR decomposition of B)
416*
417 irows = ihi + 1 - ilo
418 IF( ilv ) THEN
419 icols = n + 1 - ilo
420 ELSE
421 icols = irows
422 END IF
423 itau = 1
424 iwrk = itau + irows
425 CALL zgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
426 $ work( iwrk ), lwork+1-iwrk, ierr )
427*
428* Apply the orthogonal transformation to matrix A
429*
430 CALL zunmqr( 'L', 'C', irows, icols, irows, b( ilo, ilo ), ldb,
431 $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
432 $ lwork+1-iwrk, ierr )
433*
434* Initialize VL
435*
436 IF( ilvl ) THEN
437 CALL zlaset( 'Full', n, n, czero, cone, vl, ldvl )
438 IF( irows.GT.1 ) THEN
439 CALL zlacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
440 $ vl( ilo+1, ilo ), ldvl )
441 END IF
442 CALL zungqr( irows, irows, irows, vl( ilo, ilo ), ldvl,
443 $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
444 END IF
445*
446* Initialize VR
447*
448 IF( ilvr )
449 $ CALL zlaset( 'Full', n, n, czero, cone, vr, ldvr )
450*
451* Reduce to generalized Hessenberg form
452*
453 IF( ilv ) THEN
454*
455* Eigenvectors requested -- work on whole matrix.
456*
457 CALL zgghd3( jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb, vl,
458 $ ldvl, vr, ldvr, work( iwrk ), lwork+1-iwrk, ierr )
459 ELSE
460 CALL zgghd3( 'N', 'N', irows, 1, irows, a( ilo, ilo ), lda,
461 $ b( ilo, ilo ), ldb, vl, ldvl, vr, ldvr,
462 $ work( iwrk ), lwork+1-iwrk, ierr )
463 END IF
464*
465* Perform QZ algorithm (Compute eigenvalues, and optionally, the
466* Schur form and Schur vectors)
467*
468 iwrk = itau
469 IF( ilv ) THEN
470 chtemp = 'S'
471 ELSE
472 chtemp = 'E'
473 END IF
474 CALL zlaqz0( chtemp, jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb,
475 $ alpha, beta, vl, ldvl, vr, ldvr, work( iwrk ),
476 $ lwork+1-iwrk, rwork( irwrk ), 0, ierr )
477 IF( ierr.NE.0 ) THEN
478 IF( ierr.GT.0 .AND. ierr.LE.n ) THEN
479 info = ierr
480 ELSE IF( ierr.GT.n .AND. ierr.LE.2*n ) THEN
481 info = ierr - n
482 ELSE
483 info = n + 1
484 END IF
485 GO TO 70
486 END IF
487*
488* Compute Eigenvectors
489*
490 IF( ilv ) THEN
491 IF( ilvl ) THEN
492 IF( ilvr ) THEN
493 chtemp = 'B'
494 ELSE
495 chtemp = 'L'
496 END IF
497 ELSE
498 chtemp = 'R'
499 END IF
500*
501 CALL ztgevc( chtemp, 'B', ldumma, n, a, lda, b, ldb, vl,
502 $ ldvl,
503 $ vr, ldvr, n, in, work( iwrk ), rwork( irwrk ),
504 $ ierr )
505 IF( ierr.NE.0 ) THEN
506 info = n + 2
507 GO TO 70
508 END IF
509*
510* Undo balancing on VL and VR and normalization
511*
512 IF( ilvl ) THEN
513 CALL zggbak( 'P', 'L', n, ilo, ihi, rwork( ileft ),
514 $ rwork( iright ), n, vl, ldvl, ierr )
515 DO 30 jc = 1, n
516 temp = zero
517 DO 10 jr = 1, n
518 temp = max( temp, abs1( vl( jr, jc ) ) )
519 10 CONTINUE
520 IF( temp.LT.smlnum )
521 $ GO TO 30
522 temp = one / temp
523 DO 20 jr = 1, n
524 vl( jr, jc ) = vl( jr, jc )*temp
525 20 CONTINUE
526 30 CONTINUE
527 END IF
528 IF( ilvr ) THEN
529 CALL zggbak( 'P', 'R', n, ilo, ihi, rwork( ileft ),
530 $ rwork( iright ), n, vr, ldvr, ierr )
531 DO 60 jc = 1, n
532 temp = zero
533 DO 40 jr = 1, n
534 temp = max( temp, abs1( vr( jr, jc ) ) )
535 40 CONTINUE
536 IF( temp.LT.smlnum )
537 $ GO TO 60
538 temp = one / temp
539 DO 50 jr = 1, n
540 vr( jr, jc ) = vr( jr, jc )*temp
541 50 CONTINUE
542 60 CONTINUE
543 END IF
544 END IF
545*
546* Undo scaling if necessary
547*
548 70 CONTINUE
549*
550 IF( ilascl )
551 $ CALL zlascl( 'G', 0, 0, anrmto, anrm, n, 1, alpha, n, ierr )
552*
553 IF( ilbscl )
554 $ CALL zlascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
555*
556 work( 1 ) = dcmplx( lwkopt )
557 RETURN
558*
559* End of ZGGEV3
560*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
zgeqrf
subroutine zgeqrf(m, n, a, lda, tau, work, lwork, info)
ZGEQRF
Definition zgeqrf.f:144
zggbak
subroutine zggbak(job, side, n, ilo, ihi, lscale, rscale, m, v, ldv, info)
ZGGBAK
Definition zggbak.f:147
zggbal
subroutine zggbal(job, n, a, lda, b, ldb, ilo, ihi, lscale, rscale, work, info)
ZGGBAL
Definition zggbal.f:175
zgghd3
subroutine zgghd3(compq, compz, n, ilo, ihi, a, lda, b, ldb, q, ldq, z, ldz, work, lwork, info)
ZGGHD3
Definition zgghd3.f:226
zlacpy
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
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
zlange
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:113
zlaqz0
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:283
zlascl
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
zlaset
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
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
ztgevc
subroutine ztgevc(side, howmny, select, n, s, lds, p, ldp, vl, ldvl, vr, ldvr, mm, m, work, rwork, info)
ZTGEVC
Definition ztgevc.f:217
zungqr
subroutine zungqr(m, n, k, a, lda, tau, work, lwork, info)
ZUNGQR
Definition zungqr.f:126
zunmqr
subroutine zunmqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
ZUNMQR
Definition zunmqr.f:165
Here is the call graph for this function:
Here is the caller graph for this function: