LAPACK 3.11.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, 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:146
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: