LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ sggev3()

subroutine sggev3 ( 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 
)

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

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

Purpose:
 SGGEV3 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

          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 SLAQZ0.
                =N+2: error return from STGEVC.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 223 of file sggev3.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, LWKOPT
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, sgghd3, slaqz0, slabad,
261  $ xerbla
262 * ..
263 * .. External Functions ..
264  LOGICAL LSAME
265  REAL SLAMCH, SLANGE
266  EXTERNAL lsame, slamch, slange
267 * ..
268 * .. Intrinsic Functions ..
269  INTRINSIC abs, max, sqrt
270 * ..
271 * .. Executable Statements ..
272 *
273 * Decode the input arguments
274 *
275  IF( lsame( jobvl, 'N' ) ) THEN
276  ijobvl = 1
277  ilvl = .false.
278  ELSE IF( lsame( jobvl, 'V' ) ) THEN
279  ijobvl = 2
280  ilvl = .true.
281  ELSE
282  ijobvl = -1
283  ilvl = .false.
284  END IF
285 *
286  IF( lsame( jobvr, 'N' ) ) THEN
287  ijobvr = 1
288  ilvr = .false.
289  ELSE IF( lsame( jobvr, 'V' ) ) THEN
290  ijobvr = 2
291  ilvr = .true.
292  ELSE
293  ijobvr = -1
294  ilvr = .false.
295  END IF
296  ilv = ilvl .OR. ilvr
297 *
298 * Test the input arguments
299 *
300  info = 0
301  lquery = ( lwork.EQ.-1 )
302  IF( ijobvl.LE.0 ) THEN
303  info = -1
304  ELSE IF( ijobvr.LE.0 ) THEN
305  info = -2
306  ELSE IF( n.LT.0 ) THEN
307  info = -3
308  ELSE IF( lda.LT.max( 1, n ) ) THEN
309  info = -5
310  ELSE IF( ldb.LT.max( 1, n ) ) THEN
311  info = -7
312  ELSE IF( ldvl.LT.1 .OR. ( ilvl .AND. ldvl.LT.n ) ) THEN
313  info = -12
314  ELSE IF( ldvr.LT.1 .OR. ( ilvr .AND. ldvr.LT.n ) ) THEN
315  info = -14
316  ELSE IF( lwork.LT.max( 1, 8*n ) .AND. .NOT.lquery ) THEN
317  info = -16
318  END IF
319 *
320 * Compute workspace
321 *
322  IF( info.EQ.0 ) THEN
323  CALL sgeqrf( n, n, b, ldb, work, work, -1, ierr )
324  lwkopt = max( 1, 8*n, 3*n+int( work( 1 ) ) )
325  CALL sormqr( 'L', 'T', n, n, n, b, ldb, work, a, lda, work,
326  $ -1, ierr )
327  lwkopt = max( lwkopt, 3*n+int( work( 1 ) ) )
328  CALL sgghd3( jobvl, jobvr, n, 1, n, a, lda, b, ldb, vl, ldvl,
329  $ vr, ldvr, work, -1, ierr )
330  lwkopt = max( lwkopt, 3*n+int( work( 1 ) ) )
331  IF( ilvl ) THEN
332  CALL sorgqr( n, n, n, vl, ldvl, work, work, -1, ierr )
333  lwkopt = max( lwkopt, 3*n+int( work( 1 ) ) )
334  CALL slaqz0( 'S', jobvl, jobvr, n, 1, n, a, lda, b, ldb,
335  $ alphar, alphai, beta, vl, ldvl, vr, ldvr,
336  $ work, -1, 0, ierr )
337  lwkopt = max( lwkopt, 2*n+int( work( 1 ) ) )
338  ELSE
339  CALL slaqz0( 'E', jobvl, jobvr, n, 1, n, a, lda, b, ldb,
340  $ alphar, alphai, beta, vl, ldvl, vr, ldvr,
341  $ work, -1, 0, ierr )
342  lwkopt = max( lwkopt, 2*n+int( work( 1 ) ) )
343  END IF
344  work( 1 ) = real( lwkopt )
345 *
346  END IF
347 *
348  IF( info.NE.0 ) THEN
349  CALL xerbla( 'SGGEV3 ', -info )
350  RETURN
351  ELSE IF( lquery ) THEN
352  RETURN
353  END IF
354 *
355 * Quick return if possible
356 *
357  IF( n.EQ.0 )
358  $ RETURN
359 *
360 * Get machine constants
361 *
362  eps = slamch( 'P' )
363  smlnum = slamch( 'S' )
364  bignum = one / smlnum
365  CALL slabad( smlnum, bignum )
366  smlnum = sqrt( smlnum ) / eps
367  bignum = one / smlnum
368 *
369 * Scale A if max element outside range [SMLNUM,BIGNUM]
370 *
371  anrm = slange( 'M', n, n, a, lda, work )
372  ilascl = .false.
373  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
374  anrmto = smlnum
375  ilascl = .true.
376  ELSE IF( anrm.GT.bignum ) THEN
377  anrmto = bignum
378  ilascl = .true.
379  END IF
380  IF( ilascl )
381  $ CALL slascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
382 *
383 * Scale B if max element outside range [SMLNUM,BIGNUM]
384 *
385  bnrm = slange( 'M', n, n, b, ldb, work )
386  ilbscl = .false.
387  IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
388  bnrmto = smlnum
389  ilbscl = .true.
390  ELSE IF( bnrm.GT.bignum ) THEN
391  bnrmto = bignum
392  ilbscl = .true.
393  END IF
394  IF( ilbscl )
395  $ CALL slascl( 'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
396 *
397 * Permute the matrices A, B to isolate eigenvalues if possible
398 *
399  ileft = 1
400  iright = n + 1
401  iwrk = iright + n
402  CALL sggbal( 'P', n, a, lda, b, ldb, ilo, ihi, work( ileft ),
403  $ work( iright ), work( iwrk ), ierr )
404 *
405 * Reduce B to triangular form (QR decomposition of B)
406 *
407  irows = ihi + 1 - ilo
408  IF( ilv ) THEN
409  icols = n + 1 - ilo
410  ELSE
411  icols = irows
412  END IF
413  itau = iwrk
414  iwrk = itau + irows
415  CALL sgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
416  $ work( iwrk ), lwork+1-iwrk, ierr )
417 *
418 * Apply the orthogonal transformation to matrix A
419 *
420  CALL sormqr( 'L', 'T', irows, icols, irows, b( ilo, ilo ), ldb,
421  $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
422  $ lwork+1-iwrk, ierr )
423 *
424 * Initialize VL
425 *
426  IF( ilvl ) THEN
427  CALL slaset( 'Full', n, n, zero, one, vl, ldvl )
428  IF( irows.GT.1 ) THEN
429  CALL slacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
430  $ vl( ilo+1, ilo ), ldvl )
431  END IF
432  CALL sorgqr( irows, irows, irows, vl( ilo, ilo ), ldvl,
433  $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
434  END IF
435 *
436 * Initialize VR
437 *
438  IF( ilvr )
439  $ CALL slaset( 'Full', n, n, zero, one, vr, ldvr )
440 *
441 * Reduce to generalized Hessenberg form
442 *
443  IF( ilv ) THEN
444 *
445 * Eigenvectors requested -- work on whole matrix.
446 *
447  CALL sgghd3( jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb, vl,
448  $ ldvl, vr, ldvr, work( iwrk ), lwork+1-iwrk, ierr )
449  ELSE
450  CALL sgghd3( 'N', 'N', irows, 1, irows, a( ilo, ilo ), lda,
451  $ b( ilo, ilo ), ldb, vl, ldvl, vr, ldvr,
452  $ work( iwrk ), lwork+1-iwrk, ierr )
453  END IF
454 *
455 * Perform QZ algorithm (Compute eigenvalues, and optionally, the
456 * Schur forms and Schur vectors)
457 *
458  iwrk = itau
459  IF( ilv ) THEN
460  chtemp = 'S'
461  ELSE
462  chtemp = 'E'
463  END IF
464  CALL slaqz0( chtemp, jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb,
465  $ alphar, alphai, beta, vl, ldvl, vr, ldvr,
466  $ work( iwrk ), lwork+1-iwrk, 0, 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 *
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 *
499  IF( ilvl ) THEN
500  CALL sggbak( 'P', 'L', n, ilo, ihi, work( ileft ),
501  $ work( iright ), n, vl, ldvl, ierr )
502  DO 50 jc = 1, n
503  IF( alphai( jc ).LT.zero )
504  $ GO TO 50
505  temp = zero
506  IF( alphai( jc ).EQ.zero ) THEN
507  DO 10 jr = 1, n
508  temp = max( temp, abs( vl( jr, jc ) ) )
509  10 CONTINUE
510  ELSE
511  DO 20 jr = 1, n
512  temp = max( temp, abs( vl( jr, jc ) )+
513  $ abs( vl( jr, jc+1 ) ) )
514  20 CONTINUE
515  END IF
516  IF( temp.LT.smlnum )
517  $ GO TO 50
518  temp = one / temp
519  IF( alphai( jc ).EQ.zero ) THEN
520  DO 30 jr = 1, n
521  vl( jr, jc ) = vl( jr, jc )*temp
522  30 CONTINUE
523  ELSE
524  DO 40 jr = 1, n
525  vl( jr, jc ) = vl( jr, jc )*temp
526  vl( jr, jc+1 ) = vl( jr, jc+1 )*temp
527  40 CONTINUE
528  END IF
529  50 CONTINUE
530  END IF
531  IF( ilvr ) THEN
532  CALL sggbak( 'P', 'R', n, ilo, ihi, work( ileft ),
533  $ work( iright ), n, vr, ldvr, ierr )
534  DO 100 jc = 1, n
535  IF( alphai( jc ).LT.zero )
536  $ GO TO 100
537  temp = zero
538  IF( alphai( jc ).EQ.zero ) THEN
539  DO 60 jr = 1, n
540  temp = max( temp, abs( vr( jr, jc ) ) )
541  60 CONTINUE
542  ELSE
543  DO 70 jr = 1, n
544  temp = max( temp, abs( vr( jr, jc ) )+
545  $ abs( vr( jr, jc+1 ) ) )
546  70 CONTINUE
547  END IF
548  IF( temp.LT.smlnum )
549  $ GO TO 100
550  temp = one / temp
551  IF( alphai( jc ).EQ.zero ) THEN
552  DO 80 jr = 1, n
553  vr( jr, jc ) = vr( jr, jc )*temp
554  80 CONTINUE
555  ELSE
556  DO 90 jr = 1, n
557  vr( jr, jc ) = vr( jr, jc )*temp
558  vr( jr, jc+1 ) = vr( jr, jc+1 )*temp
559  90 CONTINUE
560  END IF
561  100 CONTINUE
562  END IF
563 *
564 * End of eigenvector calculation
565 *
566  END IF
567 *
568 * Undo scaling if necessary
569 *
570  110 CONTINUE
571 *
572  IF( ilascl ) THEN
573  CALL slascl( 'G', 0, 0, anrmto, anrm, n, 1, alphar, n, ierr )
574  CALL slascl( 'G', 0, 0, anrmto, anrm, n, 1, alphai, n, ierr )
575  END IF
576 *
577  IF( ilbscl ) THEN
578  CALL slascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
579  END IF
580 *
581  work( 1 ) = real( lwkopt )
582  RETURN
583 *
584 * End of SGGEV3
585 *
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:74
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
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
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
recursive subroutine slaqz0(WANTS, WANTQ, WANTZ, N, ILO, IHI, A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, REC, INFO)
SLAQZ0
Definition: slaqz0.f:304
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
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 stgevc(SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL, LDVL, VR, LDVR, MM, M, WORK, INFO)
STGEVC
Definition: stgevc.f:295
subroutine sgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
SGEQRF
Definition: sgeqrf.f:145
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
subroutine sgghd3(COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, LDQ, Z, LDZ, WORK, LWORK, INFO)
SGGHD3
Definition: sgghd3.f:230
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: