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