LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ sgees()

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

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

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

Purpose:
 SGEES 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.
 The leading columns of Z then form an orthonormal basis for the
 invariant subspace corresponding to the selected eigenvalues.

 A 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 complex
          eigenvalues are selected.
          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 is set to N+2 (see INFO below).
[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 has been 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 will 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; if
          JOBVS = 'V', LDVS >= N.
[out]WORK
          WORK is REAL array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) contains the optimal LWORK.
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK.  LWORK >= max(1,3*N).
          For good performance, LWORK must 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]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 matrix 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 214 of file sgees.f.

216 *
217 * -- LAPACK driver routine --
218 * -- LAPACK is a software package provided by Univ. of Tennessee, --
219 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
220 *
221 * .. Scalar Arguments ..
222  CHARACTER JOBVS, SORT
223  INTEGER INFO, LDA, LDVS, LWORK, N, SDIM
224 * ..
225 * .. Array Arguments ..
226  LOGICAL BWORK( * )
227  REAL A( LDA, * ), VS( LDVS, * ), WI( * ), WORK( * ),
228  $ WR( * )
229 * ..
230 * .. Function Arguments ..
231  LOGICAL SELECT
232  EXTERNAL SELECT
233 * ..
234 *
235 * =====================================================================
236 *
237 * .. Parameters ..
238  REAL ZERO, ONE
239  parameter( zero = 0.0e0, one = 1.0e0 )
240 * ..
241 * .. Local Scalars ..
242  LOGICAL CURSL, LASTSL, LQUERY, LST2SL, SCALEA, WANTST,
243  $ WANTVS
244  INTEGER HSWORK, I, I1, I2, IBAL, ICOND, IERR, IEVAL,
245  $ IHI, ILO, INXT, IP, ITAU, IWRK, MAXWRK, MINWRK
246  REAL ANRM, BIGNUM, CSCALE, EPS, S, SEP, SMLNUM
247 * ..
248 * .. Local Arrays ..
249  INTEGER IDUM( 1 )
250  REAL DUM( 1 )
251 * ..
252 * .. External Subroutines ..
253  EXTERNAL scopy, sgebak, sgebal, sgehrd, shseqr, slabad,
255 * ..
256 * .. External Functions ..
257  LOGICAL LSAME
258  INTEGER ILAENV
259  REAL SLAMCH, SLANGE
260  EXTERNAL lsame, ilaenv, slamch, slange
261 * ..
262 * .. Intrinsic Functions ..
263  INTRINSIC max, sqrt
264 * ..
265 * .. Executable Statements ..
266 *
267 * Test the input arguments
268 *
269  info = 0
270  lquery = ( lwork.EQ.-1 )
271  wantvs = lsame( jobvs, 'V' )
272  wantst = lsame( sort, 'S' )
273  IF( ( .NOT.wantvs ) .AND. ( .NOT.lsame( jobvs, 'N' ) ) ) THEN
274  info = -1
275  ELSE IF( ( .NOT.wantst ) .AND. ( .NOT.lsame( sort, 'N' ) ) ) THEN
276  info = -2
277  ELSE IF( n.LT.0 ) THEN
278  info = -4
279  ELSE IF( lda.LT.max( 1, n ) ) THEN
280  info = -6
281  ELSE IF( ldvs.LT.1 .OR. ( wantvs .AND. ldvs.LT.n ) ) THEN
282  info = -11
283  END IF
284 *
285 * Compute workspace
286 * (Note: Comments in the code beginning "Workspace:" describe the
287 * minimal amount of workspace needed at that point in the code,
288 * as well as the preferred amount for good performance.
289 * NB refers to the optimal block size for the immediately
290 * following subroutine, as returned by ILAENV.
291 * HSWORK refers to the workspace preferred by SHSEQR, as
292 * calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
293 * the worst case.)
294 *
295  IF( info.EQ.0 ) THEN
296  IF( n.EQ.0 ) THEN
297  minwrk = 1
298  maxwrk = 1
299  ELSE
300  maxwrk = 2*n + n*ilaenv( 1, 'SGEHRD', ' ', n, 1, n, 0 )
301  minwrk = 3*n
302 *
303  CALL shseqr( 'S', jobvs, n, 1, n, a, lda, wr, wi, vs, ldvs,
304  $ work, -1, ieval )
305  hswork = work( 1 )
306 *
307  IF( .NOT.wantvs ) THEN
308  maxwrk = max( maxwrk, n + hswork )
309  ELSE
310  maxwrk = max( maxwrk, 2*n + ( n - 1 )*ilaenv( 1,
311  $ 'SORGHR', ' ', n, 1, n, -1 ) )
312  maxwrk = max( maxwrk, n + hswork )
313  END IF
314  END IF
315  work( 1 ) = maxwrk
316 *
317  IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
318  info = -13
319  END IF
320  END IF
321 *
322  IF( info.NE.0 ) THEN
323  CALL xerbla( 'SGEES ', -info )
324  RETURN
325  ELSE IF( lquery ) THEN
326  RETURN
327  END IF
328 *
329 * Quick return if possible
330 *
331  IF( n.EQ.0 ) THEN
332  sdim = 0
333  RETURN
334  END IF
335 *
336 * Get machine constants
337 *
338  eps = slamch( 'P' )
339  smlnum = slamch( 'S' )
340  bignum = one / smlnum
341  CALL slabad( smlnum, bignum )
342  smlnum = sqrt( smlnum ) / eps
343  bignum = one / smlnum
344 *
345 * Scale A if max element outside range [SMLNUM,BIGNUM]
346 *
347  anrm = slange( 'M', n, n, a, lda, dum )
348  scalea = .false.
349  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
350  scalea = .true.
351  cscale = smlnum
352  ELSE IF( anrm.GT.bignum ) THEN
353  scalea = .true.
354  cscale = bignum
355  END IF
356  IF( scalea )
357  $ CALL slascl( 'G', 0, 0, anrm, cscale, n, n, a, lda, ierr )
358 *
359 * Permute the matrix to make it more nearly triangular
360 * (Workspace: need N)
361 *
362  ibal = 1
363  CALL sgebal( 'P', n, a, lda, ilo, ihi, work( ibal ), ierr )
364 *
365 * Reduce to upper Hessenberg form
366 * (Workspace: need 3*N, prefer 2*N+N*NB)
367 *
368  itau = n + ibal
369  iwrk = n + itau
370  CALL sgehrd( n, ilo, ihi, a, lda, work( itau ), work( iwrk ),
371  $ lwork-iwrk+1, ierr )
372 *
373  IF( wantvs ) THEN
374 *
375 * Copy Householder vectors to VS
376 *
377  CALL slacpy( 'L', n, n, a, lda, vs, ldvs )
378 *
379 * Generate orthogonal matrix in VS
380 * (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
381 *
382  CALL sorghr( n, ilo, ihi, vs, ldvs, work( itau ), work( iwrk ),
383  $ lwork-iwrk+1, ierr )
384  END IF
385 *
386  sdim = 0
387 *
388 * Perform QR iteration, accumulating Schur vectors in VS if desired
389 * (Workspace: need N+1, prefer N+HSWORK (see comments) )
390 *
391  iwrk = itau
392  CALL shseqr( 'S', jobvs, n, ilo, ihi, a, lda, wr, wi, vs, ldvs,
393  $ work( iwrk ), lwork-iwrk+1, ieval )
394  IF( ieval.GT.0 )
395  $ info = ieval
396 *
397 * Sort eigenvalues if desired
398 *
399  IF( wantst .AND. info.EQ.0 ) THEN
400  IF( scalea ) THEN
401  CALL slascl( 'G', 0, 0, cscale, anrm, n, 1, wr, n, ierr )
402  CALL slascl( 'G', 0, 0, cscale, anrm, n, 1, wi, n, ierr )
403  END IF
404  DO 10 i = 1, n
405  bwork( i ) = SELECT( wr( i ), wi( i ) )
406  10 CONTINUE
407 *
408 * Reorder eigenvalues and transform Schur vectors
409 * (Workspace: none needed)
410 *
411  CALL strsen( 'N', jobvs, bwork, n, a, lda, vs, ldvs, wr, wi,
412  $ sdim, s, sep, work( iwrk ), lwork-iwrk+1, idum, 1,
413  $ icond )
414  IF( icond.GT.0 )
415  $ info = n + icond
416  END IF
417 *
418  IF( wantvs ) THEN
419 *
420 * Undo balancing
421 * (Workspace: need N)
422 *
423  CALL sgebak( 'P', 'R', n, ilo, ihi, work( ibal ), n, vs, ldvs,
424  $ ierr )
425  END IF
426 *
427  IF( scalea ) THEN
428 *
429 * Undo scaling for the Schur form of A
430 *
431  CALL slascl( 'H', 0, 0, cscale, anrm, n, n, a, lda, ierr )
432  CALL scopy( n, a, lda+1, wr, 1 )
433  IF( cscale.EQ.smlnum ) THEN
434 *
435 * If scaling back towards underflow, adjust WI if an
436 * offdiagonal element of a 2-by-2 block in the Schur form
437 * underflows.
438 *
439  IF( ieval.GT.0 ) THEN
440  i1 = ieval + 1
441  i2 = ihi - 1
442  CALL slascl( 'G', 0, 0, cscale, anrm, ilo-1, 1, wi,
443  $ max( ilo-1, 1 ), ierr )
444  ELSE IF( wantst ) THEN
445  i1 = 1
446  i2 = n - 1
447  ELSE
448  i1 = ilo
449  i2 = ihi - 1
450  END IF
451  inxt = i1 - 1
452  DO 20 i = i1, i2
453  IF( i.LT.inxt )
454  $ GO TO 20
455  IF( wi( i ).EQ.zero ) THEN
456  inxt = i + 1
457  ELSE
458  IF( a( i+1, i ).EQ.zero ) THEN
459  wi( i ) = zero
460  wi( i+1 ) = zero
461  ELSE IF( a( i+1, i ).NE.zero .AND. a( i, i+1 ).EQ.
462  $ zero ) THEN
463  wi( i ) = zero
464  wi( i+1 ) = zero
465  IF( i.GT.1 )
466  $ CALL sswap( i-1, a( 1, i ), 1, a( 1, i+1 ), 1 )
467  IF( n.GT.i+1 )
468  $ CALL sswap( n-i-1, a( i, i+2 ), lda,
469  $ a( i+1, i+2 ), lda )
470  IF( wantvs ) THEN
471  CALL sswap( n, vs( 1, i ), 1, vs( 1, i+1 ), 1 )
472  END IF
473  a( i, i+1 ) = a( i+1, i )
474  a( i+1, i ) = zero
475  END IF
476  inxt = i + 2
477  END IF
478  20 CONTINUE
479  END IF
480 *
481 * Undo scaling for the imaginary part of the eigenvalues
482 *
483  CALL slascl( 'G', 0, 0, cscale, anrm, n-ieval, 1,
484  $ wi( ieval+1 ), max( n-ieval, 1 ), ierr )
485  END IF
486 *
487  IF( wantst .AND. info.EQ.0 ) THEN
488 *
489 * Check if reordering successful
490 *
491  lastsl = .true.
492  lst2sl = .true.
493  sdim = 0
494  ip = 0
495  DO 30 i = 1, n
496  cursl = SELECT( wr( i ), wi( i ) )
497  IF( wi( i ).EQ.zero ) THEN
498  IF( cursl )
499  $ sdim = sdim + 1
500  ip = 0
501  IF( cursl .AND. .NOT.lastsl )
502  $ info = n + 2
503  ELSE
504  IF( ip.EQ.1 ) THEN
505 *
506 * Last eigenvalue of conjugate pair
507 *
508  cursl = cursl .OR. lastsl
509  lastsl = cursl
510  IF( cursl )
511  $ sdim = sdim + 2
512  ip = -1
513  IF( cursl .AND. .NOT.lst2sl )
514  $ info = n + 2
515  ELSE
516 *
517 * First eigenvalue of conjugate pair
518 *
519  ip = 1
520  END IF
521  END IF
522  lst2sl = lastsl
523  lastsl = cursl
524  30 CONTINUE
525  END IF
526 *
527  work( 1 ) = maxwrk
528  RETURN
529 *
530 * End of SGEES
531 *
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: