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

◆ sggev()

subroutine sggev ( character jobvl,
character jobvr,
integer n,
real, dimension( lda, * ) a,
integer lda,
real, dimension( ldb, * ) b,
integer ldb,
real, dimension( * ) alphar,
real, dimension( * ) alphai,
real, dimension( * ) beta,
real, dimension( ldvl, * ) vl,
integer ldvl,
real, dimension( ldvr, * ) vr,
integer ldvr,
real, dimension( * ) work,
integer lwork,
integer info )

SGGEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices

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

Purpose:
!>
!> SGGEV computes for a pair of N-by-N real 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 eigenvector v(j) corresponding to the eigenvalue lambda(j)
!> of (A,B) satisfies
!>
!>                  A * v(j) = lambda(j) * B * v(j).
!>
!> The left eigenvector u(j) corresponding to the eigenvalue 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 REAL 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 REAL 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]ALPHAR
!>          ALPHAR is REAL array, dimension (N)
!>
[out]ALPHAI
!>          ALPHAI is REAL array, dimension (N)
!>
[out]BETA
!>          BETA is REAL array, dimension (N)
!>          On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will
!>          be the generalized eigenvalues.  If ALPHAI(j) is zero, then
!>          the j-th eigenvalue is real; if positive, then the j-th and
!>          (j+1)-st eigenvalues are a complex conjugate pair, with
!>          ALPHAI(j+1) negative.
!>
!>          Note: the quotients ALPHAR(j)/BETA(j) and ALPHAI(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, ALPHAR and ALPHAI 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 REAL array, dimension (LDVL,N)
!>          If JOBVL = 'V', the left eigenvectors u(j) are stored one
!>          after another in the columns of VL, in the same order as
!>          their eigenvalues. If the j-th eigenvalue is real, then
!>          u(j) = VL(:,j), the j-th column of VL. If the j-th and
!>          (j+1)-th eigenvalues form a complex conjugate pair, then
!>          u(j) = VL(:,j)+i*VL(:,j+1) and u(j+1) = VL(:,j)-i*VL(:,j+1).
!>          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 REAL array, dimension (LDVR,N)
!>          If JOBVR = 'V', the right eigenvectors v(j) are stored one
!>          after another in the columns of VR, in the same order as
!>          their eigenvalues. If the j-th eigenvalue is real, then
!>          v(j) = VR(:,j), the j-th column of VR. If the j-th and
!>          (j+1)-th eigenvalues form a complex conjugate pair, then
!>          v(j) = VR(:,j)+i*VR(:,j+1) and v(j+1) = VR(:,j)-i*VR(:,j+1).
!>          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 REAL 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,8*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]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 ALPHAR(j), ALPHAI(j), and BETA(j)
!>                should be correct for j=INFO+1,...,N.
!>          > N:  =N+1: other than QZ iteration failed in SHGEQZ.
!>                =N+2: error return from STGEVC.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 222 of file sggev.f.

225*
226* -- LAPACK driver routine --
227* -- LAPACK is a software package provided by Univ. of Tennessee, --
228* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
229*
230* .. Scalar Arguments ..
231 CHARACTER JOBVL, JOBVR
232 INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, N
233* ..
234* .. Array Arguments ..
235 REAL A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
236 $ B( LDB, * ), BETA( * ), VL( LDVL, * ),
237 $ VR( LDVR, * ), WORK( * )
238* ..
239*
240* =====================================================================
241*
242* .. Parameters ..
243 REAL ZERO, ONE
244 parameter( zero = 0.0e+0, one = 1.0e+0 )
245* ..
246* .. Local Scalars ..
247 LOGICAL ILASCL, ILBSCL, ILV, ILVL, ILVR, LQUERY
248 CHARACTER CHTEMP
249 INTEGER ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT, ILO,
250 $ IN, IRIGHT, IROWS, ITAU, IWRK, JC, JR, MAXWRK,
251 $ MINWRK
252 REAL ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
253 $ SMLNUM, TEMP
254* ..
255* .. Local Arrays ..
256 LOGICAL LDUMMA( 1 )
257* ..
258* .. External Subroutines ..
259 EXTERNAL sgeqrf, sggbak, sggbal, sgghrd, shgeqz,
260 $ slacpy,
261 $ slascl, slaset, sorgqr, sormqr, stgevc, xerbla
262* ..
263* .. External Functions ..
264 LOGICAL LSAME
265 INTEGER ILAENV
266 REAL SLAMCH, SLANGE, SROUNDUP_LWORK
267 EXTERNAL lsame, ilaenv, slamch,
268 $ slange, sroundup_lwork
269* ..
270* .. Intrinsic Functions ..
271 INTRINSIC abs, max, sqrt
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 = -12
316 ELSE IF( ldvr.LT.1 .OR. ( ilvr .AND. ldvr.LT.n ) ) THEN
317 info = -14
318 END IF
319*
320* Compute workspace
321* (Note: Comments in the code beginning "Workspace:" describe the
322* minimal amount of workspace needed at that point in the code,
323* as well as the preferred amount for good performance.
324* NB refers to the optimal block size for the immediately
325* following subroutine, as returned by ILAENV. The workspace is
326* computed assuming ILO = 1 and IHI = N, the worst case.)
327*
328 IF( info.EQ.0 ) THEN
329 minwrk = max( 1, 8*n )
330 maxwrk = max( 1, n*( 7 +
331 $ ilaenv( 1, 'SGEQRF', ' ', n, 1, n, 0 ) ) )
332 maxwrk = max( maxwrk, n*( 7 +
333 $ ilaenv( 1, 'SORMQR', ' ', n, 1, n, 0 ) ) )
334 IF( ilvl ) THEN
335 maxwrk = max( maxwrk, n*( 7 +
336 $ ilaenv( 1, 'SORGQR', ' ', n, 1, n, -1 ) ) )
337 END IF
338 work( 1 ) = sroundup_lwork(maxwrk)
339*
340 IF( lwork.LT.minwrk .AND. .NOT.lquery )
341 $ info = -16
342 END IF
343*
344 IF( info.NE.0 ) THEN
345 CALL xerbla( 'SGGEV ', -info )
346 RETURN
347 ELSE IF( lquery ) THEN
348 RETURN
349 END IF
350*
351* Quick return if possible
352*
353 IF( n.EQ.0 )
354 $ RETURN
355*
356* Get machine constants
357*
358 eps = slamch( 'P' )
359 smlnum = slamch( 'S' )
360 bignum = one / smlnum
361 smlnum = sqrt( smlnum ) / eps
362 bignum = one / smlnum
363*
364* Scale A if max element outside range [SMLNUM,BIGNUM]
365*
366 anrm = slange( 'M', n, n, a, lda, work )
367 ilascl = .false.
368 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
369 anrmto = smlnum
370 ilascl = .true.
371 ELSE IF( anrm.GT.bignum ) THEN
372 anrmto = bignum
373 ilascl = .true.
374 END IF
375 IF( ilascl )
376 $ CALL slascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
377*
378* Scale B if max element outside range [SMLNUM,BIGNUM]
379*
380 bnrm = slange( 'M', n, n, b, ldb, work )
381 ilbscl = .false.
382 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
383 bnrmto = smlnum
384 ilbscl = .true.
385 ELSE IF( bnrm.GT.bignum ) THEN
386 bnrmto = bignum
387 ilbscl = .true.
388 END IF
389 IF( ilbscl )
390 $ CALL slascl( 'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
391*
392* Permute the matrices A, B to isolate eigenvalues if possible
393* (Workspace: need 6*N)
394*
395 ileft = 1
396 iright = n + 1
397 iwrk = iright + n
398 CALL sggbal( 'P', n, a, lda, b, ldb, ilo, ihi, work( ileft ),
399 $ work( iright ), work( iwrk ), ierr )
400*
401* Reduce B to triangular form (QR decomposition of B)
402* (Workspace: need N, prefer N*NB)
403*
404 irows = ihi + 1 - ilo
405 IF( ilv ) THEN
406 icols = n + 1 - ilo
407 ELSE
408 icols = irows
409 END IF
410 itau = iwrk
411 iwrk = itau + irows
412 CALL sgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
413 $ work( iwrk ), lwork+1-iwrk, ierr )
414*
415* Apply the orthogonal transformation to matrix A
416* (Workspace: need N, prefer N*NB)
417*
418 CALL sormqr( 'L', 'T', irows, icols, irows, b( ilo, ilo ), ldb,
419 $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
420 $ lwork+1-iwrk, ierr )
421*
422* Initialize VL
423* (Workspace: need N, prefer N*NB)
424*
425 IF( ilvl ) THEN
426 CALL slaset( 'Full', n, n, zero, one, vl, ldvl )
427 IF( irows.GT.1 ) THEN
428 CALL slacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
429 $ vl( ilo+1, ilo ), ldvl )
430 END IF
431 CALL sorgqr( irows, irows, irows, vl( ilo, ilo ), ldvl,
432 $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
433 END IF
434*
435* Initialize VR
436*
437 IF( ilvr )
438 $ CALL slaset( 'Full', n, n, zero, one, vr, ldvr )
439*
440* Reduce to generalized Hessenberg form
441* (Workspace: none needed)
442*
443 IF( ilv ) THEN
444*
445* Eigenvectors requested -- work on whole matrix.
446*
447 CALL sgghrd( jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb, vl,
448 $ ldvl, vr, ldvr, ierr )
449 ELSE
450 CALL sgghrd( 'N', 'N', irows, 1, irows, a( ilo, ilo ), lda,
451 $ b( ilo, ilo ), ldb, vl, ldvl, vr, ldvr, ierr )
452 END IF
453*
454* Perform QZ algorithm (Compute eigenvalues, and optionally, the
455* Schur forms and Schur vectors)
456* (Workspace: need N)
457*
458 iwrk = itau
459 IF( ilv ) THEN
460 chtemp = 'S'
461 ELSE
462 chtemp = 'E'
463 END IF
464 CALL shgeqz( chtemp, jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb,
465 $ alphar, alphai, beta, vl, ldvl, vr, ldvr,
466 $ work( iwrk ), lwork+1-iwrk, ierr )
467 IF( ierr.NE.0 ) THEN
468 IF( ierr.GT.0 .AND. ierr.LE.n ) THEN
469 info = ierr
470 ELSE IF( ierr.GT.n .AND. ierr.LE.2*n ) THEN
471 info = ierr - n
472 ELSE
473 info = n + 1
474 END IF
475 GO TO 110
476 END IF
477*
478* Compute Eigenvectors
479* (Workspace: need 6*N)
480*
481 IF( ilv ) THEN
482 IF( ilvl ) THEN
483 IF( ilvr ) THEN
484 chtemp = 'B'
485 ELSE
486 chtemp = 'L'
487 END IF
488 ELSE
489 chtemp = 'R'
490 END IF
491 CALL stgevc( chtemp, 'B', ldumma, n, a, lda, b, ldb, vl,
492 $ ldvl,
493 $ vr, ldvr, n, in, work( iwrk ), ierr )
494 IF( ierr.NE.0 ) THEN
495 info = n + 2
496 GO TO 110
497 END IF
498*
499* Undo balancing on VL and VR and normalization
500* (Workspace: none needed)
501*
502 IF( ilvl ) THEN
503 CALL sggbak( 'P', 'L', n, ilo, ihi, work( ileft ),
504 $ work( iright ), n, vl, ldvl, ierr )
505 DO 50 jc = 1, n
506 IF( alphai( jc ).LT.zero )
507 $ GO TO 50
508 temp = zero
509 IF( alphai( jc ).EQ.zero ) THEN
510 DO 10 jr = 1, n
511 temp = max( temp, abs( vl( jr, jc ) ) )
512 10 CONTINUE
513 ELSE
514 DO 20 jr = 1, n
515 temp = max( temp, abs( vl( jr, jc ) )+
516 $ abs( vl( jr, jc+1 ) ) )
517 20 CONTINUE
518 END IF
519 IF( temp.LT.smlnum )
520 $ GO TO 50
521 temp = one / temp
522 IF( alphai( jc ).EQ.zero ) THEN
523 DO 30 jr = 1, n
524 vl( jr, jc ) = vl( jr, jc )*temp
525 30 CONTINUE
526 ELSE
527 DO 40 jr = 1, n
528 vl( jr, jc ) = vl( jr, jc )*temp
529 vl( jr, jc+1 ) = vl( jr, jc+1 )*temp
530 40 CONTINUE
531 END IF
532 50 CONTINUE
533 END IF
534 IF( ilvr ) THEN
535 CALL sggbak( 'P', 'R', n, ilo, ihi, work( ileft ),
536 $ work( iright ), n, vr, ldvr, ierr )
537 DO 100 jc = 1, n
538 IF( alphai( jc ).LT.zero )
539 $ GO TO 100
540 temp = zero
541 IF( alphai( jc ).EQ.zero ) THEN
542 DO 60 jr = 1, n
543 temp = max( temp, abs( vr( jr, jc ) ) )
544 60 CONTINUE
545 ELSE
546 DO 70 jr = 1, n
547 temp = max( temp, abs( vr( jr, jc ) )+
548 $ abs( vr( jr, jc+1 ) ) )
549 70 CONTINUE
550 END IF
551 IF( temp.LT.smlnum )
552 $ GO TO 100
553 temp = one / temp
554 IF( alphai( jc ).EQ.zero ) THEN
555 DO 80 jr = 1, n
556 vr( jr, jc ) = vr( jr, jc )*temp
557 80 CONTINUE
558 ELSE
559 DO 90 jr = 1, n
560 vr( jr, jc ) = vr( jr, jc )*temp
561 vr( jr, jc+1 ) = vr( jr, jc+1 )*temp
562 90 CONTINUE
563 END IF
564 100 CONTINUE
565 END IF
566*
567* End of eigenvector calculation
568*
569 END IF
570*
571* Undo scaling if necessary
572*
573 110 CONTINUE
574*
575 IF( ilascl ) THEN
576 CALL slascl( 'G', 0, 0, anrmto, anrm, n, 1, alphar, n,
577 $ ierr )
578 CALL slascl( 'G', 0, 0, anrmto, anrm, n, 1, alphai, n,
579 $ ierr )
580 END IF
581*
582 IF( ilbscl ) THEN
583 CALL slascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
584 END IF
585*
586 work( 1 ) = sroundup_lwork(maxwrk)
587 RETURN
588*
589* End of SGGEV
590*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
sgeqrf
subroutine sgeqrf(m, n, a, lda, tau, work, lwork, info)
SGEQRF
Definition sgeqrf.f:144
sggbak
subroutine sggbak(job, side, n, ilo, ihi, lscale, rscale, m, v, ldv, info)
SGGBAK
Definition sggbak.f:146
sggbal
subroutine sggbal(job, n, a, lda, b, ldb, ilo, ihi, lscale, rscale, work, info)
SGGBAL
Definition sggbal.f:175
sgghrd
subroutine sgghrd(compq, compz, n, ilo, ihi, a, lda, b, ldb, q, ldq, z, ldz, info)
SGGHRD
Definition sgghrd.f:206
shgeqz
subroutine shgeqz(job, compq, compz, n, ilo, ihi, h, ldh, t, ldt, alphar, alphai, beta, q, ldq, z, ldz, work, lwork, info)
SHGEQZ
Definition shgeqz.f:303
ilaenv
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:160
slacpy
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
Definition slacpy.f:101
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
slange
real function slange(norm, m, n, a, lda, work)
SLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition slange.f:112
slascl
subroutine slascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition slascl.f:142
slaset
subroutine slaset(uplo, m, n, alpha, beta, a, lda)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition slaset.f:108
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
sroundup_lwork
real function sroundup_lwork(lwork)
SROUNDUP_LWORK
Definition sroundup_lwork.f:59
stgevc
subroutine stgevc(side, howmny, select, n, s, lds, p, ldp, vl, ldvl, vr, ldvr, mm, m, work, info)
STGEVC
Definition stgevc.f:293
sorgqr
subroutine sorgqr(m, n, k, a, lda, tau, work, lwork, info)
SORGQR
Definition sorgqr.f:126
sormqr
subroutine sormqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
SORMQR
Definition sormqr.f:166
Here is the call graph for this function:
Here is the caller graph for this function: