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

◆ sgeesx()

subroutine sgeesx ( character  JOBVS,
character  SORT,
external  SELECT,
character  SENSE,
integer  N,
real, dimension( lda, * )  A,
integer  LDA,
integer  SDIM,
real, dimension( * )  WR,
real, dimension( * )  WI,
real, dimension( ldvs, * )  VS,
integer  LDVS,
real  RCONDE,
real  RCONDV,
real, dimension( * )  WORK,
integer  LWORK,
integer, dimension( * )  IWORK,
integer  LIWORK,
logical, dimension( * )  BWORK,
integer  INFO 
)

SGEESX computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE matrices

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

Purpose:
 SGEESX computes for an N-by-N real nonsymmetric matrix A, the
 eigenvalues, the real Schur form T, and, optionally, the matrix of
 Schur vectors Z.  This gives the Schur factorization A = Z*T*(Z**T).

 Optionally, it also orders the eigenvalues on the diagonal of the
 real Schur form so that selected eigenvalues are at the top left;
 computes a reciprocal condition number for the average of the
 selected eigenvalues (RCONDE); and computes a reciprocal condition
 number for the right invariant subspace corresponding to the
 selected eigenvalues (RCONDV).  The leading columns of Z form an
 orthonormal basis for this invariant subspace.

 For further explanation of the reciprocal condition numbers RCONDE
 and RCONDV, see Section 4.10 of the LAPACK Users' Guide (where
 these quantities are called s and sep respectively).

 A real matrix is in real Schur form if it is upper quasi-triangular
 with 1-by-1 and 2-by-2 blocks. 2-by-2 blocks will be standardized in
 the form
           [  a  b  ]
           [  c  a  ]

 where b*c < 0. The eigenvalues of such a block are a +- sqrt(bc).
Parameters
[in]JOBVS
          JOBVS is CHARACTER*1
          = 'N': Schur vectors are not computed;
          = 'V': Schur vectors are computed.
[in]SORT
          SORT is CHARACTER*1
          Specifies whether or not to order the eigenvalues on the
          diagonal of the Schur form.
          = 'N': Eigenvalues are not ordered;
          = 'S': Eigenvalues are ordered (see SELECT).
[in]SELECT
          SELECT is a LOGICAL FUNCTION of two REAL arguments
          SELECT must be declared EXTERNAL in the calling subroutine.
          If SORT = 'S', SELECT is used to select eigenvalues to sort
          to the top left of the Schur form.
          If SORT = 'N', SELECT is not referenced.
          An eigenvalue WR(j)+sqrt(-1)*WI(j) is selected if
          SELECT(WR(j),WI(j)) is true; i.e., if either one of a
          complex conjugate pair of eigenvalues is selected, then both
          are.  Note that a selected complex eigenvalue may no longer
          satisfy SELECT(WR(j),WI(j)) = .TRUE. after ordering, since
          ordering may change the value of complex eigenvalues
          (especially if the eigenvalue is ill-conditioned); in this
          case INFO may be set to N+3 (see INFO below).
[in]SENSE
          SENSE is CHARACTER*1
          Determines which reciprocal condition numbers are computed.
          = 'N': None are computed;
          = 'E': Computed for average of selected eigenvalues only;
          = 'V': Computed for selected right invariant subspace only;
          = 'B': Computed for both.
          If SENSE = 'E', 'V' or 'B', SORT must equal 'S'.
[in]N
          N is INTEGER
          The order of the matrix A. N >= 0.
[in,out]A
          A is REAL array, dimension (LDA, N)
          On entry, the N-by-N matrix A.
          On exit, A is overwritten by its real Schur form T.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[out]SDIM
          SDIM is INTEGER
          If SORT = 'N', SDIM = 0.
          If SORT = 'S', SDIM = number of eigenvalues (after sorting)
                         for which SELECT is true. (Complex conjugate
                         pairs for which SELECT is true for either
                         eigenvalue count as 2.)
[out]WR
          WR is REAL array, dimension (N)
[out]WI
          WI is REAL array, dimension (N)
          WR and WI contain the real and imaginary parts, respectively,
          of the computed eigenvalues, in the same order that they
          appear on the diagonal of the output Schur form T.  Complex
          conjugate pairs of eigenvalues appear consecutively with the
          eigenvalue having the positive imaginary part first.
[out]VS
          VS is REAL array, dimension (LDVS,N)
          If JOBVS = 'V', VS contains the orthogonal matrix Z of Schur
          vectors.
          If JOBVS = 'N', VS is not referenced.
[in]LDVS
          LDVS is INTEGER
          The leading dimension of the array VS.  LDVS >= 1, and if
          JOBVS = 'V', LDVS >= N.
[out]RCONDE
          RCONDE is REAL
          If SENSE = 'E' or 'B', RCONDE contains the reciprocal
          condition number for the average of the selected eigenvalues.
          Not referenced if SENSE = 'N' or 'V'.
[out]RCONDV
          RCONDV is REAL
          If SENSE = 'V' or 'B', RCONDV contains the reciprocal
          condition number for the selected right invariant subspace.
          Not referenced if SENSE = 'N' or 'E'.
[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,3*N).
          Also, if SENSE = 'E' or 'V' or 'B',
          LWORK >= N+2*SDIM*(N-SDIM), where SDIM is the number of
          selected eigenvalues computed by this routine.  Note that
          N+2*SDIM*(N-SDIM) <= N+N*N/2. Note also that an error is only
          returned if LWORK < max(1,3*N), but if SENSE = 'E' or 'V' or
          'B' this may not be large enough.
          For good performance, LWORK must generally be larger.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates upper bounds on the optimal sizes of the
          arrays WORK and IWORK, returns these values as the first
          entries of the WORK and IWORK arrays, and no error messages
          related to LWORK or LIWORK are issued by XERBLA.
[out]IWORK
          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
[in]LIWORK
          LIWORK is INTEGER
          The dimension of the array IWORK.
          LIWORK >= 1; if SENSE = 'V' or 'B', LIWORK >= SDIM*(N-SDIM).
          Note that SDIM*(N-SDIM) <= N*N/4. Note also that an error is
          only returned if LIWORK < 1, but if SENSE = 'V' or 'B' this
          may not be large enough.

          If LIWORK = -1, then a workspace query is assumed; the
          routine only calculates upper bounds on the optimal sizes of
          the arrays WORK and IWORK, returns these values as the first
          entries of the WORK and IWORK arrays, and no error messages
          related to LWORK or LIWORK are issued by XERBLA.
[out]BWORK
          BWORK is LOGICAL array, dimension (N)
          Not referenced if SORT = 'N'.
[out]INFO
          INFO is INTEGER
          = 0: successful exit
          < 0: if INFO = -i, the i-th argument had an illegal value.
          > 0: if INFO = i, and i is
             <= N: the QR algorithm failed to compute all the
                   eigenvalues; elements 1:ILO-1 and i+1:N of WR and WI
                   contain those eigenvalues which have converged; if
                   JOBVS = 'V', VS contains the transformation which
                   reduces A to its partially converged Schur form.
             = N+1: the eigenvalues could not be reordered because some
                   eigenvalues were too close to separate (the problem
                   is very ill-conditioned);
             = N+2: after reordering, roundoff changed values of some
                   complex eigenvalues so that leading eigenvalues in
                   the Schur form no longer satisfy SELECT=.TRUE.  This
                   could also be caused by underflow due to scaling.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 278 of file sgeesx.f.

281*
282* -- LAPACK driver routine --
283* -- LAPACK is a software package provided by Univ. of Tennessee, --
284* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
285*
286* .. Scalar Arguments ..
287 CHARACTER JOBVS, SENSE, SORT
288 INTEGER INFO, LDA, LDVS, LIWORK, LWORK, N, SDIM
289 REAL RCONDE, RCONDV
290* ..
291* .. Array Arguments ..
292 LOGICAL BWORK( * )
293 INTEGER IWORK( * )
294 REAL A( LDA, * ), VS( LDVS, * ), WI( * ), WORK( * ),
295 $ WR( * )
296* ..
297* .. Function Arguments ..
298 LOGICAL SELECT
299 EXTERNAL SELECT
300* ..
301*
302* =====================================================================
303*
304* .. Parameters ..
305 REAL ZERO, ONE
306 parameter( zero = 0.0e0, one = 1.0e0 )
307* ..
308* .. Local Scalars ..
309 LOGICAL CURSL, LASTSL, LQUERY, LST2SL, SCALEA, WANTSB,
310 $ WANTSE, WANTSN, WANTST, WANTSV, WANTVS
311 INTEGER HSWORK, I, I1, I2, IBAL, ICOND, IERR, IEVAL,
312 $ IHI, ILO, INXT, IP, ITAU, IWRK, LWRK, LIWRK,
313 $ MAXWRK, MINWRK
314 REAL ANRM, BIGNUM, CSCALE, EPS, SMLNUM
315* ..
316* .. Local Arrays ..
317 REAL DUM( 1 )
318* ..
319* .. External Subroutines ..
320 EXTERNAL scopy, sgebak, sgebal, sgehrd, shseqr, slabad,
322* ..
323* .. External Functions ..
324 LOGICAL LSAME
325 INTEGER ILAENV
326 REAL SLAMCH, SLANGE
327 EXTERNAL lsame, ilaenv, slamch, slange
328* ..
329* .. Intrinsic Functions ..
330 INTRINSIC max, sqrt
331* ..
332* .. Executable Statements ..
333*
334* Test the input arguments
335*
336 info = 0
337 wantvs = lsame( jobvs, 'V' )
338 wantst = lsame( sort, 'S' )
339 wantsn = lsame( sense, 'N' )
340 wantse = lsame( sense, 'E' )
341 wantsv = lsame( sense, 'V' )
342 wantsb = lsame( sense, 'B' )
343 lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
344*
345 IF( ( .NOT.wantvs ) .AND. ( .NOT.lsame( jobvs, 'N' ) ) ) THEN
346 info = -1
347 ELSE IF( ( .NOT.wantst ) .AND. ( .NOT.lsame( sort, 'N' ) ) ) THEN
348 info = -2
349 ELSE IF( .NOT.( wantsn .OR. wantse .OR. wantsv .OR. wantsb ) .OR.
350 $ ( .NOT.wantst .AND. .NOT.wantsn ) ) THEN
351 info = -4
352 ELSE IF( n.LT.0 ) THEN
353 info = -5
354 ELSE IF( lda.LT.max( 1, n ) ) THEN
355 info = -7
356 ELSE IF( ldvs.LT.1 .OR. ( wantvs .AND. ldvs.LT.n ) ) THEN
357 info = -12
358 END IF
359*
360* Compute workspace
361* (Note: Comments in the code beginning "RWorkspace:" describe the
362* minimal amount of real workspace needed at that point in the
363* code, as well as the preferred amount for good performance.
364* IWorkspace refers to integer workspace.
365* NB refers to the optimal block size for the immediately
366* following subroutine, as returned by ILAENV.
367* HSWORK refers to the workspace preferred by SHSEQR, as
368* calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
369* the worst case.
370* If SENSE = 'E', 'V' or 'B', then the amount of workspace needed
371* depends on SDIM, which is computed by the routine STRSEN later
372* in the code.)
373*
374 IF( info.EQ.0 ) THEN
375 liwrk = 1
376 IF( n.EQ.0 ) THEN
377 minwrk = 1
378 lwrk = 1
379 ELSE
380 maxwrk = 2*n + n*ilaenv( 1, 'SGEHRD', ' ', n, 1, n, 0 )
381 minwrk = 3*n
382*
383 CALL shseqr( 'S', jobvs, n, 1, n, a, lda, wr, wi, vs, ldvs,
384 $ work, -1, ieval )
385 hswork = int( work( 1 ) )
386*
387 IF( .NOT.wantvs ) THEN
388 maxwrk = max( maxwrk, n + hswork )
389 ELSE
390 maxwrk = max( maxwrk, 2*n + ( n - 1 )*ilaenv( 1,
391 $ 'SORGHR', ' ', n, 1, n, -1 ) )
392 maxwrk = max( maxwrk, n + hswork )
393 END IF
394 lwrk = maxwrk
395 IF( .NOT.wantsn )
396 $ lwrk = max( lwrk, n + ( n*n )/2 )
397 IF( wantsv .OR. wantsb )
398 $ liwrk = ( n*n )/4
399 END IF
400 iwork( 1 ) = liwrk
401 work( 1 ) = lwrk
402*
403 IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
404 info = -16
405 ELSE IF( liwork.LT.1 .AND. .NOT.lquery ) THEN
406 info = -18
407 END IF
408 END IF
409*
410 IF( info.NE.0 ) THEN
411 CALL xerbla( 'SGEESX', -info )
412 RETURN
413 ELSE IF( lquery ) THEN
414 RETURN
415 END IF
416*
417* Quick return if possible
418*
419 IF( n.EQ.0 ) THEN
420 sdim = 0
421 RETURN
422 END IF
423*
424* Get machine constants
425*
426 eps = slamch( 'P' )
427 smlnum = slamch( 'S' )
428 bignum = one / smlnum
429 CALL slabad( smlnum, bignum )
430 smlnum = sqrt( smlnum ) / eps
431 bignum = one / smlnum
432*
433* Scale A if max element outside range [SMLNUM,BIGNUM]
434*
435 anrm = slange( 'M', n, n, a, lda, dum )
436 scalea = .false.
437 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
438 scalea = .true.
439 cscale = smlnum
440 ELSE IF( anrm.GT.bignum ) THEN
441 scalea = .true.
442 cscale = bignum
443 END IF
444 IF( scalea )
445 $ CALL slascl( 'G', 0, 0, anrm, cscale, n, n, a, lda, ierr )
446*
447* Permute the matrix to make it more nearly triangular
448* (RWorkspace: need N)
449*
450 ibal = 1
451 CALL sgebal( 'P', n, a, lda, ilo, ihi, work( ibal ), ierr )
452*
453* Reduce to upper Hessenberg form
454* (RWorkspace: need 3*N, prefer 2*N+N*NB)
455*
456 itau = n + ibal
457 iwrk = n + itau
458 CALL sgehrd( n, ilo, ihi, a, lda, work( itau ), work( iwrk ),
459 $ lwork-iwrk+1, ierr )
460*
461 IF( wantvs ) THEN
462*
463* Copy Householder vectors to VS
464*
465 CALL slacpy( 'L', n, n, a, lda, vs, ldvs )
466*
467* Generate orthogonal matrix in VS
468* (RWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
469*
470 CALL sorghr( n, ilo, ihi, vs, ldvs, work( itau ), work( iwrk ),
471 $ lwork-iwrk+1, ierr )
472 END IF
473*
474 sdim = 0
475*
476* Perform QR iteration, accumulating Schur vectors in VS if desired
477* (RWorkspace: need N+1, prefer N+HSWORK (see comments) )
478*
479 iwrk = itau
480 CALL shseqr( 'S', jobvs, n, ilo, ihi, a, lda, wr, wi, vs, ldvs,
481 $ work( iwrk ), lwork-iwrk+1, ieval )
482 IF( ieval.GT.0 )
483 $ info = ieval
484*
485* Sort eigenvalues if desired
486*
487 IF( wantst .AND. info.EQ.0 ) THEN
488 IF( scalea ) THEN
489 CALL slascl( 'G', 0, 0, cscale, anrm, n, 1, wr, n, ierr )
490 CALL slascl( 'G', 0, 0, cscale, anrm, n, 1, wi, n, ierr )
491 END IF
492 DO 10 i = 1, n
493 bwork( i ) = SELECT( wr( i ), wi( i ) )
494 10 CONTINUE
495*
496* Reorder eigenvalues, transform Schur vectors, and compute
497* reciprocal condition numbers
498* (RWorkspace: if SENSE is not 'N', need N+2*SDIM*(N-SDIM)
499* otherwise, need N )
500* (IWorkspace: if SENSE is 'V' or 'B', need SDIM*(N-SDIM)
501* otherwise, need 0 )
502*
503 CALL strsen( sense, jobvs, bwork, n, a, lda, vs, ldvs, wr, wi,
504 $ sdim, rconde, rcondv, work( iwrk ), lwork-iwrk+1,
505 $ iwork, liwork, icond )
506 IF( .NOT.wantsn )
507 $ maxwrk = max( maxwrk, n+2*sdim*( n-sdim ) )
508 IF( icond.EQ.-15 ) THEN
509*
510* Not enough real workspace
511*
512 info = -16
513 ELSE IF( icond.EQ.-17 ) THEN
514*
515* Not enough integer workspace
516*
517 info = -18
518 ELSE IF( icond.GT.0 ) THEN
519*
520* STRSEN failed to reorder or to restore standard Schur form
521*
522 info = icond + n
523 END IF
524 END IF
525*
526 IF( wantvs ) THEN
527*
528* Undo balancing
529* (RWorkspace: need N)
530*
531 CALL sgebak( 'P', 'R', n, ilo, ihi, work( ibal ), n, vs, ldvs,
532 $ ierr )
533 END IF
534*
535 IF( scalea ) THEN
536*
537* Undo scaling for the Schur form of A
538*
539 CALL slascl( 'H', 0, 0, cscale, anrm, n, n, a, lda, ierr )
540 CALL scopy( n, a, lda+1, wr, 1 )
541 IF( ( wantsv .OR. wantsb ) .AND. info.EQ.0 ) THEN
542 dum( 1 ) = rcondv
543 CALL slascl( 'G', 0, 0, cscale, anrm, 1, 1, dum, 1, ierr )
544 rcondv = dum( 1 )
545 END IF
546 IF( cscale.EQ.smlnum ) THEN
547*
548* If scaling back towards underflow, adjust WI if an
549* offdiagonal element of a 2-by-2 block in the Schur form
550* underflows.
551*
552 IF( ieval.GT.0 ) THEN
553 i1 = ieval + 1
554 i2 = ihi - 1
555 CALL slascl( 'G', 0, 0, cscale, anrm, ilo-1, 1, wi, n,
556 $ ierr )
557 ELSE IF( wantst ) THEN
558 i1 = 1
559 i2 = n - 1
560 ELSE
561 i1 = ilo
562 i2 = ihi - 1
563 END IF
564 inxt = i1 - 1
565 DO 20 i = i1, i2
566 IF( i.LT.inxt )
567 $ GO TO 20
568 IF( wi( i ).EQ.zero ) THEN
569 inxt = i + 1
570 ELSE
571 IF( a( i+1, i ).EQ.zero ) THEN
572 wi( i ) = zero
573 wi( i+1 ) = zero
574 ELSE IF( a( i+1, i ).NE.zero .AND. a( i, i+1 ).EQ.
575 $ zero ) THEN
576 wi( i ) = zero
577 wi( i+1 ) = zero
578 IF( i.GT.1 )
579 $ CALL sswap( i-1, a( 1, i ), 1, a( 1, i+1 ), 1 )
580 IF( n.GT.i+1 )
581 $ CALL sswap( n-i-1, a( i, i+2 ), lda,
582 $ a( i+1, i+2 ), lda )
583 IF( wantvs ) THEN
584 CALL sswap( n, vs( 1, i ), 1, vs( 1, i+1 ), 1 )
585 END IF
586 a( i, i+1 ) = a( i+1, i )
587 a( i+1, i ) = zero
588 END IF
589 inxt = i + 2
590 END IF
591 20 CONTINUE
592 END IF
593 CALL slascl( 'G', 0, 0, cscale, anrm, n-ieval, 1,
594 $ wi( ieval+1 ), max( n-ieval, 1 ), ierr )
595 END IF
596*
597 IF( wantst .AND. info.EQ.0 ) THEN
598*
599* Check if reordering successful
600*
601 lastsl = .true.
602 lst2sl = .true.
603 sdim = 0
604 ip = 0
605 DO 30 i = 1, n
606 cursl = SELECT( wr( i ), wi( i ) )
607 IF( wi( i ).EQ.zero ) THEN
608 IF( cursl )
609 $ sdim = sdim + 1
610 ip = 0
611 IF( cursl .AND. .NOT.lastsl )
612 $ info = n + 2
613 ELSE
614 IF( ip.EQ.1 ) THEN
615*
616* Last eigenvalue of conjugate pair
617*
618 cursl = cursl .OR. lastsl
619 lastsl = cursl
620 IF( cursl )
621 $ sdim = sdim + 2
622 ip = -1
623 IF( cursl .AND. .NOT.lst2sl )
624 $ info = n + 2
625 ELSE
626*
627* First eigenvalue of conjugate pair
628*
629 ip = 1
630 END IF
631 END IF
632 lst2sl = lastsl
633 lastsl = cursl
634 30 CONTINUE
635 END IF
636*
637 work( 1 ) = maxwrk
638 IF( wantsv .OR. wantsb ) THEN
639 iwork( 1 ) = sdim*(n-sdim)
640 ELSE
641 iwork( 1 ) = 1
642 END IF
643*
644 RETURN
645*
646* End of SGEESX
647*
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 slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:103
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: ilaenv.f:162
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
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 sgebal(JOB, N, A, LDA, ILO, IHI, SCALE, INFO)
SGEBAL
Definition: sgebal.f:160
subroutine sgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
SGEHRD
Definition: sgehrd.f:167
subroutine sgebak(JOB, SIDE, N, ILO, IHI, SCALE, M, V, LDV, INFO)
SGEBAK
Definition: sgebak.f:130
subroutine sorghr(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
SORGHR
Definition: sorghr.f:126
subroutine strsen(JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, WR, WI, M, S, SEP, WORK, LWORK, IWORK, LIWORK, INFO)
STRSEN
Definition: strsen.f:314
subroutine shseqr(JOB, COMPZ, N, ILO, IHI, H, LDH, WR, WI, Z, LDZ, WORK, LWORK, INFO)
SHSEQR
Definition: shseqr.f:316
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:82
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82
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: