LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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.
Date
June 2016

Definition at line 283 of file sgeesx.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: