LAPACK 3.12.0
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

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