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

◆ 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
!>          The dimension of the array WORK. LWORK >= MAX(1,8*N).
!>          For good performance, LWORK should 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 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 $ LWKMIN
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,
261 $ sgghd3, slaqz0, slacpy,
262 $ slascl, slaset, sorgqr,
263 $ sormqr, stgevc
264* ..
265* .. External Functions ..
266 LOGICAL LSAME
267 REAL SLAMCH, SLANGE, SROUNDUP_LWORK
268 EXTERNAL lsame, slamch, slange,
269 $ sroundup_lwork
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 lwkmin = max( 1, 8*n )
306 IF( ijobvl.LE.0 ) THEN
307 info = -1
308 ELSE IF( ijobvr.LE.0 ) THEN
309 info = -2
310 ELSE IF( n.LT.0 ) THEN
311 info = -3
312 ELSE IF( lda.LT.max( 1, n ) ) THEN
313 info = -5
314 ELSE IF( ldb.LT.max( 1, n ) ) THEN
315 info = -7
316 ELSE IF( ldvl.LT.1 .OR. ( ilvl .AND. ldvl.LT.n ) ) THEN
317 info = -12
318 ELSE IF( ldvr.LT.1 .OR. ( ilvr .AND. ldvr.LT.n ) ) THEN
319 info = -14
320 ELSE IF( lwork.LT.lwkmin .AND. .NOT.lquery ) THEN
321 info = -16
322 END IF
323*
324* Compute workspace
325*
326 IF( info.EQ.0 ) THEN
327 CALL sgeqrf( n, n, b, ldb, work, work, -1, ierr )
328 lwkopt = max( lwkmin, 3*n+int( work( 1 ) ) )
329 CALL sormqr( 'L', 'T', n, n, n, b, ldb, work, a, lda, work,
330 $ -1, ierr )
331 lwkopt = max( lwkopt, 3*n+int( work( 1 ) ) )
332 CALL sgghd3( jobvl, jobvr, n, 1, n, a, lda,
333 $ b, ldb, vl, ldvl,
334 $ vr, ldvr, work, -1, ierr )
335 lwkopt = max( lwkopt, 3*n+int( work( 1 ) ) )
336 IF( ilvl ) THEN
337 CALL sorgqr( n, n, n, vl, ldvl, work, work, -1, ierr )
338 lwkopt = max( lwkopt, 3*n+int( work( 1 ) ) )
339 CALL slaqz0( 'S', 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 ELSE
344 CALL slaqz0( 'E', jobvl, jobvr, n, 1, n, a, lda, b, ldb,
345 $ alphar, alphai, beta, vl, ldvl, vr, ldvr,
346 $ work, -1, 0, ierr )
347 lwkopt = max( lwkopt, 2*n+int( work( 1 ) ) )
348 END IF
349 IF( n.EQ.0 ) THEN
350 work( 1 ) = 1
351 ELSE
352 work( 1 ) = sroundup_lwork( lwkopt )
353 END IF
354 END IF
355*
356 IF( info.NE.0 ) THEN
357 CALL xerbla( 'SGGEV3 ', -info )
358 RETURN
359 ELSE IF( lquery ) THEN
360 RETURN
361 END IF
362*
363* Quick return if possible
364*
365 IF( n.EQ.0 )
366 $ RETURN
367*
368* Get machine constants
369*
370 eps = slamch( 'P' )
371 smlnum = slamch( 'S' )
372 bignum = one / smlnum
373 smlnum = sqrt( smlnum ) / eps
374 bignum = one / smlnum
375*
376* Scale A if max element outside range [SMLNUM,BIGNUM]
377*
378 anrm = slange( 'M', n, n, a, lda, work )
379 ilascl = .false.
380 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
381 anrmto = smlnum
382 ilascl = .true.
383 ELSE IF( anrm.GT.bignum ) THEN
384 anrmto = bignum
385 ilascl = .true.
386 END IF
387 IF( ilascl )
388 $ CALL slascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
389*
390* Scale B if max element outside range [SMLNUM,BIGNUM]
391*
392 bnrm = slange( 'M', n, n, b, ldb, work )
393 ilbscl = .false.
394 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
395 bnrmto = smlnum
396 ilbscl = .true.
397 ELSE IF( bnrm.GT.bignum ) THEN
398 bnrmto = bignum
399 ilbscl = .true.
400 END IF
401 IF( ilbscl )
402 $ CALL slascl( 'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
403*
404* Permute the matrices A, B to isolate eigenvalues if possible
405*
406 ileft = 1
407 iright = n + 1
408 iwrk = iright + n
409 CALL sggbal( 'P', n, a, lda, b, ldb, ilo, ihi, work( ileft ),
410 $ work( iright ), work( iwrk ), ierr )
411*
412* Reduce B to triangular form (QR decomposition of B)
413*
414 irows = ihi + 1 - ilo
415 IF( ilv ) THEN
416 icols = n + 1 - ilo
417 ELSE
418 icols = irows
419 END IF
420 itau = iwrk
421 iwrk = itau + irows
422 CALL sgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
423 $ work( iwrk ), lwork+1-iwrk, ierr )
424*
425* Apply the orthogonal transformation to matrix A
426*
427 CALL sormqr( 'L', 'T', irows, icols, irows, b( ilo, ilo ), ldb,
428 $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
429 $ lwork+1-iwrk, ierr )
430*
431* Initialize VL
432*
433 IF( ilvl ) THEN
434 CALL slaset( 'Full', n, n, zero, one, vl, ldvl )
435 IF( irows.GT.1 ) THEN
436 CALL slacpy( 'L', irows-1, irows-1,
437 $ b( ilo+1, ilo ), ldb,
438 $ vl( ilo+1, ilo ), ldvl )
439 END IF
440 CALL sorgqr( irows, irows, irows, vl( ilo, ilo ), ldvl,
441 $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
442 END IF
443*
444* Initialize VR
445*
446 IF( ilvr )
447 $ CALL slaset( 'Full', n, n, zero, one, vr, ldvr )
448*
449* Reduce to generalized Hessenberg form
450*
451 IF( ilv ) THEN
452*
453* Eigenvectors requested -- work on whole matrix.
454*
455 CALL sgghd3( jobvl, jobvr, n, ilo, ihi,
456 $ a, lda,
457 $ b, ldb, vl,
458 $ ldvl, vr, ldvr,
459 $ work( iwrk ), lwork+1-iwrk, ierr )
460 ELSE
461 CALL sgghd3( 'N', 'N', irows, 1, irows,
462 $ a( ilo, ilo ), lda,
463 $ b( ilo, ilo ), ldb, vl,
464 $ ldvl, vr, ldvr,
465 $ work( iwrk ), lwork+1-iwrk, ierr )
466 END IF
467*
468* Perform QZ algorithm (Compute eigenvalues, and optionally, the
469* Schur forms and Schur vectors)
470*
471 iwrk = itau
472 IF( ilv ) THEN
473 chtemp = 'S'
474 ELSE
475 chtemp = 'E'
476 END IF
477 CALL slaqz0( chtemp, jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb,
478 $ alphar, alphai, beta, vl, ldvl, vr, ldvr,
479 $ work( iwrk ), lwork+1-iwrk, 0, ierr )
480 IF( ierr.NE.0 ) THEN
481 IF( ierr.GT.0 .AND. ierr.LE.n ) THEN
482 info = ierr
483 ELSE IF( ierr.GT.n .AND. ierr.LE.2*n ) THEN
484 info = ierr - n
485 ELSE
486 info = n + 1
487 END IF
488 GO TO 110
489 END IF
490*
491* Compute Eigenvectors
492*
493 IF( ilv ) THEN
494 IF( ilvl ) THEN
495 IF( ilvr ) THEN
496 chtemp = 'B'
497 ELSE
498 chtemp = 'L'
499 END IF
500 ELSE
501 chtemp = 'R'
502 END IF
503 CALL stgevc( chtemp, 'B', ldumma, n, a, lda, b, ldb, vl,
504 $ ldvl,
505 $ vr, ldvr, n, in, work( iwrk ), ierr )
506 IF( ierr.NE.0 ) THEN
507 info = n + 2
508 GO TO 110
509 END IF
510*
511* Undo balancing on VL and VR and normalization
512*
513 IF( ilvl ) THEN
514 CALL sggbak( 'P', 'L', n, ilo, ihi, work( ileft ),
515 $ work( iright ), n, vl, ldvl, ierr )
516 DO 50 jc = 1, n
517 IF( alphai( jc ).LT.zero )
518 $ GO TO 50
519 temp = zero
520 IF( alphai( jc ).EQ.zero ) THEN
521 DO 10 jr = 1, n
522 temp = max( temp, abs( vl( jr, jc ) ) )
523 10 CONTINUE
524 ELSE
525 DO 20 jr = 1, n
526 temp = max( temp, abs( vl( jr, jc ) )+
527 $ abs( vl( jr, jc+1 ) ) )
528 20 CONTINUE
529 END IF
530 IF( temp.LT.smlnum )
531 $ GO TO 50
532 temp = one / temp
533 IF( alphai( jc ).EQ.zero ) THEN
534 DO 30 jr = 1, n
535 vl( jr, jc ) = vl( jr, jc )*temp
536 30 CONTINUE
537 ELSE
538 DO 40 jr = 1, n
539 vl( jr, jc ) = vl( jr, jc )*temp
540 vl( jr, jc+1 ) = vl( jr, jc+1 )*temp
541 40 CONTINUE
542 END IF
543 50 CONTINUE
544 END IF
545 IF( ilvr ) THEN
546 CALL sggbak( 'P', 'R', n, ilo, ihi, work( ileft ),
547 $ work( iright ), n, vr, ldvr, ierr )
548 DO 100 jc = 1, n
549 IF( alphai( jc ).LT.zero )
550 $ GO TO 100
551 temp = zero
552 IF( alphai( jc ).EQ.zero ) THEN
553 DO 60 jr = 1, n
554 temp = max( temp, abs( vr( jr, jc ) ) )
555 60 CONTINUE
556 ELSE
557 DO 70 jr = 1, n
558 temp = max( temp, abs( vr( jr, jc ) )+
559 $ abs( vr( jr, jc+1 ) ) )
560 70 CONTINUE
561 END IF
562 IF( temp.LT.smlnum )
563 $ GO TO 100
564 temp = one / temp
565 IF( alphai( jc ).EQ.zero ) THEN
566 DO 80 jr = 1, n
567 vr( jr, jc ) = vr( jr, jc )*temp
568 80 CONTINUE
569 ELSE
570 DO 90 jr = 1, n
571 vr( jr, jc ) = vr( jr, jc )*temp
572 vr( jr, jc+1 ) = vr( jr, jc+1 )*temp
573 90 CONTINUE
574 END IF
575 100 CONTINUE
576 END IF
577*
578* End of eigenvector calculation
579*
580 END IF
581*
582* Undo scaling if necessary
583*
584 110 CONTINUE
585*
586 IF( ilascl ) THEN
587 CALL slascl( 'G', 0, 0, anrmto, anrm, n, 1, alphar, n,
588 $ ierr )
589 CALL slascl( 'G', 0, 0, anrmto, anrm, n, 1, alphai, n,
590 $ ierr )
591 END IF
592*
593 IF( ilbscl ) THEN
594 CALL slascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
595 END IF
596*
597 work( 1 ) = sroundup_lwork( lwkopt )
598 RETURN
599*
600* End of SGGEV3
601*
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
sgghd3
subroutine sgghd3(compq, compz, n, ilo, ihi, a, lda, b, ldb, q, ldq, z, ldz, work, lwork, info)
SGGHD3
Definition sgghd3.f:229
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
slaqz0
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:303
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: