LAPACK 3.12.0
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 224 of file sggev.f.

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