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

Definition at line 228 of file sggev3.f.

228 *
229 * -- LAPACK driver routine (version 3.6.0) --
230 * -- LAPACK is a software package provided by Univ. of Tennessee, --
231 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
232 * January 2015
233 *
234 * .. Scalar Arguments ..
235  CHARACTER jobvl, jobvr
236  INTEGER info, lda, ldb, ldvl, ldvr, lwork, n
237 * ..
238 * .. Array Arguments ..
239  REAL a( lda, * ), alphai( * ), alphar( * ),
240  $ b( ldb, * ), beta( * ), vl( ldvl, * ),
241  $ vr( ldvr, * ), work( * )
242 * ..
243 *
244 * =====================================================================
245 *
246 * .. Parameters ..
247  REAL zero, one
248  parameter( zero = 0.0e+0, one = 1.0e+0 )
249 * ..
250 * .. Local Scalars ..
251  LOGICAL ilascl, ilbscl, ilv, ilvl, ilvr, lquery
252  CHARACTER chtemp
253  INTEGER icols, ierr, ihi, ijobvl, ijobvr, ileft, ilo,
254  $ in, iright, irows, itau, iwrk, jc, jr, lwkopt
255  REAL anrm, anrmto, bignum, bnrm, bnrmto, eps,
256  $ smlnum, temp
257 * ..
258 * .. Local Arrays ..
259  LOGICAL ldumma( 1 )
260 * ..
261 * .. External Subroutines ..
262  EXTERNAL sgeqrf, sggbak, sggbal, sgghd3, shgeqz, slabad,
264  $ xerbla
265 * ..
266 * .. External Functions ..
267  LOGICAL lsame
268  REAL slamch, slange
269  EXTERNAL lsame, slamch, slange
270 * ..
271 * .. Intrinsic Functions ..
272  INTRINSIC abs, max, sqrt
273 * ..
274 * .. Executable Statements ..
275 *
276 * Decode the input arguments
277 *
278  IF( lsame( jobvl, 'N' ) ) THEN
279  ijobvl = 1
280  ilvl = .false.
281  ELSE IF( lsame( jobvl, 'V' ) ) THEN
282  ijobvl = 2
283  ilvl = .true.
284  ELSE
285  ijobvl = -1
286  ilvl = .false.
287  END IF
288 *
289  IF( lsame( jobvr, 'N' ) ) THEN
290  ijobvr = 1
291  ilvr = .false.
292  ELSE IF( lsame( jobvr, 'V' ) ) THEN
293  ijobvr = 2
294  ilvr = .true.
295  ELSE
296  ijobvr = -1
297  ilvr = .false.
298  END IF
299  ilv = ilvl .OR. ilvr
300 *
301 * Test the input arguments
302 *
303  info = 0
304  lquery = ( lwork.EQ.-1 )
305  IF( ijobvl.LE.0 ) THEN
306  info = -1
307  ELSE IF( ijobvr.LE.0 ) THEN
308  info = -2
309  ELSE IF( n.LT.0 ) THEN
310  info = -3
311  ELSE IF( lda.LT.max( 1, n ) ) THEN
312  info = -5
313  ELSE IF( ldb.LT.max( 1, n ) ) THEN
314  info = -7
315  ELSE IF( ldvl.LT.1 .OR. ( ilvl .AND. ldvl.LT.n ) ) THEN
316  info = -12
317  ELSE IF( ldvr.LT.1 .OR. ( ilvr .AND. ldvr.LT.n ) ) THEN
318  info = -14
319  ELSE IF( lwork.LT.max( 1, 8*n ) .AND. .NOT.lquery ) THEN
320  info = -16
321  END IF
322 *
323 * Compute workspace
324 *
325  IF( info.EQ.0 ) THEN
326  CALL sgeqrf( n, n, b, ldb, work, work, -1, ierr )
327  lwkopt = max( 1, 8*n, 3*n+int( work( 1 ) ) )
328  CALL sormqr( 'L', 'T', n, n, n, b, ldb, work, a, lda, work,
329  $ -1, ierr )
330  lwkopt = max( lwkopt, 3*n+int( work( 1 ) ) )
331  CALL sgghd3( jobvl, jobvr, n, 1, n, a, lda, b, ldb, vl, ldvl,
332  $ vr, ldvr, work, -1, ierr )
333  lwkopt = max( lwkopt, 3*n+int( work( 1 ) ) )
334  IF( ilvl ) THEN
335  CALL sorgqr( n, n, n, vl, ldvl, work, work, -1, ierr )
336  lwkopt = max( lwkopt, 3*n+int( work( 1 ) ) )
337  CALL shgeqz( 'S', jobvl, jobvr, n, 1, n, a, lda, b, ldb,
338  $ alphar, alphai, beta, vl, ldvl, vr, ldvr,
339  $ work, -1, ierr )
340  lwkopt = max( lwkopt, 2*n+int( work( 1 ) ) )
341  ELSE
342  CALL shgeqz( 'E', jobvl, jobvr, n, 1, n, a, lda, b, ldb,
343  $ alphar, alphai, beta, vl, ldvl, vr, ldvr,
344  $ work, -1, ierr )
345  lwkopt = max( lwkopt, 2*n+int( work( 1 ) ) )
346  END IF
347  work( 1 ) = REAL( lwkopt )
348 *
349  END IF
350 *
351  IF( info.NE.0 ) THEN
352  CALL xerbla( 'SGGEV3 ', -info )
353  RETURN
354  ELSE IF( lquery ) THEN
355  RETURN
356  END IF
357 *
358 * Quick return if possible
359 *
360  IF( n.EQ.0 )
361  $ RETURN
362 *
363 * Get machine constants
364 *
365  eps = slamch( 'P' )
366  smlnum = slamch( 'S' )
367  bignum = one / smlnum
368  CALL slabad( smlnum, bignum )
369  smlnum = sqrt( smlnum ) / eps
370  bignum = one / smlnum
371 *
372 * Scale A if max element outside range [SMLNUM,BIGNUM]
373 *
374  anrm = slange( 'M', n, n, a, lda, work )
375  ilascl = .false.
376  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
377  anrmto = smlnum
378  ilascl = .true.
379  ELSE IF( anrm.GT.bignum ) THEN
380  anrmto = bignum
381  ilascl = .true.
382  END IF
383  IF( ilascl )
384  $ CALL slascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
385 *
386 * Scale B if max element outside range [SMLNUM,BIGNUM]
387 *
388  bnrm = slange( 'M', n, n, b, ldb, work )
389  ilbscl = .false.
390  IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
391  bnrmto = smlnum
392  ilbscl = .true.
393  ELSE IF( bnrm.GT.bignum ) THEN
394  bnrmto = bignum
395  ilbscl = .true.
396  END IF
397  IF( ilbscl )
398  $ CALL slascl( 'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
399 *
400 * Permute the matrices A, B to isolate eigenvalues if possible
401 *
402  ileft = 1
403  iright = n + 1
404  iwrk = iright + n
405  CALL sggbal( 'P', n, a, lda, b, ldb, ilo, ihi, work( ileft ),
406  $ work( iright ), work( iwrk ), ierr )
407 *
408 * Reduce B to triangular form (QR decomposition of B)
409 *
410  irows = ihi + 1 - ilo
411  IF( ilv ) THEN
412  icols = n + 1 - ilo
413  ELSE
414  icols = irows
415  END IF
416  itau = iwrk
417  iwrk = itau + irows
418  CALL sgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
419  $ work( iwrk ), lwork+1-iwrk, ierr )
420 *
421 * Apply the orthogonal transformation to matrix A
422 *
423  CALL sormqr( 'L', 'T', irows, icols, irows, b( ilo, ilo ), ldb,
424  $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
425  $ lwork+1-iwrk, ierr )
426 *
427 * Initialize VL
428 *
429  IF( ilvl ) THEN
430  CALL slaset( 'Full', n, n, zero, one, vl, ldvl )
431  IF( irows.GT.1 ) THEN
432  CALL slacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
433  $ vl( ilo+1, ilo ), ldvl )
434  END IF
435  CALL sorgqr( irows, irows, irows, vl( ilo, ilo ), ldvl,
436  $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
437  END IF
438 *
439 * Initialize VR
440 *
441  IF( ilvr )
442  $ CALL slaset( 'Full', n, n, zero, one, vr, ldvr )
443 *
444 * Reduce to generalized Hessenberg form
445 *
446  IF( ilv ) THEN
447 *
448 * Eigenvectors requested -- work on whole matrix.
449 *
450  CALL sgghd3( jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb, vl,
451  $ ldvl, vr, ldvr, work( iwrk ), lwork+1-iwrk, ierr )
452  ELSE
453  CALL sgghd3( 'N', 'N', irows, 1, irows, a( ilo, ilo ), lda,
454  $ b( ilo, ilo ), ldb, vl, ldvl, vr, ldvr,
455  $ work( iwrk ), lwork+1-iwrk, ierr )
456  END IF
457 *
458 * Perform QZ algorithm (Compute eigenvalues, and optionally, the
459 * Schur forms and Schur vectors)
460 *
461  iwrk = itau
462  IF( ilv ) THEN
463  chtemp = 'S'
464  ELSE
465  chtemp = 'E'
466  END IF
467  CALL shgeqz( chtemp, jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb,
468  $ alphar, alphai, beta, vl, ldvl, vr, ldvr,
469  $ work( iwrk ), lwork+1-iwrk, ierr )
470  IF( ierr.NE.0 ) THEN
471  IF( ierr.GT.0 .AND. ierr.LE.n ) THEN
472  info = ierr
473  ELSE IF( ierr.GT.n .AND. ierr.LE.2*n ) THEN
474  info = ierr - n
475  ELSE
476  info = n + 1
477  END IF
478  GO TO 110
479  END IF
480 *
481 * Compute Eigenvectors
482 *
483  IF( ilv ) THEN
484  IF( ilvl ) THEN
485  IF( ilvr ) THEN
486  chtemp = 'B'
487  ELSE
488  chtemp = 'L'
489  END IF
490  ELSE
491  chtemp = 'R'
492  END IF
493  CALL stgevc( chtemp, 'B', ldumma, n, a, lda, b, ldb, vl, ldvl,
494  $ vr, ldvr, n, in, work( iwrk ), ierr )
495  IF( ierr.NE.0 ) THEN
496  info = n + 2
497  GO TO 110
498  END IF
499 *
500 * Undo balancing on VL and VR and normalization
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, ierr )
577  CALL slascl( 'G', 0, 0, anrmto, anrm, n, 1, alphai, n, ierr )
578  END IF
579 *
580  IF( ilbscl ) THEN
581  CALL slascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
582  END IF
583 *
584  work( 1 ) = REAL( lwkopt )
585  RETURN
586 *
587 * End of SGGEV3
588 *
subroutine stgevc(SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL, LDVL, VR, LDVR, MM, M, WORK, INFO)
STGEVC
Definition: stgevc.f:297
subroutine sgghd3(COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, LDQ, Z, LDZ, WORK, LWORK, INFO)
SGGHD3
Definition: sgghd3.f:232
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:306
subroutine sormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMQR
Definition: sormqr.f:170
subroutine sgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
SGEQRF
Definition: sgeqrf.f:138
subroutine sggbak(JOB, SIDE, N, ILO, IHI, LSCALE, RSCALE, M, V, LDV, INFO)
SGGBAK
Definition: sggbak.f:149
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:116
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
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:112
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:145
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:76
subroutine sorgqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
SORGQR
Definition: sorgqr.f:130
subroutine sggbal(JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE, WORK, INFO)
SGGBAL
Definition: sggbal.f:179
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:105
Here is the call graph for this function:
Here is the caller graph for this function: