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

◆ dgeesx()

subroutine dgeesx ( character  jobvs,
character  sort,
external  select,
character  sense,
integer  n,
double precision, dimension( lda, * )  a,
integer  lda,
integer  sdim,
double precision, dimension( * )  wr,
double precision, dimension( * )  wi,
double precision, dimension( ldvs, * )  vs,
integer  ldvs,
double precision  rconde,
double precision  rcondv,
double precision, dimension( * )  work,
integer  lwork,
integer, dimension( * )  iwork,
integer  liwork,
logical, dimension( * )  bwork,
integer  info 
)

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

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

Purpose:
 DGEESX 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (N)
[out]WI
          WI is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION
          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 DOUBLE PRECISION
          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 DOUBLE PRECISION 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 dgeesx.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 DOUBLE PRECISION RCONDE, RCONDV
290* ..
291* .. Array Arguments ..
292 LOGICAL BWORK( * )
293 INTEGER IWORK( * )
294 DOUBLE PRECISION 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 DOUBLE PRECISION ZERO, ONE
306 parameter( zero = 0.0d0, one = 1.0d0 )
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, LIWRK, LWRK,
313 $ MAXWRK, MINWRK
314 DOUBLE PRECISION ANRM, BIGNUM, CSCALE, EPS, SMLNUM
315* ..
316* .. Local Arrays ..
317 DOUBLE PRECISION DUM( 1 )
318* ..
319* .. External Subroutines ..
320 EXTERNAL dcopy, dgebak, dgebal, dgehrd, dhseqr, dlacpy,
322* ..
323* .. External Functions ..
324 LOGICAL LSAME
325 INTEGER ILAENV
326 DOUBLE PRECISION DLAMCH, DLANGE
327 EXTERNAL lsame, ilaenv, dlamch, dlange
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 DHSEQR, 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 DTRSEN 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, 'DGEHRD', ' ', n, 1, n, 0 )
381 minwrk = 3*n
382*
383 CALL dhseqr( '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 $ 'DORGHR', ' ', 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( 'DGEESX', -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 = dlamch( 'P' )
427 smlnum = dlamch( 'S' )
428 bignum = one / smlnum
429 smlnum = sqrt( smlnum ) / eps
430 bignum = one / smlnum
431*
432* Scale A if max element outside range [SMLNUM,BIGNUM]
433*
434 anrm = dlange( 'M', n, n, a, lda, dum )
435 scalea = .false.
436 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
437 scalea = .true.
438 cscale = smlnum
439 ELSE IF( anrm.GT.bignum ) THEN
440 scalea = .true.
441 cscale = bignum
442 END IF
443 IF( scalea )
444 $ CALL dlascl( 'G', 0, 0, anrm, cscale, n, n, a, lda, ierr )
445*
446* Permute the matrix to make it more nearly triangular
447* (RWorkspace: need N)
448*
449 ibal = 1
450 CALL dgebal( 'P', n, a, lda, ilo, ihi, work( ibal ), ierr )
451*
452* Reduce to upper Hessenberg form
453* (RWorkspace: need 3*N, prefer 2*N+N*NB)
454*
455 itau = n + ibal
456 iwrk = n + itau
457 CALL dgehrd( n, ilo, ihi, a, lda, work( itau ), work( iwrk ),
458 $ lwork-iwrk+1, ierr )
459*
460 IF( wantvs ) THEN
461*
462* Copy Householder vectors to VS
463*
464 CALL dlacpy( 'L', n, n, a, lda, vs, ldvs )
465*
466* Generate orthogonal matrix in VS
467* (RWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
468*
469 CALL dorghr( n, ilo, ihi, vs, ldvs, work( itau ), work( iwrk ),
470 $ lwork-iwrk+1, ierr )
471 END IF
472*
473 sdim = 0
474*
475* Perform QR iteration, accumulating Schur vectors in VS if desired
476* (RWorkspace: need N+1, prefer N+HSWORK (see comments) )
477*
478 iwrk = itau
479 CALL dhseqr( 'S', jobvs, n, ilo, ihi, a, lda, wr, wi, vs, ldvs,
480 $ work( iwrk ), lwork-iwrk+1, ieval )
481 IF( ieval.GT.0 )
482 $ info = ieval
483*
484* Sort eigenvalues if desired
485*
486 IF( wantst .AND. info.EQ.0 ) THEN
487 IF( scalea ) THEN
488 CALL dlascl( 'G', 0, 0, cscale, anrm, n, 1, wr, n, ierr )
489 CALL dlascl( 'G', 0, 0, cscale, anrm, n, 1, wi, n, ierr )
490 END IF
491 DO 10 i = 1, n
492 bwork( i ) = SELECT( wr( i ), wi( i ) )
493 10 CONTINUE
494*
495* Reorder eigenvalues, transform Schur vectors, and compute
496* reciprocal condition numbers
497* (RWorkspace: if SENSE is not 'N', need N+2*SDIM*(N-SDIM)
498* otherwise, need N )
499* (IWorkspace: if SENSE is 'V' or 'B', need SDIM*(N-SDIM)
500* otherwise, need 0 )
501*
502 CALL dtrsen( sense, jobvs, bwork, n, a, lda, vs, ldvs, wr, wi,
503 $ sdim, rconde, rcondv, work( iwrk ), lwork-iwrk+1,
504 $ iwork, liwork, icond )
505 IF( .NOT.wantsn )
506 $ maxwrk = max( maxwrk, n+2*sdim*( n-sdim ) )
507 IF( icond.EQ.-15 ) THEN
508*
509* Not enough real workspace
510*
511 info = -16
512 ELSE IF( icond.EQ.-17 ) THEN
513*
514* Not enough integer workspace
515*
516 info = -18
517 ELSE IF( icond.GT.0 ) THEN
518*
519* DTRSEN failed to reorder or to restore standard Schur form
520*
521 info = icond + n
522 END IF
523 END IF
524*
525 IF( wantvs ) THEN
526*
527* Undo balancing
528* (RWorkspace: need N)
529*
530 CALL dgebak( 'P', 'R', n, ilo, ihi, work( ibal ), n, vs, ldvs,
531 $ ierr )
532 END IF
533*
534 IF( scalea ) THEN
535*
536* Undo scaling for the Schur form of A
537*
538 CALL dlascl( 'H', 0, 0, cscale, anrm, n, n, a, lda, ierr )
539 CALL dcopy( n, a, lda+1, wr, 1 )
540 IF( ( wantsv .OR. wantsb ) .AND. info.EQ.0 ) THEN
541 dum( 1 ) = rcondv
542 CALL dlascl( 'G', 0, 0, cscale, anrm, 1, 1, dum, 1, ierr )
543 rcondv = dum( 1 )
544 END IF
545 IF( cscale.EQ.smlnum ) THEN
546*
547* If scaling back towards underflow, adjust WI if an
548* offdiagonal element of a 2-by-2 block in the Schur form
549* underflows.
550*
551 IF( ieval.GT.0 ) THEN
552 i1 = ieval + 1
553 i2 = ihi - 1
554 CALL dlascl( 'G', 0, 0, cscale, anrm, ilo-1, 1, wi, n,
555 $ ierr )
556 ELSE IF( wantst ) THEN
557 i1 = 1
558 i2 = n - 1
559 ELSE
560 i1 = ilo
561 i2 = ihi - 1
562 END IF
563 inxt = i1 - 1
564 DO 20 i = i1, i2
565 IF( i.LT.inxt )
566 $ GO TO 20
567 IF( wi( i ).EQ.zero ) THEN
568 inxt = i + 1
569 ELSE
570 IF( a( i+1, i ).EQ.zero ) THEN
571 wi( i ) = zero
572 wi( i+1 ) = zero
573 ELSE IF( a( i+1, i ).NE.zero .AND. a( i, i+1 ).EQ.
574 $ zero ) THEN
575 wi( i ) = zero
576 wi( i+1 ) = zero
577 IF( i.GT.1 )
578 $ CALL dswap( i-1, a( 1, i ), 1, a( 1, i+1 ), 1 )
579 IF( n.GT.i+1 )
580 $ CALL dswap( n-i-1, a( i, i+2 ), lda,
581 $ a( i+1, i+2 ), lda )
582 IF( wantvs ) THEN
583 CALL dswap( n, vs( 1, i ), 1, vs( 1, i+1 ), 1 )
584 END IF
585 a( i, i+1 ) = a( i+1, i )
586 a( i+1, i ) = zero
587 END IF
588 inxt = i + 2
589 END IF
590 20 CONTINUE
591 END IF
592 CALL dlascl( 'G', 0, 0, cscale, anrm, n-ieval, 1,
593 $ wi( ieval+1 ), max( n-ieval, 1 ), ierr )
594 END IF
595*
596 IF( wantst .AND. info.EQ.0 ) THEN
597*
598* Check if reordering successful
599*
600 lastsl = .true.
601 lst2sl = .true.
602 sdim = 0
603 ip = 0
604 DO 30 i = 1, n
605 cursl = SELECT( wr( i ), wi( i ) )
606 IF( wi( i ).EQ.zero ) THEN
607 IF( cursl )
608 $ sdim = sdim + 1
609 ip = 0
610 IF( cursl .AND. .NOT.lastsl )
611 $ info = n + 2
612 ELSE
613 IF( ip.EQ.1 ) THEN
614*
615* Last eigenvalue of conjugate pair
616*
617 cursl = cursl .OR. lastsl
618 lastsl = cursl
619 IF( cursl )
620 $ sdim = sdim + 2
621 ip = -1
622 IF( cursl .AND. .NOT.lst2sl )
623 $ info = n + 2
624 ELSE
625*
626* First eigenvalue of conjugate pair
627*
628 ip = 1
629 END IF
630 END IF
631 lst2sl = lastsl
632 lastsl = cursl
633 30 CONTINUE
634 END IF
635*
636 work( 1 ) = maxwrk
637 IF( wantsv .OR. wantsb ) THEN
638 iwork( 1 ) = max( 1, sdim*( n-sdim ) )
639 ELSE
640 iwork( 1 ) = 1
641 END IF
642*
643 RETURN
644*
645* End of DGEESX
646*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
subroutine dgebak(job, side, n, ilo, ihi, scale, m, v, ldv, info)
DGEBAK
Definition dgebak.f:130
subroutine dgebal(job, n, a, lda, ilo, ihi, scale, info)
DGEBAL
Definition dgebal.f:163
subroutine dgehrd(n, ilo, ihi, a, lda, tau, work, lwork, info)
DGEHRD
Definition dgehrd.f:167
subroutine dhseqr(job, compz, n, ilo, ihi, h, ldh, wr, wi, z, ldz, work, lwork, info)
DHSEQR
Definition dhseqr.f:316
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
Definition dlacpy.f:103
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
double precision function dlange(norm, m, n, a, lda, work)
DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition dlange.f:114
subroutine dlascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition dlascl.f:143
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
subroutine dswap(n, dx, incx, dy, incy)
DSWAP
Definition dswap.f:82
subroutine dtrsen(job, compq, select, n, t, ldt, q, ldq, wr, wi, m, s, sep, work, lwork, iwork, liwork, info)
DTRSEN
Definition dtrsen.f:313
subroutine dorghr(n, ilo, ihi, a, lda, tau, work, lwork, info)
DORGHR
Definition dorghr.f:126
Here is the call graph for this function:
Here is the caller graph for this function: