LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ zgeesx()

subroutine zgeesx ( character  JOBVS,
character  SORT,
external  SELECT,
character  SENSE,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
integer  SDIM,
complex*16, dimension( * )  W,
complex*16, dimension( ldvs, * )  VS,
integer  LDVS,
double precision  RCONDE,
double precision  RCONDV,
complex*16, dimension( * )  WORK,
integer  LWORK,
double precision, dimension( * )  RWORK,
logical, dimension( * )  BWORK,
integer  INFO 
)

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

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

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

 Optionally, it also orders the eigenvalues on the diagonal of the
 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 complex matrix is in Schur form if it is upper triangular.
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 one COMPLEX*16 argument
          SELECT must be declared EXTERNAL in the calling subroutine.
          If SORT = 'S', SELECT is used to select eigenvalues to order
          to the top left of the Schur form.
          If SORT = 'N', SELECT is not referenced.
          An eigenvalue W(j) is selected if SELECT(W(j)) is true.
[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 COMPLEX*16 array, dimension (LDA, N)
          On entry, the N-by-N matrix A.
          On exit, A is overwritten by its 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 for which
                         SELECT is true.
[out]W
          W is COMPLEX*16 array, dimension (N)
          W contains the computed eigenvalues, in the same order
          that they appear on the diagonal of the output Schur form T.
[out]VS
          VS is COMPLEX*16 array, dimension (LDVS,N)
          If JOBVS = 'V', VS contains the unitary 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 COMPLEX*16 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,2*N).
          Also, if SENSE = 'E' or 'V' or 'B', LWORK >= 2*SDIM*(N-SDIM),
          where SDIM is the number of selected eigenvalues computed by
          this routine.  Note that 2*SDIM*(N-SDIM) <= N*N/2. Note also
          that an error is only returned if LWORK < max(1,2*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 bound on the optimal size of the
          array WORK, returns this value as the first entry of the WORK
          array, and no error message related to LWORK is issued by
          XERBLA.
[out]RWORK
          RWORK is DOUBLE PRECISION array, dimension (N)
[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 W
                   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 236 of file zgeesx.f.

239 *
240 * -- LAPACK driver routine --
241 * -- LAPACK is a software package provided by Univ. of Tennessee, --
242 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
243 *
244 * .. Scalar Arguments ..
245  CHARACTER JOBVS, SENSE, SORT
246  INTEGER INFO, LDA, LDVS, LWORK, N, SDIM
247  DOUBLE PRECISION RCONDE, RCONDV
248 * ..
249 * .. Array Arguments ..
250  LOGICAL BWORK( * )
251  DOUBLE PRECISION RWORK( * )
252  COMPLEX*16 A( LDA, * ), VS( LDVS, * ), W( * ), WORK( * )
253 * ..
254 * .. Function Arguments ..
255  LOGICAL SELECT
256  EXTERNAL SELECT
257 * ..
258 *
259 * =====================================================================
260 *
261 * .. Parameters ..
262  DOUBLE PRECISION ZERO, ONE
263  parameter( zero = 0.0d0, one = 1.0d0 )
264 * ..
265 * .. Local Scalars ..
266  LOGICAL LQUERY, SCALEA, WANTSB, WANTSE, WANTSN, WANTST,
267  $ WANTSV, WANTVS
268  INTEGER HSWORK, I, IBAL, ICOND, IERR, IEVAL, IHI, ILO,
269  $ ITAU, IWRK, LWRK, MAXWRK, MINWRK
270  DOUBLE PRECISION ANRM, BIGNUM, CSCALE, EPS, SMLNUM
271 * ..
272 * .. Local Arrays ..
273  DOUBLE PRECISION DUM( 1 )
274 * ..
275 * .. External Subroutines ..
276  EXTERNAL dlabad, dlascl, xerbla, zcopy, zgebak, zgebal,
278 * ..
279 * .. External Functions ..
280  LOGICAL LSAME
281  INTEGER ILAENV
282  DOUBLE PRECISION DLAMCH, ZLANGE
283  EXTERNAL lsame, ilaenv, dlamch, zlange
284 * ..
285 * .. Intrinsic Functions ..
286  INTRINSIC max, sqrt
287 * ..
288 * .. Executable Statements ..
289 *
290 * Test the input arguments
291 *
292  info = 0
293  wantvs = lsame( jobvs, 'V' )
294  wantst = lsame( sort, 'S' )
295  wantsn = lsame( sense, 'N' )
296  wantse = lsame( sense, 'E' )
297  wantsv = lsame( sense, 'V' )
298  wantsb = lsame( sense, 'B' )
299  lquery = ( lwork.EQ.-1 )
300 *
301  IF( ( .NOT.wantvs ) .AND. ( .NOT.lsame( jobvs, 'N' ) ) ) THEN
302  info = -1
303  ELSE IF( ( .NOT.wantst ) .AND. ( .NOT.lsame( sort, 'N' ) ) ) THEN
304  info = -2
305  ELSE IF( .NOT.( wantsn .OR. wantse .OR. wantsv .OR. wantsb ) .OR.
306  $ ( .NOT.wantst .AND. .NOT.wantsn ) ) THEN
307  info = -4
308  ELSE IF( n.LT.0 ) THEN
309  info = -5
310  ELSE IF( lda.LT.max( 1, n ) ) THEN
311  info = -7
312  ELSE IF( ldvs.LT.1 .OR. ( wantvs .AND. ldvs.LT.n ) ) THEN
313  info = -11
314  END IF
315 *
316 * Compute workspace
317 * (Note: Comments in the code beginning "Workspace:" describe the
318 * minimal amount of real workspace needed at that point in the
319 * code, as well as the preferred amount for good performance.
320 * CWorkspace refers to complex workspace, and RWorkspace to real
321 * workspace. NB refers to the optimal block size for the
322 * immediately following subroutine, as returned by ILAENV.
323 * HSWORK refers to the workspace preferred by ZHSEQR, as
324 * calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
325 * the worst case.
326 * If SENSE = 'E', 'V' or 'B', then the amount of workspace needed
327 * depends on SDIM, which is computed by the routine ZTRSEN later
328 * in the code.)
329 *
330  IF( info.EQ.0 ) THEN
331  IF( n.EQ.0 ) THEN
332  minwrk = 1
333  lwrk = 1
334  ELSE
335  maxwrk = n + n*ilaenv( 1, 'ZGEHRD', ' ', n, 1, n, 0 )
336  minwrk = 2*n
337 *
338  CALL zhseqr( 'S', jobvs, n, 1, n, a, lda, w, vs, ldvs,
339  $ work, -1, ieval )
340  hswork = dble( work( 1 ) )
341 *
342  IF( .NOT.wantvs ) THEN
343  maxwrk = max( maxwrk, hswork )
344  ELSE
345  maxwrk = max( maxwrk, n + ( n - 1 )*ilaenv( 1, 'ZUNGHR',
346  $ ' ', n, 1, n, -1 ) )
347  maxwrk = max( maxwrk, hswork )
348  END IF
349  lwrk = maxwrk
350  IF( .NOT.wantsn )
351  $ lwrk = max( lwrk, ( n*n )/2 )
352  END IF
353  work( 1 ) = lwrk
354 *
355  IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
356  info = -15
357  END IF
358  END IF
359 *
360  IF( info.NE.0 ) THEN
361  CALL xerbla( 'ZGEESX', -info )
362  RETURN
363  ELSE IF( lquery ) THEN
364  RETURN
365  END IF
366 *
367 * Quick return if possible
368 *
369  IF( n.EQ.0 ) THEN
370  sdim = 0
371  RETURN
372  END IF
373 *
374 * Get machine constants
375 *
376  eps = dlamch( 'P' )
377  smlnum = dlamch( 'S' )
378  bignum = one / smlnum
379  CALL dlabad( smlnum, bignum )
380  smlnum = sqrt( smlnum ) / eps
381  bignum = one / smlnum
382 *
383 * Scale A if max element outside range [SMLNUM,BIGNUM]
384 *
385  anrm = zlange( 'M', n, n, a, lda, dum )
386  scalea = .false.
387  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
388  scalea = .true.
389  cscale = smlnum
390  ELSE IF( anrm.GT.bignum ) THEN
391  scalea = .true.
392  cscale = bignum
393  END IF
394  IF( scalea )
395  $ CALL zlascl( 'G', 0, 0, anrm, cscale, n, n, a, lda, ierr )
396 *
397 *
398 * Permute the matrix to make it more nearly triangular
399 * (CWorkspace: none)
400 * (RWorkspace: need N)
401 *
402  ibal = 1
403  CALL zgebal( 'P', n, a, lda, ilo, ihi, rwork( ibal ), ierr )
404 *
405 * Reduce to upper Hessenberg form
406 * (CWorkspace: need 2*N, prefer N+N*NB)
407 * (RWorkspace: none)
408 *
409  itau = 1
410  iwrk = n + itau
411  CALL zgehrd( n, ilo, ihi, a, lda, work( itau ), work( iwrk ),
412  $ lwork-iwrk+1, ierr )
413 *
414  IF( wantvs ) THEN
415 *
416 * Copy Householder vectors to VS
417 *
418  CALL zlacpy( 'L', n, n, a, lda, vs, ldvs )
419 *
420 * Generate unitary matrix in VS
421 * (CWorkspace: need 2*N-1, prefer N+(N-1)*NB)
422 * (RWorkspace: none)
423 *
424  CALL zunghr( n, ilo, ihi, vs, ldvs, work( itau ), work( iwrk ),
425  $ lwork-iwrk+1, ierr )
426  END IF
427 *
428  sdim = 0
429 *
430 * Perform QR iteration, accumulating Schur vectors in VS if desired
431 * (CWorkspace: need 1, prefer HSWORK (see comments) )
432 * (RWorkspace: none)
433 *
434  iwrk = itau
435  CALL zhseqr( 'S', jobvs, n, ilo, ihi, a, lda, w, vs, ldvs,
436  $ work( iwrk ), lwork-iwrk+1, ieval )
437  IF( ieval.GT.0 )
438  $ info = ieval
439 *
440 * Sort eigenvalues if desired
441 *
442  IF( wantst .AND. info.EQ.0 ) THEN
443  IF( scalea )
444  $ CALL zlascl( 'G', 0, 0, cscale, anrm, n, 1, w, n, ierr )
445  DO 10 i = 1, n
446  bwork( i ) = SELECT( w( i ) )
447  10 CONTINUE
448 *
449 * Reorder eigenvalues, transform Schur vectors, and compute
450 * reciprocal condition numbers
451 * (CWorkspace: if SENSE is not 'N', need 2*SDIM*(N-SDIM)
452 * otherwise, need none )
453 * (RWorkspace: none)
454 *
455  CALL ztrsen( sense, jobvs, bwork, n, a, lda, vs, ldvs, w, sdim,
456  $ rconde, rcondv, work( iwrk ), lwork-iwrk+1,
457  $ icond )
458  IF( .NOT.wantsn )
459  $ maxwrk = max( maxwrk, 2*sdim*( n-sdim ) )
460  IF( icond.EQ.-14 ) THEN
461 *
462 * Not enough complex workspace
463 *
464  info = -15
465  END IF
466  END IF
467 *
468  IF( wantvs ) THEN
469 *
470 * Undo balancing
471 * (CWorkspace: none)
472 * (RWorkspace: need N)
473 *
474  CALL zgebak( 'P', 'R', n, ilo, ihi, rwork( ibal ), n, vs, ldvs,
475  $ ierr )
476  END IF
477 *
478  IF( scalea ) THEN
479 *
480 * Undo scaling for the Schur form of A
481 *
482  CALL zlascl( 'U', 0, 0, cscale, anrm, n, n, a, lda, ierr )
483  CALL zcopy( n, a, lda+1, w, 1 )
484  IF( ( wantsv .OR. wantsb ) .AND. info.EQ.0 ) THEN
485  dum( 1 ) = rcondv
486  CALL dlascl( 'G', 0, 0, cscale, anrm, 1, 1, dum, 1, ierr )
487  rcondv = dum( 1 )
488  END IF
489  END IF
490 *
491  work( 1 ) = maxwrk
492  RETURN
493 *
494 * End of ZGEESX
495 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:74
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
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
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:81
double precision function zlange(NORM, M, N, A, LDA, WORK)
ZLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: zlange.f:115
subroutine zgebal(JOB, N, A, LDA, ILO, IHI, SCALE, INFO)
ZGEBAL
Definition: zgebal.f:162
subroutine zgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
ZGEHRD
Definition: zgehrd.f:167
subroutine zgebak(JOB, SIDE, N, ILO, IHI, SCALE, M, V, LDV, INFO)
ZGEBAK
Definition: zgebak.f:131
subroutine zlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
ZLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: zlascl.f:143
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:103
subroutine zhseqr(JOB, COMPZ, N, ILO, IHI, H, LDH, W, Z, LDZ, WORK, LWORK, INFO)
ZHSEQR
Definition: zhseqr.f:299
subroutine ztrsen(JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, W, M, S, SEP, WORK, LWORK, INFO)
ZTRSEN
Definition: ztrsen.f:264
subroutine zunghr(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGHR
Definition: zunghr.f:126
Here is the call graph for this function:
Here is the caller graph for this function: