LAPACK 3.11.0
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.

          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:152
Here is the call graph for this function:
Here is the caller graph for this function: