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

◆ cggev3()

subroutine cggev3 ( character  jobvl,
character  jobvr,
integer  n,
complex, dimension( lda, * )  a,
integer  lda,
complex, dimension( ldb, * )  b,
integer  ldb,
complex, dimension( * )  alpha,
complex, dimension( * )  beta,
complex, dimension( ldvl, * )  vl,
integer  ldvl,
complex, dimension( ldvr, * )  vr,
integer  ldvr,
complex, dimension( * )  work,
integer  lwork,
real, dimension( * )  rwork,
integer  info 
)

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

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

Purpose:
 CGGEV3 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 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 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 array, dimension (N)
[out]BETA
          BETA is COMPLEX 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 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 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 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 REAL 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 CHGEQZ,
                =N+2: error return from CTGEVC.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 214 of file cggev3.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 REAL RWORK( * )
227 COMPLEX A( LDA, * ), ALPHA( * ), B( LDB, * ),
228 $ BETA( * ), VL( LDVL, * ), VR( LDVR, * ),
229 $ WORK( * )
230* ..
231*
232* =====================================================================
233*
234* .. Parameters ..
235 REAL ZERO, ONE
236 parameter( zero = 0.0e0, one = 1.0e0 )
237 COMPLEX CZERO, CONE
238 parameter( czero = ( 0.0e0, 0.0e0 ),
239 $ cone = ( 1.0e0, 0.0e0 ) )
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 REAL ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
248 $ SMLNUM, TEMP
249 COMPLEX X
250* ..
251* .. Local Arrays ..
252 LOGICAL LDUMMA( 1 )
253* ..
254* .. External Subroutines ..
255 EXTERNAL cgeqrf, cggbak, cggbal, cgghd3, claqz0, clacpy,
257* ..
258* .. External Functions ..
259 LOGICAL LSAME
260 REAL CLANGE, SLAMCH
261 EXTERNAL lsame, clange, slamch
262* ..
263* .. Intrinsic Functions ..
264 INTRINSIC abs, aimag, max, real, sqrt
265* ..
266* .. Statement Functions ..
267 REAL ABS1
268* ..
269* .. Statement Function definitions ..
270 abs1( x ) = abs( real( x ) ) + abs( aimag( 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 cgeqrf( n, n, b, ldb, work, work, -1, ierr )
325 lwkopt = max( n, n+int( work( 1 ) ) )
326 CALL cunmqr( '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 cungqr( 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 cgghd3( 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 claqz0( '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 cgghd3( 'N', 'N', n, 1, n, a, lda, b, ldb, vl, ldvl,
343 $ vr, ldvr, work, -1, ierr )
344 lwkopt = max( lwkopt, n+int( work( 1 ) ) )
345 CALL claqz0( '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 ) = cmplx( lwkopt )
351 END IF
352*
353 IF( info.NE.0 ) THEN
354 CALL xerbla( 'CGGEV3 ', -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 = slamch( 'E' )*slamch( 'B' )
368 smlnum = slamch( '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 = clange( '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 clascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
386*
387* Scale B if max element outside range [SMLNUM,BIGNUM]
388*
389 bnrm = clange( '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 clascl( '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 cggbal( '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 cgeqrf( 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 cunmqr( '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 claset( 'Full', n, n, czero, cone, vl, ldvl )
432 IF( irows.GT.1 ) THEN
433 CALL clacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
434 $ vl( ilo+1, ilo ), ldvl )
435 END IF
436 CALL cungqr( 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 claset( '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 cgghd3( jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb, vl,
452 $ ldvl, vr, ldvr, work( iwrk ), lwork+1-iwrk,
453 $ ierr )
454 ELSE
455 CALL cgghd3( 'N', 'N', irows, 1, irows, a( ilo, ilo ), lda,
456 $ b( ilo, ilo ), ldb, vl, ldvl, vr, ldvr,
457 $ work( iwrk ), lwork+1-iwrk, ierr )
458 END IF
459*
460* Perform QZ algorithm (Compute eigenvalues, and optionally, the
461* Schur form and Schur vectors)
462*
463 iwrk = itau
464 IF( ilv ) THEN
465 chtemp = 'S'
466 ELSE
467 chtemp = 'E'
468 END IF
469 CALL claqz0( chtemp, jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb,
470 $ alpha, beta, vl, ldvl, vr, ldvr, work( iwrk ),
471 $ lwork+1-iwrk, rwork( irwrk ), 0, ierr )
472 IF( ierr.NE.0 ) THEN
473 IF( ierr.GT.0 .AND. ierr.LE.n ) THEN
474 info = ierr
475 ELSE IF( ierr.GT.n .AND. ierr.LE.2*n ) THEN
476 info = ierr - n
477 ELSE
478 info = n + 1
479 END IF
480 GO TO 70
481 END IF
482*
483* Compute Eigenvectors
484*
485 IF( ilv ) THEN
486 IF( ilvl ) THEN
487 IF( ilvr ) THEN
488 chtemp = 'B'
489 ELSE
490 chtemp = 'L'
491 END IF
492 ELSE
493 chtemp = 'R'
494 END IF
495*
496 CALL ctgevc( chtemp, 'B', ldumma, n, a, lda, b, ldb, vl, ldvl,
497 $ vr, ldvr, n, in, work( iwrk ), rwork( irwrk ),
498 $ ierr )
499 IF( ierr.NE.0 ) THEN
500 info = n + 2
501 GO TO 70
502 END IF
503*
504* Undo balancing on VL and VR and normalization
505*
506 IF( ilvl ) THEN
507 CALL cggbak( 'P', 'L', n, ilo, ihi, rwork( ileft ),
508 $ rwork( iright ), n, vl, ldvl, ierr )
509 DO 30 jc = 1, n
510 temp = zero
511 DO 10 jr = 1, n
512 temp = max( temp, abs1( vl( jr, jc ) ) )
513 10 CONTINUE
514 IF( temp.LT.smlnum )
515 $ GO TO 30
516 temp = one / temp
517 DO 20 jr = 1, n
518 vl( jr, jc ) = vl( jr, jc )*temp
519 20 CONTINUE
520 30 CONTINUE
521 END IF
522 IF( ilvr ) THEN
523 CALL cggbak( 'P', 'R', n, ilo, ihi, rwork( ileft ),
524 $ rwork( iright ), n, vr, ldvr, ierr )
525 DO 60 jc = 1, n
526 temp = zero
527 DO 40 jr = 1, n
528 temp = max( temp, abs1( vr( jr, jc ) ) )
529 40 CONTINUE
530 IF( temp.LT.smlnum )
531 $ GO TO 60
532 temp = one / temp
533 DO 50 jr = 1, n
534 vr( jr, jc ) = vr( jr, jc )*temp
535 50 CONTINUE
536 60 CONTINUE
537 END IF
538 END IF
539*
540* Undo scaling if necessary
541*
542 70 CONTINUE
543*
544 IF( ilascl )
545 $ CALL clascl( 'G', 0, 0, anrmto, anrm, n, 1, alpha, n, ierr )
546*
547 IF( ilbscl )
548 $ CALL clascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
549*
550 work( 1 ) = cmplx( lwkopt )
551 RETURN
552*
553* End of CGGEV3
554*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine cgeqrf(m, n, a, lda, tau, work, lwork, info)
CGEQRF
Definition cgeqrf.f:146
subroutine cggbak(job, side, n, ilo, ihi, lscale, rscale, m, v, ldv, info)
CGGBAK
Definition cggbak.f:148
subroutine cggbal(job, n, a, lda, b, ldb, ilo, ihi, lscale, rscale, work, info)
CGGBAL
Definition cggbal.f:177
subroutine cgghd3(compq, compz, n, ilo, ihi, a, lda, b, ldb, q, ldq, z, ldz, work, lwork, info)
CGGHD3
Definition cgghd3.f:231
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.
Definition clacpy.f:103
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
real function clange(norm, m, n, a, lda, work)
CLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition clange.f:115
recursive subroutine claqz0(wants, wantq, wantz, n, ilo, ihi, a, lda, b, ldb, alpha, beta, q, ldq, z, ldz, work, lwork, rwork, rec, info)
CLAQZ0
Definition claqz0.f:284
subroutine clascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
CLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition clascl.f:143
subroutine claset(uplo, m, n, alpha, beta, a, lda)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition claset.f:106
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
subroutine ctgevc(side, howmny, select, n, s, lds, p, ldp, vl, ldvl, vr, ldvr, mm, m, work, rwork, info)
CTGEVC
Definition ctgevc.f:219
subroutine cungqr(m, n, k, a, lda, tau, work, lwork, info)
CUNGQR
Definition cungqr.f:128
subroutine cunmqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
CUNMQR
Definition cunmqr.f:168
Here is the call graph for this function:
Here is the caller graph for this function: